/[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.8 - (hide annotations) (download)
Mon Mar 20 21:36:11 2006 UTC (18 years, 2 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint58e_post, checkpoint58h_post, checkpoint58f_post, checkpoint58d_post, checkpoint58c_post, checkpoint58g_post
Changes since 1.7: +4 -1 lines 
  seaice: add an EVP solver following Hunke and Dukowicz (1997) and the
  documentation of CICE. Turn on by defining SEAICE_ALLOW_EVP in
  SEAICE_OPTIONS.h and SEAICEuseEVP=.true. in data.seaice. Works only
  with SEAICE_CGRID defined.
  Use at own risk.
1 mlosch 1.8 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_ocean_stress.F,v 1.7 2006/03/17 15:53:38 mlosch Exp $
2 mlosch 1.1 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 mlosch 1.5 #include "GRID.h"
22 mlosch 1.1 #include "FFIELDS.h"
23     #include "SEAICE.h"
24     #include "SEAICE_PARAMS.h"
25 mlosch 1.5 #include "SEAICE_FFIELDS.h"
26 mlosch 1.1
27     C === Routine arguments ===
28     C myTime - Simulation time
29     C myIter - Simulation timestep number
30     C myThid - Thread no. that called this routine.
31     _RL myTime
32     INTEGER myIter
33     INTEGER myThid
34     CEndOfInterface
35    
36     #ifdef SEAICE_CGRID
37     C === Local variables ===
38     C i,j,bi,bj - Loop counters
39    
40     INTEGER i, j, bi, bj
41 mlosch 1.5 _RL SINWAT, COSWAT, SINWIN, COSWIN
42     _RL fuIce, fvIce, FX, FY
43 mlosch 1.4 _RL areaW, areaS
44 mlosch 1.1
45 mlosch 1.5 _RL press (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
46     _RL etaPlusZeta (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47     _RL zetaMinusEta(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48     _RL etaMeanZ (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49     _RL etaMeanU (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50     _RL etaMeanV (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
51     _RL dVdx (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
52     _RL dVdy (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
53     _RL dUdx (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
54     _RL dUdy (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
55    
56 mlosch 1.1 c introduce turning angle (default is zero)
57     SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
58     COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
59 mlosch 1.5 SINWIN=SIN(SEAICE_airTurnAngle*deg2rad)
60     COSWIN=COS(SEAICE_airTurnAngle*deg2rad)
61 mlosch 1.1
62     C-- Update overlap regions
63     CALL EXCH_UV_XY_RL(WINDX, WINDY, .TRUE., myThid)
64    
65     #ifndef SEAICE_EXTERNAL_FLUXES
66 mlosch 1.3 C-- Interpolate wind stress (N/m^2) from C-points of C-grid
67     C to U and V points of C-grid for forcing the ocean model.
68 mlosch 1.1 DO bj=myByLo(myThid),myByHi(myThid)
69     DO bi=myBxLo(myThid),myBxHi(myThid)
70     DO j=1,sNy
71     DO i=1,sNx
72 mlosch 1.3 fu(I,J,bi,bj)=0.5*(WINDX(I,J,bi,bj) + WINDX(I-1,J,bi,bj))
73     fv(I,J,bi,bj)=0.5*(WINDY(I,J,bi,bj) + WINDY(I,J-1,bi,bj))
74 mlosch 1.1 ENDDO
75     ENDDO
76     ENDDO
77     ENDDO
78     #endif /* ifndef SEAICE_EXTERNAL_FLUXES */
79    
80 mlosch 1.5 IF ( useHB87StressCoupling ) THEN
81     C
82     C use an intergral over ice and ocean surface layer to define
83     C surface stresses on ocean following Hibler and Bryan (1987, JPO)
84     C
85     C recompute viscosities from updated ice velocities
86     CALL SEAICE_CALC_VISCOSITIES(
87     I uIce(1-Olx,1-Oly,1,1,1), vIce(1-Olx,1-Oly,1,1,1),
88     I zMin, zMax, hEffM, press0,
89     O eta, zeta, press,
90 mlosch 1.8 #ifdef SEAICE_ALLOW_EVP
91     O seaice_div, seaice_tension, seaice_shear,
92     #endif /* SEAICE_ALLOW_EVP */
93 mlosch 1.5 I myThid )
94     C re-compute internal stresses with updated ice velocities
95     DO bj=myByLo(myThid),myByHi(myThid)
96     DO bi=myBxLo(myThid),myBxHi(myThid)
97     DO j=1-Oly+1,sNy+Oly-1
98     DO i=1-Olx+1,sNx+Olx-1
99     etaPlusZeta (I,J) = eta(I,J,bi,bj) + zeta(I,J,bi,bj)
100     zetaMinusEta(I,J) = zeta(I,J,bi,bj) - eta(I,J,bi,bj)
101     etaMeanU (I,J) =
102     & HALF*(ETA (I,J,bi,bj) + ETA (I-1,J ,bi,bj))
103     etaMeanV (I,J) =
104     & HALF*(ETA (I,J,bi,bj) + ETA (I ,J-1,bi,bj))
105     etaMeanZ (I,J) = QUART *
106     & ( eta(I ,J,bi,bj) + eta(I ,J-1,bi,bj)
107     & + eta(I-1,J,bi,bj) + eta(I-1,J-1,bi,bj) )
108     dUdx(I,J) = ( uIce(I+1,J,1,bi,bj) - uIce(I,J,1,bi,bj) )
109     & * _recip_dxF(I,J,bi,bj)
110     dUdy(I,J) = ( uIce(I,J+1,1,bi,bj) - uIce(I,J,1,bi,bj) )
111     & * _recip_dyU(I,J+1,bi,bj)
112     dVdx(I,J) = ( vIce(I+1,J,1,bi,bj) - vIce(I,J,1,bi,bj) )
113     & * _recip_dxV(I+1,J,bi,bj)
114     dVdy(I,J) = ( vIce(I,J+1,1,bi,bj) - vIce(I,J,1,bi,bj) )
115     & * _recip_dyF(I,J,bi,bj)
116     ENDDO
117     ENDDO
118     DO J = 1,sNy
119     DO I = 1,sNx
120 mlosch 1.7 C First FX = (d/dx)*sigma
121 mlosch 1.5 C + d/dx[ eta+zeta d/dx ] U
122     FX = _recip_dxC(I,J,bi,bj) *
123     & ( etaPlusZeta(I ,J) * dUdx(I ,J)
124     & - etaPlusZeta(I-1,J) * dUdx(I-1,J) )
125     C + (d/dy)[eta*(d/dy + tanphi/a)] U (also on UVRT1/2)
126     FX = FX + _recip_dyG(I,J,bi,bj) * (
127     & ( etaMeanZ(I,J+1) * dUdy(I,J+1)
128     & - etaMeanZ(I,J ) * dUdy(I,J )
129     & )
130     & - ( etaMeanZ(I,J+1)
131     & * ( uIce(I,J+1,1,bi,bj)+uIce(I,J,1,bi,bj) )
132     & - etaMeanZ(I,J )
133     & * ( uIce(I,J-1,1,bi,bj)+uIce(I,J,1,bi,bj) ) )
134     & * 0.5 _d 0 * _tanPhiAtU(I,J,bi,bj)
135     & * recip_rSphere )
136     C - 2*eta*(tanphi/a) * ( tanphi/a ) U
137     FX = FX - TWO * uIce(I,J,1,bi,bj)
138     & * etaMeanU(I,J)*recip_rSphere*recip_rSphere
139     & * _tanPhiAtU(I,J,bi,bj) * _tanPhiAtU(I,J,bi,bj)
140     C + d/dx[ (zeta-eta) dV/dy]
141     FX = FX +
142     & ( zetaMinusEta(I ,J ) * dVdy(I ,J )
143     & - zetaMinusEta(I-1,J ) * dVdy(I-1,J )
144     & ) * _recip_dxC(I,J,bi,bj)
145     C + d/dy[ eta dV/x ]
146     FX = FX + (
147     & etaMeanZ(I,J+1)
148     & * ( vIce(I ,J+1,1,bi,bj) - vIce(I-1,J+1,1,bi,bj) )
149     & * _recip_dxV(I,J+1,bi,bj)
150     & - etaMeanZ(I,J )
151     & * ( vIce(I ,J,1,bi,bj) - vIce(I-1,J,1,bi,bj) )
152     & * _recip_dxV(I,J,bi,bj)
153     & ) * _recip_dyG(I,J,bi,bj)
154     C - d/dx[ (eta+zeta) * v * (tanphi/a) ]
155     FX = FX - (
156     & etaPlusZeta(I ,J)
157     & * 0.5 * (vIce(I ,J,1,bi,bj)+vIce(I ,J+1,1,bi,bj))
158     & * 0.5 * ( _tanPhiAtU(I ,J,bi,bj)
159     & + _tanPhiAtU(I+1,J,bi,bj) )
160     & - etaPlusZeta(I-1,J) *
161     & * 0.5 * (vIce(I-1,J,1,bi,bj)+vIce(I-1,J+1,1,bi,bj))
162     & * 0.5 * ( _tanPhiAtU(I-1,J,bi,bj)
163     & + _tanPhiAtU(I ,J,bi,bj) )
164     & )* _recip_dxC(I,J,bi,bj)*recip_rSphere
165     C - 2*eta*(tanphi/a) * dV/dx
166     FX = FX
167     & -TWO * etaMeanU(I,J) * _tanPhiAtV(I,J,bi,bj)
168     & *recip_rSphere
169     & *(vIce(I ,J,1,bi,bj) + vIce(I ,J+1,1,bi,bj)
170     & - vIce(I-1,J,1,bi,bj) - vIce(I-1,J+1,1,bi,bj))
171     & * _recip_dxC(I,J,bi,bj)
172     C - (d/dx) P/2
173     FX = _maskW(I,J,1,bi,bj) * ( FX - _recip_dxC(I,J,bi,bj)
174     & * ( press(I,J,bi,bj) - press(I-1,J,bi,bj) ) )
175     C
176 mlosch 1.7 C then FY = (d/dy)*sigma
177 mlosch 1.5 C + d/dy [(eta+zeta) d/dy] V
178     FY = _recip_dyC(I,J,bi,bj) *
179     & ( dVdy(I,J ) * etaPlusZeta(I,J )
180     & - dVdy(I,J-1) * etaPlusZeta(I,J-1) )
181     C + d/dx [eta d/dx] V
182     FY = FY + _recip_dxC(I,J,bi,bj) *
183     & ( eta(I ,J,bi,bj) * dVdx(I ,J)
184     & - eta(I-1,J,bi,bj) * dVdx(I-1,J) )
185     C - d/dy [(zeta-eta) tanphi/a] V
186     FY = FY - _recip_dyC(I,J,bi,bj) * recip_rSphere * (
187     & zetaMinusEta(I,J ) * tanPhiAtU(I,J ,bi,bj)
188     & * 0.5 * ( vIce(I,J,1,bi,bj) + vIce(I,J+1,1,bi,bj))
189     & - zetaMinusEta(I,J-1) * tanPhiAtU(I,J-1,bi,bj)
190     & * 0.5 * ( vIce(I,J,1,bi,bj) + vIce(I,J-1,1,bi,bj)) )
191     C 2*eta tanphi/a ( - tanphi/a - d/dy) V
192     FY = FY - TWO*etaMeanV(I,J) * recip_rSphere
193     & * _tanPhiAtV(I,J,bi,bj) * (
194     & _tanPhiAtV(I,J,bi,bj) * recip_rSphere
195     & + _recip_dyC(I,J,bi,bj) *
196     & ( 0.5 * ( vIce(I,J,1,bi,bj) + vIce(I,J+1,1,bi,bj))
197     & - 0.5 * ( vIce(I,J,1,bi,bj) + vIce(I,J-1,1,bi,bj)) ) )
198     C + d/dy[ (zeta-eta) dU/dx ]
199     FY = FY +
200     & ( zetaMinusEta(I,J )*dUdx(I,J )
201     & - zetaMinusEta(I,J-1)*dUdx(I,J-1) )
202     & * _recip_dyC(I,J,bi,bj)
203     C + d/dx[ eta dU/dy ]
204     FY = FY + _recip_dxG(I,J,bi,bj) *
205     & ( etaMeanZ(I+1,J) * dUdy(I+1,J)
206     & - etaMeanZ(I ,J) * dUdy(I ,J) )
207     C + d/dx[ eta * (tanphi/a) * U ]
208     FY = FY + (
209     & etaMeanZ(I+1,J) * 0.5 *
210     & ( uIce(I+1,J ,1,bi,bj) * _tanPhiAtU(I+1,J ,bi,bj)
211     & + uIce(I+1,J-1,1,bi,bj) * _tanPhiAtU(I+1,J-1,bi,bj) )
212     & - etaMeanZ(I ,J) * 0.5 *
213     & ( uIce(I ,J ,1,bi,bj) * _tanPhiAtU(I ,J ,bi,bj)
214     & + uIce(I ,J-1,1,bi,bj) * _tanPhiAtU(I ,J ,bi,bj) )
215     & ) * _recip_dxG(I,J,bi,bj)*recip_rSphere
216     C + 2*eta*(tanphi/a) dU/dx
217     FY = FY +
218     & TWO * etaMeanV(I,J)*TWO * _tanPhiAtV(I,J,bi,bj)
219     & * ( uIce(I+1,J,1,bi,bj)+uIce(I+1,J-1,1,bi,bj)
220     & - uIce(I ,J,1,bi,bj)-uIce(I ,J-1,1,bi,bj) )
221     & * _recip_dxG(I,J,bi,bj) * recip_rSphere
222     C - (d/dy) P/2
223     FY = _maskS(I,J,1,bi,bj) * ( FY - _recip_dyC(I,J,bi,bj)
224     & * ( press(I,J,bi,bj) - press(I,J-1,bi,bj) ) )
225     C
226     C recompute wind stress over ice (done already in seaice_dynsolver,
227     C but not saved)
228     fuIce = 0.5 _d 0 *
229     & ( DAIRN(I ,J,bi,bj)*(
230 mlosch 1.6 & COSWIN*uWind(I ,J,bi,bj)
231     & -SIGN(SINWIN, _fCori(I ,J,bi,bj))*vWind(I ,J,bi,bj) )
232 mlosch 1.5 & + DAIRN(I-1,J,bi,bj)*(
233 mlosch 1.6 & COSWIN*uWind(I-1,J,bi,bj)
234     & -SIGN(SINWIN, _fCori(I-1,J,bi,bj))*vWind(I-1,J,bi,bj) )
235 mlosch 1.5 & )
236     fvIce = 0.5 _d 0 *
237     & ( DAIRN(I,J ,bi,bj)*(
238 mlosch 1.6 & SIGN(SINWIN, _fCori(I ,J,bi,bj))*uWind(I,J ,bi,bj)
239     & +COSWIN*vWind(I,J ,bi,bj) )
240 mlosch 1.5 & + DAIRN(I,J-1,bi,bj)*(
241 mlosch 1.6 & SIGN(SINWIN, _fCori(I,J-1,bi,bj))*uWind(I,J-1,bi,bj)
242     & +COSWIN*vWind(I,J-1,bi,bj) )
243 mlosch 1.5 & )
244     C average wind stress over ice and ocean and apply averaged wind
245     C stress and internal ice stresses to surface layer of ocean
246     areaW = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I-1,J,1,bi,bj))
247     areaS = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))
248 mlosch 1.7 fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)+areaW*fuIce + FX
249     fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)+areaS*fvIce + FY
250 mlosch 1.5 END DO
251     END DO
252     ENDDO
253     ENDDO
254     ELSE
255    
256     C-- Compute ice-affected wind stress (interpolate to U/V-points)
257     C by averaging wind stress and ice-ocean stress according to
258     C ice cover
259 mlosch 1.1 DO bj=myByLo(myThid),myByHi(myThid)
260     DO bi=myBxLo(myThid),myBxHi(myThid)
261     DO j=1,sNy
262     DO i=1,sNx
263 mlosch 1.6 fuIce=HALF*( DWATN(I,J,bi,bj)+DWATN(I,J+1,bi,bj) )*
264 mlosch 1.1 & COSWAT *
265     & ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )
266 mlosch 1.6 & - SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *
267     & ( DWATN(I ,J,bi,bj) *
268     & 0.5 _d 0*(vIce(I ,J ,1,bi,bj)-GWATY(I ,J ,bi,bj)
269     & +vIce(I ,J+1,1,bi,bj)-GWATY(I ,J+1,bi,bj))
270     & + DWATN(I-1,J,bi,bj) *
271     & 0.5 _d 0*(vIce(I-1,J ,1,bi,bj)-GWATY(I-1,J ,bi,bj)
272     & +vIce(I-1,J+1,1,bi,bj)-GWATY(I-1,J+1,bi,bj))
273 mlosch 1.1 & )
274 mlosch 1.6 fvIce=HALF*( DWATN(I,J,bi,bj)+DWATN(I+1,J,bi,bj) )*
275     & COSWAT *
276     & ( VICE(I,J,1,bi,bj)-GWATY(I,J,bi,bj) )
277     & + SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *
278     & ( DWATN(I,J ,bi,bj) *
279     & 0.5 _d 0*(uIce(I ,J ,1,bi,bj)-GWATX(I ,J ,bi,bj)
280     & +uIce(I+1,J ,1,bi,bj)-GWATX(I+1,J ,bi,bj))
281     & + DWATN(I,J-1,bi,bj) *
282     & 0.5 _d 0*(uIce(I ,J-1,1,bi,bj)-GWATX(I ,J-1,bi,bj)
283     & +uIce(I+1,J-1,1,bi,bj)-GWATX(I+1,J-1,bi,bj))
284 mlosch 1.1 & )
285 mlosch 1.4 areaW = 0.5 _d 0 * (AREA(I,J,1,bi,bj) + AREA(I-1,J,1,bi,bj))
286     areaS = 0.5 _d 0 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))
287     fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)+areaW*fuIce
288     fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)+areaS*fvIce
289 mlosch 1.1 ENDDO
290     ENDDO
291     ENDDO
292     ENDDO
293 mlosch 1.5 ENDIF
294 mlosch 1.1 CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
295 mlosch 1.3
296 mlosch 1.1 #endif /* not SEAICE_CGRID */
297    
298     RETURN
299     END
  ViewVC Help
Powered by ViewVC 1.1.22