/[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.13 - (hide annotations) (download)
Fri Apr 20 18:29:58 2007 UTC (17 years, 1 month ago) by mlosch
Branch: MAIN
Changes since 1.12: +13 -9 lines
fix bug in evp solver, while doing that
  - change (improve?) discretization of stress computations for both
    cgrid-lsr and evp solver
  - add a new routine to compute strain rates
  - changes lab_sea results

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

  ViewVC Help
Powered by ViewVC 1.1.22