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

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

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

Revision 1.29 - (hide annotations) (download)
Mon Jun 6 14:55:08 2011 UTC (13 years, 3 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62z, checkpoint63g, checkpoint63, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63a, checkpoint63b, checkpoint63c
Changes since 1.28: +3 -3 lines 
refine debugLevel criteria when printing messages
1 jmc 1.29 C $Header: /u/gcmpack/MITgcm/pkg/seaice/lsr.F,v 1.28 2010/03/16 00:23:18 jmc Exp $
2 edhill 1.8 C $Name: $
3 heimbach 1.2
4 mlosch 1.21 C for an alternative discretization of d/dx[ (zeta-eta) dV/dy]
5     C and d/dy[ (zeta-eta) dU/dx] uncomment this option
6     C#define SEAICE_TEST
7    
8 heimbach 1.2 #include "SEAICE_OPTIONS.h"
9    
10     CStartOfInterface
11 heimbach 1.7 SUBROUTINE lsr( ilcall, myThid )
12 heimbach 1.2 C /==========================================================\
13     C | SUBROUTINE lsr |
14 dimitri 1.6 C | o Solve ice momentum equation with an LSR dynamics solver|
15     C | (see Zhang and Hibler, JGR, 102, 8691-8702, 1997 |
16     C | and Zhang and Rothrock, MWR, 131, 845- 861, 2003) |
17     C | Written by Jinlun Zhang, PSC/UW, Feb-2001 |
18     C | zhang@apl.washington.edu |
19 heimbach 1.2 C |==========================================================|
20     C \==========================================================/
21     IMPLICIT NONE
22    
23     C === Global variables ===
24     #include "SIZE.h"
25     #include "EEPARAMS.h"
26     #include "PARAMS.h"
27 mlosch 1.16 #include "GRID.h"
28 heimbach 1.2 #include "SEAICE.h"
29     #include "SEAICE_PARAMS.h"
30 mlosch 1.16 C#include "SEAICE_GRID.h"
31 heimbach 1.2
32 heimbach 1.7 #ifdef ALLOW_AUTODIFF_TAMC
33     # include "tamc.h"
34     #endif
35 heimbach 1.2
36     C === Routine arguments ===
37     C myThid - Thread no. that called this routine.
38 heimbach 1.7 INTEGER ilcall
39 heimbach 1.2 INTEGER myThid
40     CEndOfInterface
41    
42 mlosch 1.18 #ifndef SEAICE_CGRID
43 dimitri 1.5 #ifdef SEAICE_ALLOW_DYNAMICS
44 heimbach 1.2
45     C === Local variables ===
46 dimitri 1.4 C i,j,bi,bj - Loop counters
47 heimbach 1.2
48 dimitri 1.6 INTEGER i, j, m, bi, bj, j1, j2, im, jm
49 mlosch 1.23 INTEGER ICOUNT1, ICOUNT2
50 heimbach 1.7 INTEGER phexit
51    
52 dimitri 1.6 _RL WFAU, WFAV, WFAU1, WFAV1, WFAU2, WFAV2
53     _RL AA1, AA2, AA3, AA4, AA5, AA6, S1, S2, S1A, S2A
54    
55 mlosch 1.16 _RL AU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
56     _RL BU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
57     _RL CU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
58     _RL AV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
59     _RL BV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
60     _RL CV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
61     _RL UERR (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
62     _RL FXY (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
63    
64     _RL URT(1-Olx:sNx+Olx), CUU(1-Olx:sNx+Olx)
65     _RL VRT(1-Oly:sNy+Oly), CVV(1-Oly:sNy+Oly)
66    
67     _RL etaPlusZeta (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
68     _RL zetaMinusEta(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
69     _RL ETAMEAN (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
70     _RL ZETAMEAN (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
71    
72     _RL UVRT1 (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
73     _RL UVRT2 (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
74    
75 mlosch 1.26 _RL UTMP (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
76 mlosch 1.27 _RL VTMP (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
77 mlosch 1.26
78 mlosch 1.21 _RL dVdx (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
79     _RL dVdy (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
80     _RL dUdx (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
81     _RL dUdy (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
82     #ifdef SEAICE_TEST
83 mlosch 1.26 _RL uz (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
84     _RL vz (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
85 mlosch 1.21 #endif
86 dimitri 1.6
87 heimbach 1.2 C SET SOME VALUES
88 dimitri 1.6 WFAU1=0.95 _d 0
89     WFAV1=0.95 _d 0
90     WFAU2=ZERO
91     WFAV2=ZERO
92    
93     S1A=0.80 _d 0
94     S2A=0.80 _d 0
95     WFAU=WFAU1
96     WFAV=WFAV1
97    
98     ICOUNT1=SOLV_MAX_ITERS
99     ICOUNT2=SOLV_MAX_ITERS
100 heimbach 1.2
101     C SOLVE FOR UICE
102    
103 heimbach 1.7 #ifdef ALLOW_AUTODIFF_TAMC
104 jmc 1.28 cph That is an important one! Note, that
105 heimbach 1.7 cph * lsr is called twice, thus the icall index
106     cph * this storing is still outside the iteration loop
107 heimbach 1.25 CADJ STORE uice = comlev1_dynsol,
108 heimbach 1.7 CADJ & key = ikey_dynamics + (ilcall-1)*nchklev_1
109 heimbach 1.25 CADJ STORE vice = comlev1_dynsol,
110 heimbach 1.7 CADJ & key = ikey_dynamics + (ilcall-1)*nchklev_1
111     #endif /* ALLOW_AUTODIFF_TAMC */
112    
113 heimbach 1.2 DO bj=myByLo(myThid),myByHi(myThid)
114     DO bi=myBxLo(myThid),myBxHi(myThid)
115 mlosch 1.16 DO j=1-Oly,sNy+Oly
116     DO i=1-Olx,sNx+Olx
117 heimbach 1.2 FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
118 mlosch 1.26 & +AMASS(I,J,bi,bj)/SEAICE_deltaTdyn*UICENM1(I,J,bi,bj)
119 heimbach 1.2 FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)
120 mlosch 1.26 & +AMASS(I,J,bi,bj)/SEAICE_deltaTdyn*VICENM1(I,J,bi,bj)
121 heimbach 1.2 FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)*UVM(I,J,bi,bj)
122     FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)*UVM(I,J,bi,bj)
123 mlosch 1.16 etaPlusZeta(I,J,bi,bj) = ETA(I,J,bi,bj)+ZETA(I,J,bi,bj)
124     zetaMinusEta(I,J,bi,bj) = ZETA(I,J,bi,bj)-ETA(I,J,bi,bj)
125     ENDDO
126     ENDDO
127     DO j=1-Oly+1,sNy+Oly
128     DO i=1-Olx+1,sNx+Olx
129     ETAMEAN(I,J,bi,bj) =QUART*(
130     & ETA(I,J-1,bi,bj) + ETA(I-1,J-1,bi,bj)
131     & +ETA(I,J ,bi,bj) + ETA(I-1,J ,bi,bj))
132     ZETAMEAN(I,J,bi,bj)=QUART*(
133     & ZETA(I,J-1,bi,bj) + ZETA(I-1,J-1,bi,bj)
134     & +ZETA(I,J ,bi,bj) + ZETA(I-1,J ,bi,bj))
135 heimbach 1.2 ENDDO
136     ENDDO
137     ENDDO
138     ENDDO
139    
140     DO bj=myByLo(myThid),myByHi(myThid)
141     DO bi=myBxLo(myThid),myBxHi(myThid)
142    
143     DO J=1,sNy
144     DO I=1,sNx
145 mlosch 1.16 AA1=( etaPlusZeta(I ,J-1,bi,bj) * _recip_dxF(I ,J-1,bi,bj)
146     & +etaPlusZeta(I ,J ,bi,bj) * _recip_dxF(I ,J ,bi,bj)
147     & )*0.5 _d 0 * _recip_dxV(I,J,bi,bj) * UVM(I,J,bi,bj)
148     AA2=( etaPlusZeta(I-1,J-1,bi,bj) * _recip_dxF(I-1,J-1,bi,bj)
149     & +etaPlusZeta(I-1,J ,bi,bj) * _recip_dxF(I-1,J ,bi,bj)
150     & )*0.5 _d 0 * _recip_dxV(I,J,bi,bj) * UVM(I,J,bi,bj)
151     AA3= 0.5 _d 0 *(ETA(I-1,J ,bi,bj)+ETA(I,J ,bi,bj))
152     AA4= 0.5 _d 0 *(ETA(I-1,J-1,bi,bj)+ETA(I,J-1,bi,bj))
153 jmc 1.22 AA5= -(AA3-AA4) * _tanPhiAtV(I,J,bi,bj)
154 mlosch 1.17 & * _recip_dyU(I,J,bi,bj)*recip_rSphere
155     AA6=TWO*ETAMEAN(I,J,bi,bj) *recip_rSphere*recip_rSphere
156 mlosch 1.16 & * _tanPhiAtV(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj)
157     AU(I,J,bi,bj)=-AA2
158     CU(I,J,bi,bj)=-AA1
159     BU(I,J,bi,bj)=(ONE-UVM(I,J,bi,bj))
160     & - AU(I,J,bi,bj) - CU(I,J,bi,bj)
161     & + ((AA3+AA4)*_recip_dyU(I,J,bi,bj)*_recip_dyU(I,J,bi,bj)
162     & + AA5 + AA6
163     & + AMASS(I,J,bi,bj)/SEAICE_deltaTdyn
164     & + DRAGS(I,J,bi,bj)
165     & )*UVM(I,J,bi,bj)
166 heimbach 1.2 END DO
167     END DO
168    
169     DO J=1,sNy
170 dimitri 1.6 AU(1,J,bi,bj)=ZERO
171     CU(sNx,J,bi,bj)=ZERO
172     CU(1,J,bi,bj)=CU(1,J,bi,bj)/BU(1,J,bi,bj)
173 heimbach 1.2 END DO
174    
175 mlosch 1.16 C now set up right-hand side
176     DO J=1-Oly,sNy+Oly-1
177     DO I=1-Olx,sNx+Olx-1
178 mlosch 1.21 dVdy(I,J) = 0.5 _d 0 * (
179     & ( VICEC(I+1,J+1,bi,bj) - VICEC(I+1,J ,bi,bj) )
180     & * _recip_dyG(I+1,J,bi,bj)
181     & +(VICEC(I ,J+1,bi,bj) - VICEC(I ,J ,bi,bj) )
182     & * _recip_dyG(I, J,bi,bj) )
183     dVdx(I,J) = 0.5 _d 0 * (
184     & ( VICEC(I+1,J+1,bi,bj) - VICEC(I ,J+1,bi,bj) )
185     & * _recip_dxG(I,J+1,bi,bj)
186     & +(VICEC(I+1,J ,bi,bj) - VICEC(I ,J ,bi,bj) )
187     & * _recip_dxG(I,J, bi,bj) )
188     ENDDO
189     ENDDO
190     #ifdef SEAICE_TEST
191     DO j=1-Oly,sNy+Oly-1
192     DO i=1-Olx,sNx+Olx-1
193 mlosch 1.26 vz(i,j) = quart * (
194 mlosch 1.21 & vicec(i,j,bi,bj) + vicec(i+1,j,bi,bj) )
195 mlosch 1.26 vz(i,j)= vz(i,j) + quart * (
196 mlosch 1.21 & vicec(i,j+1,bi,bj) + vicec(i+1,j+1,bi,bj) )
197 mlosch 1.16 ENDDO
198     ENDDO
199 mlosch 1.21 #endif
200 heimbach 1.2 DO J=1,sNy
201     DO I=1,sNx
202 dimitri 1.6 FXY(I,J,bi,bj)=DRAGA(I,J,bi,bj)*VICEC(I,J,bi,bj)
203 mlosch 1.16 & +FORCEX(I,J,bi,bj)
204 mlosch 1.21 #ifdef SEAICE_TEST
205 jmc 1.22 & + ( 0.5 _d 0 *
206 mlosch 1.21 & (zetaMinusEta(i,j,bi,bj)+zetaMinusEta(i,j-1,bi,bj))
207 mlosch 1.26 & *(vz(i,j)-vz(i,j-1)) * _recip_dyC(i,j,bi,bj)
208 jmc 1.22 & - 0.5 _d 0 *
209 mlosch 1.21 & (zetaMinusEta(i-1,j,bi,bj)+zetaMinusEta(i-1,j-1,bi,bj))
210 mlosch 1.26 & *(vz(i-1,j)-vz(i-1,j-1)) * _recip_dyC(i-1,j,bi,bj)
211 mlosch 1.21 & ) * _recip_dxV(i,j,bi,bj)
212     #else
213     & + ( zetaMinusEta(I ,J ,bi,bj) * dVdy(I ,J )
214     & + zetaMinusEta(I ,J-1,bi,bj) * dVdy(I ,J-1)
215     & - zetaMinusEta(I-1,J ,bi,bj) * dVdy(I-1,J )
216     & - zetaMinusEta(I-1,J-1,bi,bj) * dVdy(I-1,J-1)
217 jmc 1.22 & )* 0.5 _d 0 * _recip_dxV(I,J,bi,bj)
218 mlosch 1.21 #endif
219 mlosch 1.16 &
220 mlosch 1.21 & + ( ETA (I ,J ,bi,bj) * dVdx(I ,J )
221     & + ETA (I-1,J ,bi,bj) * dVdx(I-1,J )
222     & - ETA (I ,J-1,bi,bj) * dVdx(I ,J-1)
223     & - ETA (I-1,J-1,bi,bj) * dVdx(I-1,J-1)
224 mlosch 1.16 & ) * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)
225     &
226     & -(etaPlusZeta(I ,J ,bi,bj)+etaPlusZeta(I ,J-1,bi,bj)
227     & -etaPlusZeta(I-1,J-1,bi,bj)-etaPlusZeta(I-1,J ,bi,bj))
228     & * VICEC(I,J,bi,bj)
229     & * _tanPhiAtV(I,J,bi,bj)
230 mlosch 1.17 & * 0.5 _d 0 * _recip_dxV(I,J,bi,bj)*recip_rSphere
231 mlosch 1.16 &
232     & -(ETAMEAN(I,J,bi,bj)+ZETAMEAN(I,J,bi,bj))
233     & *(VICEC(I+1,J,bi,bj) - VICEC(I-1,J,bi,bj))
234     & * _tanPhiAtV(I,J,bi,bj)
235     & * 1.0 _d 0 /( _dxG(I,J,bi,bj) + _dxG(I-1,J,bi,bj) )
236 mlosch 1.17 & *recip_rSphere
237 mlosch 1.16 &
238     & -ETAMEAN(I,J,bi,bj)
239     & *(VICEC(I+1,J,bi,bj) - VICEC(I-1,J,bi,bj))
240     & *TWO* _tanPhiAtV(I,J,bi,bj)
241     & * 1.0 _d 0 /( _dxG(I,J,bi,bj) + _dxG(I-1,J,bi,bj) )
242 mlosch 1.17 & *recip_rSphere
243 mlosch 1.16
244     UVRT1(I,J,bi,bj)=
245     & 0.5 _d 0 * (ETA(I-1,J-1,bi,bj)+ETA(I,J-1,bi,bj))
246     & * _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj)
247     & - ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J-1,bi,bj)
248 mlosch 1.17 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
249 mlosch 1.16 & + TWO*ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj)
250 mlosch 1.17 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
251 mlosch 1.16 UVRT2(I,J,bi,bj)=
252     & 0.5 _d 0 * (ETA(I-1,J,bi,bj)+ETA(I,J,bi,bj))
253     & * _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj)
254     & + ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J+1,bi,bj)
255 mlosch 1.17 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
256 mlosch 1.16 & - TWO*ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj)
257 mlosch 1.17 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
258 heimbach 1.2 END DO
259     END DO
260    
261     ENDDO
262     ENDDO
263    
264     C NOW DO ITERATION
265 dimitri 1.6
266 heimbach 1.7 cph--- iteration starts here
267     cph--- need to kick out goto
268     phexit = -1
269    
270     C ITERATION START -----------------------------------------------------
271     #ifdef ALLOW_AUTODIFF_TAMC
272     CADJ LOOP = iteration uice
273     #endif /* ALLOW_AUTODIFF_TAMC */
274 jmc 1.22 DO M=1, solv_max_iters
275     IF ( phexit .EQ. -1 ) THEN
276 heimbach 1.7
277 heimbach 1.2 DO bj=myByLo(myThid),myByHi(myThid)
278     DO bi=myBxLo(myThid),myBxHi(myThid)
279     C NOW SET U(3)=U(1)
280     DO J=1,sNy
281     DO I=1,sNx
282 mlosch 1.27 UTMP(I,J,bi,bj)=UICE(I,J,bi,bj)
283 heimbach 1.2 END DO
284     END DO
285    
286 jmc 1.22 DO J=1,sNy
287 heimbach 1.2 DO I=1,sNx
288 dimitri 1.6 IF(I.EQ.1) THEN
289 mlosch 1.16 AA2=(etaPlusZeta(I-1,J-1,bi,bj) * _recip_dxF(I-1,J-1,bi,bj)
290     & +etaPlusZeta(I-1,J ,bi,bj) * _recip_dxF(I-1,J ,bi,bj)
291     & )*0.5 _d 0 * _recip_dxV(I,J,bi,bj)
292 mlosch 1.27 AA3=AA2*UICE(I-1,J,bi,bj)*UVM(I,J,bi,bj)
293 dimitri 1.6 ELSE IF(I.EQ.sNx) THEN
294 mlosch 1.16 AA1=(etaPlusZeta(I ,J-1,bi,bj) * _recip_dxF(I ,J-1,bi,bj)
295     & +etaPlusZeta(I ,J ,bi,bj) * _recip_dxF(I ,J ,bi,bj)
296     & )*0.5 _d 0 * _recip_dxV(I,J,bi,bj)
297 mlosch 1.27 AA3=AA1*UICE(I+1,J,bi,bj)*UVM(I,J,bi,bj)
298 dimitri 1.6 ELSE
299     AA3=ZERO
300     END IF
301     URT(I)=FXY(I,J,bi,bj)+AA3
302 mlosch 1.27 & +UVRT1(I,J,bi,bj)*UICE(I,J-1,bi,bj)
303     & +UVRT2(I,J,bi,bj)*UICE(I,J+1,bi,bj)
304 heimbach 1.2 URT(I)=URT(I)*UVM(I,J,bi,bj)
305     END DO
306    
307     DO I=1,sNx
308 dimitri 1.6 CUU(I)=CU(I,J,bi,bj)
309 heimbach 1.2 END DO
310 dimitri 1.6 URT(1)=URT(1)/BU(1,J,bi,bj)
311 heimbach 1.2 DO I=2,sNx
312     IM=I-1
313 dimitri 1.6 CUU(I)=CUU(I)/(BU(I,J,bi,bj)-AU(I,J,bi,bj)*CUU(IM))
314     URT(I)=(URT(I)-AU(I,J,bi,bj)*URT(IM))
315     & /(BU(I,J,bi,bj)-AU(I,J,bi,bj)*CUU(IM))
316 heimbach 1.2 END DO
317     DO I=1,sNx-1
318     J1=sNx-I
319     J2=J1+1
320     URT(J1)=URT(J1)-CUU(J1)*URT(J2)
321     END DO
322     DO I=1,sNx
323 mlosch 1.27 UICE(I,J,bi,bj)=UTMP(I,J,bi,bj)
324 mlosch 1.26 & +WFAU*(URT(I)-UTMP(I,J,bi,bj))
325 heimbach 1.2 END DO
326    
327 jmc 1.22 END DO
328 heimbach 1.2
329     ENDDO
330     ENDDO
331    
332 dimitri 1.6 IF(MOD(M,SOLV_NCHECK).EQ.0) THEN
333 jmc 1.22 S1=ZERO
334 dimitri 1.6 DO bj=myByLo(myThid),myByHi(myThid)
335     DO bi=myBxLo(myThid),myBxHi(myThid)
336     DO J=1,sNy
337     DO I=1,sNx
338 mlosch 1.27 UERR(I,J,bi,bj)=(UICE(I,J,bi,bj)-UTMP(I,J,bi,bj))
339 heimbach 1.2 & *UVM(I,J,bi,bj)
340 dimitri 1.6 S1=MAX(ABS(UERR(I,J,bi,bj)),S1)
341     END DO
342 heimbach 1.2 END DO
343 dimitri 1.6 ENDDO
344 heimbach 1.2 ENDDO
345 jmc 1.24 _GLOBAL_MAX_RL( S1, myThid )
346 dimitri 1.6 C SAFEGUARD AGAINST BAD FORCING ETC
347     IF(M.GT.1.AND.S1.GT.S1A) WFAU=WFAU2
348     S1A=S1
349     IF(S1.LT.LSR_ERROR) THEN
350     ICOUNT1=M
351 heimbach 1.7 phexit = 1
352 dimitri 1.6 END IF
353     END IF
354 mlosch 1.27 _EXCH_XY_RL( UICE, myThid )
355 dimitri 1.6
356 jmc 1.22 ENDIF
357     ENDDO
358 heimbach 1.7 C ITERATION END -----------------------------------------------------
359    
360 jmc 1.29 IF ( debugLevel .GE. debLevC ) THEN
361 dimitri 1.11 _BEGIN_MASTER( myThid )
362     write(*,'(A,I6,1P2E22.14)')' U lsr iters, error = ',ICOUNT1,S1
363 edhill 1.14 _END_MASTER( myThid )
364 dimitri 1.11 ENDIF
365 heimbach 1.2
366     C NOW FOR VICE
367     DO bj=myByLo(myThid),myByHi(myThid)
368     DO bi=myBxLo(myThid),myBxHi(myThid)
369    
370     DO J=1,sNy
371     DO I=1,sNx
372 mlosch 1.16 AA1=0.5 _d 0 * _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj)
373     & * (etaPlusZeta(I-1,J ,bi,bj) + etaPlusZeta(I,J ,bi,bj))
374     AA2=0.5 _d 0 * _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj)
375     & * (etaPlusZeta(I-1,J-1,bi,bj) + etaPlusZeta(I,J-1,bi,bj))
376     AA3= (ETA(I ,J-1,bi,bj) * _recip_dxV(I,J,bi,bj)
377     & +ETA(I ,J ,bi,bj) * _recip_dxV(I,J,bi,bj)
378     & )* 0.5 _d 0 * _recip_dxV(I,J,bi,bj)
379     AA4= (ETA(I-1,J-1,bi,bj)+ETA(I-1,J,bi,bj))*0.5 _d 0
380     & *_recip_dxV(I,J,bi,bj) * _recip_dxV(I,J,bi,bj)
381     AA5=(zetaMinusEta(I-1,J ,bi,bj) + zetaMinusEta(I,J ,bi,bj)
382     & -zetaMinusEta(I-1,J-1,bi,bj) - zetaMinusEta(I,J-1,bi,bj)
383 jmc 1.22 & )* _tanPhiAtV(I,J,bi,bj)
384 mlosch 1.17 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
385 mlosch 1.16
386 mlosch 1.17 AA6=TWO*ETAMEAN(I,J,bi,bj) * recip_rSphere*recip_rSphere
387 mlosch 1.16 & * _tanPhiAtV(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj)
388    
389     AV(I,J,bi,bj)=(
390     & - AA2
391     & - (ZETAMEAN(I,J,bi,bj)-ETAMEAN(I,J,bi,bj))
392 jmc 1.22 & * _tanPhiAtV(I,J-1,bi,bj)
393 mlosch 1.17 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
394 mlosch 1.16 & -ETAMEAN(I,J,bi,bj)*TWO* _tanPhiAtV(I,J,bi,bj)
395 mlosch 1.17 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
396 mlosch 1.16 & )*UVM(I,J,bi,bj)
397     CV(I,J,bi,bj)=(
398 jmc 1.22 & -AA1
399 mlosch 1.16 & +(ZETAMEAN(I,J,bi,bj)-ETAMEAN(I,J,bi,bj))
400 jmc 1.22 & * _tanPhiAtV(I,J+1,bi,bj)
401 mlosch 1.17 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
402 mlosch 1.16 & +ETAMEAN(I,J,bi,bj)*TWO* _tanPhiAtV(I,J,bi,bj)
403 mlosch 1.17 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
404 mlosch 1.16 & )*UVM(I,J,bi,bj)
405     BV(I,J,bi,bj)= (ONE-UVM(I,J,bi,bj))
406     & +( (AA1+AA2) + (AA3+AA4) + AA5 + AA6
407     & +AMASS(I,J,bi,bj)/SEAICE_deltaTdyn+DRAGS(I,J,bi,bj))
408     & *UVM(I,J,bi,bj)
409 heimbach 1.2 END DO
410     END DO
411    
412     DO I=1,sNx
413 dimitri 1.6 AV(I,1,bi,bj)=ZERO
414     CV(I,sNy,bi,bj)=ZERO
415     CV(I,1,bi,bj)=CV(I,1,bi,bj)/BV(I,1,bi,bj)
416 heimbach 1.2 END DO
417    
418 mlosch 1.16 C now set up right-hand-side
419     DO J=1-Oly,sNy+Oly-1
420     DO I=1-Olx,sNx+Olx-1
421 mlosch 1.21 dUdx(I,J) = 0.5 _d 0 * (
422     & ( UICEC(I+1,J+1,bi,bj) - UICEC(I ,J+1,bi,bj) )
423     & * _recip_dxG(I,J+1,bi,bj)
424     & +(UICEC(I+1,J ,bi,bj) - UICEC(I ,J ,bi,bj) )
425     & * _recip_dxG(I,J ,bi,bj) )
426     dUdy(I,J) = 0.5 _d 0 * (
427     & ( UICEC(I+1,J+1,bi,bj) - UICEC(I+1,J ,bi,bj) )
428     & * _recip_dyG(I+1,J,bi,bj)
429     & +(UICEC(I ,J+1,bi,bj) - UICEC(I ,J ,bi,bj) )
430     & * _recip_dyG(I, J,bi,bj) )
431 mlosch 1.16 ENDDO
432     ENDDO
433 mlosch 1.21 #ifdef SEAICE_TEST
434 mlosch 1.20 DO j=1-Oly,sNy+Oly-1
435     DO i=1-Olx,sNx+Olx-1
436 mlosch 1.26 uz(i,j) = quart * (
437 mlosch 1.21 & uicec(i,j,bi,bj) + uicec(i+1,j,bi,bj) )
438 mlosch 1.26 uz(i,j)= uz(i,j) + quart * (
439 mlosch 1.20 & uicec(i,j+1,bi,bj) + uicec(i+1,j+1,bi,bj) )
440     ENDDO
441     ENDDO
442 mlosch 1.21 #endif
443 heimbach 1.2 DO J=1,sNy
444     DO I=1,sNx
445 dimitri 1.6 FXY(I,J,bi,bj)=-DRAGA(I,J,bi,bj)*UICEC(I,J,bi,bj)
446 mlosch 1.16 & +FORCEY(I,J,bi,bj)
447 mlosch 1.20 &
448 mlosch 1.21 #ifdef SEAICE_TEST
449 jmc 1.22 & + ( 0.5 _d 0 *
450 mlosch 1.20 & (zetaMinusEta(i,j,bi,bj)+zetaMinusEta(i-1,j,bi,bj))
451 mlosch 1.26 & *(uz(i,j)-uz(i-1,j)) * _recip_dxC(i,j,bi,bj)
452 jmc 1.22 & - 0.5 _d 0 *
453 mlosch 1.20 & (zetaMinusEta(i,j-1,bi,bj)+zetaMinusEta(i-1,j-1,bi,bj))
454 mlosch 1.26 & *(uz(i,j-1)-uz(i-1,j-1)) * _recip_dxC(i,j-1,bi,bj)
455 mlosch 1.20 & ) * _recip_dyU(i,j,bi,bj)
456 mlosch 1.21 #else
457     & + ( zetaMinusEta(I ,J ,bi,bj) * dUdx(I ,J )
458     & + zetaMinusEta(I-1,J ,bi,bj) * dUdx(I-1,J )
459     & - zetaMinusEta(I ,J-1,bi,bj) * dUdx(I ,J-1)
460     & - zetaMinusEta(I-1,J-1,bi,bj) * dUdx(I-1,J-1)
461     & )* 0.5 _d 0 * _recip_dyU(I,J,bi,bj)
462     #endif
463 jmc 1.22 &
464 mlosch 1.21 & + ( ETA (I ,J ,bi,bj) * dUdy(I ,J )
465     & + ETA (I ,J-1,bi,bj) * dUdy(I ,J-1)
466     & - ETA (I-1,J ,bi,bj) * dUdy(I-1,J )
467     & - ETA (I-1,J-1,bi,bj) * dUdy(I-1,J-1)
468     & )*0.5 _d 0* _recip_dxV(I,J,bi,bj)
469 mlosch 1.16 &
470     & +(ETA(I ,J ,bi,bj) + ETA(I ,J-1,bi,bj)
471     & -ETA(I-1,J-1,bi,bj) - ETA(I-1,J ,bi,bj))
472     & * UICEC(I,J,bi,bj)
473     & * _tanPhiAtV(I,J,bi,bj)
474 mlosch 1.17 & * 0.5 _d 0 * _recip_dxV(I,J,bi,bj)*recip_rSphere
475 mlosch 1.16 & +ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj)
476     & *(UICEC(I+1,J,bi,bj)-UICEC(I-1,J,bi,bj))
477 mlosch 1.17 & * 0.5 _d 0 * _recip_dxV(I,J,bi,bj)*recip_rSphere
478 mlosch 1.16 &
479     & +ETAMEAN(I,J,bi,bj)*TWO * _tanPhiAtV(I,J,bi,bj)
480     & *(UICEC(I+1,J,bi,bj)-UICEC(I-1,J,bi,bj))
481     & * 1. _d 0 /( _dxG(I,J,bi,bj) + _dxG(I-1,J,bi,bj))
482 mlosch 1.17 & *recip_rSphere
483 mlosch 1.16 UVRT1(I,J,bi,bj)= 0.5 _d 0 * (
484     & ETA(I-1,J-1,bi,bj) * _recip_dxV(I,J,bi,bj)
485     & +ETA(I-1,J ,bi,bj) * _recip_dxV(I,J,bi,bj)
486     & ) * _recip_dxV(I,J,bi,bj)
487     UVRT2(I,J,bi,bj)= 0.5 _d 0 * (
488     & ETA(I ,J-1,bi,bj) * _recip_dxV(I,J,bi,bj)
489     & +ETA(I ,J ,bi,bj) * _recip_dxV(I,J,bi,bj)
490     & ) * _recip_dxV(I,J,bi,bj)
491 heimbach 1.2
492     END DO
493     END DO
494    
495     ENDDO
496     ENDDO
497    
498     C NOW DO ITERATION
499 dimitri 1.6
500 heimbach 1.7 cph--- iteration starts here
501     cph--- need to kick out goto
502     phexit = -1
503    
504     C ITERATION START -----------------------------------------------------
505     #ifdef ALLOW_AUTODIFF_TAMC
506     CADJ LOOP = iteration vice
507     #endif /* ALLOW_AUTODIFF_TAMC */
508 jmc 1.22 DO M=1, solv_max_iters
509     IF ( phexit .EQ. -1 ) THEN
510 heimbach 1.7
511 heimbach 1.2 C NOW SET U(3)=U(1)
512     DO bj=myByLo(myThid),myByHi(myThid)
513     DO bi=myBxLo(myThid),myBxHi(myThid)
514    
515     DO J=1,sNy
516     DO I=1,sNx
517 mlosch 1.27 VTMP(I,J,bi,bj)=VICE(I,J,bi,bj)
518 heimbach 1.2 END DO
519     END DO
520    
521 heimbach 1.7 DO I=1,sNx
522 heimbach 1.2 DO J=1,sNy
523 dimitri 1.6 IF(J.EQ.1) THEN
524 mlosch 1.16 AA2= _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj)
525 jmc 1.22 & * 0.5 _d 0 *(
526 mlosch 1.16 & etaPlusZeta(I-1,J-1,bi,bj) + etaPlusZeta(I,J-1,bi,bj)
527     & )
528     AA3=( AA2
529     & +( ZETAMEAN(I,J,bi,bj)-ETAMEAN(I,J,bi,bj) )
530     & * _tanPhiAtV(I,J-1,bi,bj)
531 jmc 1.22 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
532 mlosch 1.16 & + ETAMEAN(I,J,bi,bj)*TWO* _tanPhiAtV(I,J,bi,bj)
533 mlosch 1.17 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere )
534 mlosch 1.27 & *VICE(I,J-1,bi,bj)*UVM(I,J,bi,bj)
535 dimitri 1.6 ELSE IF(J.EQ.sNy) THEN
536 mlosch 1.16 AA1= _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj)
537     & * 0.5 _d 0 * (
538     & etaPlusZeta(I-1,J,bi,bj) + etaPlusZeta(I,J,bi,bj)
539     & )
540     AA3=( AA1
541     & -( ZETAMEAN(I,J,bi,bj)-ETAMEAN(I,J,bi,bj))
542     & * _tanPhiAtV(I,J+1,bi,bj)
543 jmc 1.22 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
544 mlosch 1.16 & - ETAMEAN(I,J,bi,bj)*TWO* _tanPhiAtV(I,J,bi,bj)
545 mlosch 1.17 & * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere )
546 mlosch 1.27 & *VICE(I,J+1,bi,bj)*UVM(I,J,bi,bj)
547 dimitri 1.6 ELSE
548     AA3=ZERO
549     END IF
550    
551 mlosch 1.27 VRT(J)=FXY(I,J,bi,bj)+AA3+UVRT1(I,J,bi,bj)*VICE(I-1,J,bi,bj)
552     & +UVRT2(I,J,bi,bj)*VICE(I+1,J,bi,bj)
553 heimbach 1.2 VRT(J)=VRT(J)*UVM(I,J,bi,bj)
554     END DO
555    
556     DO J=1,sNy
557 dimitri 1.6 CVV(J)=CV(I,J,bi,bj)
558 heimbach 1.2 END DO
559 dimitri 1.6 VRT(1)=VRT(1)/BV(I,1,bi,bj)
560 heimbach 1.2 DO J=2,sNy
561     JM=J-1
562 dimitri 1.6 CVV(J)=CVV(J)/(BV(I,J,bi,bj)-AV(I,J,bi,bj)*CVV(JM))
563     VRT(J)=(VRT(J)-AV(I,J,bi,bj)*VRT(JM))
564     & /(BV(I,J,bi,bj)-AV(I,J,bi,bj)*CVV(JM))
565 heimbach 1.2 END DO
566     DO J=1,sNy-1
567     J1=sNy-J
568     J2=J1+1
569     VRT(J1)=VRT(J1)-CVV(J1)*VRT(J2)
570     END DO
571     DO J=1,sNy
572 mlosch 1.27 VICE(I,J,bi,bj)=VTMP(I,J,bi,bj)
573 mlosch 1.26 & +WFAV*(VRT(J)-VTMP(I,J,bi,bj))
574 heimbach 1.2 END DO
575 heimbach 1.7 ENDDO
576 heimbach 1.2
577     ENDDO
578     ENDDO
579    
580 dimitri 1.6 IF(MOD(M,SOLV_NCHECK).EQ.0) THEN
581 jmc 1.22 S2=ZERO
582 dimitri 1.6 DO bj=myByLo(myThid),myByHi(myThid)
583     DO bi=myBxLo(myThid),myBxHi(myThid)
584     DO J=1,sNy
585     DO I=1,sNx
586 mlosch 1.27 UERR(I,J,bi,bj)=(VICE(I,J,bi,bj)-VTMP(I,J,bi,bj))
587 heimbach 1.2 & *UVM(I,J,bi,bj)
588 dimitri 1.6 S2=MAX(ABS(UERR(I,J,bi,bj)),S2)
589     END DO
590 heimbach 1.2 END DO
591 dimitri 1.6 ENDDO
592 heimbach 1.2 ENDDO
593 jmc 1.24 _GLOBAL_MAX_RL( S2, myThid )
594 dimitri 1.6 C SAFEGUARD AGAINST BAD FORCING ETC
595     IF(M.GT.1.AND.S2.GT.S2A) WFAV=WFAV2
596     S2A=S2
597     IF(S2.LT.LSR_ERROR) THEN
598     ICOUNT2=M
599 heimbach 1.7 phexit = 1
600 dimitri 1.6 END IF
601     END IF
602    
603 mlosch 1.27 _EXCH_XY_RL( VICE, myThid )
604 dimitri 1.6
605 jmc 1.22 ENDIF
606     ENDDO
607 heimbach 1.7 C ITERATION END -----------------------------------------------------
608    
609 jmc 1.29 IF ( debugLevel .GE. debLevC ) THEN
610 dimitri 1.11 _BEGIN_MASTER( myThid )
611     write(*,'(A,I6,1P2E22.14)')' V lsr iters, error = ',ICOUNT2,S2
612 edhill 1.14 _END_MASTER( myThid )
613 dimitri 1.11 ENDIF
614 heimbach 1.2
615     C NOW END
616     C NOW MAKE COROLIS TERM IMPLICIT
617     DO bj=myByLo(myThid),myByHi(myThid)
618     DO bi=myBxLo(myThid),myBxHi(myThid)
619     DO J=1,sNy
620     DO I=1,sNx
621 mlosch 1.27 UICE(I,J,bi,bj)=UICE(I,J,bi,bj)*UVM(I,J,bi,bj)
622     VICE(I,J,bi,bj)=VICE(I,J,bi,bj)*UVM(I,J,bi,bj)
623 heimbach 1.2 END DO
624     END DO
625     ENDDO
626     ENDDO
627 mlosch 1.27 CALL EXCH_UV_XY_RL( UICE, VICE,.TRUE.,myThid)
628 heimbach 1.2
629 dimitri 1.5 #endif /* SEAICE_ALLOW_DYNAMICS */
630 mlosch 1.18 #endif /* SEAICE_CGRID */
631 heimbach 1.2
632     RETURN
633     END
  ViewVC Help
Powered by ViewVC 1.1.22