/[MITgcm]/MITgcm/pkg/seaice/seaice_ocean_stress.F
ViewVC logotype

Annotation of /MITgcm/pkg/seaice/seaice_ocean_stress.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.1 - (hide annotations) (download)
Mon Mar 6 13:17:37 2006 UTC (18 years, 2 months ago) by mlosch
Branch: MAIN 
 - add c-grid version of the dynamic solver:
   + three new routines that are the c-grid substitute lsr, dynsolver
     and ostres
   + put a few fields that were local to dynsolver into global common
     blocks, so that I can move the computation of stresses etc into
     seaice_lsr (saves coding but may break the adjoint; Patrick, I am
     sorry!).
 - replace more hardwired parameters by runtime parameters
 - add ice masks that mask the rhs of the implicit solvers where there
   is no ice (commented out in seaice_dynsolver, because i am not sure
   if this works properly), eventually this should replace the clipping
   of ice velocities in seaice_dynsolver to +/-40cm/s.
1 mlosch 1.1 C $Header: $
2     C $Name: $
3    
4     #include "SEAICE_OPTIONS.h"
5    
6     CStartOfInterface
7     SUBROUTINE SEAICE_OCEAN_STRESS(
8     I myTime, myIter, myThid )
9     C /==========================================================\
10     C | SUBROUTINE SEAICE_OCEAN_STRESS |
11     C | o Calculate ocean surface stresses |
12     C | - C-grid version |
13     C |==========================================================|
14     C \==========================================================/
15     IMPLICIT NONE
16    
17     C === Global variables ===
18     #include "SIZE.h"
19     #include "EEPARAMS.h"
20     #include "PARAMS.h"
21     #include "FFIELDS.h"
22     #include "SEAICE.h"
23     #include "SEAICE_PARAMS.h"
24    
25     C === Routine arguments ===
26     C myTime - Simulation time
27     C myIter - Simulation timestep number
28     C myThid - Thread no. that called this routine.
29     _RL myTime
30     INTEGER myIter
31     INTEGER myThid
32     CML _RL COR_ICE (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
33     CEndOfInterface
34    
35     #ifdef SEAICE_CGRID
36     C === Local variables ===
37     C i,j,bi,bj - Loop counters
38    
39     INTEGER i, j, bi, bj
40     _RL SINWIN, COSWIN, SINWAT, COSWAT
41     #ifdef SEAICE_TEST_ICE_STRESS_1
42     _RL fuIce, fvIce
43     #endif
44    
45     c introduce turning angle (default is zero)
46     SINWIN=SIN(SEAICE_airTurnAngle*deg2rad)
47     COSWIN=COS(SEAICE_airTurnAngle*deg2rad)
48     SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
49     COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
50    
51     CML#ifdef SEAICE_ORIGINAL_BAD_ICE_STRESS
52     CMLC-- Following formulation is problematic and is no longer used.
53     CML#ifdef SEAICE_ALLOW_DYNAMICS
54     CML IF ( SEAICEuseDYNAMICS ) THEN
55     CMLC-- Compute ice-affected wind stress
56     CML DO bj=myByLo(myThid),myByHi(myThid)
57     CML DO bi=myBxLo(myThid),myBxHi(myThid)
58     CML DO j=1,sNy
59     CML DO i=1,sNx
60     CML WINDX(I,J,bi,bj)=DWATN(I,J,bi,bj)
61     CML & *(COSWAT*(GWATX(I,J,bi,bj)-UICE(I,J,1,bi,bj))
62     CML & -SINWAT*(GWATY(I,J,bi,bj)-VICEC(I,J,bi,bj)))
63     CML WINDY(I,J,bi,bj)=DWATN(I,J,bi,bj)
64     CML & *(SINWAT*(GWATX(I,J,bi,bj)-UICEC(I,J,bi,bj))
65     CML & +COSWAT*(GWATY(I,J,bi,bj)-VICE(I,J,1,bi,bj)))
66     CML WINDX(I,J,bi,bj)=WINDX(I,J,bi,bj)-( COR_ICE(I,J,bi,bj)
67     CML & *GWATY(I,J,bi,bj)-COR_ICE(I,J,bi,bj)*VICEC(I,J,bi,bj))
68     CML WINDY(I,J,bi,bj)=WINDY(I,J,bi,bj)-(-COR_ICE(I,J,bi,bj)
69     CML & *GWATX(I,J,bi,bj)+COR_ICE(I,J,bi,bj)*UICEC(I,J,bi,bj))
70     CML WINDX(I,J,bi,bj)=WINDX(I,J,bi,bj)-(UICE(I,J,1,bi,bj)
71     CML & -UICE(I,J,3,bi,bj))*AMASS(I,J,bi,bj)/SEAICE_DT*TWO
72     CML WINDY(I,J,bi,bj)=WINDY(I,J,bi,bj)-(VICE(I,J,1,bi,bj)
73     CML & -VICE(I,J,3,bi,bj))*AMASS(I,J,bi,bj)/SEAICE_DT*TWO
74     CML ENDDO
75     CML ENDDO
76     CML ENDDO
77     CML ENDDO
78     CML DO bj=myByLo(myThid),myByHi(myThid)
79     CML DO bi=myBxLo(myThid),myBxHi(myThid)
80     CML DO j=1,sNy
81     CML DO i=1,sNx
82     CML WINDX(I,J,bi,bj)=-WINDX(I,J,bi,bj)
83     CML WINDY(I,J,bi,bj)=-WINDY(I,J,bi,bj)
84     CML ENDDO
85     CML ENDDO
86     CML ENDDO
87     CML ENDDO
88     CML ENDIF
89     CML#endif /* SEAICE_ALLOW_DYNAMICS */
90     CML#endif /* SEAICE_ORIGINAL_BAD_ICE_STRESS */
91    
92     C-- Update overlap regions
93     CALL EXCH_UV_XY_RL(WINDX, WINDY, .TRUE., myThid)
94    
95     #ifndef SEAICE_EXTERNAL_FLUXES
96     C-- Interpolate wind stress (N/m^2) from South-West B-grid
97     C to South-West C-grid for forcing ocean model.
98     DO bj=myByLo(myThid),myByHi(myThid)
99     DO bi=myBxLo(myThid),myBxHi(myThid)
100     DO j=1,sNy
101     DO i=1,sNx
102     fu(I,J,bi,bj)=WINDX(I,J,bi,bj)
103     fv(I,J,bi,bj)=WINDY(I,J,bi,bj)
104     ENDDO
105     ENDDO
106     ENDDO
107     ENDDO
108     CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
109     #endif /* ifndef SEAICE_EXTERNAL_FLUXES */
110    
111     #ifdef SEAICE_TEST_ICE_STRESS_1
112     C-- Compute ice-affected wind stress (interpolate to U/V-points)
113     DO bj=myByLo(myThid),myByHi(myThid)
114     DO bi=myBxLo(myThid),myBxHi(myThid)
115     DO j=1,sNy
116     DO i=1,sNx
117     fuIce=HALF*( DWATN(I,J,bi,bj)+DWATN(I,J+1,bi,bj) )*(
118     & COSWAT *
119     & ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )
120     & - SINWAT* 0.5 _d 0 * (
121     & 0.5 _d 0*(vIce(I ,J ,1,bi,bj)-GWATY(I ,J ,bi,bj)
122     & +vIce(I-1,J ,1,bi,bj)-GWATY(I-1,J ,bi,bj))
123     & +0.5 _d 0*(vIce(I ,J+1,1,bi,bj)-GWATY(I ,J+1,bi,bj)
124     & +vIce(I-1,J+1,1,bi,bj)-GWATY(I-1,J+1,bi,bj)) )
125     & )
126     fvIce=HALF*( DWATN(I,J,bi,bj)+DWATN(I+1,J,bi,bj) )*(
127     & SINWAT *
128     & ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )
129     & + COSWAT * 0.5 _d 0 * (
130     & 0.5 _d 0*(uIce(I ,J ,1,bi,bj)-GWATY(I ,J ,bi,bj)
131     & +uIce(I+1,J ,1,bi,bj)-GWATX(I+1,J ,bi,bj))
132     & +0.5 _d 0*(uIce(I ,J-1,1,bi,bj)-GWATY(I ,J-1,bi,bj)
133     & +uIce(I+1,J-1,1,bi,bj)-GWATX(I+1,J-1,bi,bj)) )
134     & )
135     fu(I,J,bi,bj)=(ONE-AREA(I,J,1,bi,bj))*fu(I,J,bi,bj)+
136     & AREA(I,J,1,bi,bj)*fuIce
137     fv(I,J,bi,bj)=(ONE-AREA(I,J,1,bi,bj))*fv(I,J,bi,bj)+
138     & AREA(I,J,1,bi,bj)*fvIce
139     ENDDO
140     ENDDO
141     ENDDO
142     ENDDO
143     CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
144     #endif /* SEAICE_TEST_ICE_STRESS_1 */
145     #endif /* not SEAICE_CGRID */
146    
147     RETURN
148     END
  ViewVC Help
Powered by ViewVC 1.1.22