/[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.6 - (hide annotations) (download)
Fri May 23 20:19:16 2003 UTC (21 years, 4 months ago) by dimitri
Branch: MAIN
CVS Tags: checkpoint50g_post, checkpoint50h_post, checkpoint50i_post
Changes since 1.5: +287 -242 lines 
checkpoint50g_post
o merged with release1_p17 (pkg/seaice and verification/lab_sea)
  - added SEAICE_MULTILEVEL for 8-category sea-ice thermodynamics
  - LSR sea-ice dynamic solver moved to SouthWest B-grid location and
    made the default because of faster convergence than ADI
1 heimbach 1.2 C $Header:
2    
3     #include "SEAICE_OPTIONS.h"
4    
5     CStartOfInterface
6     SUBROUTINE lsr( myThid )
7     C /==========================================================\
8     C | SUBROUTINE lsr |
9 dimitri 1.6 C | o Solve ice momentum equation with an LSR dynamics solver|
10     C | (see Zhang and Hibler, JGR, 102, 8691-8702, 1997 |
11     C | and Zhang and Rothrock, MWR, 131, 845- 861, 2003) |
12     C | Written by Jinlun Zhang, PSC/UW, Feb-2001 |
13     C | zhang@apl.washington.edu |
14 heimbach 1.2 C |==========================================================|
15     C \==========================================================/
16     IMPLICIT NONE
17    
18     C === Global variables ===
19     #include "SIZE.h"
20     #include "EEPARAMS.h"
21     #include "PARAMS.h"
22     #include "SEAICE.h"
23     #include "SEAICE_PARAMS.h"
24     #include "SEAICE_GRID.h"
25    
26    
27     C === Routine arguments ===
28     C myThid - Thread no. that called this routine.
29     INTEGER myThid
30     CEndOfInterface
31    
32     #ifdef ALLOW_SEAICE
33 dimitri 1.5 #ifdef SEAICE_ALLOW_DYNAMICS
34 heimbach 1.2
35     C === Local variables ===
36 dimitri 1.4 C i,j,bi,bj - Loop counters
37 heimbach 1.2
38 dimitri 1.6 INTEGER i, j, m, bi, bj, j1, j2, im, jm
39     INTEGER ICOUNT1, ICOUNT2, SOLV_MAX_ITERS, SOLV_NCHECK
40     _RL RADIUS, RADIUS2
41     _RL WFAU, WFAV, WFAU1, WFAV1, WFAU2, WFAV2
42     _RL AA1, AA2, AA3, AA4, AA5, AA6, S1, S2, S1A, S2A
43    
44     _RL AU (1:sNx,1:sNy,nSx,nSy), BU (1:sNx,1:sNy,nSx,nSy)
45     &, CU (1:sNx,1:sNy,nSx,nSy)
46     _RL AV (1:sNx,1:sNy,nSx,nSy), BV (1:sNx,1:sNy,nSx,nSy)
47     &, CV (1:sNx,1:sNy,nSx,nSy)
48 heimbach 1.2 _RL UERR (1:sNx,1:sNy,nSx,nSy)
49 dimitri 1.6 _RL FXY (1:sNx, 1:sNy,nSx,nSy)
50 heimbach 1.2
51     _RL URT(1:sNx), VRT(1:sNy), CUU(1:sNx), CVV(1:sNy)
52    
53 dimitri 1.6 _RL ETAMEAN (1:sNx,1:sNy,nSx,nSy)
54     _RL ZETAMEAN (1:sNx,1:sNy,nSx,nSy)
55     _RL DELXY (1:sNx,1:sNy,nSx,nSy)
56     _RL DELXR (1:sNx,1:sNy,nSx,nSy)
57     _RL DELYR (1:sNx,1:sNy,nSx,nSy)
58     _RL DELX2 (1:sNx,1:sNy,nSx,nSy)
59     _RL DELY2 (1:sNx,1:sNy,nSx,nSy)
60    
61     _RL UVRT1 (1:sNx,1:sNy,nSx,nSy)
62     _RL UVRT2 (1:sNx,1:sNy,nSx,nSy)
63    
64 heimbach 1.2 C SET SOME VALUES
65 dimitri 1.3 RADIUS=6370. _d 3
66 dimitri 1.6 RADIUS2=ONE/(RADIUS*RADIUS)
67     WFAU1=0.95 _d 0
68     WFAV1=0.95 _d 0
69     WFAU2=ZERO
70     WFAV2=ZERO
71    
72     S1A=0.80 _d 0
73     S2A=0.80 _d 0
74     WFAU=WFAU1
75     WFAV=WFAV1
76    
77     SOLV_MAX_ITERS=1500
78     SOLV_NCHECK=2
79    
80     ICOUNT1=SOLV_MAX_ITERS
81     ICOUNT2=SOLV_MAX_ITERS
82 heimbach 1.2
83     C SOLVE FOR UICE
84    
85     DO bj=myByLo(myThid),myByHi(myThid)
86     DO bi=myBxLo(myThid),myBxHi(myThid)
87     DO j=1,sNy
88     DO i=1,sNx
89     FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
90 dimitri 1.6 & +AMASS(I,J,bi,bj)/DELTAT*UICE(I,J,2,bi,bj)
91 heimbach 1.2 FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)
92 dimitri 1.6 & +AMASS(I,J,bi,bj)/DELTAT*VICE(I,J,2,bi,bj)
93 heimbach 1.2 FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)*UVM(I,J,bi,bj)
94     FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)*UVM(I,J,bi,bj)
95 dimitri 1.6 ETAMEAN(I,J,bi,bj) =QUART*(ETA(I,J-1,bi,bj)+ETA(I-1,J-1,bi,bj)
96     & +ETA(I,J,bi,bj)+ETA(I-1,J,bi,bj))
97     ZETAMEAN(I,J,bi,bj)=QUART*(ZETA(I,J-1,bi,bj)
98     & +ZETA(I-1,J-1,bi,bj)
99     & +ZETA(I,J,bi,bj)+ZETA(I-1,J,bi,bj))
100     DELX2(I,J,bi,bj)=HALF/(DXUICE(I,J,bi,bj)*DXUICE(I,J,bi,bj))
101     DELY2(I,J,bi,bj)=HALF/(DYUICE(I,J,bi,bj)*DYUICE(I,J,bi,bj))
102     DELXY(I,J,bi,bj)=HALF/(DXUICE(I,J,bi,bj)*DYUICE(I,J,bi,bj))
103     DELXR(I,J,bi,bj)=HALF/(DXUICE(I,J,bi,bj)*RADIUS)
104     DELYR(I,J,bi,bj)=HALF/(DYUICE(I,J,bi,bj)*RADIUS)
105 heimbach 1.2 ENDDO
106     ENDDO
107     ENDDO
108     ENDDO
109    
110     DO bj=myByLo(myThid),myByHi(myThid)
111     DO bi=myBxLo(myThid),myBxHi(myThid)
112    
113     DO J=1,sNy
114     DO I=1,sNx
115 dimitri 1.6 AA1=((ETA(I,J-1,bi,bj)+ZETA(I,J-1,bi,bj))
116     & /CSTICE(I-1,J-1,bi,bj)
117     & +(ETA(I,J,bi,bj)+ZETA(I,J,bi,bj))
118     & /CSTICE(I-1,J,bi,bj))/CSUICE(I,J,bi,bj)
119     AA2=((ETA(I-1,J-1,bi,bj)+ZETA(I-1,J-1,bi,bj))
120     & /CSTICE(I-1,J-1,bi,bj)
121     & +(ETA(I-1,J,bi,bj)+ZETA(I-1,J,bi,bj))
122     & /CSTICE(I-1,J,bi,bj))/CSUICE(I,J,bi,bj)
123     AA3=ETA(I-1,J,bi,bj)+ETA(I,J,bi,bj)
124     AA4=ETA(I-1,J-1,bi,bj)+ETA(I,J-1,bi,bj)
125     AA5=-(ETA(I-1,J,bi,bj)+ETA(I,J,bi,bj)-ETA(I-1,J-1,bi,bj)
126     & -ETA(I,J-1,bi,bj))*TNGICE(I,J,bi,bj)
127     AA6=TWO*ETAMEAN(I,J,bi,bj)*TNGICE(I,J,bi,bj)*TNGICE(I,J,bi,bj)
128     AU(I,J,bi,bj)=-AA2*DELX2(I,J,bi,bj)*UVM(I,J,bi,bj)
129     BU(I,J,bi,bj)=((AA1+AA2)*DELX2(I,J,bi,bj)+(AA3+AA4)
130     & *DELY2(I,J,bi,bj)
131     & +AA5*DELYR(I,J,bi,bj)+AA6*RADIUS2
132 heimbach 1.2 & +AMASS(I,J,bi,bj)/DELTAT+DRAGS(I,J,bi,bj))
133 dimitri 1.3 & *UVM(I,J,bi,bj)+(ONE-UVM(I,J,bi,bj))
134 dimitri 1.6 CU(I,J,bi,bj)=-AA1*DELX2(I,J,bi,bj)*UVM(I,J,bi,bj)
135 heimbach 1.2 END DO
136     END DO
137    
138     DO J=1,sNy
139 dimitri 1.6 AU(1,J,bi,bj)=ZERO
140     CU(sNx,J,bi,bj)=ZERO
141     CU(1,J,bi,bj)=CU(1,J,bi,bj)/BU(1,J,bi,bj)
142 heimbach 1.2 END DO
143    
144     DO J=1,sNy
145     DO I=1,sNx
146 dimitri 1.6 FXY(I,J,bi,bj)=DRAGA(I,J,bi,bj)*VICEC(I,J,bi,bj)
147 heimbach 1.2 3 +FORCEX(I,J,bi,bj)
148 dimitri 1.6 3 +HALF*(ZETA(I,J,bi,bj)*(VICEC(I+1,J+1,bi,bj)
149 heimbach 1.2 3 +VICEC(I,J+1,bi,bj)
150     3 -VICEC(I+1,J,bi,bj)-VICEC(I,J,bi,bj))
151 dimitri 1.6 3 +ZETA(I,J-1,bi,bj)*(VICEC(I+1,J,bi,bj)
152 heimbach 1.2 3 +VICEC(I,J,bi,bj)-VICEC(I+1,J-1,bi,bj)
153 dimitri 1.6 3 -VICEC(I,J-1,bi,bj))+ZETA(I-1,J,bi,bj)
154 heimbach 1.2 3 *(VICEC(I,J,bi,bj)+VICEC(I-1,J,bi,bj)
155     3 -VICEC(I,J+1,bi,bj)-VICEC(I-1,J+1,bi,bj))
156 dimitri 1.6 3 +ZETA(I-1,J-1,bi,bj)*(VICEC(I,J-1,bi,bj)
157 heimbach 1.2 3 +VICEC(I-1,J-1,bi,bj)-VICEC(I,J,bi,bj)
158 dimitri 1.6 3 -VICEC(I-1,J,bi,bj)))*DELXY(I,J,bi,bj)
159     3 /CSUICE(I,J,bi,bj)
160 heimbach 1.2 3
161 dimitri 1.6 4 -HALF*(ETA(I,J,bi,bj)*(VICEC(I+1,J+1,bi,bj)
162 heimbach 1.2 4 +VICEC(I,J+1,bi,bj)
163     4 -VICEC(I+1,J,bi,bj)-VICEC(I,J,bi,bj))
164 dimitri 1.6 4 +ETA(I,J-1,bi,bj)*(VICEC(I+1,J,bi,bj)
165 heimbach 1.2 4 +VICEC(I,J,bi,bj)-VICEC(I+1,J-1,bi,bj)
166 dimitri 1.6 4 -VICEC(I,J-1,bi,bj))+ETA(I-1,J,bi,bj)
167 heimbach 1.2 4 *(VICEC(I,J,bi,bj)+VICEC(I-1,J,bi,bj)
168     4 -VICEC(I,J+1,bi,bj)-VICEC(I-1,J+1,bi,bj))
169 dimitri 1.6 4 +ETA(I-1,J-1,bi,bj)*(VICEC(I,J-1,bi,bj)
170 heimbach 1.2 4 +VICEC(I-1,J-1,bi,bj)-VICEC(I,J,bi,bj)
171 dimitri 1.6 4 -VICEC(I-1,J,bi,bj)))*DELXY(I,J,bi,bj)
172     4 /CSUICE(I,J,bi,bj)
173 heimbach 1.2 4
174 dimitri 1.6 5 +HALF*(ETA(I,J,bi,bj)/CSTICE(I-1,J,bi,bj)
175 heimbach 1.2 5 *(VICEC(I+1,J+1,bi,bj)+VICEC(I+1,J,bi,bj)
176     5 -VICEC(I,J+1,bi,bj)-VICEC(I,J,bi,bj))
177 dimitri 1.6 5 +ETA(I-1,J,bi,bj)/CSTICE(I-1,J,bi,bj)
178 heimbach 1.2 5 *(VICEC(I,J+1,bi,bj)
179     5 +VICEC(I,J,bi,bj)-VICEC(I-1,J+1,bi,bj)
180     5 -VICEC(I-1,J,bi,bj))
181 dimitri 1.6 5 +ETA(I,J-1,bi,bj)/CSTICE(I-1,J-1,bi,bj)
182 heimbach 1.2 5 *(VICEC(I,J,bi,bj)+VICEC(I,J-1,bi,bj)
183     5 -VICEC(I+1,J,bi,bj)-VICEC(I+1,J-1,bi,bj))
184 dimitri 1.6 5 +ETA(I-1,J-1,bi,bj)/CSTICE(I-1,J-1,bi,bj)
185     5 *(VICEC(I-1,J,bi,bj)
186 heimbach 1.2 5 +VICEC(I-1,J-1,bi,bj)-VICEC(I,J,bi,bj)
187 dimitri 1.6 5 -VICEC(I,J-1,bi,bj)))*DELXY(I,J,bi,bj)
188 heimbach 1.2 5
189 dimitri 1.6 6 -((ZETA(I,J,bi,bj)+ZETA(I,J-1,bi,bj)
190     6 -ZETA(I-1,J-1,bi,bj)-ZETA(I-1,J,bi,bj))
191     6 +(ETA(I,J,bi,bj)+ETA(I,J-1,bi,bj)-ETA(I-1,J-1,bi,bj)
192     6 -ETA(I-1,J,bi,bj)))
193     6 *TNGICE(I,J,bi,bj)*VICEC(I,J,bi,bj)*DELXR(I,J,bi,bj)
194 heimbach 1.2 6 /CSUICE(I,J,bi,bj)
195 dimitri 1.6 6 -(ETAMEAN(I,J,bi,bj)+ZETAMEAN(I,J,bi,bj))
196     6 *TNGICE(I,J,bi,bj)
197 heimbach 1.2 6 *(VICEC(I+1,J,bi,bj)-VICEC(I-1,J,bi,bj))
198 dimitri 1.6 6 *DELXR(I,J,bi,bj)/CSUICE(I,J,bi,bj)
199 heimbach 1.2 6
200 dimitri 1.6 7 -ETAMEAN(I,J,bi,bj)*TWO*TNGICE(I,J,bi,bj)
201     7 *(VICEC(I+1,J,bi,bj)
202     7 -VICEC(I-1,J,bi,bj))*DELXR(I,J,bi,bj)
203     7 /CSUICE(I,J,bi,bj)
204     UVRT1(I,J,bi,bj)=(ETA(I-1,J-1,bi,bj)+ETA(I,J-1,bi,bj))
205     & *DELY2(I,J,bi,bj)
206     & -ETAMEAN(I,J,bi,bj)*DELYR(I,J,bi,bj)
207     & *TNGICE(I,J-1,bi,bj)
208     & +ETAMEAN(I,J,bi,bj)*DELYR(I,J,bi,bj)
209     & *TWO*TNGICE(I,J,bi,bj)
210     UVRT2(I,J,bi,bj)=(ETA(I-1,J,bi,bj)+ETA(I,J,bi,bj))
211     & *DELY2(I,J,bi,bj)
212     & +ETAMEAN(I,J,bi,bj)*DELYR(I,J,bi,bj)
213     & *TNGICE(I,J+1,bi,bj)
214     & -ETAMEAN(I,J,bi,bj)*DELYR(I,J,bi,bj)
215     & *TWO*TNGICE(I,J,bi,bj)
216 heimbach 1.2 END DO
217     END DO
218    
219     ENDDO
220     ENDDO
221    
222     C NOW DO ITERATION
223     100 CONTINUE
224 dimitri 1.6
225     DO 8000 M=1, solv_max_iters
226 heimbach 1.2
227     DO bj=myByLo(myThid),myByHi(myThid)
228     DO bi=myBxLo(myThid),myBxHi(myThid)
229     C NOW SET U(3)=U(1)
230     DO J=1,sNy
231     DO I=1,sNx
232     UICE(I,J,3,bi,bj)=UICE(I,J,1,bi,bj)
233     END DO
234     END DO
235    
236     DO 1200 J=1,sNy
237     DO I=1,sNx
238 dimitri 1.6 IF(I.EQ.1) THEN
239     AA2=((ETA(I-1,J-1,bi,bj)+ZETA(I-1,J-1,bi,bj))
240     & /CSTICE(I-1,J-1,bi,bj)
241     & +(ETA(I-1,J,bi,bj)+ZETA(I-1,J,bi,bj))
242     & /CSTICE(I-1,J,bi,bj))/CSUICE(I,J,bi,bj)
243     AA3=AA2*DELX2(I,J,bi,bj)*UICE(I-1,J,1,bi,bj)*UVM(I,J,bi,bj)
244     ELSE IF(I.EQ.sNx) THEN
245     AA1=((ETA(I,J-1,bi,bj)+ZETA(I,J-1,bi,bj))
246     & /CSTICE(I-1,J-1,bi,bj)
247     & +(ETA(I,J,bi,bj)+ZETA(I,J,bi,bj))
248     & /CSTICE(I-1,J,bi,bj))/CSUICE(I,J,bi,bj)
249     AA3=AA1*DELX2(I,J,bi,bj)*UICE(I+1,J,1,bi,bj)*UVM(I,J,bi,bj)
250     ELSE
251     AA3=ZERO
252     END IF
253     URT(I)=FXY(I,J,bi,bj)+AA3
254     & +UVRT1(I,J,bi,bj)*UICE(I,J-1,1,bi,bj)
255     & +UVRT2(I,J,bi,bj)*UICE(I,J+1,1,bi,bj)
256 heimbach 1.2 URT(I)=URT(I)*UVM(I,J,bi,bj)
257     END DO
258    
259     DO I=1,sNx
260 dimitri 1.6 CUU(I)=CU(I,J,bi,bj)
261 heimbach 1.2 END DO
262 dimitri 1.6 URT(1)=URT(1)/BU(1,J,bi,bj)
263 heimbach 1.2 DO I=2,sNx
264     IM=I-1
265 dimitri 1.6 CUU(I)=CUU(I)/(BU(I,J,bi,bj)-AU(I,J,bi,bj)*CUU(IM))
266     URT(I)=(URT(I)-AU(I,J,bi,bj)*URT(IM))
267     & /(BU(I,J,bi,bj)-AU(I,J,bi,bj)*CUU(IM))
268 heimbach 1.2 END DO
269     DO I=1,sNx-1
270     J1=sNx-I
271     J2=J1+1
272     URT(J1)=URT(J1)-CUU(J1)*URT(J2)
273     END DO
274     DO I=1,sNx
275     UICE(I,J,1,bi,bj)=UICE(I,J,3,bi,bj)
276     & +WFAU*(URT(I)-UICE(I,J,3,bi,bj))
277     END DO
278    
279     1200 CONTINUE
280    
281     ENDDO
282     ENDDO
283    
284 dimitri 1.6 IF(MOD(M,SOLV_NCHECK).EQ.0) THEN
285     DO bj=myByLo(myThid),myByHi(myThid)
286     DO bi=myBxLo(myThid),myBxHi(myThid)
287     S1=ZERO
288     DO J=1,sNy
289     DO I=1,sNx
290     UERR(I,J,bi,bj)=(UICE(I,J,1,bi,bj)-UICE(I,J,3,bi,bj))
291 heimbach 1.2 & *UVM(I,J,bi,bj)
292 dimitri 1.6 S1=MAX(ABS(UERR(I,J,bi,bj)),S1)
293     END DO
294 heimbach 1.2 END DO
295 dimitri 1.6 _GLOBAL_MAX_R8( S1, myThid )
296     ENDDO
297 heimbach 1.2 ENDDO
298 dimitri 1.6 C SAFEGUARD AGAINST BAD FORCING ETC
299     IF(M.GT.1.AND.S1.GT.S1A) WFAU=WFAU2
300     S1A=S1
301     IF(S1.LT.LSR_ERROR) THEN
302     ICOUNT1=M
303     GO TO 8001
304     END IF
305     END IF
306 dimitri 1.5 CALL SEAICE_EXCH ( UICE, myThid )
307 dimitri 1.6
308     8000 CONTINUE
309     8001 CONTINUE
310     write(*,'(A,I6,1P2E22.14)')' U lsr iters, error = ',ICOUNT1,S1
311 heimbach 1.2
312     C NOW FOR VICE
313     DO bj=myByLo(myThid),myByHi(myThid)
314     DO bi=myBxLo(myThid),myBxHi(myThid)
315    
316     DO J=1,sNy
317     DO I=1,sNx
318 dimitri 1.6 AA1=ETA(I-1,J,bi,bj)+ZETA(I-1,J,bi,bj)+ETA(I,J,bi,bj)
319     & +ZETA(I,J,bi,bj)
320     AA2=ETA(I-1,J-1,bi,bj)+ZETA(I-1,J-1,bi,bj)+ETA(I,J-1,bi,bj)
321     & +ZETA(I,J-1,bi,bj)
322     AA3=(ETA(I,J-1,bi,bj)/CSUICE(I,J,bi,bj)+ETA(I,J,bi,bj)
323 heimbach 1.2 & /CSUICE(I,J,bi,bj))/CSUICE(I,J,bi,bj)
324 dimitri 1.6 AA4=(ETA(I-1,J-1,bi,bj)/CSUICE(I,J,bi,bj)+ETA(I-1,J,bi,bj)
325 heimbach 1.2 & /CSUICE(I,J,bi,bj))/CSUICE(I,J,bi,bj)
326 dimitri 1.6 AA5=((ZETA(I-1,J,bi,bj)-ETA(I-1,J,bi,bj))
327     & +(ZETA(I,J,bi,bj)-ETA(I,J,bi,bj))
328     & -(ZETA(I-1,J-1,bi,bj)-ETA(I-1,J-1,bi,bj))
329     & -(ZETA(I,J-1,bi,bj)
330     & -ETA(I,J-1,bi,bj)))*TNGICE(I,J,bi,bj)
331     AA6=TWO*ETAMEAN(I,J,bi,bj)*TNGICE(I,J,bi,bj)*TNGICE(I,J,bi,bj)
332    
333     AV(I,J,bi,bj)=(-AA2*DELY2(I,J,bi,bj)
334     & -(ZETAMEAN(I,J,bi,bj)-ETAMEAN(I,J,bi,bj))
335     & *TNGICE(I,J-1,bi,bj)*DELYR(I,J,bi,bj)
336     & -ETAMEAN(I,J,bi,bj)*TWO*TNGICE(I,J,bi,bj)*DELYR(I,J,bi,bj))
337     & *UVM(I,J,bi,bj)
338     BV(I,J,bi,bj)=((AA1+AA2)*DELY2(I,J,bi,bj)
339     & +(AA3+AA4)*DELX2(I,J,bi,bj)
340     & +AA5*DELYR(I,J,bi,bj)+AA6*RADIUS2
341     & +AMASS(I,J,bi,bj)/DELTAT+DRAGS(I,J,bi,bj))
342     & *UVM(I,J,bi,bj)+(ONE-UVM(I,J,bi,bj))
343     CV(I,J,bi,bj)=(-AA1*DELY2(I,J,bi,bj)
344     & +(ZETAMEAN(I,J,bi,bj)-ETAMEAN(I,J,bi,bj))
345     & *TNGICE(I,J+1,bi,bj)*DELYR(I,J,bi,bj)
346     & +ETAMEAN(I,J,bi,bj)*TWO*TNGICE(I,J,bi,bj)*DELYR(I,J,bi,bj))
347     & *UVM(I,J,bi,bj)
348 heimbach 1.2 END DO
349     END DO
350    
351     DO I=1,sNx
352 dimitri 1.6 AV(I,1,bi,bj)=ZERO
353     CV(I,sNy,bi,bj)=ZERO
354     CV(I,1,bi,bj)=CV(I,1,bi,bj)/BV(I,1,bi,bj)
355 heimbach 1.2 END DO
356    
357     DO J=1,sNy
358     DO I=1,sNx
359 dimitri 1.6 FXY(I,J,bi,bj)=-DRAGA(I,J,bi,bj)*UICEC(I,J,bi,bj)
360 heimbach 1.2 2 +FORCEY(I,J,bi,bj)
361 dimitri 1.3 3 +(HALF*(UICEC(I+1,J,bi,bj)-UICEC(I-1,J,bi,bj))
362 dimitri 1.6 3 *(ZETA(I-1,J,bi,bj)+ZETA(I,J,bi,bj)
363     3 -ZETA(I-1,J-1,bi,bj)-ZETA(I,J-1,bi,bj))*DELXY(I,J,bi,bj)
364     3 /CSUICE(I,J,bi,bj)+HALF*ZETAMEAN(I,J,bi,bj)
365     3 *((UICEC(I+1,J+1,bi,bj)
366 heimbach 1.2 3 -UICEC(I-1,J+1,bi,bj))/CSUICE(I,J+1,bi,bj)
367     3 -(UICEC(I+1,J-1,bi,bj)-UICEC(I-1,J-1,bi,bj))
368 dimitri 1.6 3 /CSUICE(I,J-1,bi,bj))*DELXY(I,J,bi,bj))
369 heimbach 1.2 3
370 dimitri 1.3 4 -(HALF*(UICEC(I+1,J,bi,bj)-UICEC(I-1,J,bi,bj))
371 dimitri 1.6 4 *(ETA(I-1,J,bi,bj)+ETA(I,J,bi,bj)
372     4 -ETA(I-1,J-1,bi,bj)-ETA(I,J-1,bi,bj))*DELXY(I,J,bi,bj)
373     4 /CSUICE(I,J,bi,bj)
374     4 +HALF*ETAMEAN(I,J,bi,bj)*((UICEC(I+1,J+1,bi,bj)
375 heimbach 1.2 4 -UICEC(I-1,J+1,bi,bj))/CSUICE(I,J+1,bi,bj)
376     4 -(UICEC(I+1,J-1,bi,bj)-UICEC(I-1,J-1,bi,bj))
377 dimitri 1.6 4 /CSUICE(I,J-1,bi,bj))*DELXY(I,J,bi,bj))
378 heimbach 1.2 4
379 dimitri 1.6 5 +HALF*(ETA(I,J,bi,bj)*(UICEC(I+1,J+1,bi,bj)
380 heimbach 1.2 5 +UICEC(I,J+1,bi,bj)
381 dimitri 1.6 5 -UICEC(I+1,J,bi,bj)-UICEC(I,J,bi,bj))+ETA(I,J-1,bi,bj)
382 heimbach 1.2 5 *(UICEC(I+1,J,bi,bj)
383     5 +UICEC(I,J,bi,bj)-UICEC(I+1,J-1,bi,bj)
384 dimitri 1.6 5 -UICEC(I,J-1,bi,bj))+ETA(I-1,J,bi,bj)
385 heimbach 1.2 5 *(UICEC(I,J,bi,bj)+UICEC(I-1,J,bi,bj)-UICEC(I,J+1,bi,bj)
386     5 -UICEC(I-1,J+1,bi,bj))
387 dimitri 1.6 5 +ETA(I-1,J-1,bi,bj)*(UICEC(I,J-1,bi,bj)
388     5 +UICEC(I-1,J-1,bi,bj)
389 heimbach 1.2 5 -UICEC(I,J,bi,bj)
390 dimitri 1.6 5 -UICEC(I-1,J,bi,bj)))*DELXY(I,J,bi,bj)/CSUICE(I,J,bi,bj)
391 heimbach 1.2 5
392 dimitri 1.6 6 +(ETA(I,J,bi,bj)+ETA(I,J-1,bi,bj)-ETA(I-1,J-1,bi,bj)
393     6 -ETA(I-1,J,bi,bj))
394 heimbach 1.2 6 *TNGICE(I,J,bi,bj)*UICEC(I,J,bi,bj)
395 dimitri 1.6 6 *DELXR(I,J,bi,bj)/CSUICE(I,J,bi,bj)
396     6 +ETAMEAN(I,J,bi,bj)*TNGICE(I,J,bi,bj)*(UICEC(I+1,J,bi,bj)
397     6 -UICEC(I-1,J,bi,bj))*DELXR(I,J,bi,bj)/CSUICE(I,J,bi,bj)
398 heimbach 1.2 6
399 dimitri 1.6 7 +ETAMEAN(I,J,bi,bj)*TWO*TNGICE(I,J,bi,bj)
400     7 *(UICEC(I+1,J,bi,bj)
401     7 -UICEC(I-1,J,bi,bj))*DELXR(I,J,bi,bj)/CSUICE(I,J,bi,bj)
402     UVRT1(I,J,bi,bj)=(ETA(I-1,J-1,bi,bj)/CSUICE(I,J,bi,bj)
403     & +ETA(I-1,J,bi,bj)/CSUICE(I,J,bi,bj))
404     & *DELX2(I,J,bi,bj)/CSUICE(I,J,bi,bj)
405     UVRT2(I,J,bi,bj)=(ETA(I,J-1,bi,bj)/CSUICE(I,J,bi,bj)
406     & +ETA(I,J,bi,bj)/CSUICE(I,J,bi,bj))
407     & *DELX2(I,J,bi,bj)/CSUICE(I,J,bi,bj)
408 heimbach 1.2
409     END DO
410     END DO
411    
412     ENDDO
413     ENDDO
414    
415     C NOW DO ITERATION
416     300 CONTINUE
417 dimitri 1.6
418     DO 9000 M=1, solv_max_iters
419 heimbach 1.2 C NOW SET U(3)=U(1)
420     DO bj=myByLo(myThid),myByHi(myThid)
421     DO bi=myBxLo(myThid),myBxHi(myThid)
422    
423     DO J=1,sNy
424     DO I=1,sNx
425     VICE(I,J,3,bi,bj)=VICE(I,J,1,bi,bj)
426     END DO
427     END DO
428    
429     DO 1300 I=1,sNx
430     DO J=1,sNy
431 dimitri 1.6 IF(J.EQ.1) THEN
432     AA2=ETA(I-1,J-1,bi,bj)+ZETA(I-1,J-1,bi,bj)
433     & +ETA(I,J-1,bi,bj)+ZETA(I,J-1,bi,bj)
434     AA3=(AA2*DELY2(I,J,bi,bj)
435     & +(ZETAMEAN(I,J,bi,bj)-ETAMEAN(I,J,bi,bj))
436     & *TNGICE(I,J-1,bi,bj)*DELYR(I,J,bi,bj)
437     & +ETAMEAN(I,J,bi,bj)*TWO*TNGICE(I,J,bi,bj)
438     & *DELYR(I,J,bi,bj))
439     & *VICE(I,J-1,1,bi,bj)*UVM(I,J,bi,bj)
440     ELSE IF(J.EQ.sNy) THEN
441     AA1=ETA(I-1,J,bi,bj)+ZETA(I-1,J,bi,bj)+ETA(I,J,bi,bj)
442     & +ZETA(I,J,bi,bj)
443     AA3=(AA1*DELY2(I,J,bi,bj)
444     & -(ZETAMEAN(I,J,bi,bj)-ETAMEAN(I,J,bi,bj))
445     & *TNGICE(I,J+1,bi,bj)*DELYR(I,J,bi,bj)
446     & -ETAMEAN(I,J,bi,bj)*TWO*TNGICE(I,J,bi,bj)
447     & *DELYR(I,J,bi,bj))
448     & *VICE(I,J+1,1,bi,bj)*UVM(I,J,bi,bj)
449     ELSE
450     AA3=ZERO
451     END IF
452    
453     VRT(J)=FXY(I,J,bi,bj)+AA3+UVRT1(I,J,bi,bj)*VICE(I-1,J,1,bi,bj)
454     & +UVRT2(I,J,bi,bj)*VICE(I+1,J,1,bi,bj)
455 heimbach 1.2 VRT(J)=VRT(J)*UVM(I,J,bi,bj)
456     END DO
457    
458     DO J=1,sNy
459 dimitri 1.6 CVV(J)=CV(I,J,bi,bj)
460 heimbach 1.2 END DO
461 dimitri 1.6 VRT(1)=VRT(1)/BV(I,1,bi,bj)
462 heimbach 1.2 DO J=2,sNy
463     JM=J-1
464 dimitri 1.6 CVV(J)=CVV(J)/(BV(I,J,bi,bj)-AV(I,J,bi,bj)*CVV(JM))
465     VRT(J)=(VRT(J)-AV(I,J,bi,bj)*VRT(JM))
466     & /(BV(I,J,bi,bj)-AV(I,J,bi,bj)*CVV(JM))
467 heimbach 1.2 END DO
468     DO J=1,sNy-1
469     J1=sNy-J
470     J2=J1+1
471     VRT(J1)=VRT(J1)-CVV(J1)*VRT(J2)
472     END DO
473     DO J=1,sNy
474     VICE(I,J,1,bi,bj)=VICE(I,J,3,bi,bj)
475     & +WFAV*(VRT(J)-VICE(I,J,3,bi,bj))
476     END DO
477     1300 CONTINUE
478    
479     ENDDO
480     ENDDO
481    
482 dimitri 1.6 IF(MOD(M,SOLV_NCHECK).EQ.0) THEN
483     DO bj=myByLo(myThid),myByHi(myThid)
484     DO bi=myBxLo(myThid),myBxHi(myThid)
485     S2=ZERO
486     DO J=1,sNy
487     DO I=1,sNx
488     UERR(I,J,bi,bj)=(VICE(I,J,1,bi,bj)-VICE(I,J,3,bi,bj))
489 heimbach 1.2 & *UVM(I,J,bi,bj)
490 dimitri 1.6 S2=MAX(ABS(UERR(I,J,bi,bj)),S2)
491     END DO
492 heimbach 1.2 END DO
493 dimitri 1.6 _GLOBAL_MAX_R8( S2, myThid )
494     ENDDO
495 heimbach 1.2 ENDDO
496 dimitri 1.6 C SAFEGUARD AGAINST BAD FORCING ETC
497     IF(M.GT.1.AND.S2.GT.S2A) WFAV=WFAV2
498     S2A=S2
499     IF(S2.LT.LSR_ERROR) THEN
500     ICOUNT2=M
501     GO TO 9001
502     END IF
503     END IF
504    
505     CALL SEAICE_EXCH ( VICE, myThid )
506    
507     9000 CONTINUE
508     9001 CONTINUE
509     write(*,'(A,I6,1P2E22.14)')' V lsr iters, error = ',ICOUNT2,S2
510 heimbach 1.2
511     C NOW END
512     C NOW MAKE COROLIS TERM IMPLICIT
513     DO bj=myByLo(myThid),myByHi(myThid)
514     DO bi=myBxLo(myThid),myBxHi(myThid)
515     DO J=1,sNy
516     DO I=1,sNx
517     UICE(I,J,1,bi,bj)=UICE(I,J,1,bi,bj)*UVM(I,J,bi,bj)
518     VICE(I,J,1,bi,bj)=VICE(I,J,1,bi,bj)*UVM(I,J,bi,bj)
519     END DO
520     END DO
521     ENDDO
522     ENDDO
523 dimitri 1.6 CALL SEAICE_EXCH ( UICE, myThid )
524     CALL SEAICE_EXCH ( VICE, myThid )
525 heimbach 1.2
526 dimitri 1.5 #endif /* SEAICE_ALLOW_DYNAMICS */
527     #endif /* ALLOW_SEAICE */
528 heimbach 1.2
529     RETURN
530     END
  ViewVC Help
Powered by ViewVC 1.1.22