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

Contents of /MITgcm/pkg/seaice/seaice_lsr.F

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


Revision 1.3 - (show annotations) (download)
Fri Mar 10 16:06:40 2006 UTC (18 years, 2 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint58b_post
Changes since 1.2: +2 -2 lines
o another bug fix: wrong index for grid spacing

1 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_lsr.F,v 1.2 2006/03/10 13:59:40 mlosch Exp $
2 C $Name: $
3
4 #include "SEAICE_OPTIONS.h"
5
6 CStartOfInterface
7 SUBROUTINE SEAICE_LSR( ilcall, myThid )
8 C /==========================================================\
9 C | SUBROUTINE SEAICE_LSR |
10 C | o Solve ice momentum equation with an LSR dynamics solver|
11 C | (see Zhang and Hibler, JGR, 102, 8691-8702, 1997 |
12 C | and Zhang and Rothrock, MWR, 131, 845- 861, 2003) |
13 C | Written by Jinlun Zhang, PSC/UW, Feb-2001 |
14 C | zhang@apl.washington.edu |
15 C |==========================================================|
16 C | C-grid version by Martin Losch |
17 C \==========================================================/
18 IMPLICIT NONE
19
20 C === Global variables ===
21 #include "SIZE.h"
22 #include "EEPARAMS.h"
23 #include "PARAMS.h"
24 #include "GRID.h"
25 #include "SEAICE.h"
26 #include "SEAICE_PARAMS.h"
27
28 #ifdef ALLOW_AUTODIFF_TAMC
29 # include "tamc.h"
30 #endif
31
32 C === Routine arguments ===
33 C myThid - Thread no. that called this routine.
34 INTEGER ilcall
35 INTEGER myThid
36 CEndOfInterface
37
38 #ifdef SEAICE_CGRID
39 #ifdef SEAICE_ALLOW_DYNAMICS
40
41 C === Local variables ===
42 C i,j,bi,bj - Loop counters
43
44 INTEGER i, j, m, bi, bj, j1, j2, im, jm
45 INTEGER ICOUNT1, ICOUNT2, SOLV_MAX_ITERS, SOLV_NCHECK
46 INTEGER phexit
47
48 _RL WFAU, WFAV, WFAU1, WFAV1, WFAU2, WFAV2
49 _RL AA1, AA2, AA3, AA4, AA5, AA6, AA7, S1, S2, S1A, S2A
50
51 _RL AU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
52 _RL BU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
53 _RL CU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
54 _RL AV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
55 _RL BV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
56 _RL CV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
57 _RL UERR (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
58 _RL FXY (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
59
60 _RL URT(1-Olx:sNx+Olx), CUU(1-Olx:sNx+Olx)
61 _RL VRT(1-Oly:sNy+Oly), CVV(1-Oly:sNy+Oly)
62
63 _RL etaPlusZeta (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
64 _RL zetaMinusEta(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
65 _RL etaMeanZ (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
66 _RL etaMeanU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
67 _RL etaMeanV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
68
69 _RL UVRT1 (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
70 _RL UVRT2 (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
71
72 _RL dVdy (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
73 _RL dUdx (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
74 _RL dUdy (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
75
76 _RL SINWAT, COSWAT
77 _RL ECCEN, ECM2, DELT1, DELT2
78 _RL TEMPVAR
79
80 _RL PRESS (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
81 _RL E11 (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
82 _RL E22 (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
83 _RL E12 (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
84
85 C-- FIRST SET UP BASIC CONSTANTS
86 ECCEN=SEAICE_eccen
87 ECM2=ONE/(ECCEN**2)
88
89 C-- introduce turning angles
90 SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
91 COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
92
93 C SET SOME VALUES
94 WFAU1=0.95 _d 0
95 WFAV1=0.95 _d 0
96 WFAU2=ZERO
97 WFAV2=ZERO
98
99 S1A=0.80 _d 0
100 S2A=0.80 _d 0
101 WFAU=WFAU1
102 WFAV=WFAV1
103
104 SOLV_MAX_ITERS=1500
105 SOLV_NCHECK=2
106
107 ICOUNT1=SOLV_MAX_ITERS
108 ICOUNT2=SOLV_MAX_ITERS
109
110 #ifdef ALLOW_AUTODIFF_TAMC
111 cph That's an important one! Note, that
112 cph * lsr is called twice, thus the icall index
113 cph * this storing is still outside the iteration loop
114 CADJ STORE uice = comlev1_lsr,
115 CADJ & key = ikey_dynamics + (ilcall-1)*nchklev_1
116 CADJ STORE vice = comlev1_lsr,
117 CADJ & key = ikey_dynamics + (ilcall-1)*nchklev_1
118 #endif /* ALLOW_AUTODIFF_TAMC */
119
120 DO bj=myByLo(myThid),myByHi(myThid)
121 DO bi=myBxLo(myThid),myBxHi(myThid)
122 DO j=1,sNy
123 DO i=1,sNx
124 C NOW EVALUATE STRAIN RATES
125 E11(I,J,bi,bj)= _recip_dxF(I,J,bi,bj) *
126 & (uIceC(I,J,bi,bj)-uIceC(I-1,J ,bi,bj))
127 & -HALF*
128 & (vIceC(I,J,bi,bj)+vIceC(I ,J+1,bi,bj))
129 & * _tanPhiAtU(I,J,bi,bj)*recip_rSphere
130 E22(I,J,bi,bj)= _recip_dyF(I,J,bi,bj) *
131 & (vIceC(I,J,bi,bj)+vIceC(I,J-1,bi,bj))
132 E12(I,J,bi,bj)=HALF*(
133 & (uIceC(I,J+1,bi,bj)+uIceC(I+1,J+1,bi,bj)
134 & -uIceC(I,J-1,bi,bj)-uIceC(I+1,J-1,bi,bj))
135 & * 1. _d 0 / (dyC(I,J,bi,bj) + dyC(I,J-1,bi,bj))
136 & +
137 & (vIceC(I+1,J+1,bi,bj)+vIceC(I+1,J,bi,bj)
138 & -vIceC(I-1,J+1,bi,bj)-vIceC(I-1,J,bi,bj))
139 & * 1. _d 0 / (dxC(I,J,bi,bj) + dxC(I-1,J,bi,bj))
140 & +HALF*
141 & (uIceC(I, J, bi,bj)+uIceC(I+1,J, bi,bj))
142 & * _tanPhiAtU(I,J,bi,bj)*recip_rSphere)
143 C NOW EVALUATE VISCOSITIES
144 DELT1=(E11(I,J,bi,bj)**2+E22(I,J,bi,bj)**2)*(ONE+ECM2)
145 & +4. _d 0*ECM2*E12(I,J,bi,bj)**2
146 & +TWO*E11(I,J,bi,bj)*E22(I,J,bi,bj)*(ONE-ECM2)
147 IF ( DELT1 .LE. SEAICE_EPS_SQ ) THEN
148 DELT2=SEAICE_EPS
149 ELSE
150 DELT2=SQRT(DELT1)
151 ENDIF
152 ZETA(I,J,bi,bj)=HALF*PRESS0(I,J,bi,bj)/DELT2
153 C NOW PUT MIN AND MAX VISCOSITIES IN
154 ZETA(I,J,bi,bj)=MIN(ZMAX(I,J,bi,bj),ZETA(I,J,bi,bj))
155 ZETA(I,J,bi,bj)=MAX(ZMIN(I,J,bi,bj),ZETA(I,J,bi,bj))
156 C NOW SET VISCOSITIES TO ZERO AT HEFFMFLOW PTS
157 ZETA(I,J,bi,bj)=ZETA(I,J,bi,bj)*HEFFM(I,J,bi,bj)
158 ETA(I,J,bi,bj)=ECM2*ZETA(I,J,bi,bj)
159 PRESS(I,J,bi,bj)=TWO*ZETA(I,J,bi,bj)*DELT2
160 ENDDO
161 ENDDO
162 ENDDO
163 ENDDO
164
165 C-- Update overlap regions
166 _EXCH_XY_R8(ETA, myThid)
167 _EXCH_XY_R8(ZETA, myThid)
168 _EXCH_XY_R8(PRESS, myThid)
169
170 DO bj=myByLo(myThid),myByHi(myThid)
171 DO bi=myBxLo(myThid),myBxHi(myThid)
172 DO j=1,sNy
173 DO i=1,sNx
174 C NOW SET UP NON-LINEAR WATER DRAG, FORCEX, FORCEY
175 TEMPVAR = QUART*(
176 & (uIceC(I ,J,bi,bj)-GWATX(I ,J,bi,bj)
177 & +uIceC(I+1,J,bi,bj)-GWATX(I+1,J,bi,bj))**2
178 & +(vIceC(I,J ,bi,bj)-GWATY(I,J ,bi,bj)
179 & +vIceC(I,J+1,bi,bj)-GWATY(I,J+1,bi,bj))**2)
180 IF ( TEMPVAR .LE. (QUART/SEAICE_waterDrag)**2 ) THEN
181 DWATN(I,J,bi,bj)=QUART
182 ELSE
183 DWATN(I,J,bi,bj)=SEAICE_waterDrag*SQRT(TEMPVAR)
184 ENDIF
185 C NOW SET UP SYMMETTRIC DRAG
186 DRAGS(I,J,bi,bj)=DWATN(I,J,bi,bj)*COSWAT
187 C NOW SET UP ANTI SYMMETTRIC DRAG PLUS CORIOLIS
188 CML DRAGA(I,J,bi,bj)=DWATN(I,J,bi,bj)*SINWAT+COR_ICE(I,J,bi,bj)
189 DRAGA(I,J,bi,bj)=DWATN(I,J,bi,bj)*SINWAT
190 & + seaiceMassC(I,J,bi,bj) * _fCori(I,J,bi,bj)
191 C NOW ADD IN CURRENT FORCE ( remember to average to correct velocity points )
192 FORCEX(I,J,bi,bj)=FORCEX0(I,J,bi,bj)+DWATN(I,J,bi,bj) *
193 & (COSWAT * GWATX(I,J,bi,bj) - SINWAT* 0.5 _d 0 * (
194 & 0.5 _d 0 * (GWATY(I,J ,bi,bj)+GWATY(I-1,J ,bi,bj))
195 & +0.5 _d 0 * (GWATY(I,J+1,bi,bj)+GWATY(I-1,J+1,bi,bj)) )
196 & )
197 FORCEY(I,J,bi,bj)=FORCEY0(I,J,bi,bj)+DWATN(I,J,bi,bj) *
198 & (SINWAT * GWATX(I,J,bi,bj) + COSWAT * 0.5 _d 0 * (
199 & 0.5 _d 0 * (GWATY(I,J ,bi,bj) + GWATX(I+1,J ,bi,bj))
200 & +0.5 _d 0 * (GWATY(I,J-1,bi,bj) + GWATX(I+1,J-1,bi,bj)) )
201 & )
202 C NOW CALCULATE PRESSURE FORCE AND ADD TO EXTERNAL FORCE
203 FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
204 & - _recip_dxC(I,J,bi,bj)
205 & *(PRESS(I, J,bi,bj) - PRESS(I-1,J,bi,bj))
206 FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)
207 & - _recip_dyC(I,J,bi,bj)
208 & *(PRESS(I,J, bi,bj) - PRESS(I,J-1,bi,bj))
209 ENDDO
210 ENDDO
211 ENDDO
212 ENDDO
213
214 DO bj=myByLo(myThid),myByHi(myThid)
215 DO bi=myBxLo(myThid),myBxHi(myThid)
216 DO j=1-Oly,sNy+Oly
217 DO i=1-Olx,sNx+Olx
218 FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
219 & +seaiceMassU(I,J,bi,bj)/SEAICE_deltaTdyn
220 & *uIce(I,J,2,bi,bj)
221 FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)
222 & +seaiceMassV(I,J,bi,bj)/SEAICE_deltaTdyn
223 & *vIce(I,J,2,bi,bj)
224 FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)* _maskW(I,J,1,bi,bj)
225 FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)* _maskS(I,J,1,bi,bj)
226 etaPlusZeta (I,J,bi,bj) = ETA (I,J,bi,bj)+ZETA(I,J,bi,bj)
227 zetaMinusEta(I,J,bi,bj) = ZETA(I,J,bi,bj)-ETA (I,J,bi,bj)
228 ENDDO
229 ENDDO
230 DO j=1-Oly+1,sNy+Oly
231 DO i=1-Olx+1,sNx+Olx
232 etaMeanU (I,J,bi,bj) =
233 & HALF*(ETA (I,J,bi,bj) + ETA (I-1,J ,bi,bj))
234 etaMeanV (I,J,bi,bj) =
235 & HALF*(ETA (I,J,bi,bj) + ETA (I ,J-1,bi,bj))
236 ENDDO
237 ENDDO
238 DO j=1-Oly+1,sNy+Oly
239 DO i=1-Olx+1,sNx+Olx
240 etaMeanZ (I,J,bi,bj) =
241 & HALF * ( etaMeanU(I,J,bi,bj) + etaMeanU(I,J-1,bi,bj) )
242 ENDDO
243 ENDDO
244 ENDDO
245 ENDDO
246
247 C SOLVE FOR uIce
248 DO bj=myByLo(myThid),myByHi(myThid)
249 DO bi=myBxLo(myThid),myBxHi(myThid)
250
251 DO J=1,sNy
252 DO I=1,sNx
253 C coefficients of uIce(I,J)
254 C (d/dx)[(eta+zeta)*d/dx)] U
255 AA1 = etaPlusZeta(I ,J,bi,bj)
256 & * _recip_dxF(I ,J,bi,bj)
257 & * _recip_dxC(I ,J,bi,bj) * _maskW(I,J,1,bi,bj)
258 AA2 = etaPlusZeta(I-1,J,bi,bj)
259 & * _recip_dxF(I-1,J,bi,bj)
260 & * _recip_dxC(I ,J,bi,bj) * _maskW(I,J,1,bi,bj)
261 C (d/dy)[eta*(d/dy + tanphi/a)] U (also on UVRT1/2)
262 AA3= ( etaMeanZ(I,J+1,bi,bj) * _recip_dyU(I,J+1,bi,bj)
263 & + etaMeanZ(I,J ,bi,bj) * _recip_dyU(I,J ,bi,bj)
264 & ) * _recip_dyG(I,J,bi,bj)
265 & - (etaMeanZ(I,J+1,bi,bj) - etaMeanZ(I,J,bi,bj))
266 & * 0.5 _d 0 * _tanPhiAtU(I,J,bi,bj)
267 & * recip_rSphere * _recip_dyG(I,J,bi,bj)
268 C 2*eta*(tanphi/a) * ( tanphi/a ) U
269 AA6=TWO*etaMeanU(I,J,bi,bj)*recip_rSphere*recip_rSphere
270 & * _tanPhiAtU(I,J,bi,bj) * _tanPhiAtU(I,J,bi,bj)
271 AU(I,J,bi,bj)=-AA2
272 CU(I,J,bi,bj)=-AA1
273 BU(I,J,bi,bj)=(ONE - _maskW(I,J,1,bi,bj))
274 & - AU(I,J,bi,bj) - CU(I,J,bi,bj)
275 & + ( AA3 + AA6
276 & + seaiceMassU(I,J,bi,bj)/SEAICE_deltaTdyn
277 & + 0.5*(DRAGS(I,J,bi,bj)+DRAGS(I-1,J,bi,bj))
278 & )* _maskW(I,J,1,bi,bj)
279 C coefficients of uIce(I,J-1)
280 UVRT1(I,J,bi,bj)=
281 & etaMeanZ(I,J ,bi,bj) * _recip_dyG(I,J ,bi,bj) * (
282 & _recip_dyU(I,J ,bi,bj)
283 & - _tanPhiAtU(I,J ,bi,bj) * 0.5 _d 0 * recip_rSphere )
284 & + TWO*etaMeanU(I,J,bi,bj) * _tanPhiAtU(I,J,bi,bj)
285 & * 1.0 _d 0 / ( _dyU(I,J,bi,bj) + _dyU(I,J+1,bi,bj) )
286 & *recip_rSphere
287 C coefficients of uIce(I,J+1)
288 UVRT2(I,J,bi,bj)=
289 & etaMeanZ(I,J+1,bi,bj) * _recip_dyG(I,J ,bi,bj) * (
290 & _recip_dyU(I,J+1,bi,bj)
291 & + _tanPhiAtU(I,J+1,bi,bj) * 0.5 _d 0 * recip_rSphere )
292 & - TWO*etaMeanU(I,J,bi,bj) * _tanPhiAtU(I,J,bi,bj)
293 & * 1.0 _d 0 / ( _dyU(I,J,bi,bj) + _dyU(I,J+1,bi,bj) )
294 & *recip_rSphere
295 END DO
296 END DO
297
298 DO J=1,sNy
299 AU(1,J,bi,bj)=ZERO
300 CU(sNx,J,bi,bj)=ZERO
301 CU(1,J,bi,bj)=CU(1,J,bi,bj)/BU(1,J,bi,bj)
302 END DO
303
304 C now set up right-hand side
305 DO J=1-Oly,sNy+Oly-1
306 DO I=1-Olx,sNx+Olx
307 dVdy(I,J) = ( vIceC(I,J+1,bi,bj) - vIceC(I,J,bi,bj) )
308 & * _recip_dyF(I,J,bi,bj)
309 ENDDO
310 ENDDO
311 DO J=1,sNy
312 DO I=1,sNx
313 C coriolis, "off-diagonal" drag terms and other forcing
314 FXY(I,J,bi,bj)=
315 & 0.5*( DRAGA( i ,j,bi,bj)
316 & *0.5*( vIceC( i ,j,bi,bj)+vIceC( i ,j+1,bi,bj) )
317 & + DRAGA(i-1,j,bi,bj)
318 & *0.5*( vIceC(i-1,j,bi,bj)+vIceC(i-1,j+1,bi,bj) ) )
319 & +FORCEX(I,J,bi,bj)
320 C + d/dx[ (zeta-eta) dV/dy]
321 FXY(I,J,bi,bj)=FXY(I,J,bi,bj) +
322 & ( zetaMinusEta(I ,J ,bi,bj) * dVdy(I ,J )
323 & - zetaMinusEta(I-1,J ,bi,bj) * dVdy(I-1,J )
324 & ) * _recip_dxC(I,J,bi,bj)
325 C + d/dy[ eta dV/x ]
326 FXY(I,J,bi,bj)=FXY(I,J,bi,bj) + (
327 & etaMeanZ(I,J+1,bi,bj)
328 & * ( vIceC(I ,J+1,bi,bj) - vIceC(I-1,J+1,bi,bj) )
329 & * _recip_dxV(I,J+1,bi,bj)
330 & - etaMeanZ(I,J,bi,bj)
331 & * ( vIceC(I ,J,bi,bj) - vIceC(I-1,J,bi,bj) )
332 & * _recip_dxV(I,J,bi,bj)
333 & ) * _recip_dyG(I,J,bi,bj)
334 C - d/dx[ (eta+zeta) * v * (tanphi/a) ]
335 FXY(I,J,bi,bj)=FXY(I,J,bi,bj) - (
336 & etaPlusZeta(I ,J ,bi,bj)
337 & * 0.5 _d 0 * (vIceC(I ,J,bi,bj)+vIceC(I ,J+1,bi,bj))
338 & * 0.5 _d 0 * ( _tanPhiAtU(I ,J,bi,bj)
339 & + _tanPhiAtU(I+1,J,bi,bj) )
340 & - etaPlusZeta(I-1,J,bi,bj) *
341 & * 0.5 _d 0 * (vIceC(I-1,J,bi,bj)+vIceC(I-1,J+1,bi,bj))
342 & * 0.5 _d 0 * ( _tanPhiAtU(I-1,J,bi,bj)
343 & + _tanPhiAtU(I ,J,bi,bj) )
344 & )* _recip_dxC(I,J,bi,bj)*recip_rSphere
345 C - 2*eta*(tanphi/a) * dV/dx
346 FXY(I,J,bi,bj)=FXY(I,J,bi,bj) -
347 & -TWO * etaMeanU(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj)
348 & *recip_rSphere
349 & *(vIceC(I ,J,bi,bj) + vIceC(I ,J+1,bi,bj)
350 & -vIceC(I-1,J,bi,bj) - vIceC(I-1,J+1,bi,bj))
351 & * _recip_dxC(I,J,bi,bj)
352 END DO
353 END DO
354
355 ENDDO
356 ENDDO
357
358 C NOW DO ITERATION
359 100 CONTINUE
360
361 cph--- iteration starts here
362 cph--- need to kick out goto
363 phexit = -1
364
365 C ITERATION START -----------------------------------------------------
366 #ifdef ALLOW_AUTODIFF_TAMC
367 CADJ LOOP = iteration uice
368 #endif /* ALLOW_AUTODIFF_TAMC */
369
370 DO 8000 M=1, solv_max_iters
371 cph(
372 IF ( phexit .EQ. -1 ) THEN
373 cph)
374 DO bj=myByLo(myThid),myByHi(myThid)
375 DO bi=myBxLo(myThid),myBxHi(myThid)
376 C NOW SET U(3)=U(1)
377 DO J=1,sNy
378 DO I=1,sNx
379 uIce(I,J,3,bi,bj)=uIce(I,J,1,bi,bj)
380 END DO
381 END DO
382
383 DO 1200 J=1,sNy
384 DO I=1,sNx
385 IF(I.EQ.1) THEN
386 AA2 = etaPlusZeta(I-1,J,bi,bj)
387 & * _recip_dxF(I-1,J,bi,bj)
388 & * _recip_dxC(I ,J,bi,bj)
389 AA3=AA2*uIce(I-1,J,1,bi,bj)* _maskW(I,J,1,bi,bj)
390 ELSE IF(I.EQ.sNx) THEN
391 AA1 = etaPlusZeta(I ,J,bi,bj)
392 & * _recip_dxF(I ,J,bi,bj)
393 & * _recip_dxC(I ,J,bi,bj)
394 AA3=AA1*uIce(I+1,J,1,bi,bj) * _maskW(I,J,1,bi,bj)
395 ELSE
396 AA3=ZERO
397 END IF
398 URT(I)=FXY(I,J,bi,bj)+AA3
399 & +UVRT1(I,J,bi,bj)*uIce(I,J-1,1,bi,bj)
400 & +UVRT2(I,J,bi,bj)*uIce(I,J+1,1,bi,bj)
401 URT(I)=URT(I)* _maskW(I,J,1,bi,bj) * seaiceMaskU(I,J,bi,bj)
402 END DO
403
404 DO I=1,sNx
405 CUU(I)=CU(I,J,bi,bj)
406 END DO
407 URT(1)=URT(1)/BU(1,J,bi,bj)
408 DO I=2,sNx
409 IM=I-1
410 CUU(I)=CUU(I)/(BU(I,J,bi,bj)-AU(I,J,bi,bj)*CUU(IM))
411 URT(I)=(URT(I)-AU(I,J,bi,bj)*URT(IM))
412 & /(BU(I,J,bi,bj)-AU(I,J,bi,bj)*CUU(IM))
413 END DO
414 DO I=1,sNx-1
415 J1=sNx-I
416 J2=J1+1
417 URT(J1)=URT(J1)-CUU(J1)*URT(J2)
418 END DO
419 DO I=1,sNx
420 uIce(I,J,1,bi,bj)=uIce(I,J,3,bi,bj)
421 & +WFAU*(URT(I)-uIce(I,J,3,bi,bj))
422 END DO
423
424 1200 CONTINUE
425
426 ENDDO
427 ENDDO
428
429 IF(MOD(M,SOLV_NCHECK).EQ.0) THEN
430 DO bj=myByLo(myThid),myByHi(myThid)
431 DO bi=myBxLo(myThid),myBxHi(myThid)
432 S1=ZERO
433 DO J=1,sNy
434 DO I=1,sNx
435 UERR(I,J,bi,bj)=(uIce(I,J,1,bi,bj)-uIce(I,J,3,bi,bj))
436 & * _maskW(I,J,1,bi,bj)
437 S1=MAX(ABS(UERR(I,J,bi,bj)),S1)
438 END DO
439 END DO
440 _GLOBAL_MAX_R8( S1, myThid )
441 ENDDO
442 ENDDO
443 C SAFEGUARD AGAINST BAD FORCING ETC
444 IF(M.GT.1.AND.S1.GT.S1A) WFAU=WFAU2
445 S1A=S1
446 IF(S1.LT.LSR_ERROR) THEN
447 ICOUNT1=M
448 cph(
449 cph GO TO 8001
450 phexit = 1
451 cph)
452 END IF
453 END IF
454 CALL SEAICE_EXCH_UV ( uIce, vIce, myThid )
455
456 cph(
457 END IF
458 cph)
459
460 8000 CONTINUE
461 cph 8001 CONTINUE
462 C ITERATION END -----------------------------------------------------
463
464 IF ( debugLevel .GE. debLevB ) THEN
465 _BEGIN_MASTER( myThid )
466 write(*,'(A,I6,1P2E22.14)')' U lsr iters, error = ',ICOUNT1,S1
467 _END_MASTER( myThid )
468 ENDIF
469
470 C NOW FOR vIce
471 DO bj=myByLo(myThid),myByHi(myThid)
472 DO bi=myBxLo(myThid),myBxHi(myThid)
473
474 DO J=1,sNy
475 DO I=1,sNx
476 C coefficients for VICE(I,J)
477 C d/dy [(eta+zeta) d/dy] V
478 AA1= etaPlusZeta(I,J ,bi,bj)
479 & *_recip_dyF(I,J ,bi,bj) * _recip_dyC(I,J,bi,bj)
480 AA2= etaPlusZeta(I,J-1,bi,bj)
481 & * _recip_dyF(I,J-1,bi,bj) * _recip_dyC(I,J,bi,bj)
482 C d/dx [eta d/dx] V
483 AA3= etaMeanZ(I+1,J,bi,bj)
484 & * _recip_dxG(I,J,bi,bj) * _recip_dxV(I+1,J,bi,bj)
485 AA4= etaMeanZ(I ,J,bi,bj)
486 & *_recip_dxG(I,J,bi,bj) * _recip_dxV(I ,J,bi,bj)
487 C d/dy [(zeta-eta) tanphi/a] V
488 AA5= zetaMinusEta(I,J ,bi,bj) * tanPhiAtU(I,J ,bi,bj)
489 & * _recip_dyC(I,J,bi,bj)*recip_rSphere * 0.5 _d 0
490 AA6= zetaMinusEta(I,J-1,bi,bj) * tanPhiAtU(I,J-1,bi,bj)
491 & * _recip_dyC(I,J,bi,bj)*recip_rSphere * 0.5 _d 0
492 C 2*eta tanphi/a (tanphi/a - d/dy) V
493 AA7=TWO*etaMeanV(I,J,bi,bj) * recip_rSphere
494 & * _tanPhiAtV(I,J,bi,bj)
495 C* _tanPhiAtV(I,J,bi,bj)*recip_rSphere
496
497 AV(I,J,bi,bj)=(
498 & - AA2
499 & - AA6
500 & - AA7*1.0 _d 0 / ( _dyF(I,J,bi,bj) + _dyF(I,J-1,bi,bj) )
501 & )* _maskS(I,J,1,bi,bj)
502 CV(I,J,bi,bj)=(
503 & -AA1
504 & + AA5
505 & + AA7*1.0 _d 0 / ( _dyF(I,J,bi,bj) + _dyF(I,J-1,bi,bj) )
506 & )* _maskS(I,J,1,bi,bj)
507 BV(I,J,bi,bj)= (ONE- _maskS(I,J,1,bi,bj))
508 & +( (AA1+AA2) + (AA3+AA4) + (AA5-AA6)
509 & + AA7 * _tanPhiAtV(I,J,bi,bj)*recip_rSphere
510 & + seaiceMassV(I,J,bi,bj)/SEAICE_deltaTdyn
511 & + 0.5 _d 0 * ( DRAGS(I,J,bi,bj) + DRAGS(I,J-1,bi,bj) )
512 & )* _maskS(I,J,1,bi,bj)
513 C coefficients for V(I-1,J)
514 UVRT1(I,J,bi,bj)= AA4
515 C coefficients for V(I+1,J)
516 UVRT2(I,J,bi,bj)= AA3
517 END DO
518 END DO
519
520 DO I=1,sNx
521 AV(I,1,bi,bj)=ZERO
522 CV(I,sNy,bi,bj)=ZERO
523 CV(I,1,bi,bj)=CV(I,1,bi,bj)/BV(I,1,bi,bj)
524 END DO
525
526 C now set up right-hand-side
527 DO J=1-Oly,sNy+Oly-1
528 DO I=1-Olx,sNx+Olx-1
529 dUdx(I,J) = ( uIceC(I+1,J,bi,bj) - uIceC(I,J,bi,bj) )
530 & * _recip_dxF(I,J,bi,bj)
531 dUdy(I,J) = ( uIceC(I,J+1,bi,bj) - uIceC(I,J,bi,bj) )
532 & * _recip_dyU(I,J+1,bi,bj)
533 ENDDO
534 ENDDO
535 DO J=1,sNy
536 DO I=1,sNx
537 C coriols, "off-diagonal" drag terms and other foring
538 FXY(I,J,bi,bj)=
539 & -0.5*( _fCori(i, j ,bi,bj)
540 & *0.5*( uIceC(i ,j ,bi,bj)+uIceC(i+1, j,bi,bj) )
541 & + _fCori(i,j-1,bi,bj)
542 & *0.5*( uIceC(i ,j-1,bi,bj)+uIceC(i+1,j-1,bi,bj) ) )
543 & + FORCEY(I,J,bi,bj)
544 C + d/dy[ (zeta-eta) dU/dx ]
545 FXY(I,J,bi,bj)=FXY(I,J,bi,bj) +
546 & ( zetaMinusEta(I,J ,bi,bj)*dUdx(I,J )
547 & - zetaMinusEta(I,J-1,bi,bj)*dUdx(I,J-1) )
548 & * _recip_dyC(I,J,bi,bj)
549 C + d/dx[ eta dU/dy ]
550 FXY(I,J,bi,bj)=FXY(I,J,bi,bj) +
551 & ( etaMeanZ(I+1,J ,bi,bj) * dUdy(I+1,J)
552 & - etaMeanZ(I ,J ,bi,bj) * dUdy(I ,J))
553 & * _recip_dxG(I,J,bi,bj)
554 C + d/dx[ eta * (tanphi/a) * U ]
555 FXY(I,J,bi,bj)=FXY(I,J,bi,bj) + (
556 & etaMeanZ(I+1,J,bi,bj) * 0.5 *
557 & ( uIceC(I+1,J ,bi,bj) * _tanPhiAtU(I+1,J ,bi,bj)
558 & + uIceC(I+1,J-1,bi,bj) * _tanPhiAtU(I+1,J-1,bi,bj) )
559 & - etaMeanZ(I ,J,bi,bj) * 0.5 *
560 & ( uIceC(I ,J ,bi,bj) * _tanPhiAtU(I ,J ,bi,bj)
561 & + uIceC(I ,J-1,bi,bj) * _tanPhiAtU(I ,J ,bi,bj) )
562 & ) * _recip_dxG(I,J,bi,bj)*recip_rSphere
563 C + 2*eta*(tanphi/a) dU/dx
564 FXY(I,J,bi,bj)=FXY(I,J,bi,bj) +
565 & TWO * etaMeanV(I,J,bi,bj)*TWO * _tanPhiAtV(I,J,bi,bj)
566 & * ( uIceC(I+1,J,bi,bj)+uIceC(I+1,J-1,bi,bj)
567 & - uIceC(I ,J,bi,bj)-uIceC(I ,J-1,bi,bj) )
568 & * _recip_dxG(I,J,bi,bj)
569 & *recip_rSphere
570 END DO
571 END DO
572
573 ENDDO
574 ENDDO
575
576 C NOW DO ITERATION
577 300 CONTINUE
578
579 cph--- iteration starts here
580 cph--- need to kick out goto
581 phexit = -1
582
583 C ITERATION START -----------------------------------------------------
584 #ifdef ALLOW_AUTODIFF_TAMC
585 CADJ LOOP = iteration vice
586 #endif /* ALLOW_AUTODIFF_TAMC */
587
588 DO 9000 M=1, solv_max_iters
589 cph(
590 IF ( phexit .EQ. -1 ) THEN
591 cph)
592 C NOW SET U(3)=U(1)
593 DO bj=myByLo(myThid),myByHi(myThid)
594 DO bi=myBxLo(myThid),myBxHi(myThid)
595
596 DO J=1,sNy
597 DO I=1,sNx
598 vIce(I,J,3,bi,bj)=vIce(I,J,1,bi,bj)
599 END DO
600 END DO
601
602 DO I=1,sNx
603 DO J=1,sNy
604 IF(J.EQ.1) THEN
605 AA2= etaPlusZeta(I,J-1,bi,bj)
606 & * _recip_dyF(I,J-1,bi,bj) * _recip_dyC(I,J,bi,bj)
607 AA3=( AA2
608 & + zetaMinusEta(I,J-1,bi,bj) * tanPhiAtU(I,J-1,bi,bj)
609 & * _recip_dyC(I,J,bi,bj)*recip_rSphere
610 & + TWO*etaMeanV(I,J,bi,bj) * recip_rSphere
611 & * _tanPhiAtV(I,J,bi,bj)
612 & *1.0 _d 0 / ( _dyF(I,J,bi,bj) + _dyF(I,J-1,bi,bj) )
613 & ) * vIce(I,J-1,1,bi,bj) * _maskS(I,J,1,bi,bj)
614 ELSE IF(J.EQ.sNy) THEN
615 AA1= etaPlusZeta(I,J ,bi,bj)
616 & *_recip_dyF(I,J ,bi,bj) * _recip_dyC(I,J,bi,bj)
617 AA3=( AA1
618 & - zetaMinusEta(I,J ,bi,bj) * tanPhiAtU(I,J ,bi,bj)
619 & * _recip_dyC(I,J,bi,bj)*recip_rSphere
620 & - TWO*etaMeanV(I,J,bi,bj) * recip_rSphere
621 & * _tanPhiAtV(I,J,bi,bj)
622 & *1.0 _d 0 / ( _dyF(I,J,bi,bj) + _dyF(I,J-1,bi,bj) )
623 & ) * vIce(I,J+1,1,bi,bj) * _maskS(I,J,1,bi,bj)
624 ELSE
625 AA3=ZERO
626 END IF
627
628 VRT(J)=FXY(I,J,bi,bj)+AA3+UVRT1(I,J,bi,bj)*vIce(I-1,J,1,bi,bj)
629 & +UVRT2(I,J,bi,bj)*vIce(I+1,J,1,bi,bj)
630 VRT(J)=VRT(J)* _maskS(I,J,1,bi,bj) * seaiceMaskV(I,J,bi,bj)
631 END DO
632
633 DO J=1,sNy
634 CVV(J)=CV(I,J,bi,bj)
635 END DO
636 VRT(1)=VRT(1)/BV(I,1,bi,bj)
637 DO J=2,sNy
638 JM=J-1
639 CVV(J)=CVV(J)/(BV(I,J,bi,bj)-AV(I,J,bi,bj)*CVV(JM))
640 VRT(J)=(VRT(J)-AV(I,J,bi,bj)*VRT(JM))
641 & /(BV(I,J,bi,bj)-AV(I,J,bi,bj)*CVV(JM))
642 END DO
643 DO J=1,sNy-1
644 J1=sNy-J
645 J2=J1+1
646 VRT(J1)=VRT(J1)-CVV(J1)*VRT(J2)
647 END DO
648 DO J=1,sNy
649 vIce(I,J,1,bi,bj)=vIce(I,J,3,bi,bj)
650 & +WFAV*(VRT(J)-vIce(I,J,3,bi,bj))
651 END DO
652 ENDDO
653
654 ENDDO
655 ENDDO
656
657 IF(MOD(M,SOLV_NCHECK).EQ.0) THEN
658 DO bj=myByLo(myThid),myByHi(myThid)
659 DO bi=myBxLo(myThid),myBxHi(myThid)
660 S2=ZERO
661 DO J=1,sNy
662 DO I=1,sNx
663 UERR(I,J,bi,bj)=(vIce(I,J,1,bi,bj)-vIce(I,J,3,bi,bj))
664 & * _maskS(I,J,1,bi,bj)
665 S2=MAX(ABS(UERR(I,J,bi,bj)),S2)
666 END DO
667 END DO
668 _GLOBAL_MAX_R8( S2, myThid )
669 ENDDO
670 ENDDO
671 C SAFEGUARD AGAINST BAD FORCING ETC
672 IF(M.GT.1.AND.S2.GT.S2A) WFAV=WFAV2
673 S2A=S2
674 IF(S2.LT.LSR_ERROR) THEN
675 ICOUNT2=M
676 cph(
677 cph GO TO 9001
678 phexit = 1
679 cph)
680 END IF
681 END IF
682
683 CALL SEAICE_EXCH_UV ( uIce, vIce, myThid )
684
685 cph(
686 END IF
687 cph)
688
689 9000 CONTINUE
690 cph 9001 CONTINUE
691 C ITERATION END -----------------------------------------------------
692
693 IF ( debugLevel .GE. debLevB ) THEN
694 _BEGIN_MASTER( myThid )
695 write(*,'(A,I6,1P2E22.14)')' V lsr iters, error = ',ICOUNT2,S2
696 _END_MASTER( myThid )
697 ENDIF
698
699 C NOW END
700 C NOW MAKE COROLIS TERM IMPLICIT
701 DO bj=myByLo(myThid),myByHi(myThid)
702 DO bi=myBxLo(myThid),myBxHi(myThid)
703 DO J=1,sNy
704 DO I=1,sNx
705 uIce(I,J,1,bi,bj)=uIce(I,J,1,bi,bj)* _maskW(I,J,1,bi,bj)
706 vIce(I,J,1,bi,bj)=vIce(I,J,1,bi,bj)* _maskS(I,J,1,bi,bj)
707 END DO
708 END DO
709 ENDDO
710 ENDDO
711 CALL SEAICE_EXCH_UV ( uIce, vIce, myThid )
712
713 #endif /* SEAICE_ALLOW_DYNAMICS */
714 #endif /* SEAICE_CGRID */
715
716 RETURN
717 END

  ViewVC Help
Powered by ViewVC 1.1.22