/[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.5 - (show annotations) (download)
Tue Feb 18 05:33:55 2003 UTC (21 years, 2 months ago) by dimitri
Branch: MAIN
CVS Tags: checkpoint50c_post, checkpoint50c_pre, checkpoint48i_post, checkpoint50, checkpoint50d_post, checkpoint50b_pre, checkpoint48f_post, checkpoint48h_post, checkpoint50f_post, checkpoint50a_post, checkpoint50f_pre, checkpoint50e_pre, checkpoint50e_post, checkpoint50d_pre, checkpoint49, checkpoint48g_post, checkpoint50b_post
Changes since 1.4: +7 -9 lines 
Merging from release1_p12:
o Modifications for using pkg/exf with pkg/seaice
  - improved description of the various forcing configurations
  - added basic radiation bulk formulae to pkg/exf
  - units/sign fix for evap computation in exf_getffields.F
  - updated verification/global_with_exf/results/output.txt
o Added pkg/sbo for computing IERS Special Bureau for the Oceans
  (SBO) core products, including oceanic mass, center-of-mass,
  angular, and bottom pressure (see pkg/sbo/README.sbo).
o Lower bound for viscosity/diffusivity in pkg/kpp/kpp_routines.F
  to avoid negative values in shallow regions.
  - updated verification/natl_box/results/output.txt
  - updated verification/lab_sea/results/output.txt
o MPI gather, scatter: eesupp/src/gather_2d.F and scatter_2d.F
o Added useSingleCpuIO option (see PARAMS.h).
o Updated useSingleCpuIO option in mdsio_writefield.F to
  work with multi-field files, e.g., for single-file pickup.
o pkg/seaice:
  - bug fix in growth.F: QNET for no shortwave case
  - added HeffFile for specifying initial sea-ice thickness
  - changed SEAICE_EXTERNAL_FLUXES wind stress implementation
o Added missing /* */ to CPP comments in pkg/seaice, pkg/exf,
  kpp_transport_t.F, forward_step.F, and the_main_loop.F
o pkg/seaice:
  - adjoint-friendly modifications
  - added a SEAICE_WRITE_PICKUP at end of the_model_main.F
1 C $Header:
2
3 #include "SEAICE_OPTIONS.h"
4
5 CStartOfInterface
6 SUBROUTINE lsr( myThid )
7 C /==========================================================\
8 C | SUBROUTINE lsr |
9 C | o Solve ice momentum equation with an LSR method |
10 C | (see Zhang and Hibler, JGR, 102, 8691-8702, 1997) |
11 C |==========================================================|
12 C \==========================================================/
13 IMPLICIT NONE
14
15 C === Global variables ===
16 #include "SIZE.h"
17 #include "EEPARAMS.h"
18 #include "PARAMS.h"
19 #include "SEAICE.h"
20 #include "SEAICE_PARAMS.h"
21 #include "SEAICE_GRID.h"
22
23
24 C === Routine arguments ===
25 C myThid - Thread no. that called this routine.
26 INTEGER myThid
27 CEndOfInterface
28
29 #ifdef ALLOW_SEAICE
30 #ifdef SEAICE_ALLOW_DYNAMICS
31
32 C === Local variables ===
33 C i,j,bi,bj - Loop counters
34
35 INTEGER i, j, bi, bj, j1, j2, im, jm, icou
36 _RL RADIUS, DELXY, DELXR, DELYR, DELX2, DELY2, RADIUS2
37 _RL ETAMEAN, ZETAMEAN, WFAU, WFAV
38 _RL AA1, AA2, AA3, AA4, AA5, AA6, S1, S2
39
40 _RL AU (1:sNx,1:sNy), BU (1:sNx,1:sNy), CU (1:sNx,1:sNy)
41 _RL AV (1:sNx,1:sNy), BV (1:sNx,1:sNy), CV (1:sNx,1:sNy)
42 _RL UERR (1:sNx,1:sNy,nSx,nSy)
43 _RL FXY (1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
44
45 _RL URT(1:sNx), VRT(1:sNy), CUU(1:sNx), CVV(1:sNy)
46
47 C SET SOME VALUES
48 RADIUS=6370. _d 3
49 ICOUNT=0
50 ICOU=0
51 WFAU=0.95 _d 0
52 WFAV=0.95 _d 0
53
54 C SOLVE FOR UICE
55
56 DO bj=myByLo(myThid),myByHi(myThid)
57 DO bi=myBxLo(myThid),myBxHi(myThid)
58 DO j=1,sNy
59 DO i=1,sNx
60 FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
61 & +AMASS(I,J,bi,bj)/DELTAT*UICE(I,J,2,bi,bj)
62 FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)
63 & +AMASS(I,J,bi,bj)/DELTAT*VICE(I,J,2,bi,bj)
64 FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)*UVM(I,J,bi,bj)
65 FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)*UVM(I,J,bi,bj)
66 ENDDO
67 ENDDO
68 ENDDO
69 ENDDO
70
71 DO bj=myByLo(myThid),myByHi(myThid)
72 DO bi=myBxLo(myThid),myBxHi(myThid)
73
74 DO J=1,sNy
75 DO I=1,sNx
76 DELX2=HALF/(DXUICE(I,J,bi,bj)*DXUICE(I,J,bi,bj))
77 DELY2=HALF/(DYUICE(I,J,bi,bj)*DYUICE(I,J,bi,bj))
78 DELYR=HALF/(DYUICE(I,J,bi,bj)*RADIUS)
79 RADIUS2=ONE/(RADIUS*RADIUS)
80 ETAMEAN=QUART*(ETA(I,J+1,bi,bj)+ETA(I+1,J+1,bi,bj)
81 & +ETA(I,J,bi,bj)+ETA(I+1,J,bi,bj))
82 AA1=((ETA(I+1,J,bi,bj)+ZETA(I+1,J,bi,bj))/CSTICE(I,J,bi,bj)
83 & +(ETA(I+1,J+1,bi,bj)+ZETA(I+1,J+1,bi,bj))
84 & /CSTICE(I,J+1,bi,bj))/CSUICE(I,J,bi,bj)
85 AA2=((ETA(I,J,bi,bj)+ZETA(I,J,bi,bj))/CSTICE(I,J,bi,bj)
86 & +(ETA(I,J+1,bi,bj)+ZETA(I,J+1,bi,bj))
87 & /CSTICE(I,J+1,bi,bj))/CSUICE(I,J,bi,bj)
88 AA3=ETA(I,J+1,bi,bj)+ETA(I+1,J+1,bi,bj)
89 AA4=ETA(I,J,bi,bj)+ETA(I+1,J,bi,bj)
90 AA5=-(ETA(I,J+1,bi,bj)+ETA(I+1,J+1,bi,bj)-ETA(I,J,bi,bj)
91 & -ETA(I+1,J,bi,bj))*TNGICE(I,J,bi,bj)
92 AA6=TWO*ETAMEAN*TNGICE(I,J,bi,bj)*TNGICE(I,J,bi,bj)
93 AU(I,J)=-AA2*DELX2*UVM(I,J,bi,bj)
94 BU(I,J)=((AA1+AA2)*DELX2+(AA3+AA4)*DELY2+AA5*DELYR+AA6*RADIUS2
95 & +AMASS(I,J,bi,bj)/DELTAT+DRAGS(I,J,bi,bj))
96 & *UVM(I,J,bi,bj)+(ONE-UVM(I,J,bi,bj))
97 CU(I,J)=-AA1*DELX2*UVM(I,J,bi,bj)
98 END DO
99 END DO
100
101 DO J=1,sNy
102 AU(1,J)=ZERO
103 CU(sNx,J)=ZERO
104 CU(1,J)=CU(1,J)/BU(1,J)
105 END DO
106
107 DO J=1,sNy
108 DO I=1,sNx
109 DELXY=HALF/(DXUICE(I,J,bi,bj)*DYUICE(I,J,bi,bj))
110 DELXR=HALF/(DXUICE(I,J,bi,bj)*RADIUS)
111 DELX2=HALF/(DXUICE(I,J,bi,bj)*DXUICE(I,J,bi,bj))
112 ETAMEAN=QUART*(ETA(I,J+1,bi,bj)+ETA(I+1,J+1,bi,bj)
113 & +ETA(I,J,bi,bj)+ETA(I+1,J,bi,bj))
114 ZETAMEAN=QUART*(ZETA(I,J+1,bi,bj)+ZETA(I+1,J+1,bi,bj)
115 & +ZETA(I,J,bi,bj)+ZETA(I+1,J,bi,bj))
116 AA1=((ETA(I+1,J,bi,bj)+ZETA(I+1,J,bi,bj))/CSUICE(I,J,bi,bj)
117 & +(ETA(I+1,J+1,bi,bj)+ZETA(I+1,J+1,bi,bj))
118 & /CSUICE(I,J,bi,bj))/CSUICE(I,J,bi,bj)
119 AA2=((ETA(I,J,bi,bj)+ZETA(I,J,bi,bj))/CSUICE(I,J,bi,bj)
120 & +(ETA(I,J+1,bi,bj)+ZETA(I,J+1,bi,bj))
121 & /CSUICE(I,J,bi,bj))/CSUICE(I,J,bi,bj)
122
123 IF(I.EQ.1) THEN
124 AA3=AA2*DELX2*UICE(I-1,J,1,bi,bj)*UVM(I,J,bi,bj)
125 ELSE IF(I.EQ.sNx) THEN
126 AA3=AA1*DELX2*UICE(I+1,J,1,bi,bj)*UVM(I,J,bi,bj)
127 ELSE
128 AA3=ZERO
129 END IF
130
131 FXY(I,J)=AA3+DRAGA(I,J,bi,bj)*VICEC(I,J,bi,bj)
132 3 +FORCEX(I,J,bi,bj)
133 3 +HALF*(ZETA(I+1,J+1,bi,bj)*(VICEC(I+1,J+1,bi,bj)
134 3 +VICEC(I,J+1,bi,bj)
135 3 -VICEC(I+1,J,bi,bj)-VICEC(I,J,bi,bj))
136 3 +ZETA(I+1,J,bi,bj)*(VICEC(I+1,J,bi,bj)
137 3 +VICEC(I,J,bi,bj)-VICEC(I+1,J-1,bi,bj)
138 3 -VICEC(I,J-1,bi,bj))+ZETA(I,J+1,bi,bj)
139 3 *(VICEC(I,J,bi,bj)+VICEC(I-1,J,bi,bj)
140 3 -VICEC(I,J+1,bi,bj)-VICEC(I-1,J+1,bi,bj))
141 3 +ZETA(I,J,bi,bj)*(VICEC(I,J-1,bi,bj)
142 3 +VICEC(I-1,J-1,bi,bj)-VICEC(I,J,bi,bj)
143 3 -VICEC(I-1,J,bi,bj)))*DELXY/CSUICE(I,J,bi,bj)
144 3
145 4 -HALF*(ETA(I+1,J+1,bi,bj)*(VICEC(I+1,J+1,bi,bj)
146 4 +VICEC(I,J+1,bi,bj)
147 4 -VICEC(I+1,J,bi,bj)-VICEC(I,J,bi,bj))
148 4 +ETA(I+1,J,bi,bj)*(VICEC(I+1,J,bi,bj)
149 4 +VICEC(I,J,bi,bj)-VICEC(I+1,J-1,bi,bj)
150 4 -VICEC(I,J-1,bi,bj))+ETA(I,J+1,bi,bj)
151 4 *(VICEC(I,J,bi,bj)+VICEC(I-1,J,bi,bj)
152 4 -VICEC(I,J+1,bi,bj)-VICEC(I-1,J+1,bi,bj))
153 4 +ETA(I,J,bi,bj)*(VICEC(I,J-1,bi,bj)
154 4 +VICEC(I-1,J-1,bi,bj)-VICEC(I,J,bi,bj)
155 4 -VICEC(I-1,J,bi,bj)))*DELXY/CSUICE(I,J,bi,bj)
156 4
157 5 +HALF*(ETA(I+1,J+1,bi,bj)/CSTICE(I,J+1,bi,bj)
158 5 *(VICEC(I+1,J+1,bi,bj)+VICEC(I+1,J,bi,bj)
159 5 -VICEC(I,J+1,bi,bj)-VICEC(I,J,bi,bj))
160 5 +ETA(I,J+1,bi,bj)/CSTICE(I,J+1,bi,bj)
161 5 *(VICEC(I,J+1,bi,bj)
162 5 +VICEC(I,J,bi,bj)-VICEC(I-1,J+1,bi,bj)
163 5 -VICEC(I-1,J,bi,bj))
164 5 +ETA(I+1,J,bi,bj)/CSTICE(I,J,bi,bj)
165 5 *(VICEC(I,J,bi,bj)+VICEC(I,J-1,bi,bj)
166 5 -VICEC(I+1,J,bi,bj)-VICEC(I+1,J-1,bi,bj))
167 5 +ETA(I,J,bi,bj)/CSTICE(I,J,bi,bj)*(VICEC(I-1,J,bi,bj)
168 5 +VICEC(I-1,J-1,bi,bj)-VICEC(I,J,bi,bj)
169 5 -VICEC(I,J-1,bi,bj)))*DELXY
170 5
171 6 -((ZETA(I+1,J+1,bi,bj)+ZETA(I+1,J,bi,bj)
172 6 -ZETA(I,J,bi,bj)-ZETA(I,J+1,bi,bj))
173 6 +(ETA(I+1,J+1,bi,bj)+ETA(I+1,J,bi,bj)-ETA(I,J,bi,bj)
174 6 -ETA(I,J+1,bi,bj)))
175 6 *TNGICE(I,J,bi,bj)*VICEC(I,J,bi,bj)*DELXR
176 6 /CSUICE(I,J,bi,bj)
177 6 -(ETAMEAN+ZETAMEAN)*TNGICE(I,J,bi,bj)
178 6 *(VICEC(I+1,J,bi,bj)-VICEC(I-1,J,bi,bj))
179 6 *DELXR/CSUICE(I,J,bi,bj)
180 6
181 7 -ETAMEAN*TWO*TNGICE(I,J,bi,bj)*(VICEC(I+1,J,bi,bj)
182 7 -VICEC(I-1,J,bi,bj))*DELXR/CSUICE(I,J,bi,bj)
183 END DO
184 END DO
185
186 ENDDO
187 ENDDO
188
189 C NOW DO ITERATION
190 100 CONTINUE
191
192 DO bj=myByLo(myThid),myByHi(myThid)
193 DO bi=myBxLo(myThid),myBxHi(myThid)
194 C NOW SET U(3)=U(1)
195 DO J=1,sNy
196 DO I=1,sNx
197 UICE(I,J,3,bi,bj)=UICE(I,J,1,bi,bj)
198 END DO
199 END DO
200
201 DO 1200 J=1,sNy
202 DO I=1,sNx
203 DELY2=HALF/(DYUICE(I,J,bi,bj)*DYUICE(I,J,bi,bj))
204 DELYR=HALF/(DYUICE(I,J,bi,bj)*RADIUS)
205 ETAMEAN=QUART*(ETA(I,J+1,bi,bj)+ETA(I+1,J+1,bi,bj)
206 & +ETA(I,J,bi,bj)+ETA(I+1,J,bi,bj))
207 URT(I)=FXY(I,J)
208 1 +(ETA(I,J+1,bi,bj)+ETA(I+1,J+1,bi,bj))
209 1 *UICE(I,J+1,1,bi,bj)*DELY2
210 2 +(ETA(I,J,bi,bj)+ETA(I+1,J,bi,bj))
211 2 *UICE(I,J-1,1,bi,bj)*DELY2
212 3 +ETAMEAN*DELYR*(UICE(I,J+1,1,bi,bj)*TNGICE(I,J+1,bi,bj)
213 3 -UICE(I,J-1,1,bi,bj)*TNGICE(I,J-1,bi,bj))
214 4 -ETAMEAN*DELYR*TWO*TNGICE(I,J,bi,bj)*(UICE(I,J+1,1,bi,bj)
215 4 -UICE(I,J-1,1,bi,bj))
216 URT(I)=URT(I)*UVM(I,J,bi,bj)
217 END DO
218
219 DO I=1,sNx
220 CUU(I)=CU(I,J)
221 END DO
222 URT(1)=URT(1)/BU(1,J)
223 DO I=2,sNx
224 IM=I-1
225 CUU(I)=CUU(I)/(BU(I,J)-AU(I,J)*CUU(IM))
226 URT(I)=(URT(I)-AU(I,J)*URT(IM))/(BU(I,J)-AU(I,J)*CUU(IM))
227 END DO
228 DO I=1,sNx-1
229 J1=sNx-I
230 J2=J1+1
231 URT(J1)=URT(J1)-CUU(J1)*URT(J2)
232 END DO
233 DO I=1,sNx
234 UICE(I,J,1,bi,bj)=UICE(I,J,3,bi,bj)
235 & +WFAU*(URT(I)-UICE(I,J,3,bi,bj))
236 END DO
237
238 1200 CONTINUE
239
240 ENDDO
241 ENDDO
242
243 ICOUNT=ICOUNT+1
244 IF(ICOUNT.GT.1500) GO TO 201
245 202 CONTINUE
246 C NOW CHECK MAX ERROR
247 C The following loop has to be global operation
248 S1=ZERO
249 C FORM ERROR MATRIX
250 DO bj=myByLo(myThid),myByHi(myThid)
251 DO bi=myBxLo(myThid),myBxHi(myThid)
252 DO J=1,sNy
253 DO I=1,sNx
254 UERR(I,J,bi,bj)=(UICE(I,J,1,bi,bj)-UICE(I,J,3,bi,bj))
255 & *UVM(I,J,bi,bj)
256 S1=MAX(ABS(UERR(I,J,bi,bj)),S1)
257 END DO
258 END DO
259 ENDDO
260 ENDDO
261 _GLOBAL_MAX_R8( S1, myThid )
262 C NOW FIND ERROR
263 IF(S1.LT.LSR_ERROR) GO TO 200
264 CALL SEAICE_EXCH ( UICE, myThid )
265 GO TO 100
266 201 CONTINUE
267 PRINT 11
268 200 CONTINUE
269 write(*,'(A,I6,1P2E22.14)')' lsr iters, error = ',ICOUNT,S1
270
271 C NOW FOR VICE
272 DO bj=myByLo(myThid),myByHi(myThid)
273 DO bi=myBxLo(myThid),myBxHi(myThid)
274
275 DO J=1,sNy
276 DO I=1,sNx
277 DELX2=HALF/(DXUICE(I,J,bi,bj)*DXUICE(I,J,bi,bj))
278 DELY2=HALF/(DYUICE(I,J,bi,bj)*DYUICE(I,J,bi,bj))
279 DELYR=HALF/(DYUICE(I,J,bi,bj)*RADIUS)
280 RADIUS2=ONE/(RADIUS*RADIUS)
281 ETAMEAN=QUART*(ETA(I,J+1,bi,bj)+ETA(I+1,J+1,bi,bj)
282 & +ETA(I,J,bi,bj)+ETA(I+1,J,bi,bj))
283 ZETAMEAN=QUART*(ZETA(I,J+1,bi,bj)+ZETA(I+1,J+1,bi,bj)
284 & +ZETA(I,J,bi,bj)+ZETA(I+1,J,bi,bj))
285 AA1=ETA(I,J+1,bi,bj)+ZETA(I,J+1,bi,bj)+ETA(I+1,J+1,bi,bj)
286 & +ZETA(I+1,J+1,bi,bj)
287 AA2=ETA(I,J,bi,bj)+ZETA(I,J,bi,bj)+ETA(I+1,J,bi,bj)
288 & +ZETA(I+1,J,bi,bj)
289 AA3=(ETA(I+1,J,bi,bj)/CSUICE(I,J,bi,bj)+ETA(I+1,J+1,bi,bj)
290 & /CSUICE(I,J,bi,bj))/CSUICE(I,J,bi,bj)
291 AA4=(ETA(I,J,bi,bj)/CSUICE(I,J,bi,bj)+ETA(I,J+1,bi,bj)
292 & /CSUICE(I,J,bi,bj))/CSUICE(I,J,bi,bj)
293 AA5=((ZETA(I,J+1,bi,bj)-ETA(I,J+1,bi,bj))
294 & +(ZETA(I+1,J+1,bi,bj)-ETA(I+1,J+1,bi,bj))
295 & -(ZETA(I,J,bi,bj)-ETA(I,J,bi,bj))-(ZETA(I+1,J,bi,bj)
296 & -ETA(I+1,J,bi,bj)))*TNGICE(I,J,bi,bj)
297 AA6=TWO*ETAMEAN*TNGICE(I,J,bi,bj)*TNGICE(I,J,bi,bj)
298
299 AV(I,J)=(-AA2*DELY2
300 & -(ZETAMEAN-ETAMEAN)*TNGICE(I,J-1,bi,bj)*DELYR
301 & -ETAMEAN*TWO*TNGICE(I,J,bi,bj)*DELYR)*UVM(I,J,bi,bj)
302 BV(I,J)=((AA1+AA2)*DELY2+(AA3+AA4)*DELX2+AA5*DELYR+AA6*RADIUS2
303 & +AMASS(I,J,bi,bj)/DELTAT+DRAGS(I,J,bi,bj))
304 & *UVM(I,J,bi,bj)+(ONE-UVM(I,J,bi,bj))
305 CV(I,J)=(-AA1*DELY2
306 & +(ZETAMEAN-ETAMEAN)*TNGICE(I,J+1,bi,bj)*DELYR
307 & +ETAMEAN*TWO*TNGICE(I,J,bi,bj)*DELYR)*UVM(I,J,bi,bj)
308 END DO
309 END DO
310
311 DO I=1,sNx
312 AV(I,1)=ZERO
313 CV(I,sNy)=ZERO
314 CV(I,1)=CV(I,1)/BV(I,1)
315 END DO
316
317 DO J=1,sNy
318 DO I=1,sNx
319 DELXY=HALF/(DXUICE(I,J,bi,bj)*DYUICE(I,J,bi,bj))
320 DELXR=HALF/(DXUICE(I,J,bi,bj)*RADIUS)
321 DELY2=HALF/(DYUICE(I,J,bi,bj)*DYUICE(I,J,bi,bj))
322 DELYR=HALF/(DYUICE(I,J,bi,bj)*RADIUS)
323 ETAMEAN=QUART*(ETA(I,J+1,bi,bj)+ETA(I+1,J+1,bi,bj)
324 & +ETA(I,J,bi,bj)+ETA(I+1,J,bi,bj))
325 ZETAMEAN=QUART*(ZETA(I,J+1,bi,bj)+ZETA(I+1,J+1,bi,bj)
326 & +ZETA(I,J,bi,bj)+ZETA(I+1,J,bi,bj))
327 AA1=ETA(I,J+1,bi,bj)+ZETA(I,J+1,bi,bj)+ETA(I+1,J+1,bi,bj)
328 & +ZETA(I+1,J+1,bi,bj)
329 AA2=ETA(I,J,bi,bj)+ZETA(I,J,bi,bj)+ETA(I+1,J,bi,bj)
330 & +ZETA(I+1,J,bi,bj)
331
332 IF(J.EQ.1) THEN
333 AA3=(AA2*DELY2+(ZETAMEAN-ETAMEAN)*TNGICE(I,J-1,bi,bj)*DELYR
334 & +ETAMEAN*TWO*TNGICE(I,J,bi,bj)*DELYR)
335 & *VICE(I,J-1,1,bi,bj)*UVM(I,J,bi,bj)
336 ELSE IF(J.EQ.sNy) THEN
337 AA3=(AA1*DELY2-(ZETAMEAN-ETAMEAN)*TNGICE(I,J+1,bi,bj)*DELYR
338 & -ETAMEAN*TWO*TNGICE(I,J,bi,bj)*DELYR)
339 & *VICE(I,J+1,1,bi,bj)*UVM(I,J,bi,bj)
340 ELSE
341 AA3=ZERO
342 END IF
343
344 FXY(I,J)=AA3-DRAGA(I,J,bi,bj)*UICEC(I,J,bi,bj)
345 2 +FORCEY(I,J,bi,bj)
346 3 +(HALF*(UICEC(I+1,J,bi,bj)-UICEC(I-1,J,bi,bj))
347 3 *(ZETA(I,J+1,bi,bj)+ZETA(I+1,J+1,bi,bj)
348 3 -ZETA(I,J,bi,bj)-ZETA(I+1,J,bi,bj))*DELXY
349 3 /CSUICE(I,J,bi,bj)+HALF*ZETAMEAN*((UICEC(I+1,J+1,bi,bj)
350 3 -UICEC(I-1,J+1,bi,bj))/CSUICE(I,J+1,bi,bj)
351 3 -(UICEC(I+1,J-1,bi,bj)-UICEC(I-1,J-1,bi,bj))
352 3 /CSUICE(I,J-1,bi,bj))*DELXY)
353 3
354 4 -(HALF*(UICEC(I+1,J,bi,bj)-UICEC(I-1,J,bi,bj))
355 4 *(ETA(I,J+1,bi,bj)+ETA(I+1,J+1,bi,bj)
356 4 -ETA(I,J,bi,bj)-ETA(I+1,J,bi,bj))*DELXY/CSUICE(I,J,bi,bj)
357 4 +HALF*ETAMEAN*((UICEC(I+1,J+1,bi,bj)
358 4 -UICEC(I-1,J+1,bi,bj))/CSUICE(I,J+1,bi,bj)
359 4 -(UICEC(I+1,J-1,bi,bj)-UICEC(I-1,J-1,bi,bj))
360 4 /CSUICE(I,J-1,bi,bj))*DELXY)
361 4
362 5 +HALF*(ETA(I+1,J+1,bi,bj)*(UICEC(I+1,J+1,bi,bj)
363 5 +UICEC(I,J+1,bi,bj)
364 5 -UICEC(I+1,J,bi,bj)-UICEC(I,J,bi,bj))+ETA(I+1,J,bi,bj)
365 5 *(UICEC(I+1,J,bi,bj)
366 5 +UICEC(I,J,bi,bj)-UICEC(I+1,J-1,bi,bj)
367 5 -UICEC(I,J-1,bi,bj))+ETA(I,J+1,bi,bj)
368 5 *(UICEC(I,J,bi,bj)+UICEC(I-1,J,bi,bj)-UICEC(I,J+1,bi,bj)
369 5 -UICEC(I-1,J+1,bi,bj))
370 5 +ETA(I,J,bi,bj)*(UICEC(I,J-1,bi,bj)+UICEC(I-1,J-1,bi,bj)
371 5 -UICEC(I,J,bi,bj)
372 5 -UICEC(I-1,J,bi,bj)))*DELXY/CSUICE(I,J,bi,bj)
373 5
374 6 +(ETA(I+1,J+1,bi,bj)+ETA(I+1,J,bi,bj)-ETA(I,J,bi,bj)
375 6 -ETA(I,J+1,bi,bj))
376 6 *TNGICE(I,J,bi,bj)*UICEC(I,J,bi,bj)
377 6 *DELXR/CSUICE(I,J,bi,bj)
378 6 +ETAMEAN*TNGICE(I,J,bi,bj)*(UICEC(I+1,J,bi,bj)
379 6 -UICEC(I-1,J,bi,bj))*DELXR/CSUICE(I,J,bi,bj)
380 6
381 7 +ETAMEAN*TWO*TNGICE(I,J,bi,bj)*(UICEC(I+1,J,bi,bj)
382 7 -UICEC(I-1,J,bi,bj))*DELXR/CSUICE(I,J,bi,bj)
383
384 END DO
385 END DO
386
387 ENDDO
388 ENDDO
389
390 C NOW DO ITERATION
391 300 CONTINUE
392 C NOW SET U(3)=U(1)
393 DO bj=myByLo(myThid),myByHi(myThid)
394 DO bi=myBxLo(myThid),myBxHi(myThid)
395
396 DO J=1,sNy
397 DO I=1,sNx
398 VICE(I,J,3,bi,bj)=VICE(I,J,1,bi,bj)
399 END DO
400 END DO
401
402 DO 1300 I=1,sNx
403 DO J=1,sNy
404 DELX2=HALF/(DXUICE(I,J,bi,bj)*DXUICE(I,J,bi,bj))
405 VRT(J)=FXY(I,J)
406 6 +((ETA(I+1,J,bi,bj)/CSUICE(I,J,bi,bj)+ETA(I+1,J+1,bi,bj)
407 6 /CSUICE(I,J,bi,bj))*VICE(I+1,J,1,bi,bj)*DELX2
408 7 +(ETA(I,J,bi,bj)/CSUICE(I,J,bi,bj)+ETA(I,J+1,bi,bj)
409 7 /CSUICE(I,J,bi,bj))*VICE(I-1,J,1,bi,bj)*DELX2)
410 7 /CSUICE(I,J,bi,bj)
411 VRT(J)=VRT(J)*UVM(I,J,bi,bj)
412 END DO
413
414 DO J=1,sNy
415 CVV(J)=CV(I,J)
416 END DO
417 VRT(1)=VRT(1)/BV(I,1)
418 DO J=2,sNy
419 JM=J-1
420 CVV(J)=CVV(J)/(BV(I,J)-AV(I,J)*CVV(JM))
421 VRT(J)=(VRT(J)-AV(I,J)*VRT(JM))/(BV(I,J)-AV(I,J)*CVV(JM))
422 END DO
423 DO J=1,sNy-1
424 J1=sNy-J
425 J2=J1+1
426 VRT(J1)=VRT(J1)-CVV(J1)*VRT(J2)
427 END DO
428 DO J=1,sNy
429 VICE(I,J,1,bi,bj)=VICE(I,J,3,bi,bj)
430 & +WFAV*(VRT(J)-VICE(I,J,3,bi,bj))
431 END DO
432 1300 CONTINUE
433
434 ENDDO
435 ENDDO
436
437 ICOU=ICOU+1
438 IF(ICOU.GT.1500) GO TO 301
439 C NOW CHECK MAX ERROR
440 C The following loop has to be global operation
441 S2=ZERO
442 C FORM ERROR MATRIX
443 DO bj=myByLo(myThid),myByHi(myThid)
444 DO bi=myBxLo(myThid),myBxHi(myThid)
445 DO J=1,sNy
446 DO I=1,sNx
447 UERR(I,J,bi,bj)=(VICE(I,J,1,bi,bj)-VICE(I,J,3,bi,bj))
448 & *UVM(I,J,bi,bj)
449 S2=MAX(ABS(UERR(I,J,bi,bj)),S2)
450 END DO
451 END DO
452 ENDDO
453 ENDDO
454 _GLOBAL_MAX_R8( S2, myThid )
455 C NOW FIND ERROR
456 IF(S2.LT.LSR_ERROR) GO TO 400
457 CALL SEAICE_EXCH( VICE, myThid )
458 GO TO 300
459 301 CONTINUE
460 PRINT 121
461 11 FORMAT(1X,'NO CONVERGENCE FOR U AFTER 1500 ITERATIONS')
462 121 FORMAT(1X,'NO CONVERGENCE FOR V AFTER 1500 ITERATIONS')
463 400 CONTINUE
464 write(*,'(A,I6,1P2E22.14)')' lsr iters, error = ',ICOU,S2
465
466 C NOW END
467 C NOW MAKE COROLIS TERM IMPLICIT
468 DO bj=myByLo(myThid),myByHi(myThid)
469 DO bi=myBxLo(myThid),myBxHi(myThid)
470 DO J=1,sNy
471 DO I=1,sNx
472 UICE(I,J,1,bi,bj)=UICE(I,J,1,bi,bj)*UVM(I,J,bi,bj)
473 VICE(I,J,1,bi,bj)=VICE(I,J,1,bi,bj)*UVM(I,J,bi,bj)
474 END DO
475 END DO
476 ENDDO
477 ENDDO
478 CALL SEAICE_EXCH( UICE, myThid )
479 CALL SEAICE_EXCH( VICE, myThid )
480
481 #endif /* SEAICE_ALLOW_DYNAMICS */
482 #endif /* ALLOW_SEAICE */
483
484 RETURN
485 END
  ViewVC Help
Powered by ViewVC 1.1.22