1 |
C $Header: /u/gcmpack/MITgcm/pkg/seaice/lsr.F,v 1.18 2006/03/06 13:17:37 mlosch Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "SEAICE_OPTIONS.h" |
5 |
|
6 |
CStartOfInterface |
7 |
SUBROUTINE lsr( ilcall, myThid ) |
8 |
C /==========================================================\ |
9 |
C | SUBROUTINE 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 \==========================================================/ |
17 |
IMPLICIT NONE |
18 |
|
19 |
C === Global variables === |
20 |
#include "SIZE.h" |
21 |
#include "EEPARAMS.h" |
22 |
#include "PARAMS.h" |
23 |
#include "GRID.h" |
24 |
#include "SEAICE.h" |
25 |
#include "SEAICE_PARAMS.h" |
26 |
C#include "SEAICE_GRID.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 |
#ifndef 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, 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 ETAMEAN (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy) |
66 |
_RL ZETAMEAN (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy) |
67 |
|
68 |
_RL UVRT1 (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy) |
69 |
_RL UVRT2 (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy) |
70 |
|
71 |
_RL vDiffX (1-Olx:sNx+Olx,1-Oly:sNy+Oly) |
72 |
_RL vDiffY (1-Olx:sNx+Olx,1-Oly:sNy+Oly) |
73 |
_RL uDiffY (1-Olx:sNx+Olx,1-Oly:sNy+Oly) |
74 |
|
75 |
C SET SOME VALUES |
76 |
WFAU1=0.95 _d 0 |
77 |
WFAV1=0.95 _d 0 |
78 |
WFAU2=ZERO |
79 |
WFAV2=ZERO |
80 |
|
81 |
S1A=0.80 _d 0 |
82 |
S2A=0.80 _d 0 |
83 |
WFAU=WFAU1 |
84 |
WFAV=WFAV1 |
85 |
|
86 |
SOLV_MAX_ITERS=1500 |
87 |
SOLV_NCHECK=2 |
88 |
|
89 |
ICOUNT1=SOLV_MAX_ITERS |
90 |
ICOUNT2=SOLV_MAX_ITERS |
91 |
|
92 |
C SOLVE FOR UICE |
93 |
|
94 |
#ifdef ALLOW_AUTODIFF_TAMC |
95 |
cph That's an important one! Note, that |
96 |
cph * lsr is called twice, thus the icall index |
97 |
cph * this storing is still outside the iteration loop |
98 |
CADJ STORE uice = comlev1_lsr, |
99 |
CADJ & key = ikey_dynamics + (ilcall-1)*nchklev_1 |
100 |
CADJ STORE vice = comlev1_lsr, |
101 |
CADJ & key = ikey_dynamics + (ilcall-1)*nchklev_1 |
102 |
#endif /* ALLOW_AUTODIFF_TAMC */ |
103 |
|
104 |
DO bj=myByLo(myThid),myByHi(myThid) |
105 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
106 |
DO j=1-Oly,sNy+Oly |
107 |
DO i=1-Olx,sNx+Olx |
108 |
FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj) |
109 |
& +AMASS(I,J,bi,bj)/SEAICE_deltaTdyn*UICE(I,J,2,bi,bj) |
110 |
FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj) |
111 |
& +AMASS(I,J,bi,bj)/SEAICE_deltaTdyn*VICE(I,J,2,bi,bj) |
112 |
FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)*UVM(I,J,bi,bj) |
113 |
FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)*UVM(I,J,bi,bj) |
114 |
etaPlusZeta(I,J,bi,bj) = ETA(I,J,bi,bj)+ZETA(I,J,bi,bj) |
115 |
zetaMinusEta(I,J,bi,bj) = ZETA(I,J,bi,bj)-ETA(I,J,bi,bj) |
116 |
ENDDO |
117 |
ENDDO |
118 |
DO j=1-Oly+1,sNy+Oly |
119 |
DO i=1-Olx+1,sNx+Olx |
120 |
ETAMEAN(I,J,bi,bj) =QUART*( |
121 |
& ETA(I,J-1,bi,bj) + ETA(I-1,J-1,bi,bj) |
122 |
& +ETA(I,J ,bi,bj) + ETA(I-1,J ,bi,bj)) |
123 |
ZETAMEAN(I,J,bi,bj)=QUART*( |
124 |
& ZETA(I,J-1,bi,bj) + ZETA(I-1,J-1,bi,bj) |
125 |
& +ZETA(I,J ,bi,bj) + ZETA(I-1,J ,bi,bj)) |
126 |
ENDDO |
127 |
ENDDO |
128 |
ENDDO |
129 |
ENDDO |
130 |
|
131 |
DO bj=myByLo(myThid),myByHi(myThid) |
132 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
133 |
|
134 |
DO J=1,sNy |
135 |
DO I=1,sNx |
136 |
AA1=( etaPlusZeta(I ,J-1,bi,bj) * _recip_dxF(I ,J-1,bi,bj) |
137 |
& +etaPlusZeta(I ,J ,bi,bj) * _recip_dxF(I ,J ,bi,bj) |
138 |
& )*0.5 _d 0 * _recip_dxV(I,J,bi,bj) * UVM(I,J,bi,bj) |
139 |
AA2=( etaPlusZeta(I-1,J-1,bi,bj) * _recip_dxF(I-1,J-1,bi,bj) |
140 |
& +etaPlusZeta(I-1,J ,bi,bj) * _recip_dxF(I-1,J ,bi,bj) |
141 |
& )*0.5 _d 0 * _recip_dxV(I,J,bi,bj) * UVM(I,J,bi,bj) |
142 |
AA3= 0.5 _d 0 *(ETA(I-1,J ,bi,bj)+ETA(I,J ,bi,bj)) |
143 |
AA4= 0.5 _d 0 *(ETA(I-1,J-1,bi,bj)+ETA(I,J-1,bi,bj)) |
144 |
AA5= -(AA3-AA4) * _tanPhiAtV(I,J,bi,bj) |
145 |
& * _recip_dyU(I,J,bi,bj)*recip_rSphere |
146 |
AA6=TWO*ETAMEAN(I,J,bi,bj) *recip_rSphere*recip_rSphere |
147 |
& * _tanPhiAtV(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj) |
148 |
AU(I,J,bi,bj)=-AA2 |
149 |
CU(I,J,bi,bj)=-AA1 |
150 |
BU(I,J,bi,bj)=(ONE-UVM(I,J,bi,bj)) |
151 |
& - AU(I,J,bi,bj) - CU(I,J,bi,bj) |
152 |
& + ((AA3+AA4)*_recip_dyU(I,J,bi,bj)*_recip_dyU(I,J,bi,bj) |
153 |
& + AA5 + AA6 |
154 |
& + AMASS(I,J,bi,bj)/SEAICE_deltaTdyn |
155 |
& + DRAGS(I,J,bi,bj) |
156 |
& )*UVM(I,J,bi,bj) |
157 |
END DO |
158 |
END DO |
159 |
|
160 |
DO J=1,sNy |
161 |
AU(1,J,bi,bj)=ZERO |
162 |
CU(sNx,J,bi,bj)=ZERO |
163 |
CU(1,J,bi,bj)=CU(1,J,bi,bj)/BU(1,J,bi,bj) |
164 |
END DO |
165 |
|
166 |
C now set up right-hand side |
167 |
DO J=1-Oly,sNy+Oly-1 |
168 |
DO I=1-Olx,sNx+Olx-1 |
169 |
vDiffY(I,J) = 0.5 _d 0 * ( |
170 |
& VICEC(I+1,J+1,bi,bj) + VICEC(I ,J+1,bi,bj) |
171 |
& -VICEC(I+1,J ,bi,bj) - VICEC(I ,J ,bi,bj) ) |
172 |
vDiffX(I,J) = 0.5 _d 0 * ( |
173 |
& VICEC(I+1,J+1,bi,bj) + VICEC(I+1,J ,bi,bj) |
174 |
& -VICEC(I ,J+1,bi,bj) - VICEC(I ,J ,bi,bj) ) |
175 |
ENDDO |
176 |
ENDDO |
177 |
DO J=1,sNy |
178 |
DO I=1,sNx |
179 |
FXY(I,J,bi,bj)=DRAGA(I,J,bi,bj)*VICEC(I,J,bi,bj) |
180 |
& +FORCEX(I,J,bi,bj) |
181 |
& + ( zetaMinusEta(I ,J ,bi,bj) * vDiffY(I ,J ) |
182 |
& + zetaMinusEta(I ,J-1,bi,bj) * vDiffY(I ,J-1) |
183 |
& - zetaMinusEta(I-1,J ,bi,bj) * vDiffY(I-1,J ) |
184 |
& - zetaMinusEta(I-1,J-1,bi,bj) * vDiffY(I-1,J-1) |
185 |
& )* 0.5 _d 0 * |
186 |
& _recip_dxV(I,J,bi,bj) * _recip_dyU(I,J,bi,bj) |
187 |
& |
188 |
& + ( ETA (I ,J ,bi,bj) * vDiffX(I ,J ) |
189 |
& * _recip_dxF(I ,J,bi,bj) |
190 |
& + ETA (I-1,J ,bi,bj) * vDiffX(I-1,J ) |
191 |
& * _recip_dxF(I-1,J,bi,bj) |
192 |
& - ETA (I ,J-1,bi,bj) * vDiffX(I ,J-1) |
193 |
& * _recip_dxF(I ,J-1,bi,bj) |
194 |
& - ETA (I-1,J-1,bi,bj) * vDiffX(I-1,J-1) |
195 |
& * _recip_dxF(I-1,J-1,bi,bj) |
196 |
& ) * 0.5 _d 0 * _recip_dyU(I,J,bi,bj) |
197 |
& |
198 |
& -(etaPlusZeta(I ,J ,bi,bj)+etaPlusZeta(I ,J-1,bi,bj) |
199 |
& -etaPlusZeta(I-1,J-1,bi,bj)-etaPlusZeta(I-1,J ,bi,bj)) |
200 |
& * VICEC(I,J,bi,bj) |
201 |
& * _tanPhiAtV(I,J,bi,bj) |
202 |
& * 0.5 _d 0 * _recip_dxV(I,J,bi,bj)*recip_rSphere |
203 |
& |
204 |
& -(ETAMEAN(I,J,bi,bj)+ZETAMEAN(I,J,bi,bj)) |
205 |
& *(VICEC(I+1,J,bi,bj) - VICEC(I-1,J,bi,bj)) |
206 |
& * _tanPhiAtV(I,J,bi,bj) |
207 |
& * 1.0 _d 0 /( _dxG(I,J,bi,bj) + _dxG(I-1,J,bi,bj) ) |
208 |
& *recip_rSphere |
209 |
& |
210 |
& -ETAMEAN(I,J,bi,bj) |
211 |
& *(VICEC(I+1,J,bi,bj) - VICEC(I-1,J,bi,bj)) |
212 |
& *TWO* _tanPhiAtV(I,J,bi,bj) |
213 |
& * 1.0 _d 0 /( _dxG(I,J,bi,bj) + _dxG(I-1,J,bi,bj) ) |
214 |
& *recip_rSphere |
215 |
|
216 |
UVRT1(I,J,bi,bj)= |
217 |
& 0.5 _d 0 * (ETA(I-1,J-1,bi,bj)+ETA(I,J-1,bi,bj)) |
218 |
& * _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj) |
219 |
& - ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J-1,bi,bj) |
220 |
& * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere |
221 |
& + TWO*ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj) |
222 |
& * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere |
223 |
UVRT2(I,J,bi,bj)= |
224 |
& 0.5 _d 0 * (ETA(I-1,J,bi,bj)+ETA(I,J,bi,bj)) |
225 |
& * _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj) |
226 |
& + ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J+1,bi,bj) |
227 |
& * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere |
228 |
& - TWO*ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj) |
229 |
& * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere |
230 |
END DO |
231 |
END DO |
232 |
|
233 |
ENDDO |
234 |
ENDDO |
235 |
|
236 |
C NOW DO ITERATION |
237 |
100 CONTINUE |
238 |
|
239 |
cph--- iteration starts here |
240 |
cph--- need to kick out goto |
241 |
phexit = -1 |
242 |
|
243 |
C ITERATION START ----------------------------------------------------- |
244 |
#ifdef ALLOW_AUTODIFF_TAMC |
245 |
CADJ LOOP = iteration uice |
246 |
#endif /* ALLOW_AUTODIFF_TAMC */ |
247 |
|
248 |
DO 8000 M=1, solv_max_iters |
249 |
cph( |
250 |
IF ( phexit .EQ. -1 ) THEN |
251 |
cph) |
252 |
DO bj=myByLo(myThid),myByHi(myThid) |
253 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
254 |
C NOW SET U(3)=U(1) |
255 |
DO J=1,sNy |
256 |
DO I=1,sNx |
257 |
UICE(I,J,3,bi,bj)=UICE(I,J,1,bi,bj) |
258 |
END DO |
259 |
END DO |
260 |
|
261 |
DO 1200 J=1,sNy |
262 |
DO I=1,sNx |
263 |
IF(I.EQ.1) THEN |
264 |
AA2=(etaPlusZeta(I-1,J-1,bi,bj) * _recip_dxF(I-1,J-1,bi,bj) |
265 |
& +etaPlusZeta(I-1,J ,bi,bj) * _recip_dxF(I-1,J ,bi,bj) |
266 |
& )*0.5 _d 0 * _recip_dxV(I,J,bi,bj) |
267 |
AA3=AA2*UICE(I-1,J,1,bi,bj)*UVM(I,J,bi,bj) |
268 |
ELSE IF(I.EQ.sNx) THEN |
269 |
AA1=(etaPlusZeta(I ,J-1,bi,bj) * _recip_dxF(I ,J-1,bi,bj) |
270 |
& +etaPlusZeta(I ,J ,bi,bj) * _recip_dxF(I ,J ,bi,bj) |
271 |
& )*0.5 _d 0 * _recip_dxV(I,J,bi,bj) |
272 |
AA3=AA1*UICE(I+1,J,1,bi,bj)*UVM(I,J,bi,bj) |
273 |
ELSE |
274 |
AA3=ZERO |
275 |
END IF |
276 |
URT(I)=FXY(I,J,bi,bj)+AA3 |
277 |
& +UVRT1(I,J,bi,bj)*UICE(I,J-1,1,bi,bj) |
278 |
& +UVRT2(I,J,bi,bj)*UICE(I,J+1,1,bi,bj) |
279 |
URT(I)=URT(I)*UVM(I,J,bi,bj) |
280 |
END DO |
281 |
|
282 |
DO I=1,sNx |
283 |
CUU(I)=CU(I,J,bi,bj) |
284 |
END DO |
285 |
URT(1)=URT(1)/BU(1,J,bi,bj) |
286 |
DO I=2,sNx |
287 |
IM=I-1 |
288 |
CUU(I)=CUU(I)/(BU(I,J,bi,bj)-AU(I,J,bi,bj)*CUU(IM)) |
289 |
URT(I)=(URT(I)-AU(I,J,bi,bj)*URT(IM)) |
290 |
& /(BU(I,J,bi,bj)-AU(I,J,bi,bj)*CUU(IM)) |
291 |
END DO |
292 |
DO I=1,sNx-1 |
293 |
J1=sNx-I |
294 |
J2=J1+1 |
295 |
URT(J1)=URT(J1)-CUU(J1)*URT(J2) |
296 |
END DO |
297 |
DO I=1,sNx |
298 |
UICE(I,J,1,bi,bj)=UICE(I,J,3,bi,bj) |
299 |
& +WFAU*(URT(I)-UICE(I,J,3,bi,bj)) |
300 |
END DO |
301 |
|
302 |
1200 CONTINUE |
303 |
|
304 |
ENDDO |
305 |
ENDDO |
306 |
|
307 |
IF(MOD(M,SOLV_NCHECK).EQ.0) THEN |
308 |
DO bj=myByLo(myThid),myByHi(myThid) |
309 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
310 |
S1=ZERO |
311 |
DO J=1,sNy |
312 |
DO I=1,sNx |
313 |
UERR(I,J,bi,bj)=(UICE(I,J,1,bi,bj)-UICE(I,J,3,bi,bj)) |
314 |
& *UVM(I,J,bi,bj) |
315 |
S1=MAX(ABS(UERR(I,J,bi,bj)),S1) |
316 |
END DO |
317 |
END DO |
318 |
_GLOBAL_MAX_R8( S1, myThid ) |
319 |
ENDDO |
320 |
ENDDO |
321 |
C SAFEGUARD AGAINST BAD FORCING ETC |
322 |
IF(M.GT.1.AND.S1.GT.S1A) WFAU=WFAU2 |
323 |
S1A=S1 |
324 |
IF(S1.LT.LSR_ERROR) THEN |
325 |
ICOUNT1=M |
326 |
cph( |
327 |
cph GO TO 8001 |
328 |
phexit = 1 |
329 |
cph) |
330 |
END IF |
331 |
END IF |
332 |
CALL EXCH_3D_RL( UICE, 3, myThid ) |
333 |
|
334 |
cph( |
335 |
END IF |
336 |
cph) |
337 |
|
338 |
8000 CONTINUE |
339 |
cph 8001 CONTINUE |
340 |
C ITERATION END ----------------------------------------------------- |
341 |
|
342 |
IF ( debugLevel .GE. debLevB ) THEN |
343 |
_BEGIN_MASTER( myThid ) |
344 |
write(*,'(A,I6,1P2E22.14)')' U lsr iters, error = ',ICOUNT1,S1 |
345 |
_END_MASTER( myThid ) |
346 |
ENDIF |
347 |
|
348 |
C NOW FOR VICE |
349 |
DO bj=myByLo(myThid),myByHi(myThid) |
350 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
351 |
|
352 |
DO J=1,sNy |
353 |
DO I=1,sNx |
354 |
AA1=0.5 _d 0 * _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj) |
355 |
& * (etaPlusZeta(I-1,J ,bi,bj) + etaPlusZeta(I,J ,bi,bj)) |
356 |
AA2=0.5 _d 0 * _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj) |
357 |
& * (etaPlusZeta(I-1,J-1,bi,bj) + etaPlusZeta(I,J-1,bi,bj)) |
358 |
AA3= (ETA(I ,J-1,bi,bj) * _recip_dxV(I,J,bi,bj) |
359 |
& +ETA(I ,J ,bi,bj) * _recip_dxV(I,J,bi,bj) |
360 |
& )* 0.5 _d 0 * _recip_dxV(I,J,bi,bj) |
361 |
AA4= (ETA(I-1,J-1,bi,bj)+ETA(I-1,J,bi,bj))*0.5 _d 0 |
362 |
& *_recip_dxV(I,J,bi,bj) * _recip_dxV(I,J,bi,bj) |
363 |
AA5=(zetaMinusEta(I-1,J ,bi,bj) + zetaMinusEta(I,J ,bi,bj) |
364 |
& -zetaMinusEta(I-1,J-1,bi,bj) - zetaMinusEta(I,J-1,bi,bj) |
365 |
& )* _tanPhiAtV(I,J,bi,bj) |
366 |
& * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere |
367 |
|
368 |
AA6=TWO*ETAMEAN(I,J,bi,bj) * recip_rSphere*recip_rSphere |
369 |
& * _tanPhiAtV(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj) |
370 |
|
371 |
AV(I,J,bi,bj)=( |
372 |
& - AA2 |
373 |
& - (ZETAMEAN(I,J,bi,bj)-ETAMEAN(I,J,bi,bj)) |
374 |
& * _tanPhiAtV(I,J-1,bi,bj) |
375 |
& * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere |
376 |
& -ETAMEAN(I,J,bi,bj)*TWO* _tanPhiAtV(I,J,bi,bj) |
377 |
& * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere |
378 |
& )*UVM(I,J,bi,bj) |
379 |
CV(I,J,bi,bj)=( |
380 |
& -AA1 |
381 |
& +(ZETAMEAN(I,J,bi,bj)-ETAMEAN(I,J,bi,bj)) |
382 |
& * _tanPhiAtV(I,J+1,bi,bj) |
383 |
& * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere |
384 |
& +ETAMEAN(I,J,bi,bj)*TWO* _tanPhiAtV(I,J,bi,bj) |
385 |
& * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere |
386 |
& )*UVM(I,J,bi,bj) |
387 |
BV(I,J,bi,bj)= (ONE-UVM(I,J,bi,bj)) |
388 |
& +( (AA1+AA2) + (AA3+AA4) + AA5 + AA6 |
389 |
& +AMASS(I,J,bi,bj)/SEAICE_deltaTdyn+DRAGS(I,J,bi,bj)) |
390 |
& *UVM(I,J,bi,bj) |
391 |
END DO |
392 |
END DO |
393 |
|
394 |
DO I=1,sNx |
395 |
AV(I,1,bi,bj)=ZERO |
396 |
CV(I,sNy,bi,bj)=ZERO |
397 |
CV(I,1,bi,bj)=CV(I,1,bi,bj)/BV(I,1,bi,bj) |
398 |
END DO |
399 |
|
400 |
C now set up right-hand-side |
401 |
DO J=1-Oly,sNy+Oly-1 |
402 |
DO I=1-Olx,sNx+Olx-1 |
403 |
uDiffY(I,J) = 0.5 _d 0 * ( |
404 |
& UICEC(I+1,J+1,bi,bj) + UICEC(I ,J+1,bi,bj) |
405 |
5 -UICEC(I+1,J ,bi,bj) - UICEC(I ,J ,bi,bj) ) |
406 |
ENDDO |
407 |
ENDDO |
408 |
DO J=1,sNy |
409 |
DO I=1,sNx |
410 |
FXY(I,J,bi,bj)=-DRAGA(I,J,bi,bj)*UICEC(I,J,bi,bj) |
411 |
& +FORCEY(I,J,bi,bj) |
412 |
& |
413 |
& + HALF*(UICEC(I+1,J,bi,bj)-UICEC(I-1,J,bi,bj)) |
414 |
& *(zetaMinusEta(I-1,J ,bi,bj)+zetaMinusEta(I,J ,bi,bj) |
415 |
& -zetaMinusEta(I-1,J-1,bi,bj)-zetaMinusEta(I,J-1,bi,bj)) |
416 |
& * 1. _d 0 /( _dxG(I,J,bi,bj) + _dxG(I-1,J,bi,bj)) |
417 |
& * _recip_dyU(I,J,bi,bj) |
418 |
& |
419 |
& +HALF*(ZETAMEAN(I,J,bi,bj) - ETAMEAN(I,J,bi,bj)) |
420 |
& *((UICEC(I+1,J+1,bi,bj) - UICEC(I-1,J+1,bi,bj)) |
421 |
& * 2. _d 0 /( _dxG(I,J+1,bi,bj) + _dxG(I-1,J+1,bi,bj)) |
422 |
& -(UICEC(I+1,J-1,bi,bj)-UICEC(I-1,J-1,bi,bj)) |
423 |
& * 2. _d 0 /( _dxG(I,J-1,bi,bj) + _dxG(I-1,J-1,bi,bj))) |
424 |
& * 0.5 _d 0 * _recip_dyU(I,J,bi,bj) |
425 |
& |
426 |
& +(ETA(I ,J ,bi,bj) * uDiffY(I ,J ) |
427 |
& +ETA(I ,J-1,bi,bj) * uDiffY(I ,J-1) |
428 |
& -ETA(I-1,J ,bi,bj) * uDiffY(I-1,J ) |
429 |
& -ETA(I-1,J-1,bi,bj) * uDiffY(I-1,J-1) |
430 |
& )*0.5 _d 0* _recip_dxV(I,J,bi,bj)* _recip_dyU(I,J,bi,bj) |
431 |
& |
432 |
& +(ETA(I ,J ,bi,bj) + ETA(I ,J-1,bi,bj) |
433 |
& -ETA(I-1,J-1,bi,bj) - ETA(I-1,J ,bi,bj)) |
434 |
& * UICEC(I,J,bi,bj) |
435 |
& * _tanPhiAtV(I,J,bi,bj) |
436 |
& * 0.5 _d 0 * _recip_dxV(I,J,bi,bj)*recip_rSphere |
437 |
& +ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj) |
438 |
& *(UICEC(I+1,J,bi,bj)-UICEC(I-1,J,bi,bj)) |
439 |
& * 0.5 _d 0 * _recip_dxV(I,J,bi,bj)*recip_rSphere |
440 |
& |
441 |
& +ETAMEAN(I,J,bi,bj)*TWO * _tanPhiAtV(I,J,bi,bj) |
442 |
& *(UICEC(I+1,J,bi,bj)-UICEC(I-1,J,bi,bj)) |
443 |
& * 1. _d 0 /( _dxG(I,J,bi,bj) + _dxG(I-1,J,bi,bj)) |
444 |
& *recip_rSphere |
445 |
UVRT1(I,J,bi,bj)= 0.5 _d 0 * ( |
446 |
& ETA(I-1,J-1,bi,bj) * _recip_dxV(I,J,bi,bj) |
447 |
& +ETA(I-1,J ,bi,bj) * _recip_dxV(I,J,bi,bj) |
448 |
& ) * _recip_dxV(I,J,bi,bj) |
449 |
UVRT2(I,J,bi,bj)= 0.5 _d 0 * ( |
450 |
& ETA(I ,J-1,bi,bj) * _recip_dxV(I,J,bi,bj) |
451 |
& +ETA(I ,J ,bi,bj) * _recip_dxV(I,J,bi,bj) |
452 |
& ) * _recip_dxV(I,J,bi,bj) |
453 |
|
454 |
END DO |
455 |
END DO |
456 |
|
457 |
ENDDO |
458 |
ENDDO |
459 |
|
460 |
C NOW DO ITERATION |
461 |
300 CONTINUE |
462 |
|
463 |
cph--- iteration starts here |
464 |
cph--- need to kick out goto |
465 |
phexit = -1 |
466 |
|
467 |
C ITERATION START ----------------------------------------------------- |
468 |
#ifdef ALLOW_AUTODIFF_TAMC |
469 |
CADJ LOOP = iteration vice |
470 |
#endif /* ALLOW_AUTODIFF_TAMC */ |
471 |
|
472 |
DO 9000 M=1, solv_max_iters |
473 |
cph( |
474 |
IF ( phexit .EQ. -1 ) THEN |
475 |
cph) |
476 |
C NOW SET U(3)=U(1) |
477 |
DO bj=myByLo(myThid),myByHi(myThid) |
478 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
479 |
|
480 |
DO J=1,sNy |
481 |
DO I=1,sNx |
482 |
VICE(I,J,3,bi,bj)=VICE(I,J,1,bi,bj) |
483 |
END DO |
484 |
END DO |
485 |
|
486 |
DO I=1,sNx |
487 |
DO J=1,sNy |
488 |
IF(J.EQ.1) THEN |
489 |
AA2= _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj) |
490 |
& * 0.5 _d 0 *( |
491 |
& etaPlusZeta(I-1,J-1,bi,bj) + etaPlusZeta(I,J-1,bi,bj) |
492 |
& ) |
493 |
AA3=( AA2 |
494 |
& +( ZETAMEAN(I,J,bi,bj)-ETAMEAN(I,J,bi,bj) ) |
495 |
& * _tanPhiAtV(I,J-1,bi,bj) |
496 |
& * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere |
497 |
& + ETAMEAN(I,J,bi,bj)*TWO* _tanPhiAtV(I,J,bi,bj) |
498 |
& * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere ) |
499 |
& *VICE(I,J-1,1,bi,bj)*UVM(I,J,bi,bj) |
500 |
ELSE IF(J.EQ.sNy) THEN |
501 |
AA1= _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj) |
502 |
& * 0.5 _d 0 * ( |
503 |
& etaPlusZeta(I-1,J,bi,bj) + etaPlusZeta(I,J,bi,bj) |
504 |
& ) |
505 |
AA3=( AA1 |
506 |
& -( ZETAMEAN(I,J,bi,bj)-ETAMEAN(I,J,bi,bj)) |
507 |
& * _tanPhiAtV(I,J+1,bi,bj) |
508 |
& * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere |
509 |
& - ETAMEAN(I,J,bi,bj)*TWO* _tanPhiAtV(I,J,bi,bj) |
510 |
& * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere ) |
511 |
& *VICE(I,J+1,1,bi,bj)*UVM(I,J,bi,bj) |
512 |
ELSE |
513 |
AA3=ZERO |
514 |
END IF |
515 |
|
516 |
VRT(J)=FXY(I,J,bi,bj)+AA3+UVRT1(I,J,bi,bj)*VICE(I-1,J,1,bi,bj) |
517 |
& +UVRT2(I,J,bi,bj)*VICE(I+1,J,1,bi,bj) |
518 |
VRT(J)=VRT(J)*UVM(I,J,bi,bj) |
519 |
END DO |
520 |
|
521 |
DO J=1,sNy |
522 |
CVV(J)=CV(I,J,bi,bj) |
523 |
END DO |
524 |
VRT(1)=VRT(1)/BV(I,1,bi,bj) |
525 |
DO J=2,sNy |
526 |
JM=J-1 |
527 |
CVV(J)=CVV(J)/(BV(I,J,bi,bj)-AV(I,J,bi,bj)*CVV(JM)) |
528 |
VRT(J)=(VRT(J)-AV(I,J,bi,bj)*VRT(JM)) |
529 |
& /(BV(I,J,bi,bj)-AV(I,J,bi,bj)*CVV(JM)) |
530 |
END DO |
531 |
DO J=1,sNy-1 |
532 |
J1=sNy-J |
533 |
J2=J1+1 |
534 |
VRT(J1)=VRT(J1)-CVV(J1)*VRT(J2) |
535 |
END DO |
536 |
DO J=1,sNy |
537 |
VICE(I,J,1,bi,bj)=VICE(I,J,3,bi,bj) |
538 |
& +WFAV*(VRT(J)-VICE(I,J,3,bi,bj)) |
539 |
END DO |
540 |
ENDDO |
541 |
|
542 |
ENDDO |
543 |
ENDDO |
544 |
|
545 |
IF(MOD(M,SOLV_NCHECK).EQ.0) THEN |
546 |
DO bj=myByLo(myThid),myByHi(myThid) |
547 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
548 |
S2=ZERO |
549 |
DO J=1,sNy |
550 |
DO I=1,sNx |
551 |
UERR(I,J,bi,bj)=(VICE(I,J,1,bi,bj)-VICE(I,J,3,bi,bj)) |
552 |
& *UVM(I,J,bi,bj) |
553 |
S2=MAX(ABS(UERR(I,J,bi,bj)),S2) |
554 |
END DO |
555 |
END DO |
556 |
_GLOBAL_MAX_R8( S2, myThid ) |
557 |
ENDDO |
558 |
ENDDO |
559 |
C SAFEGUARD AGAINST BAD FORCING ETC |
560 |
IF(M.GT.1.AND.S2.GT.S2A) WFAV=WFAV2 |
561 |
S2A=S2 |
562 |
IF(S2.LT.LSR_ERROR) THEN |
563 |
ICOUNT2=M |
564 |
cph( |
565 |
cph GO TO 9001 |
566 |
phexit = 1 |
567 |
cph) |
568 |
END IF |
569 |
END IF |
570 |
|
571 |
CALL EXCH_3D_RL( VICE, 3, myThid ) |
572 |
|
573 |
cph( |
574 |
END IF |
575 |
cph) |
576 |
|
577 |
9000 CONTINUE |
578 |
cph 9001 CONTINUE |
579 |
C ITERATION END ----------------------------------------------------- |
580 |
|
581 |
IF ( debugLevel .GE. debLevB ) THEN |
582 |
_BEGIN_MASTER( myThid ) |
583 |
write(*,'(A,I6,1P2E22.14)')' V lsr iters, error = ',ICOUNT2,S2 |
584 |
_END_MASTER( myThid ) |
585 |
ENDIF |
586 |
|
587 |
C NOW END |
588 |
C NOW MAKE COROLIS TERM IMPLICIT |
589 |
DO bj=myByLo(myThid),myByHi(myThid) |
590 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
591 |
DO J=1,sNy |
592 |
DO I=1,sNx |
593 |
UICE(I,J,1,bi,bj)=UICE(I,J,1,bi,bj)*UVM(I,J,bi,bj) |
594 |
VICE(I,J,1,bi,bj)=VICE(I,J,1,bi,bj)*UVM(I,J,bi,bj) |
595 |
END DO |
596 |
END DO |
597 |
ENDDO |
598 |
ENDDO |
599 |
CALL EXCH_UV_3D_RL( UICE, VICE, .TRUE., 3, myThid ) |
600 |
|
601 |
#endif /* SEAICE_ALLOW_DYNAMICS */ |
602 |
#endif /* SEAICE_CGRID */ |
603 |
|
604 |
RETURN |
605 |
END |