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

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

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

Revision 1.34 - (show annotations) (download)
Fri Jan 11 20:18:57 2013 UTC (11 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65d, checkpoint65e, checkpoint65, checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64c, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f
Changes since 1.33: +2 -2 lines 
change type of local variable SINWAT from _RL to _RS (has been done already
 in seaice_evp.F, seaice_lsr.F & seaice_ocean_stress.F) to avoid mixed-type
 arguments in sign function (since fCori is _RS).
A better way: keep _RL SINWAT and use: SINWAT*SIGN(oneRS,fCori)
1 C $Header: /u/gcmpack/MITgcm/pkg/seaice/lsr.F,v 1.33 2012/10/18 10:06:42 mlosch Exp $
2 C $Name: $
3
4 #ifndef SEAICE_LSRBNEW
5 C for an alternative discretization of d/dx[ (zeta-eta) dV/dy]
6 C and d/dy[ (zeta-eta) dU/dx] uncomment this option
7 C#define SEAICE_TEST
8 #endif
9
10 #include "SEAICE_OPTIONS.h"
11
12 CStartOfInterface
13 SUBROUTINE lsr( ilcall, myThid )
14 C *==========================================================*
15 C | SUBROUTINE lsr |
16 C | o Solve ice momentum equation with an LSR dynamics solver|
17 C | (see Zhang and Hibler, JGR, 102, 8691-8702, 1997 |
18 C | and Zhang and Rothrock, MWR, 131, 845- 861, 2003) |
19 C | Written by Jinlun Zhang, PSC/UW, Feb-2001 |
20 C | zhang@apl.washington.edu |
21 C *==========================================================*
22 IMPLICIT NONE
23
24 C === Global variables ===
25 #include "SIZE.h"
26 #include "EEPARAMS.h"
27 #include "PARAMS.h"
28 #include "GRID.h"
29 #include "SEAICE_SIZE.h"
30 #include "SEAICE_PARAMS.h"
31 #include "SEAICE.h"
32 C#include "SEAICE_GRID.h"
33
34 #ifdef ALLOW_AUTODIFF_TAMC
35 # include "tamc.h"
36 #endif
37
38 C === Routine arguments ===
39 C myThid - Thread no. that called this routine.
40 INTEGER ilcall
41 INTEGER myThid
42 CEndOfInterface
43
44 #ifndef SEAICE_CGRID
45 #ifdef SEAICE_ALLOW_DYNAMICS
46
47 C === Local variables ===
48 C i,j,bi,bj - Loop counters
49
50 INTEGER i, j, m, bi, bj, j1, j2, im, jm
51 INTEGER ICOUNT1, ICOUNT2
52
53 _RL WFAU, WFAV, WFAU1, WFAV1, WFAU2, WFAV2
54 _RL AA3, S1, S2, S1A, S2A
55
56 _RL e11loc, e22loc, e12loc
57
58 _RL COSWAT
59 _RS SINWAT
60 _RL ECM2, DELT1, DELT2
61 _RL TEMPVAR
62
63 C diagonals of coefficient matrices
64 _RL AU (1:sNx,1:sNy,nSx,nSy)
65 _RL BU (1:sNx,1:sNy,nSx,nSy)
66 _RL CU (1:sNx,1:sNy,nSx,nSy)
67 _RL AV (1:sNx,1:sNy,nSx,nSy)
68 _RL BV (1:sNx,1:sNy,nSx,nSy)
69 _RL CV (1:sNx,1:sNy,nSx,nSy)
70
71 C coefficients for lateral points, u(j+/-1)
72 _RL uRt1(1:sNx,1:sNy,nSx,nSy)
73 _RL uRt2(1:sNx,1:sNy,nSx,nSy)
74 C coefficients for lateral points, v(i+/-1)
75 _RL vRt1(1:sNx,1:sNy,nSx,nSy)
76 _RL vRt2(1:sNx,1:sNy,nSx,nSy)
77 C RHS
78 _RL rhsU (1:sNx,1:sNy,nSx,nSy)
79 _RL rhsV (1:sNx,1:sNy,nSx,nSy)
80 C symmetric and anti-symmetric drag coefficients
81 _RL DRAGS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
82 _RL DRAGA (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
83 #ifdef SEAICE_LSRBNEW
84 LOGICAL doIterate4u, doIterate4v
85 CHARACTER*(MAX_LEN_MBUF) msgBuf
86 C coefficients of ice velocities in coefficient matrix
87 C for both U and V-equation
88 C XX: double derivative in X
89 C YY: double derivative in Y
90 C XM: metric term with derivative in X
91 C YM: metric term with derivative in Y
92 _RL UXX (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93 _RL UYY (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94 _RL UXM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95 _RL UYM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96 _RL VXX (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97 _RL VYY (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98 _RL VXM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99 _RL VYM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100 C abbreviations
101 _RL etaU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102 _RL zetaU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103 _RL etaV (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104 _RL zetaV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105 C contribution of sigma on righ hand side
106 _RL sig11(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107 _RL sig22(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108 _RL sig12(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109 _RL sig21(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110
111 C auxillary fields
112 _RL CUU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113 _RL CVV (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114 _RL URT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115 _RL VRT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116 #else
117 _RL URT(1-OLx:sNx+OLx), CUU(1-OLx:sNx+OLx)
118 _RL VRT(1-OLy:sNy+OLy), CVV(1-OLy:sNy+OLy)
119
120 _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
121 _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
122 _RL ETAMEAN (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
123 _RL ZETAMEAN (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
124
125 _RL dVdx (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126 _RL dVdy (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
127 _RL dUdx (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128 _RL dUdy (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
129 #ifdef SEAICE_TEST
130 _RL uz (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
131 _RL vz (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
132 #endif
133
134 INTEGER phexit
135
136 _RL AA1, AA2, AA4, AA5, AA6
137
138 #endif /* SEAICE_LSRBNEW */
139
140 _RL UERR
141
142 C abbreviations
143 _RL ucLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
144 _RL vcLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
145 _RL uTmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
146 _RL vTmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
147
148 C SET SOME VALUES
149 WFAU1=0.95 _d 0
150 WFAV1=0.95 _d 0
151 WFAU2=ZERO
152 WFAV2=ZERO
153
154 S1A=0.80 _d 0
155 S2A=0.80 _d 0
156 WFAU=WFAU1
157 WFAV=WFAV1
158
159 ICOUNT1=SOLV_MAX_ITERS
160 ICOUNT2=SOLV_MAX_ITERS
161
162 ECM2= 1. _d 0/(SEAICE_eccen**2)
163
164 SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
165 COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
166
167 DO bj=myByLo(myThid),myByHi(myThid)
168 DO bi=myBxLo(myThid),myBxHi(myThid)
169 DO j=1-OLy,sNy+OLy-1
170 DO i=1-OLx,sNx+OLx-1
171 ucLoc(I,J,bi,bj) = 0.25 _d 0 * (
172 & + uIceC(I+1,J,bi,bj) + uIceC(I+1,J+1,bi,bj)
173 & + uIceC(I ,J,bi,bj) + uIceC(I, J+1,bi,bj) )
174 vcLoc(I,J,bi,bj) = 0.25 _d 0 * (
175 & + vIceC(I+1,J,bi,bj) + vIceC(I+1,J+1,bi,bj)
176 & + vIceC(I ,J,bi,bj) + vIceC(I, J+1,bi,bj) )
177 ENDDO
178 ENDDO
179 DO J=1-OLy,sNy+OLy-1
180 DO I=1-OLy,sNx+OLy-1
181 e11loc = 0.5 _d 0 * _recip_dxF(I,J,bi,bj) *
182 & (uIceC(I+1,J+1,bi,bj)+uIceC(I+1,J,bi,bj)
183 & -uIceC(I, J+1,bi,bj)-uIceC(I, J,bi,bj))
184 & + vcLoc(I,J,bi,bj) * k2AtC(I,J,bi,bj)
185 e22loc = 0.5 _d 0 * _recip_dyF(I,J,bi,bj) *
186 & (vIceC(I+1,J+1,bi,bj)+vIceC(I,J+1,bi,bj)
187 & -vIceC(I+1,J, bi,bj)-vIceC(I,J, bi,bj))
188 & + ucLoc(I,J,bi,bj) * k1AtC(I,J,bi,bj)
189 e12loc = 0.5 _d 0*(
190 & 0.5 _d 0 * _recip_dyF(I,J,bi,bj) *
191 & (uIceC(I+1,J+1,bi,bj)+uIceC(I,J+1,bi,bj)
192 & -uIceC(I+1,J, bi,bj)-uIceC(I,J, bi,bj))
193 & + 0.5 _d 0 * _recip_dxF(I,J,bi,bj) *
194 & (vIceC(I+1,J+1,bi,bj)+vIceC(I+1,J,bi,bj)
195 & -vIceC(I, J+1,bi,bj)-vIceC(I, J,bi,bj))
196 & - vcLoc(I,J,bi,bj)*k1AtC(I,J,bi,bj)
197 & - ucLoc(I,J,bi,bj)*k2AtC(I,J,bi,bj) )
198 C NOW EVALUATE VISCOSITIES
199 DELT1=(e11loc**2+e22loc**2)*(1. _d 0+ECM2)
200 & +4.0 _d 0*ECM2*e12loc**2
201 & +2.0 _d 0*e11loc*e22loc*(1. _d 0-ECM2)
202 IF ( DELT1 .LE. SEAICE_EPS_SQ ) THEN
203 DELT2=SEAICE_EPS
204 ELSE
205 DELT2=SQRT(DELT1)
206 ENDIF
207 ZETA(I,J,bi,bj) = 0.5 _d 0*PRESS0(I,J,bi,bj)/DELT2
208 C NOW PUT MIN AND MAX VISCOSITIES IN
209 ZETA(I,J,bi,bj) = MIN(ZMAX(I,J,bi,bj),ZETA(I,J,bi,bj))
210 ZETA(I,J,bi,bj) = MAX(ZMIN(I,J,bi,bj),ZETA(I,J,bi,bj))
211 C NOW SET VISCOSITIES TO ZERO AT HEFFMFLOW PTS
212 ZETA(I,J,bi,bj) = ZETA(I,J,bi,bj)*HEFFM(I,J,bi,bj)
213 ETA(I,J,bi,bj) = ECM2*ZETA(I,J,bi,bj)
214 PRESS(I,J,bi,bj) = 2.0 _d 0*ZETA(I,J,bi,bj)*DELT2
215 ENDDO
216 ENDDO
217 DO j=1,sNy
218 DO i=1,sNx
219 C NOW SET UP NON-LINEAR WATER DRAG, FORCEX, FORCEY
220 TEMPVAR=(uIce(I,J,bi,bj)-GWATX(I,J,bi,bj))**2
221 & +(vIce(I,J,bi,bj)-GWATY(I,J,bi,bj))**2
222 IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN
223 IF ( TEMPVAR .LE. (QUART/SEAICE_waterDrag_south)**2 ) THEN
224 DWATN(I,J,bi,bj)=QUART
225 ELSE
226 DWATN(I,J,bi,bj)=SEAICE_waterDrag_south*SQRT(TEMPVAR)
227 ENDIF
228 ELSE
229 IF ( TEMPVAR .LE. (QUART/SEAICE_waterDrag)**2 ) THEN
230 DWATN(I,J,bi,bj)=QUART
231 ELSE
232 DWATN(I,J,bi,bj)=SEAICE_waterDrag*SQRT(TEMPVAR)
233 ENDIF
234 ENDIF
235 C NOW SET UP SYMMETTRIC DRAG
236 DRAGS(I,J,bi,bj)=DWATN(I,J,bi,bj)*COSWAT
237 C NOW SET UP ANTI SYMMETTRIC DRAG PLUS CORIOLIS
238 DRAGA(I,J,bi,bj)=DWATN(I,J,bi,bj)
239 & *SIGN(SINWAT, _fCori(I,J,bi,bj))
240 & + AMASS(I,J,bi,bj) * _fCoriG(I,J,bi,bj)
241 C NOW ADD IN CURRENT FORCE
242 FORCEX(I,J,bi,bj)=FORCEX0(I,J,bi,bj)+DWATN(I,J,bi,bj)
243 & *(COSWAT*GWATX(I,J,bi,bj)
244 & -SIGN(SINWAT, _fCori(I,J,bi,bj))*GWATY(I,J,bi,bj))
245 FORCEY(I,J,bi,bj)=FORCEY0(I,J,bi,bj)+DWATN(I,J,bi,bj)
246 & *(SIGN(SINWAT, _fCori(I,J,bi,bj))*GWATX(I,J,bi,bj)
247 & +COSWAT*GWATY(I,J,bi,bj))
248 #ifndef SEAICE_LSRBNEW
249 C NOW CALCULATE PRESSURE FORCE AND ADD TO EXTERNAL FORCE
250 C only for old code, in the new code the pressure force
251 C is added to the rhs stress tensor terms
252 FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
253 & -QUART * _recip_dxV(I,J,bi,bj)
254 & *(PRESS(I, J,bi,bj) + PRESS(I, J-1,bi,bj)
255 & - PRESS(I-1,J,bi,bj) - PRESS(I-1,J-1,bi,bj))
256 FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)
257 & -QUART * _recip_dyU(I,J,bi,bj)
258 & *(PRESS(I,J, bi,bj) + PRESS(I-1,J, bi,bj)
259 & - PRESS(I,J-1,bi,bj) - PRESS(I-1,J-1,bi,bj))
260 #endif /* not SEAICE_LSRBNEW */
261 FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
262 & +AMASS(I,J,bi,bj)/SEAICE_deltaTdyn*uIceNm1(I,J,bi,bj)
263 FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)
264 & +AMASS(I,J,bi,bj)/SEAICE_deltaTdyn*vIceNm1(I,J,bi,bj)
265 FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)*UVM(I,J,bi,bj)
266 FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)*UVM(I,J,bi,bj)
267 ENDDO
268 ENDDO
269 ENDDO
270 ENDDO
271
272 #ifdef SEAICE_LSRBNEW
273 DO bj=myByLo(myThid),myByHi(myThid)
274 DO bi=myBxLo(myThid),myBxHi(myThid)
275 C coefficients for matrices
276 DO j=1-OLy,sNy+OLy-1
277 DO i=1-OLx+1,sNx+OLx-1
278 etaU(I,J) = 0.5 _d 0 * (
279 & + eta(I,J,bi,bj) + eta(I-1,J,bi,bj) )
280 zetaU(I,J)= 0.5 _d 0 *(
281 & + zeta(I,J,bi,bj) + zeta(I-1,J,bi,bj) )
282 ENDDO
283 ENDDO
284 DO j=1-OLy+1,sNy+OLy-1
285 DO i=1-OLx,sNx+OLx-1
286 etaV(I,J) = 0.5 _d 0 * (
287 & + eta(I,J,bi,bj) + eta(I,J-1,bi,bj) )
288 zetaV(I,J)= 0.5 _d 0 *(
289 & + zeta(I,J,bi,bj) + zeta(I,J-1,bi,bj) )
290 ENDDO
291 ENDDO
292 C coefficients of uIce(I,J) and vIce(I,J) belonging to ...
293 DO J=1,sNy
294 DO I=0,sNx
295 C ... d/dx (eta+zeta) d/dx u
296 UXX(I,J) = _dyC(I,J,bi,bj) * (zetaV(I,J)+etaV(I,J))
297 & * _recip_dxG(I,J,bi,bj)
298 C ... d/dx (zeta-eta) k1 u
299 UXM(I,J) = _dyC(I,J,bi,bj) * (zetaV(I,J)-etaV(I,J))
300 & * k1AtV(I,J,bi,bj) * 0.5 _d 0
301 ENDDO
302 ENDDO
303 DO J=0,sNy
304 DO I=1,sNx
305 C ... d/dy eta d/dy u
306 UYY(I,J) = _dxC(I,J,bi,bj) * etaU(I,J)
307 & * _recip_dyG(I,J,bi,bj)
308 C ... d/dy eta k2 u
309 UYM(I,J) = _dxC(I,J,bi,bj) * etaU(I,J)
310 & * k2AtU(I,J,bi,bj) * 0.5 _d 0
311 ENDDO
312 ENDDO
313 DO J=1,sNy
314 DO I=0,sNx
315 C ... d/dx eta dv/dx
316 VXX(I,J) = _dyC(I,J,bi,bj) * etaV(I,J)
317 & * _recip_dxG(I,J,bi,bj)
318 C ... d/dx eta k1 v
319 VXM(I,J) = _dyC(I,J,bi,bj) * etaV(I,J)
320 & * k1AtV(I,J,bi,bj) * 0.5 _d 0
321 ENDDO
322 ENDDO
323 DO J=0,sNy
324 DO I=1,sNx
325 C ... d/dy eta+zeta dv/dy
326 VYY(I,J) = _dxC(I,J,bi,bj) * (zetaU(I,J)+etaU(I,J))
327 & * _recip_dyG(I,J,bi,bj)
328 C ... d/dy (zeta-eta) k2 v
329 VYM(I,J) = _dxC(I,J,bi,bj) * (zetaU(I,J)-etaU(I,J))
330 & * k2AtU(I,J,bi,bj) * 0.5 _d 0
331 ENDDO
332 ENDDO
333
334 C assemble coefficient matrix of uIce
335 C beware of sign convention: because this
336 C is the left hand side we calculate -grad(sigma), but the coefficients
337 C of U(I,J+/-1) are counted on the right hand side
338 DO J=1,sNy
339 DO I=1,sNx
340 C coefficients for uIce(I-1,J)
341 AU(I,J,bi,bj)= ( - UXX(I-1,J) + UXM(I-1,J) )
342 & * UVM(I,J,bi,bj)
343 C coefficients for uIce(I+1,J)
344 CU(I,J,bi,bj)= ( - UXX(I ,J) - UXM(I ,J) )
345 & * UVM(I,J,bi,bj)
346 C coefficients for uIce(I,J)
347 BU(I,J,bi,bj)=(ONE - UVM(I,J,bi,bj)) +
348 & ( UXX(I-1,J) + UXX(I,J) + UYY(I,J-1) + UYY(I,J)
349 & + UXM(I-1,J) - UXM(I,J) - UYM(I,J-1) + UYM(I,J)
350 & ) * UVM(I,J,bi,bj)
351 C coefficients of uIce(I,J-1)
352 uRt1(I,J,bi,bj)= UYY(I,J-1) + UYM(I,J-1)
353 C coefficients of uIce(I,J+1)
354 uRt2(I,J,bi,bj)= UYY(I,J ) - UYM(I,J )
355 ENDDO
356 ENDDO
357
358 C now we need to normalize everything by the grid cell area
359 DO J=1,sNy
360 DO I=1,sNx
361 AU(I,J,bi,bj) = AU(I,J,bi,bj) * recip_rAz(I,J,bi,bj)
362 CU(I,J,bi,bj) = CU(I,J,bi,bj) * recip_rAz(I,J,bi,bj)
363 C here we need add in the contribution from the time derivative
364 C and the symmetric drag term; must be done after normalizing
365 BU(I,J,bi,bj) = BU(I,J,bi,bj) * recip_rAz(I,J,bi,bj)
366 & + UVM(I,J,bi,bj) *
367 & ( AMASS(I,J,bi,bj)/SEAICE_deltaTdyn
368 & + DRAGS(I,J,bi,bj) )
369 uRt1(I,J,bi,bj) = uRt1(I,J,bi,bj) * recip_rAz(I,J,bi,bj)
370 uRt2(I,J,bi,bj) = uRt2(I,J,bi,bj) * recip_rAz(I,J,bi,bj)
371 ENDDO
372 ENDDO
373
374 C assemble coefficient matrix of vIce
375 C beware of sign convention: because this
376 C is the left hand side we calculate -grad(sigma), but the coefficients
377 C of V(I,J+/-1) are counted on the right hand side
378 DO J=1,sNy
379 DO I=1,sNx
380 C coefficients for vIce(I,J-1)
381 AV(I,J,bi,bj)=( - VYY(I,J-1) + VYM(I,J-1)
382 & ) * UVM(I,J,bi,bj)
383 C coefficients for vIce(I,J+1)
384 CV(I,J,bi,bj)=( - VYY(I,J ) - VYM(I,J )
385 & ) * UVM(I,J,bi,bj)
386 C coefficients for vIce(I,J)
387 BV(I,J,bi,bj)= (ONE - UVM(I,J,bi,bj)) +
388 & ( VXX(I,J) + VXX(I-1,J) + VYY(I,J) + VYY(I,J-1)
389 & + VXM(I,J) - VXM(I-1,J) - VYM(I,J) + VYM(I,J-1)
390 & ) * UVM(I,J,bi,bj)
391 C coefficients for V(I-1,J)
392 vRt1(I,J,bi,bj) = VXX(I-1,J) + VXM(I-1,J)
393 C coefficients for V(I+1,J)
394 vRt2(I,J,bi,bj) = VXX(I ,J) - VXM(I ,J)
395 ENDDO
396 ENDDO
397
398 C now we need to normalize everything by the grid cell area
399 DO J=1,sNy
400 DO I=1,sNx
401 AV(I,J,bi,bj) = AV(I,J,bi,bj) * recip_rAz(I,J,bi,bj)
402 CV(I,J,bi,bj) = CV(I,J,bi,bj) * recip_rAz(I,J,bi,bj)
403 C here we need add in the contribution from the time derivative
404 C and the symmetric drag term; must be done after normalizing
405 BV(I,J,bi,bj) = BV(I,J,bi,bj) * recip_rAz(I,J,bi,bj)
406 & + UVM(I,J,bi,bj) *
407 & ( AMASS(I,J,bi,bj)/SEAICE_deltaTdyn
408 & + DRAGS(I,J,bi,bj) )
409 vRt1(I,J,bi,bj) = vRt1(I,J,bi,bj) * recip_rAz(I,J,bi,bj)
410 vRt2(I,J,bi,bj) = vRt2(I,J,bi,bj) * recip_rAz(I,J,bi,bj)
411 ENDDO
412 ENDDO
413
414 C set up right-hand sides
415 DO j=1,sNy
416 DO i=0,sNx
417 C at V-point
418 sig11(I,J) =
419 & (zetaV(I,J)-etaV(I,J))
420 & * (vcLoc(I,J,bi,bj)-vcLoc(I,J-1,bi,bj))
421 & * _recip_dyC(I,J,bi,bj)
422 & + (zetaV(I,J)+etaV(I,J))*k2atV(I,J,bi,bj)
423 & * 0.5 _d 0 * (vIceC(I,J,bi,bj)+vIceC(I+1,J,bi,bj))
424 & - 0.5 _d 0 * (PRESS(I,J,bi,bj)+PRESS(I,J-1,bi,bj))
425 sig12(I,J) = etaV(I,J)
426 & * (ucLoc(I,J,bi,bj)-ucLoc(I,J-1,bi,bj))
427 & * _recip_dyC(I,J,bi,bj)
428 & - etaV(I,J) * k2AtV(I,J,bi,bj)
429 & * 0.5 _d 0 * (uIceC(I,J,bi,bj)+uIceC(I+1,J,bi,bj))
430 ENDDO
431 ENDDO
432 DO j=0,sNy
433 DO i=1,sNx
434 C at U-point
435 sig22(I,J) =
436 & (zetaU(I,J)-etaU(I,J))
437 & * (ucLoc(I,J,bi,bj)-ucLoc(I-1,J,bi,bj))
438 & * _recip_dxC(I,J,bi,bj)
439 & + (zetaU(I,J)+etaU(I,J))*k2atU(I,J,bi,bj)
440 & * 0.5 _d 0 * (uIceC(I,J,bi,bj)+uIceC(I,J+1,bi,bj))
441 & - 0.5 _d 0 * (PRESS(I,J,bi,bj)+PRESS(I-1,J,bi,bj))
442 sig21(I,J) = etaU(I,J)
443 & * (vcLoc(I,J,bi,bj)-vcLoc(I-1,J,bi,bj))
444 & * _recip_dxC(I,J,bi,bj)
445 & - etaU(I,J) * k1AtU(I,J,bi,bj)
446 & * 0.5 _d 0 * (vIceC(I,J,bi,bj)+vIceC(I,J+1,bi,bj))
447 ENDDO
448 ENDDO
449 DO j=1,sNy
450 DO i=1,sNx
451 rhsU(I,J,bi,bj) = DRAGA(I,J,bi,bj)*vIceC(I,J,bi,bj)
452 & +FORCEX(I,J,bi,bj)
453 & + ( _dyC(I, J, bi,bj) * sig11(I, J )
454 & - _dyC(I-1,J, bi,bj) * sig11(I-1,J )
455 & + _dxC(I, J, bi,bj) * sig21(I, J )
456 & - _dxC(I, J-1,bi,bj) * sig21(I, J-1) )
457 & * recip_rAz(I,J,bi,bj) * UVM(I,J,bi,bj)
458 rhsV(I,J,bi,bj) = - DRAGA(I,J,bi,bj)*uIceC(I,J,bi,bj)
459 & +FORCEY(I,J,bi,bj)
460 & + ( _dyC(I, J, bi,bj) * sig12(I, J )
461 & - _dyC(I-1,J, bi,bj) * sig12(I-1,J )
462 & + _dxC(I, J, bi,bj) * sig22(I, J )
463 & - _dxC(I, J-1,bi,bj) * sig22(I, J-1) )
464 & * recip_rAz(I,J,bi,bj) * UVM(I,J,bi,bj)
465 ENDDO
466 ENDDO
467 C bi/bj-loops
468 ENDDO
469 ENDDO
470
471 C NOW DO ITERATION
472
473 doIterate4u = .TRUE.
474 doIterate4v = .TRUE.
475
476 C ITERATION START -----------------------------------------------------
477
478 DO m = 1, SOLV_MAX_ITERS
479 IF ( doIterate4u .OR. doIterate4v ) THEN
480
481 IF ( useCubedSphereExchange ) THEN
482 doIterate4u = .TRUE.
483 doIterate4v = .TRUE.
484 ENDIF
485
486 DO bj=myByLo(myThid),myByHi(myThid)
487 DO bi=myBxLo(myThid),myBxHi(myThid)
488
489 C-jmc: get less TAF warnings when always (no if doIterate) saving uIce,vIce:
490 C save uIce prior to iteration, NOW SET U(3)=U(1)
491 DO j=1-OLy,sNy+OLy
492 DO i=1-OLx,sNx+OLx
493 uTmp(I,J,bi,bj)=uIce(I,J,bi,bj)
494 ENDDO
495 ENDDO
496 C save vIce prior to iteration, NOW SET V(3)=V(1)
497 DO j=1-OLy,sNy+OLy
498 DO i=1-OLx,sNx+OLx
499 vTmp(I,J,bi,bj)=vIce(I,J,bi,bj)
500 ENDDO
501 ENDDO
502
503 IF ( doIterate4u ) THEN
504 C Solve for uIce :
505 DO J=1,sNy
506 DO I=1,sNx
507 AA3 = ZERO
508 IF (I.EQ.1) AA3 = AA3 - AU(I,J,bi,bj)*uIce(I-1,J,bi,bj)
509 IF (I.EQ.sNx) AA3 = AA3 - CU(I,J,bi,bj)*uIce(I+1,J,bi,bj)
510
511 URT(I,J)=rhsU(I,J,bi,bj)
512 & + AA3
513 #ifdef SEAICE_VECTORIZE_LSR
514 & + uRt1(I,J,bi,bj)*uTmp(I,J-1,bi,bj)
515 & + uRt2(I,J,bi,bj)*uTmp(I,J+1,bi,bj)
516 #else
517 & + uRt1(I,J,bi,bj)*uIce(I,J-1,bi,bj)
518 & + uRt2(I,J,bi,bj)*uIce(I,J+1,bi,bj)
519 #endif /* SEAICE_VECTORIZE_LSR */
520 URT(I,J)=URT(I,J)* UVM(I,J,bi,bj)
521 ENDDO
522
523 DO I=1,sNx
524 CUU(I,J)=CU(I,J,bi,bj)
525 ENDDO
526 CUU(1,J)=CUU(1,J)/BU(1,J,bi,bj)
527 URT(1,J)=URT(1,J)/BU(1,J,bi,bj)
528 #ifdef SEAICE_VECTORIZE_LSR
529 ENDDO
530 C start a new loop with reversed order to support automatic vectorization
531 DO I=2,sNx
532 IM=I-1
533 DO J=1,sNy
534 #else /* do not SEAICE_VECTORIZE_LSR */
535 DO I=2,sNx
536 IM=I-1
537 #endif /* SEAICE_VECTORIZE_LSR */
538 CUU(I,J)=CUU(I,J)/(BU(I,J,bi,bj)-AU(I,J,bi,bj)*CUU(IM,J))
539 URT(I,J)=(URT(I,J)-AU(I,J,bi,bj)*URT(IM,J))
540 & /(BU(I,J,bi,bj)-AU(I,J,bi,bj)*CUU(IM,J))
541 ENDDO
542 #ifdef SEAICE_VECTORIZE_LSR
543 ENDDO
544 C go back to original order
545 DO J=1,sNy
546 #endif /* SEAICE_VECTORIZE_LSR */
547 DO I=1,sNx-1
548 J1=sNx-I
549 J2=J1+1
550 URT(J1,J)=URT(J1,J)-CUU(J1,J)*URT(J2,J)
551 ENDDO
552 DO I=1,sNx
553 uIce(I,J,bi,bj)=uTmp(I,J,bi,bj)
554 & +WFAU*(URT(I,J)-uTmp(I,J,bi,bj))
555 ENDDO
556 ENDDO
557 C-- end doIterate4u
558 ENDIF
559
560 IF ( doIterate4v ) THEN
561 C Solve for vIce
562 DO I=1,sNx
563 DO J=1,sNy
564 AA3 = ZERO
565 IF (J.EQ.1) AA3 = AA3 - AV(I,J,bi,bj)*vIce(I,J-1,bi,bj)
566 IF (J.EQ.sNy) AA3 = AA3 - CV(I,J,bi,bj)*vIce(I,J+1,bi,bj)
567
568 VRT(I,J)=rhsV(I,J,bi,bj)
569 & + AA3
570 #ifdef SEAICE_VECTORIZE_LSR
571 & + vRt1(I,J,bi,bj)*vTmp(I-1,J,bi,bj)
572 & + vRt2(I,J,bi,bj)*vTmp(I+1,J,bi,bj)
573 #else
574 & + vRt1(I,J,bi,bj)*vIce(I-1,J,bi,bj)
575 & + vRt2(I,J,bi,bj)*vIce(I+1,J,bi,bj)
576 #endif /* SEAICE_VECTORIZE_LSR */
577 VRT(I,J)=VRT(I,J)* UVM(I,J,bi,bj)
578 ENDDO
579
580 DO J=1,sNy
581 CVV(I,J)=CV(I,J,bi,bj)
582 ENDDO
583 CVV(I,1)=CVV(I,1)/BV(I,1,bi,bj)
584 VRT(I,1)=VRT(I,1)/BV(I,1,bi,bj)
585 DO J=2,sNy
586 JM=J-1
587 CVV(I,J)=CVV(I,J)/(BV(I,J,bi,bj)-AV(I,J,bi,bj)*CVV(I,JM))
588 VRT(I,J)=(VRT(I,J)-AV(I,J,bi,bj)*VRT(I,JM))
589 & /(BV(I,J,bi,bj)-AV(I,J,bi,bj)*CVV(I,JM))
590 ENDDO
591 DO J=1,sNy-1
592 J1=sNy-J
593 J2=J1+1
594 VRT(I,J1)=VRT(I,J1)-CVV(I,J1)*VRT(I,J2)
595 ENDDO
596 DO J=1,sNy
597 vIce(I,J,bi,bj)=vTmp(I,J,bi,bj)
598 & +WFAV*(VRT(I,J)-vTmp(I,J,bi,bj))
599 ENDDO
600 ENDDO
601 C-- end doIterate4v
602 ENDIF
603
604 C end bi,bj-loops
605 ENDDO
606 ENDDO
607
608 IF ( doIterate4u.AND.MOD(m,SOLV_NCHECK).EQ.0) THEN
609 S1=ZERO
610 DO bj=myByLo(myThid),myByHi(myThid)
611 DO bi=myBxLo(myThid),myBxHi(myThid)
612 DO J=1,sNy
613 DO I=1,sNx
614 UERR=(uIce(I,J,bi,bj)-uTmp(I,J,bi,bj))
615 & * UVM(I,J,bi,bj)
616 S1=MAX(ABS(UERR),S1)
617 ENDDO
618 ENDDO
619 ENDDO
620 ENDDO
621 _GLOBAL_MAX_RL( S1, myThid )
622 c WRITE(standardMessageUnit,'(A,2I6,1P4E16.9)')
623 c & ' U iters,error,WF = ',ilcall,M,S1,S1A,WFAU
624 C SAFEGUARD AGAINST BAD FORCING ETC
625 IF(m.GT.1.AND.S1.GT.S1A) WFAU=WFAU2
626 S1A=S1
627 IF(S1.LT.LSR_ERROR) THEN
628 ICOUNT1=m
629 doIterate4u = .FALSE.
630 ENDIF
631 ENDIF
632
633 IF ( doIterate4v.AND.MOD(m,SOLV_NCHECK).EQ.0) THEN
634 S2=ZERO
635 DO bj=myByLo(myThid),myByHi(myThid)
636 DO bi=myBxLo(myThid),myBxHi(myThid)
637 DO J=1,sNy
638 DO I=1,sNx
639 UERR=(vIce(I,J,bi,bj)-vTmp(I,J,bi,bj))
640 & * UVM(I,J,bi,bj)
641 S2=MAX(ABS(UERR),S2)
642 ENDDO
643 ENDDO
644 ENDDO
645 ENDDO
646 _GLOBAL_MAX_RL( S2, myThid )
647 C SAFEGUARD AGAINST BAD FORCING ETC
648 IF(m.GT.1.AND.S2.GT.S2A) WFAV=WFAV2
649 S2A=S2
650 IF(S2.LT.LSR_ERROR) THEN
651 ICOUNT2=m
652 doIterate4v = .FALSE.
653 ENDIF
654 ENDIF
655
656 CALL EXCH_UV_XY_RL( uIce, vIce,.TRUE.,myThid)
657
658 C-- end doIterate4u or doIterate4v
659 ENDIF
660 ENDDO
661 C ITERATION END -----------------------------------------------------
662 #ifdef ALLOW_DEBUG
663 IF ( debugLevel .GE. debLevD ) THEN
664 WRITE(msgBuf,'(A,I3,A)')
665 & 'Uice post iter (SEAICE_LSR', MOD(ilcall,1000), ')'
666 CALL DEBUG_STATS_RL( 1, uIce, msgBuf, myThid )
667 WRITE(msgBuf,'(A,I3,A)')
668 & 'Vice post iter (SEAICE_LSR', MOD(ilcall,1000), ')'
669 CALL DEBUG_STATS_RL( 1, vIce, msgBuf, myThid )
670 ENDIF
671 #endif /* ALLOW_DEBUG */
672 IF ( doIterate4u .OR. doIterate4v ) THEN
673 WRITE(msgBuf,'(2A,I10,A,I4,A)') '** WARNING ** SEAICE_LSR ',
674 & '(it=', -9999, ',', ilcall,') did not converge :'
675 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
676 & SQUEEZE_RIGHT, myThid )
677 WRITE(msgBuf,'(2(A,I6,0PF6.3,1PE14.6))')
678 & ' nIt,wFU,dU=', ICOUNT1, WFAU, S1,
679 & ' ; nIt,wFV,dV=', ICOUNT2, WFAV, S2
680 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
681 & SQUEEZE_RIGHT, myThid )
682 ENDIF
683
684 #else /* SEAICE_LSRBNEW */
685 C SOLVE FOR UICE
686
687 #ifdef ALLOW_AUTODIFF_TAMC
688 cph That is an important one! Note, that
689 cph * lsr is called twice, thus the icall index
690 cph * this storing is still outside the iteration loop
691 CADJ STORE uice = comlev1_dynsol,
692 CADJ & key = ikey_dynamics + (ilcall-1)*nchklev_1
693 CADJ STORE vice = comlev1_dynsol,
694 CADJ & key = ikey_dynamics + (ilcall-1)*nchklev_1
695 #endif /* ALLOW_AUTODIFF_TAMC */
696
697 DO bj=myByLo(myThid),myByHi(myThid)
698 DO bi=myBxLo(myThid),myBxHi(myThid)
699 DO j=1-OLy,sNy+OLy
700 DO i=1-OLx,sNx+OLx
701 etaPlusZeta(I,J,bi,bj) = ETA(I,J,bi,bj)+ZETA(I,J,bi,bj)
702 zetaMinusEta(I,J,bi,bj) = ZETA(I,J,bi,bj)-ETA(I,J,bi,bj)
703 ENDDO
704 ENDDO
705 DO j=1-OLy+1,sNy+OLy
706 DO i=1-OLx+1,sNx+OLx
707 ETAMEAN(I,J,bi,bj) =QUART*(
708 & ETA(I,J-1,bi,bj) + ETA(I-1,J-1,bi,bj)
709 & +ETA(I,J ,bi,bj) + ETA(I-1,J ,bi,bj))
710 ZETAMEAN(I,J,bi,bj)=QUART*(
711 & ZETA(I,J-1,bi,bj) + ZETA(I-1,J-1,bi,bj)
712 & +ZETA(I,J ,bi,bj) + ZETA(I-1,J ,bi,bj))
713 ENDDO
714 ENDDO
715 ENDDO
716 ENDDO
717
718 DO bj=myByLo(myThid),myByHi(myThid)
719 DO bi=myBxLo(myThid),myBxHi(myThid)
720
721 DO J=1,sNy
722 DO I=1,sNx
723 AA1=( etaPlusZeta(I ,J-1,bi,bj) * _recip_dxF(I ,J-1,bi,bj)
724 & +etaPlusZeta(I ,J ,bi,bj) * _recip_dxF(I ,J ,bi,bj)
725 & )*0.5 _d 0 * _recip_dxV(I,J,bi,bj) * UVM(I,J,bi,bj)
726 AA2=( etaPlusZeta(I-1,J-1,bi,bj) * _recip_dxF(I-1,J-1,bi,bj)
727 & +etaPlusZeta(I-1,J ,bi,bj) * _recip_dxF(I-1,J ,bi,bj)
728 & )*0.5 _d 0 * _recip_dxV(I,J,bi,bj) * UVM(I,J,bi,bj)
729 AA3= 0.5 _d 0 *(ETA(I-1,J ,bi,bj)+ETA(I,J ,bi,bj))
730 AA4= 0.5 _d 0 *(ETA(I-1,J-1,bi,bj)+ETA(I,J-1,bi,bj))
731 AA5= -(AA3-AA4) * _tanPhiAtV(I,J,bi,bj)
732 & * _recip_dyU(I,J,bi,bj)*recip_rSphere
733 AA6=TWO*ETAMEAN(I,J,bi,bj) *recip_rSphere*recip_rSphere
734 & * _tanPhiAtV(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj)
735 AU(I,J,bi,bj)=-AA2
736 CU(I,J,bi,bj)=-AA1
737 BU(I,J,bi,bj)=(ONE-UVM(I,J,bi,bj))
738 & - AU(I,J,bi,bj) - CU(I,J,bi,bj)
739 & + ((AA3+AA4)*_recip_dyU(I,J,bi,bj)*_recip_dyU(I,J,bi,bj)
740 & + AA5 + AA6
741 & + AMASS(I,J,bi,bj)/SEAICE_deltaTdyn
742 & + DRAGS(I,J,bi,bj)
743 & )*UVM(I,J,bi,bj)
744 END DO
745 END DO
746
747 DO J=1,sNy
748 AU(1,J,bi,bj)=ZERO
749 CU(sNx,J,bi,bj)=ZERO
750 CU(1,J,bi,bj)=CU(1,J,bi,bj)/BU(1,J,bi,bj)
751 END DO
752
753 C now set up right-hand side
754 DO J=1-OLy,sNy+OLy-1
755 DO I=1-OLx,sNx+OLx-1
756 dVdy(I,J) = 0.5 _d 0 * (
757 & ( VICEC(I+1,J+1,bi,bj) - VICEC(I+1,J ,bi,bj) )
758 & * _recip_dyG(I+1,J,bi,bj)
759 & +(VICEC(I ,J+1,bi,bj) - VICEC(I ,J ,bi,bj) )
760 & * _recip_dyG(I, J,bi,bj) )
761 dVdx(I,J) = 0.5 _d 0 * (
762 & ( VICEC(I+1,J+1,bi,bj) - VICEC(I ,J+1,bi,bj) )
763 & * _recip_dxG(I,J+1,bi,bj)
764 & +(VICEC(I+1,J ,bi,bj) - VICEC(I ,J ,bi,bj) )
765 & * _recip_dxG(I,J, bi,bj) )
766 ENDDO
767 ENDDO
768 #ifdef SEAICE_TEST
769 DO j=1-OLy,sNy+OLy-1
770 DO i=1-OLx,sNx+OLx-1
771 vz(i,j) = quart * (
772 & vicec(i,j,bi,bj) + vicec(i+1,j,bi,bj) )
773 vz(i,j)= vz(i,j) + quart * (
774 & vicec(i,j+1,bi,bj) + vicec(i+1,j+1,bi,bj) )
775 ENDDO
776 ENDDO
777 #endif
778 DO J=1,sNy
779 DO I=1,sNx
780 rhsU(I,J,bi,bj)=DRAGA(I,J,bi,bj)*VICEC(I,J,bi,bj)
781 & +FORCEX(I,J,bi,bj)
782 #ifdef SEAICE_TEST
783 & + ( 0.5 _d 0 *
784 & (zetaMinusEta(i,j,bi,bj)+zetaMinusEta(i,j-1,bi,bj))
785 & *(vz(i,j)-vz(i,j-1)) * _recip_dyC(i,j,bi,bj)
786 & - 0.5 _d 0 *
787 & (zetaMinusEta(i-1,j,bi,bj)+zetaMinusEta(i-1,j-1,bi,bj))
788 & *(vz(i-1,j)-vz(i-1,j-1)) * _recip_dyC(i-1,j,bi,bj)
789 & ) * _recip_dxV(i,j,bi,bj)
790 #else
791 & + ( zetaMinusEta(I ,J ,bi,bj) * dVdy(I ,J )
792 & + zetaMinusEta(I ,J-1,bi,bj) * dVdy(I ,J-1)
793 & - zetaMinusEta(I-1,J ,bi,bj) * dVdy(I-1,J )
794 & - zetaMinusEta(I-1,J-1,bi,bj) * dVdy(I-1,J-1)
795 & )* 0.5 _d 0 * _recip_dxV(I,J,bi,bj)
796 #endif
797 &
798 & + ( ETA (I ,J ,bi,bj) * dVdx(I ,J )
799 & + ETA (I-1,J ,bi,bj) * dVdx(I-1,J )
800 & - ETA (I ,J-1,bi,bj) * dVdx(I ,J-1)
801 & - ETA (I-1,J-1,bi,bj) * dVdx(I-1,J-1)
802 & ) * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)
803 &
804 & -(etaPlusZeta(I ,J ,bi,bj)+etaPlusZeta(I ,J-1,bi,bj)
805 & -etaPlusZeta(I-1,J-1,bi,bj)-etaPlusZeta(I-1,J ,bi,bj))
806 & * VICEC(I,J,bi,bj)
807 & * _tanPhiAtV(I,J,bi,bj)
808 & * 0.5 _d 0 * _recip_dxV(I,J,bi,bj)*recip_rSphere
809 &
810 & -(ETAMEAN(I,J,bi,bj)+ZETAMEAN(I,J,bi,bj))
811 & *(VICEC(I+1,J,bi,bj) - VICEC(I-1,J,bi,bj))
812 & * _tanPhiAtV(I,J,bi,bj)
813 & * 1.0 _d 0 /( _dxG(I,J,bi,bj) + _dxG(I-1,J,bi,bj) )
814 & *recip_rSphere
815 &
816 & -ETAMEAN(I,J,bi,bj)
817 & *(VICEC(I+1,J,bi,bj) - VICEC(I-1,J,bi,bj))
818 & *TWO* _tanPhiAtV(I,J,bi,bj)
819 & * 1.0 _d 0 /( _dxG(I,J,bi,bj) + _dxG(I-1,J,bi,bj) )
820 & *recip_rSphere
821
822 URT1(I,J,bi,bj)=
823 & 0.5 _d 0 * (ETA(I-1,J-1,bi,bj)+ETA(I,J-1,bi,bj))
824 & * _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj)
825 & - ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J-1,bi,bj)
826 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
827 & + TWO*ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj)
828 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
829 URT2(I,J,bi,bj)=
830 & 0.5 _d 0 * (ETA(I-1,J,bi,bj)+ETA(I,J,bi,bj))
831 & * _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj)
832 & + ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J+1,bi,bj)
833 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
834 & - TWO*ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj)
835 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
836 END DO
837 END DO
838
839 ENDDO
840 ENDDO
841
842 C NOW DO ITERATION
843
844 cph--- iteration starts here
845 cph--- need to kick out goto
846 phexit = -1
847
848 C ITERATION START -----------------------------------------------------
849 #ifdef ALLOW_AUTODIFF_TAMC
850 CADJ LOOP = iteration uice
851 #endif /* ALLOW_AUTODIFF_TAMC */
852 DO M=1, solv_max_iters
853 IF ( phexit .EQ. -1 ) THEN
854
855 DO bj=myByLo(myThid),myByHi(myThid)
856 DO bi=myBxLo(myThid),myBxHi(myThid)
857 C NOW SET U(3)=U(1)
858 DO J=1,sNy
859 DO I=1,sNx
860 UTMP(I,J,bi,bj)=UICE(I,J,bi,bj)
861 END DO
862 END DO
863
864 DO J=1,sNy
865 DO I=1,sNx
866 IF(I.EQ.1) THEN
867 AA2=(etaPlusZeta(I-1,J-1,bi,bj) * _recip_dxF(I-1,J-1,bi,bj)
868 & +etaPlusZeta(I-1,J ,bi,bj) * _recip_dxF(I-1,J ,bi,bj)
869 & )*0.5 _d 0 * _recip_dxV(I,J,bi,bj)
870 AA3=AA2*UICE(I-1,J,bi,bj)*UVM(I,J,bi,bj)
871 ELSE IF(I.EQ.sNx) THEN
872 AA1=(etaPlusZeta(I ,J-1,bi,bj) * _recip_dxF(I ,J-1,bi,bj)
873 & +etaPlusZeta(I ,J ,bi,bj) * _recip_dxF(I ,J ,bi,bj)
874 & )*0.5 _d 0 * _recip_dxV(I,J,bi,bj)
875 AA3=AA1*UICE(I+1,J,bi,bj)*UVM(I,J,bi,bj)
876 ELSE
877 AA3=ZERO
878 END IF
879 URT(I)=rhsU(I,J,bi,bj)+AA3
880 & +URT1(I,J,bi,bj)*UICE(I,J-1,bi,bj)
881 & +URT2(I,J,bi,bj)*UICE(I,J+1,bi,bj)
882 URT(I)=URT(I)*UVM(I,J,bi,bj)
883 END DO
884
885 DO I=1,sNx
886 CUU(I)=CU(I,J,bi,bj)
887 END DO
888 URT(1)=URT(1)/BU(1,J,bi,bj)
889 DO I=2,sNx
890 IM=I-1
891 CUU(I)=CUU(I)/(BU(I,J,bi,bj)-AU(I,J,bi,bj)*CUU(IM))
892 URT(I)=(URT(I)-AU(I,J,bi,bj)*URT(IM))
893 & /(BU(I,J,bi,bj)-AU(I,J,bi,bj)*CUU(IM))
894 END DO
895 DO I=1,sNx-1
896 J1=sNx-I
897 J2=J1+1
898 URT(J1)=URT(J1)-CUU(J1)*URT(J2)
899 END DO
900 DO I=1,sNx
901 UICE(I,J,bi,bj)=UTMP(I,J,bi,bj)
902 & +WFAU*(URT(I)-UTMP(I,J,bi,bj))
903 END DO
904
905 END DO
906
907 ENDDO
908 ENDDO
909
910 IF(MOD(M,SOLV_NCHECK).EQ.0) THEN
911 S1=ZERO
912 DO bj=myByLo(myThid),myByHi(myThid)
913 DO bi=myBxLo(myThid),myBxHi(myThid)
914 DO J=1,sNy
915 DO I=1,sNx
916 UERR=(UICE(I,J,bi,bj)-UTMP(I,J,bi,bj))
917 & *UVM(I,J,bi,bj)
918 S1=MAX(ABS(UERR),S1)
919 END DO
920 END DO
921 ENDDO
922 ENDDO
923 _GLOBAL_MAX_RL( S1, myThid )
924 C SAFEGUARD AGAINST BAD FORCING ETC
925 IF(M.GT.1.AND.S1.GT.S1A) WFAU=WFAU2
926 S1A=S1
927 IF(S1.LT.LSR_ERROR) THEN
928 ICOUNT1=M
929 phexit = 1
930 END IF
931 END IF
932 _EXCH_XY_RL( UICE, myThid )
933
934 ENDIF
935 ENDDO
936 C ITERATION END -----------------------------------------------------
937
938 IF ( debugLevel .GE. debLevC ) THEN
939 _BEGIN_MASTER( myThid )
940 write(*,'(A,I6,1P2E22.14)')' U lsr iters, error = ',ICOUNT1,S1
941 _END_MASTER( myThid )
942 ENDIF
943
944 C NOW FOR VICE
945 DO bj=myByLo(myThid),myByHi(myThid)
946 DO bi=myBxLo(myThid),myBxHi(myThid)
947
948 DO J=1,sNy
949 DO I=1,sNx
950 AA1=0.5 _d 0 * _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj)
951 & * (etaPlusZeta(I-1,J ,bi,bj) + etaPlusZeta(I,J ,bi,bj))
952 AA2=0.5 _d 0 * _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj)
953 & * (etaPlusZeta(I-1,J-1,bi,bj) + etaPlusZeta(I,J-1,bi,bj))
954 AA3= (ETA(I ,J-1,bi,bj) * _recip_dxV(I,J,bi,bj)
955 & +ETA(I ,J ,bi,bj) * _recip_dxV(I,J,bi,bj)
956 & )* 0.5 _d 0 * _recip_dxV(I,J,bi,bj)
957 AA4= (ETA(I-1,J-1,bi,bj)+ETA(I-1,J,bi,bj))*0.5 _d 0
958 & *_recip_dxV(I,J,bi,bj) * _recip_dxV(I,J,bi,bj)
959 AA5=(zetaMinusEta(I-1,J ,bi,bj) + zetaMinusEta(I,J ,bi,bj)
960 & -zetaMinusEta(I-1,J-1,bi,bj) - zetaMinusEta(I,J-1,bi,bj)
961 & )* _tanPhiAtV(I,J,bi,bj)
962 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
963
964 AA6=TWO*ETAMEAN(I,J,bi,bj) * recip_rSphere*recip_rSphere
965 & * _tanPhiAtV(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj)
966
967 AV(I,J,bi,bj)=(
968 & - AA2
969 & - (ZETAMEAN(I,J,bi,bj)-ETAMEAN(I,J,bi,bj))
970 & * _tanPhiAtV(I,J-1,bi,bj)
971 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
972 & -ETAMEAN(I,J,bi,bj)*TWO* _tanPhiAtV(I,J,bi,bj)
973 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
974 & )*UVM(I,J,bi,bj)
975 CV(I,J,bi,bj)=(
976 & -AA1
977 & +(ZETAMEAN(I,J,bi,bj)-ETAMEAN(I,J,bi,bj))
978 & * _tanPhiAtV(I,J+1,bi,bj)
979 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
980 & +ETAMEAN(I,J,bi,bj)*TWO* _tanPhiAtV(I,J,bi,bj)
981 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
982 & )*UVM(I,J,bi,bj)
983 BV(I,J,bi,bj)= (ONE-UVM(I,J,bi,bj))
984 & +( (AA1+AA2) + (AA3+AA4) + AA5 + AA6
985 & +AMASS(I,J,bi,bj)/SEAICE_deltaTdyn+DRAGS(I,J,bi,bj))
986 & *UVM(I,J,bi,bj)
987 END DO
988 END DO
989
990 DO I=1,sNx
991 AV(I,1,bi,bj)=ZERO
992 CV(I,sNy,bi,bj)=ZERO
993 CV(I,1,bi,bj)=CV(I,1,bi,bj)/BV(I,1,bi,bj)
994 END DO
995
996 C now set up right-hand-side
997 DO J=1-OLy,sNy+OLy-1
998 DO I=1-OLx,sNx+OLx-1
999 dUdx(I,J) = 0.5 _d 0 * (
1000 & ( UICEC(I+1,J+1,bi,bj) - UICEC(I ,J+1,bi,bj) )
1001 & * _recip_dxG(I,J+1,bi,bj)
1002 & +(UICEC(I+1,J ,bi,bj) - UICEC(I ,J ,bi,bj) )
1003 & * _recip_dxG(I,J ,bi,bj) )
1004 dUdy(I,J) = 0.5 _d 0 * (
1005 & ( UICEC(I+1,J+1,bi,bj) - UICEC(I+1,J ,bi,bj) )
1006 & * _recip_dyG(I+1,J,bi,bj)
1007 & +(UICEC(I ,J+1,bi,bj) - UICEC(I ,J ,bi,bj) )
1008 & * _recip_dyG(I, J,bi,bj) )
1009 ENDDO
1010 ENDDO
1011 #ifdef SEAICE_TEST
1012 DO j=1-OLy,sNy+OLy-1
1013 DO i=1-OLx,sNx+OLx-1
1014 uz(i,j) = quart * (
1015 & uicec(i,j,bi,bj) + uicec(i+1,j,bi,bj) )
1016 uz(i,j)= uz(i,j) + quart * (
1017 & uicec(i,j+1,bi,bj) + uicec(i+1,j+1,bi,bj) )
1018 ENDDO
1019 ENDDO
1020 #endif
1021 DO J=1,sNy
1022 DO I=1,sNx
1023 rhsV(I,J,bi,bj)=-DRAGA(I,J,bi,bj)*UICEC(I,J,bi,bj)
1024 & +FORCEY(I,J,bi,bj)
1025 &
1026 #ifdef SEAICE_TEST
1027 & + ( 0.5 _d 0 *
1028 & (zetaMinusEta(i,j,bi,bj)+zetaMinusEta(i-1,j,bi,bj))
1029 & *(uz(i,j)-uz(i-1,j)) * _recip_dxC(i,j,bi,bj)
1030 & - 0.5 _d 0 *
1031 & (zetaMinusEta(i,j-1,bi,bj)+zetaMinusEta(i-1,j-1,bi,bj))
1032 & *(uz(i,j-1)-uz(i-1,j-1)) * _recip_dxC(i,j-1,bi,bj)
1033 & ) * _recip_dyU(i,j,bi,bj)
1034 #else
1035 & + ( zetaMinusEta(I ,J ,bi,bj) * dUdx(I ,J )
1036 & + zetaMinusEta(I-1,J ,bi,bj) * dUdx(I-1,J )
1037 & - zetaMinusEta(I ,J-1,bi,bj) * dUdx(I ,J-1)
1038 & - zetaMinusEta(I-1,J-1,bi,bj) * dUdx(I-1,J-1)
1039 & )* 0.5 _d 0 * _recip_dyU(I,J,bi,bj)
1040 #endif
1041 &
1042 & + ( ETA (I ,J ,bi,bj) * dUdy(I ,J )
1043 & + ETA (I ,J-1,bi,bj) * dUdy(I ,J-1)
1044 & - ETA (I-1,J ,bi,bj) * dUdy(I-1,J )
1045 & - ETA (I-1,J-1,bi,bj) * dUdy(I-1,J-1)
1046 & )*0.5 _d 0* _recip_dxV(I,J,bi,bj)
1047 &
1048 & +(ETA(I ,J ,bi,bj) + ETA(I ,J-1,bi,bj)
1049 & -ETA(I-1,J-1,bi,bj) - ETA(I-1,J ,bi,bj))
1050 & * UICEC(I,J,bi,bj)
1051 & * _tanPhiAtV(I,J,bi,bj)
1052 & * 0.5 _d 0 * _recip_dxV(I,J,bi,bj)*recip_rSphere
1053 & +ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj)
1054 & *(UICEC(I+1,J,bi,bj)-UICEC(I-1,J,bi,bj))
1055 & * 0.5 _d 0 * _recip_dxV(I,J,bi,bj)*recip_rSphere
1056 &
1057 & +ETAMEAN(I,J,bi,bj)*TWO * _tanPhiAtV(I,J,bi,bj)
1058 & *(UICEC(I+1,J,bi,bj)-UICEC(I-1,J,bi,bj))
1059 & * 1. _d 0 /( _dxG(I,J,bi,bj) + _dxG(I-1,J,bi,bj))
1060 & *recip_rSphere
1061 VRT1(I,J,bi,bj)= 0.5 _d 0 * (
1062 & ETA(I-1,J-1,bi,bj) * _recip_dxV(I,J,bi,bj)
1063 & +ETA(I-1,J ,bi,bj) * _recip_dxV(I,J,bi,bj)
1064 & ) * _recip_dxV(I,J,bi,bj)
1065 VRT2(I,J,bi,bj)= 0.5 _d 0 * (
1066 & ETA(I ,J-1,bi,bj) * _recip_dxV(I,J,bi,bj)
1067 & +ETA(I ,J ,bi,bj) * _recip_dxV(I,J,bi,bj)
1068 & ) * _recip_dxV(I,J,bi,bj)
1069
1070 END DO
1071 END DO
1072
1073 ENDDO
1074 ENDDO
1075
1076 C NOW DO ITERATION
1077
1078 cph--- iteration starts here
1079 cph--- need to kick out goto
1080 phexit = -1
1081
1082 C ITERATION START -----------------------------------------------------
1083 #ifdef ALLOW_AUTODIFF_TAMC
1084 CADJ LOOP = iteration vice
1085 #endif /* ALLOW_AUTODIFF_TAMC */
1086 DO M=1, solv_max_iters
1087 IF ( phexit .EQ. -1 ) THEN
1088
1089 C NOW SET U(3)=U(1)
1090 DO bj=myByLo(myThid),myByHi(myThid)
1091 DO bi=myBxLo(myThid),myBxHi(myThid)
1092
1093 DO J=1,sNy
1094 DO I=1,sNx
1095 VTMP(I,J,bi,bj)=VICE(I,J,bi,bj)
1096 END DO
1097 END DO
1098
1099 DO I=1,sNx
1100 DO J=1,sNy
1101 IF(J.EQ.1) THEN
1102 AA2= _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj)
1103 & * 0.5 _d 0 *(
1104 & etaPlusZeta(I-1,J-1,bi,bj) + etaPlusZeta(I,J-1,bi,bj)
1105 & )
1106 AA3=( AA2
1107 & +( ZETAMEAN(I,J,bi,bj)-ETAMEAN(I,J,bi,bj) )
1108 & * _tanPhiAtV(I,J-1,bi,bj)
1109 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
1110 & + ETAMEAN(I,J,bi,bj)*TWO* _tanPhiAtV(I,J,bi,bj)
1111 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere )
1112 & *VICE(I,J-1,bi,bj)*UVM(I,J,bi,bj)
1113 ELSE IF(J.EQ.sNy) THEN
1114 AA1= _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj)
1115 & * 0.5 _d 0 * (
1116 & etaPlusZeta(I-1,J,bi,bj) + etaPlusZeta(I,J,bi,bj)
1117 & )
1118 AA3=( AA1
1119 & -( ZETAMEAN(I,J,bi,bj)-ETAMEAN(I,J,bi,bj))
1120 & * _tanPhiAtV(I,J+1,bi,bj)
1121 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
1122 & - ETAMEAN(I,J,bi,bj)*TWO* _tanPhiAtV(I,J,bi,bj)
1123 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere )
1124 & *VICE(I,J+1,bi,bj)*UVM(I,J,bi,bj)
1125 ELSE
1126 AA3=ZERO
1127 END IF
1128
1129 VRT(J)=rhsV(I,J,bi,bj)+AA3+VRT1(I,J,bi,bj)*VICE(I-1,J,bi,bj)
1130 & +VRT2(I,J,bi,bj)*VICE(I+1,J,bi,bj)
1131 VRT(J)=VRT(J)*UVM(I,J,bi,bj)
1132 END DO
1133
1134 DO J=1,sNy
1135 CVV(J)=CV(I,J,bi,bj)
1136 END DO
1137 VRT(1)=VRT(1)/BV(I,1,bi,bj)
1138 DO J=2,sNy
1139 JM=J-1
1140 CVV(J)=CVV(J)/(BV(I,J,bi,bj)-AV(I,J,bi,bj)*CVV(JM))
1141 VRT(J)=(VRT(J)-AV(I,J,bi,bj)*VRT(JM))
1142 & /(BV(I,J,bi,bj)-AV(I,J,bi,bj)*CVV(JM))
1143 END DO
1144 DO J=1,sNy-1
1145 J1=sNy-J
1146 J2=J1+1
1147 VRT(J1)=VRT(J1)-CVV(J1)*VRT(J2)
1148 END DO
1149 DO J=1,sNy
1150 VICE(I,J,bi,bj)=VTMP(I,J,bi,bj)
1151 & +WFAV*(VRT(J)-VTMP(I,J,bi,bj))
1152 END DO
1153 ENDDO
1154
1155 ENDDO
1156 ENDDO
1157
1158 IF(MOD(M,SOLV_NCHECK).EQ.0) THEN
1159 S2=ZERO
1160 DO bj=myByLo(myThid),myByHi(myThid)
1161 DO bi=myBxLo(myThid),myBxHi(myThid)
1162 DO J=1,sNy
1163 DO I=1,sNx
1164 UERR=(VICE(I,J,bi,bj)-VTMP(I,J,bi,bj))
1165 & *UVM(I,J,bi,bj)
1166 S2=MAX(ABS(UERR),S2)
1167 END DO
1168 END DO
1169 ENDDO
1170 ENDDO
1171 _GLOBAL_MAX_RL( S2, myThid )
1172 C SAFEGUARD AGAINST BAD FORCING ETC
1173 IF(M.GT.1.AND.S2.GT.S2A) WFAV=WFAV2
1174 S2A=S2
1175 IF(S2.LT.LSR_ERROR) THEN
1176 ICOUNT2=M
1177 phexit = 1
1178 END IF
1179 END IF
1180
1181 _EXCH_XY_RL( VICE, myThid )
1182
1183 ENDIF
1184 ENDDO
1185 C ITERATION END -----------------------------------------------------
1186
1187 IF ( debugLevel .GE. debLevC ) THEN
1188 _BEGIN_MASTER( myThid )
1189 write(*,'(A,I6,1P2E22.14)')' V lsr iters, error = ',ICOUNT2,S2
1190 _END_MASTER( myThid )
1191 ENDIF
1192
1193 #endif /* SEAICE_LSRBNEW */
1194
1195 C NOW END
1196 C NOW MAKE COROLIS TERM IMPLICIT
1197 DO bj=myByLo(myThid),myByHi(myThid)
1198 DO bi=myBxLo(myThid),myBxHi(myThid)
1199 DO J=1-OLy,sNy+OLy
1200 DO I=1-OLx,sNx+OLx
1201 UICE(I,J,bi,bj)=UICE(I,J,bi,bj)*UVM(I,J,bi,bj)
1202 VICE(I,J,bi,bj)=VICE(I,J,bi,bj)*UVM(I,J,bi,bj)
1203 END DO
1204 END DO
1205 ENDDO
1206 ENDDO
1207
1208 #endif /* SEAICE_ALLOW_DYNAMICS */
1209 #endif /* SEAICE_CGRID */
1210
1211 RETURN
1212 END
  ViewVC Help
Powered by ViewVC 1.1.22