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

Contents of /MITgcm/pkg/seaice/dynsolver.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, 3 months ago) by dimitri
Branch: MAIN
CVS Tags: checkpoint48i_post, checkpoint50, checkpoint50b_pre, checkpoint48f_post, checkpoint48h_post, checkpoint50a_post, checkpoint49, checkpoint48g_post, checkpoint50b_post
Changes since 1.4: +8 -27 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 dynsolver( myTime, myIter, myThid )
7 C /==========================================================\
8 C | SUBROUTINE dynsolver |
9 C | o Ice dynamics solver |
10 C |==========================================================|
11 C \==========================================================/
12 IMPLICIT NONE
13
14 C === Global variables ===
15 #include "SIZE.h"
16 #include "EEPARAMS.h"
17 #include "PARAMS.h"
18 #include "FFIELDS.h"
19 #include "SEAICE.h"
20 #include "SEAICE_GRID.h"
21 #include "SEAICE_PARAMS.h"
22 #include "SEAICE_FFIELDS.h"
23
24 #ifdef ALLOW_AUTODIFF_TAMC
25 # include "tamc.h"
26 #endif
27
28 C === Routine arguments ===
29 C myTime - Simulation time
30 C myIter - Simulation timestep number
31 C myThid - Thread no. that called this routine.
32 _RL myTime
33 INTEGER myIter
34 INTEGER myThid
35 CEndOfInterface
36
37 #ifdef ALLOW_SEAICE
38
39 C === Local variables ===
40 C i,j,bi,bj - Loop counters
41
42 INTEGER i, j, bi, bj, kii
43 _RL DWAT, DAIR, RHOICE, RHOAIR, SINWIN, COSWIN, SINWAT, COSWAT
44 _RL GRAV, ECCEN, ECM2, GMIN, RADIUS, DELT1, DELT2, PSTAR, AAA
45
46 _RL PRESS (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
47 _RL DAIRN (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
48 _RL DWATN (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
49 _RL FORCEX0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
50 _RL FORCEY0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
51 _RL E11 (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
52 _RL E22 (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
53 _RL E12 (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
54 _RL COR_ICE (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
55 _RL ZMAX (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
56 _RL ZMIN (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
57
58 C-- FIRST SET UP BASIC CONSTANTS
59 DWAT=0.59 _d 0
60 DAIR=0.01462 _d 0
61 RHOICE=0.91 _d +03
62 RHOAIR=1.3 _d 0
63 GRAV=9.832 _d 0
64 ECCEN=TWO
65 ECM2=ONE/(ECCEN**2)
66 GMIN=1.0 _d -20
67 RADIUS=6370. _d 3
68 PSTAR=SEAICE_strength
69
70 C-- 25 DEG GIVES SIN EQUAL TO 0.4226
71 SINWIN=0.4226 _d 0
72 COSWIN=0.9063 _d 0
73 SINWAT=0.4226 _d 0
74 COSWAT=0.9063 _d 0
75
76 C-- Do not introduce turning angle
77 SINWIN=ZERO
78 COSWIN=ONE
79 SINWAT=ZERO
80 COSWAT=ONE
81
82 C-- NOW SET UP MASS PER UNIT AREA AND CORIOLIS TERM
83 DO bj=myByLo(myThid),myByHi(myThid)
84 DO bi=myBxLo(myThid),myBxHi(myThid)
85 DO j=1,sNy
86 DO i=1,sNx
87 AMASS(I,J,bi,bj)=RHOICE*QUART*(HEFF(i,j,1,bi,bj)
88 & +HEFF(i+1,j,1,bi,bj)
89 & +HEFF(i,j+1,1,bi,bj)
90 & +HEFF(i+1,j+1,1,bi,bj))
91 COR_ICE(I,J,bi,bj)=AMASS(I,J,bi,bj)
92 & *TWO*OMEGA*SINEICE(I,J,bi,bj)
93 ENDDO
94 ENDDO
95 ENDDO
96 ENDDO
97
98 C-- NOW SET UP FORCING FIELDS
99
100 IF (SEAICEwindOnCgrid) THEN
101 C-- Wind stress computed here from wind on South-West C-grid
102 DO bj=myByLo(myThid),myByHi(myThid)
103 DO bi=myBxLo(myThid),myBxHi(myThid)
104 DO j=1,sNy
105 DO i=1,sNx
106 AAA=SQRT((HALF*(UWIND(I+1,J,bi,bj)+UWIND(I+1,J+1,bi,bj)))**2
107 & +(HALF*(VWIND(I,J+1,bi,bj)+VWIND(I+1,J+1,bi,bj)))**2)
108 DAIRN(I,J,bi,bj)=RHOAIR*SEAICE_drag*
109 & (2.70 _d 0+0.142 _d 0*AAA+0.0764 _d 0*AAA*AAA)
110 FORCEX(I,J,bi,bj)=DAIRN(I,J,bi,bj)*
111 & (COSWIN*HALF*(UWIND(I+1,J,bi,bj)+UWIND(I+1,J+1,bi,bj))
112 & -SINWIN*HALF*(VWIND(I,J+1,bi,bj)+VWIND(I+1,J+1,bi,bj)))
113 FORCEY(I,J,bi,bj)=DAIRN(I,J,bi,bj)*
114 & (SINWIN*HALF*(UWIND(I+1,J,bi,bj)+UWIND(I+1,J+1,bi,bj))
115 & +COSWIN*HALF*(VWIND(I,J+1,bi,bj)+VWIND(I+1,J+1,bi,bj)))
116 ENDDO
117 ENDDO
118 ENDDO
119 ENDDO
120
121 ELSE
122 C-- Wind stress computed here from wind on South-West B-grid
123 DO bj=myByLo(myThid),myByHi(myThid)
124 DO bi=myBxLo(myThid),myBxHi(myThid)
125 DO j=1,sNy
126 DO i=1,sNx
127 AAA=SQRT(UWIND(I,J,bi,bj)**2+VWIND(I,J,bi,bj)**2)
128 DAIRN(I,J,bi,bj)=RHOAIR*SEAICE_drag*
129 & (2.70 _d 0+0.142 _d 0*AAA+0.0764 _d 0*AAA*AAA)
130 FORCEX(I,J,bi,bj)=DAIRN(I,J,bi,bj)*(COSWIN*UWIND(I,J,bi,bj)
131 & -SINWIN*VWIND(I,J,bi,bj))
132 FORCEY(I,J,bi,bj)=DAIRN(I,J,bi,bj)*(SINWIN*UWIND(I,J,bi,bj)
133 & +COSWIN*VWIND(I,J,bi,bj))
134 ENDDO
135 ENDDO
136 ENDDO
137 ENDDO
138 ENDIF
139
140 DO bj=myByLo(myThid),myByHi(myThid)
141 DO bi=myBxLo(myThid),myBxHi(myThid)
142 DO j=1,sNy
143 DO i=1,sNx
144 C-- STORE WIND ONLY STRESS
145 WINDX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
146 WINDY(I,J,bi,bj)=FORCEY(I,J,bi,bj)
147 C-- NOW ADD IN TILT
148 FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
149 & -COR_ICE(I,J,bi,bj)*GWATY(I,J,bi,bj)
150 FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)
151 & +COR_ICE(I,J,bi,bj)*GWATX(I,J,bi,bj)
152 C-- NOW SET UP ICE PRESSURE AND VISCOSITIES
153 PRESS(I,J,bi,bj)=PSTAR*HEFF(I,J,1,bi,bj)
154 & *EXP(-20.0 _d 0*(ONE-AREA(I,J,1,bi,bj)))
155 ZMAX(I,J,bi,bj)=(5.0 _d +12/(2.0 _d +04))*PRESS(I,J,bi,bj)
156 ZMIN(I,J,bi,bj)=4.0 _d +08
157 PRESS(I,J,bi,bj)=PRESS(I,J,bi,bj)*HEFFM(I,J,bi,bj)
158 ENDDO
159 ENDDO
160 ENDDO
161 ENDDO
162
163 #ifdef SEAICE_ALLOW_DYNAMICS
164 IF ( SEAICEuseDYNAMICS ) THEN
165
166 C-- Update overlap regions for PRESS
167 _EXCH_XY_R8(PRESS, myThid)
168
169 DO bj=myByLo(myThid),myByHi(myThid)
170 DO bi=myBxLo(myThid),myBxHi(myThid)
171 DO j=1,sNy
172 DO i=1,sNx
173 C NOW CALCULATE PRESSURE FORCE AND ADD TO EXTERNAL FORCE
174 FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
175 & -(QUART/(DXUICE(I,J,bi,bj)*CSUICE(I,J,bi,bj)))
176 & *(PRESS(I+1,J,bi,bj)+PRESS(I+1,J+1,bi,bj)
177 & -PRESS(I,J,bi,bj)-PRESS(I,J+1,bi,bj))
178 FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)-QUART/DYUICE(I,J,bi,bj)
179 & *(PRESS(I,J+1,bi,bj)+PRESS(I+1,J+1,bi,bj)
180 & -PRESS(I,J,bi,bj)-PRESS(I+1,J,bi,bj))
181 C NOW KEEP FORCEX0
182 FORCEX0(I,J,bi,bj)=FORCEX(I,J,bi,bj)
183 FORCEY0(I,J,bi,bj)=FORCEY(I,J,bi,bj)
184 ENDDO
185 ENDDO
186 ENDDO
187 ENDDO
188
189 C DO PSEUDO-TIMESTEPS TO OBTAIN AN ACCURATE VISCOUS-PLASTIC SOLUTION
190 C 10 PSEUDO-TIMESTEPS OR MORE ARE SUGGESTED FOR HIGH-RESOLUTION (~10KM)
191 C 1 PSEUDO-TIMESTEP CAN BE USED FOR LOW-RESOLUTION GLOBAL MODELING
192 C NPSEUDO is now set in data.seaice input file
193 C TIMESTEP FOR PSEUDO-TIMESTEPPING
194 SEAICE_DT = DELTAT/NPSEUDO
195
196 crg what about DWAIN,DRAGS,DRAGA,ETA,ZETA
197
198 crg later c$taf loop = iteration uice,vice
199
200 DO 5000 KII=1,NPSEUDO
201
202 c$taf store uice,vice = comlev1_seaice_ds,
203 c$taf& key = kii + (ikey_dynamics-1)*NPSEUDO
204
205 C NOW DO PREDICTOR TIME STEP
206 DO bj=myByLo(myThid),myByHi(myThid)
207 DO bi=myBxLo(myThid),myBxHi(myThid)
208 DO j=1-OLy,sNy+OLy
209 DO i=1-OLx,sNx+OLx
210 UICE(I,J,2,bi,bj)=UICE(I,J,1,bi,bj)
211 VICE(I,J,2,bi,bj)=VICE(I,J,1,bi,bj)
212 UICEC(I,J,bi,bj)=UICE(I,J,1,bi,bj)
213 VICEC(I,J,bi,bj)=VICE(I,J,1,bi,bj)
214 ENDDO
215 ENDDO
216 ENDDO
217 ENDDO
218
219 DO bj=myByLo(myThid),myByHi(myThid)
220 DO bi=myBxLo(myThid),myBxHi(myThid)
221 DO j=1,sNy
222 DO i=1,sNx
223 C NOW SET UP NON-LINEAR WATER DRAG, FORCEX, FORCEY
224 DWATN(I,J,bi,bj)=SEAICE_waterDrag*SQRT((UICE(I,J,1,bi,bj)
225 & -GWATX(I,J,bi,bj))**2
226 & +(VICE(I,J,1,bi,bj)-GWATY(I,J,bi,bj))**2)
227 DWATN(I,J,bi,bj)=MAX(DWATN(I,J,bi,bj),QUART)
228 C NOW SET UP SYMMETTRIC DRAG
229 DRAGS(I,J,bi,bj)=DWATN(I,J,bi,bj)*COSWAT
230 C NOW SET UP ANTI SYMMETTRIC DRAG PLUS CORIOLIS
231 DRAGA(I,J,bi,bj)=DWATN(I,J,bi,bj)*SINWAT+COR_ICE(I,J,bi,bj)
232 C NOW ADD IN CURRENT FORCE
233 FORCEX(I,J,bi,bj)=FORCEX0(I,J,bi,bj)+DWATN(I,J,bi,bj)
234 & *(COSWAT*GWATX(I,J,bi,bj)
235 & -SINWAT*GWATY(I,J,bi,bj))
236 FORCEY(I,J,bi,bj)=FORCEY0(I,J,bi,bj)+DWATN(I,J,bi,bj)
237 & *(SINWAT*GWATX(I,J,bi,bj)
238 & +COSWAT*GWATY(I,J,bi,bj))
239 ENDDO
240 ENDDO
241 ENDDO
242 ENDDO
243
244 DO bj=myByLo(myThid),myByHi(myThid)
245 DO bi=myBxLo(myThid),myBxHi(myThid)
246 DO j=1,sNy
247 DO i=1,sNx
248 C NOW EVALUATE STRAIN RATES
249 E11(I,J,bi,bj)=HALF/(DXTICE(I,J,bi,bj)*CSTICE(I,J,bi,bj))
250 & *(UICE(I,J,1,bi,bj)+UICE(I,J-1,1,bi,bj)
251 & -UICE(I-1,J,1,bi,bj)-UICE(I-1,J-1,1,bi,bj))
252 & -QUART*(VICE(I,J,1,bi,bj)+VICE(I-1,J,1,bi,bj)
253 & +VICE(I-1,J-1,1,bi,bj)+VICE(I,J-1,1,bi,bj))
254 & *TNGTICE(I,J,bi,bj)/RADIUS
255 E22(I,J,bi,bj)=HALF/DYTICE(I,J,bi,bj)
256 & *(VICE(I,J,1,bi,bj)+VICE(I-1,J,1,bi,bj)
257 & -VICE(I,J-1,1,bi,bj)-VICE(I-1,J-1,1,bi,bj))
258 E12(I,J,bi,bj)=HALF*(HALF/DYTICE(I,J,bi,bj)
259 & *(UICE(I,J,1,bi,bj)+UICE(I-1,J,1,bi,bj)
260 & -UICE(I,J-1,1,bi,bj)-UICE(I-1,J-1,1,bi,bj))
261 & +HALF/(DXTICE(I,J,bi,bj)*CSTICE(I,J,bi,bj))
262 & *(VICE(I,J,1,bi,bj)+VICE(I,J-1,1,bi,bj)
263 & -VICE(I-1,J,1,bi,bj)-VICE(I-1,J-1,1,bi,bj))
264 & +QUART*(UICE(I,J,1,bi,bj)+UICE(I-1,J,1,bi,bj)
265 & +UICE(I-1,J-1,1,bi,bj)+UICE(I,J-1,1,bi,bj))
266 & *TNGTICE(I,J,bi,bj)/RADIUS)
267 C NOW EVALUATE VISCOSITIES
268 DELT1=(E11(I,J,bi,bj)**2+E22(I,J,bi,bj)**2)*(ONE+ECM2)
269 & +4.0 _d 0*ECM2*E12(I,J,bi,bj)**2
270 1 +TWO*E11(I,J,bi,bj)*E22(I,J,bi,bj)*(ONE-ECM2)
271 DELT2=SQRT(DELT1)
272 DELT2=MAX(GMIN,DELT2)
273 ZETA(I,J,bi,bj)=HALF*PRESS(I,J,bi,bj)/DELT2
274 C NOW PUT MIN AND MAX VISCOSITIES IN
275 ZETA(I,J,bi,bj)=MIN(ZMAX(I,J,bi,bj),ZETA(I,J,bi,bj))
276 ZETA(I,J,bi,bj)=MAX(ZMIN(I,J,bi,bj),ZETA(I,J,bi,bj))
277 C NOW SET VISCOSITIES TO ZERO AT HEFFMFLOW PTS
278 ZETA(I,J,bi,bj)=ZETA(I,J,bi,bj)*HEFFM(I,J,bi,bj)
279 ETA(I,J,bi,bj)=ECM2*ZETA(I,J,bi,bj)
280 ENDDO
281 ENDDO
282 ENDDO
283 ENDDO
284
285 C-- Update overlap regions
286 _EXCH_XY_R8(ETA, myThid)
287 _EXCH_XY_R8(ZETA, myThid)
288
289 C NOW ADI SCHEME (ZHANG-J/ROTHROCK 1999,bi,bj)
290 IF ( SEAICEuseLSR ) THEN
291 CALL LSR( myThid )
292 ELSE
293 CALL ADI( myThid )
294 ENDIF
295
296 C NOW DO MODIFIED EULER STEP
297 DO bj=myByLo(myThid),myByHi(myThid)
298 DO bi=myBxLo(myThid),myBxHi(myThid)
299 DO j=1-OLy,sNy+OLy
300 DO i=1-OLx,sNx+OLx
301 UICE(I,J,1,bi,bj)=HALF*(UICE(I,J,1,bi,bj)+UICE(I,J,2,bi,bj))
302 VICE(I,J,1,bi,bj)=HALF*(VICE(I,J,1,bi,bj)+VICE(I,J,2,bi,bj))
303 UICEC(I,J,bi,bj)=UICE(I,J,1,bi,bj)
304 VICEC(I,J,bi,bj)=VICE(I,J,1,bi,bj)
305 ENDDO
306 ENDDO
307 ENDDO
308 ENDDO
309
310 DO bj=myByLo(myThid),myByHi(myThid)
311 DO bi=myBxLo(myThid),myBxHi(myThid)
312 DO j=1,sNy
313 DO i=1,sNx
314 C NOW SET UP NON-LINEAR WATER DRAG
315 DWATN(I,J,bi,bj)=SEAICE_waterDrag*SQRT((UICE(I,J,1,bi,bj)
316 & -GWATX(I,J,bi,bj))**2
317 & +(VICE(I,J,1,bi,bj)-GWATY(I,J,bi,bj))**2)
318 DWATN(I,J,bi,bj)=MAX(DWATN(I,J,bi,bj),QUART)
319 C NOW SET UP SYMMETTRIC DRAG
320 DRAGS(I,J,bi,bj)=DWATN(I,J,bi,bj)*COSWAT
321 C NOW SET UP ANTI SYMMETTRIC DRAG PLUS CORIOLIS
322 DRAGA(I,J,bi,bj)=DWATN(I,J,bi,bj)*SINWAT+COR_ICE(I,J,bi,bj)
323 C NOW ADD IN CURRENT FORCE
324 FORCEX(I,J,bi,bj)=FORCEX0(I,J,bi,bj)+DWATN(I,J,bi,bj)
325 & *(COSWAT*GWATX(I,J,bi,bj)
326 & -SINWAT*GWATY(I,J,bi,bj))
327 FORCEY(I,J,bi,bj)=FORCEY0(I,J,bi,bj)+DWATN(I,J,bi,bj)
328 & *(SINWAT*GWATX(I,J,bi,bj)
329 & +COSWAT*GWATY(I,J,bi,bj))
330 ENDDO
331 ENDDO
332 ENDDO
333 ENDDO
334
335 DO bj=myByLo(myThid),myByHi(myThid)
336 DO bi=myBxLo(myThid),myBxHi(myThid)
337 DO j=1,sNy
338 DO i=1,sNx
339 C NOW EVALUATE STRAIN RATES
340 E11(I,J,bi,bj)=HALF/(DXTICE(I,J,bi,bj)*CSTICE(I,J,bi,bj))
341 & *(UICE(I,J,1,bi,bj)+UICE(I,J-1,1,bi,bj)
342 & -UICE(I-1,J,1,bi,bj)-UICE(I-1,J-1,1,bi,bj))
343 & -QUART*(VICE(I,J,1,bi,bj)+VICE(I-1,J,1,bi,bj)
344 & +VICE(I-1,J-1,1,bi,bj)+VICE(I,J-1,1,bi,bj))
345 & *TNGTICE(I,J,bi,bj)/RADIUS
346 E22(I,J,bi,bj)=HALF/DYTICE(I,J,bi,bj)
347 & *(VICE(I,J,1,bi,bj)+VICE(I-1,J,1,bi,bj)
348 & -VICE(I,J-1,1,bi,bj)-VICE(I-1,J-1,1,bi,bj))
349 E12(I,J,bi,bj)=HALF*(HALF/DYTICE(I,J,bi,bj)
350 & *(UICE(I,J,1,bi,bj)+UICE(I-1,J,1,bi,bj)
351 & -UICE(I,J-1,1,bi,bj)-UICE(I-1,J-1,1,bi,bj))
352 & +HALF/(DXTICE(I,J,bi,bj)*CSTICE(I,J,bi,bj))
353 & *(VICE(I,J,1,bi,bj)+VICE(I,J-1,1,bi,bj)
354 & -VICE(I-1,J,1,bi,bj)-VICE(I-1,J-1,1,bi,bj))
355 & +QUART*(UICE(I,J,1,bi,bj)+UICE(I-1,J,1,bi,bj)
356 & +UICE(I-1,J-1,1,bi,bj)+UICE(I,J-1,1,bi,bj))
357 & *TNGTICE(I,J,bi,bj)/RADIUS)
358 C NOW EVALUATE VISCOSITIES
359 DELT1=(E11(I,J,bi,bj)**2+E22(I,J,bi,bj)**2)*(ONE+ECM2)
360 & +4. _d 0*ECM2*E12(I,J,bi,bj)**2
361 1 +TWO*E11(I,J,bi,bj)*E22(I,J,bi,bj)*(ONE-ECM2)
362 DELT2=SQRT(DELT1)
363 DELT2=MAX(GMIN,DELT2)
364 ZETA(I,J,bi,bj)=HALF*PRESS(I,J,bi,bj)/DELT2
365 C NOW PUT MIN AND MAX VISCOSITIES IN
366 ZETA(I,J,bi,bj)=MIN(ZMAX(I,J,bi,bj),ZETA(I,J,bi,bj))
367 ZETA(I,J,bi,bj)=MAX(ZMIN(I,J,bi,bj),ZETA(I,J,bi,bj))
368 C NOW SET VISCOSITIES TO ZERO AT HEFFMFLOW PTS
369 ZETA(I,J,bi,bj)=ZETA(I,J,bi,bj)*HEFFM(I,J,bi,bj)
370 ETA(I,J,bi,bj)=ECM2*ZETA(I,J,bi,bj)
371 ENDDO
372 ENDDO
373 ENDDO
374 ENDDO
375
376 C-- Update overlap regions
377 _EXCH_XY_R8(ETA, myThid)
378 _EXCH_XY_R8(ZETA, myThid)
379
380 C GET READY FOR SECOND CALL OF ADI
381 DO bj=myByLo(myThid),myByHi(myThid)
382 DO bi=myBxLo(myThid),myBxHi(myThid)
383 DO j=1-OLy,sNy+OLy
384 DO i=1-OLx,sNx+OLx
385 UICE(I,J,2,bi,bj)=UICEC(I,J,bi,bj)
386 VICE(I,J,2,bi,bj)=VICEC(I,J,bi,bj)
387 ENDDO
388 ENDDO
389 ENDDO
390 ENDDO
391
392 C NOW ADI SCHEME (ZHANG-J/ROTHROCK 1999)
393 IF ( SEAICEuseLSR ) THEN
394 CALL LSR( myThid )
395 ELSE
396 CALL ADI( myThid )
397 ENDIF
398
399 5000 CONTINUE
400
401 c$taf store uice,vice = comlev1, key=ikey_dynamics
402
403 ENDIF
404 #endif /* SEAICE_ALLOW_DYNAMICS */
405
406 C Calculate ocean surface stress
407 CALL OSTRES ( DWATN, COR_ICE, myThid )
408
409 #ifdef SEAICE_ALLOW_DYNAMICS
410 IF ( SEAICEuseDYNAMICS ) THEN
411
412 c Put a cap on ice velocity
413 c limit velocity to 0.40 m s-1 to avoid potential CFL violations
414 c in open water areas (drift of zero thickness ice)
415 DO bj=myByLo(myThid),myByHi(myThid)
416 DO bi=myBxLo(myThid),myBxHi(myThid)
417 DO j=1-OLy,sNy+OLy
418 DO i=1-OLx,sNx+OLx
419 #ifdef SEAICE_DEBUG
420 c write(*,'(2i4,2i2,f7.1,7f12.3)')
421 c & i,j,bi,bj,UVM(I,J,bi,bj),amass(i,j,bi,bj)
422 c & ,gwatx(I,J,bi,bj),gwaty(i,j,bi,bj)
423 c & ,forcex(I,J,bi,bj),forcey(i,j,bi,bj)
424 c & ,uice(i,j,1,bi,bj)
425 c & ,vice(i,j,1,bi,bj)
426 #endif /* SEAICE_DEBUG */
427 UICE(i,j,1,bi,bj)=min(UICE(i,j,1,bi,bj),0.40 _d +00)
428 VICE(i,j,1,bi,bj)=min(VICE(i,j,1,bi,bj),0.40 _d +00)
429 UICE(i,j,1,bi,bj)=max(UICE(i,j,1,bi,bj),-0.40 _d +00)
430 VICE(i,j,1,bi,bj)=max(VICE(i,j,1,bi,bj),-0.40 _d +00)
431 ENDDO
432 ENDDO
433 ENDDO
434 ENDDO
435
436 ENDIF
437 #endif /* SEAICE_ALLOW_DYNAMICS */
438
439 #endif /* ALLOW_SEAICE */
440
441 RETURN
442 END
  ViewVC Help
Powered by ViewVC 1.1.22