/[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.3 - (show annotations) (download)
Thu Dec 5 08:43:02 2002 UTC (21 years, 7 months ago) by dimitri
Branch: MAIN
CVS Tags: checkpoint47e_post, checkpoint47c_post, checkpoint47d_pre, checkpoint47d_post, branch-exfmods-tag, checkpoint47b_post
Branch point for: branch-exfmods-curt
Changes since 1.2: +77 -66 lines 
checkpoint47b_post
Merging from release1_p9:
o pkg/seaice
  - removed GOTO's and added taf directives
  - double precision constants to reduce the g77 (Linux)
    to F77 (SGI) differences reported in release1_p8
o tools/genmake
  - added SGI options
o verification/testscript
  - updated to that of checkpoint47a_post
o verification/global_ocean.90x40x15/input/eedata
  - modified for SGI f77 compatibility
o verification/lab_sea
  - added description of sea-ice model
  - added missing matlab routines
  - added test of thermodynamics parallelization
Modified Files:
   doc/tag-index pkg/seaice/SEAICE_FFIELDS.h
   pkg/seaice/SEAICE_PARAMS.h pkg/seaice/adi.F
   pkg/seaice/advect.F pkg/seaice/budget.F pkg/seaice/diffus.F
   pkg/seaice/dynsolver.F pkg/seaice/groatb.F pkg/seaice/growth.F
   pkg/seaice/lsr.F pkg/seaice/ostres.F
   pkg/seaice/seaice_do_diags.F pkg/seaice/seaice_get_forcing.F
   pkg/seaice/seaice_init.F pkg/seaice/seaice_model.F
   pkg/seaice/seaice_readparms.F tools/genmake
   verification/global_ocean.90x40x15/input/eedata
   verification/lab_sea/README
   verification/lab_sea/matlab/lookat_exp1.m
   verification/lab_sea/matlab/lookat_exp2.m
   verification/lab_sea/matlab/lookat_exp3.m
   verification/lab_sea/matlab/lookat_exp4.m
   verification/lab_sea/matlab/lookat_exp5.m
   verification/lab_sea/matlab/lookat_exp6.m
   verification/lab_sea/results/AREAtave.0000000010.data
   verification/lab_sea/results/HEFFtave.0000000010.data
   verification/lab_sea/results/UICEtave.0000000010.data
   verification/lab_sea/results/VICEtave.0000000010.data
   verification/lab_sea/results/output.txt
Added Files:
   verification/lab_sea/seaice.ps
   verification/lab_sea/matlab/lookat_exp7.m
   verification/lab_sea/matlab/mmax.m
   verification/lab_sea/matlab/mypcolor.m
   verification/lab_sea/matlab/myquiver.m
   verification/lab_sea/matlab/readbin.m
   verification/lab_sea/matlab/wysiwyg.m
Removed Files:
   verification/lab_sea/code/KPP_OPTIONS.h
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 /* ALLOW_AUTODIFF_TAMC */
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,k,bi,bj - Loop counters
41
42 INTEGER i, j, k, bi, bj, kii
43 _RL DWAT, DAIR, RHOICE, RHOAIR, SINWIN, COSWIN, SINWAT, COSWAT
44 _RL GRAV, ECCEN, ECM2, GMIN, RADIUS, DELT1, DELT2, PSTAR
45 _RL AAA
46
47 _RL PRESS (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
48 _RL DAIRN (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
49 _RL DWATN (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
50 _RL FORCEX0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
51 _RL FORCEY0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
52 _RL E11 (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
53 _RL E22 (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
54 _RL E12 (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
55 _RL COR_ICE (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
56 _RL ZMAX (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
57 _RL ZMIN (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
58
59 C FIRST SET UP BASIC CONSTANTS
60 DWAT=0.59 _d 0
61 DAIR=0.01462 _d 0
62 RHOICE=0.91 _d +03
63 RHOAIR=1.3 _d 0
64 C 25 DEG GIVES SIN EQUAL TO 0.4226
65 SINWIN=0.4226 _d 0
66 COSWIN=0.9063 _d 0
67 SINWAT=0.4226 _d 0
68 COSWAT=0.9063 _d 0
69 c do not introduce truning angle
70 SINWIN=ZERO
71 COSWIN=ONE
72 SINWAT=ZERO
73 COSWAT=ONE
74
75 GRAV=9.832 _d 0
76
77 ECCEN=TWO
78 ECM2=ONE/(ECCEN**2)
79 GMIN=1.0 _d -20
80 RADIUS=6370. _d 3
81
82 PSTAR=SEAICE_strength
83
84 DO bj=myByLo(myThid),myByHi(myThid)
85 DO bi=myBxLo(myThid),myBxHi(myThid)
86 DO j=1,sNy
87 DO i=1,sNx
88
89 C NOW SET UP MASS PER UNIT AREA AND CORIOLIS TERM
90 AMASS(I,J,bi,bj)=RHOICE*QUART*(HEFF(i,j,1,bi,bj)
91 & +HEFF(i+1,j,1,bi,bj)
92 & +HEFF(i,j+1,1,bi,bj)
93 & +HEFF(i+1,j+1,1,bi,bj))
94 COR_ICE(I,J,bi,bj)=AMASS(I,J,bi,bj)
95 & *TWO*OMEGA*SINEICE(I,J,bi,bj)
96 C NOW SET UP FORCING FIELD
97 C FIRST WIND
98 C IF WIND ON B GRID
99 AAA=SQRT(GAIRX(I,J,bi,bj)**2+GAIRY(I,J,bi,bj)**2)
100 C IF WIND ON C GRID
101 C AAA=SQRT((HALF*(GAIRX(I+1,J,bi,bj)+GAIRX(I+1,J+1,bi,bj)))**2
102 C & +(HALF*(GAIRY(I,J+1,bi,bj)+GAIRY(I+1,J+1,bi,bj)))**2)
103
104 DAIRN(I,J,bi,bj)=RHOAIR*SEAICE_drag*
105 & (2.70 _d 0+0.142 _d 0*AAA+0.0764 _d 0*AAA*AAA)
106
107 C IF WIND ON B GRID
108 FORCEX(I,J,bi,bj)=DAIRN(I,J,bi,bj)*(COSWIN*GAIRX(I,J,bi,bj)
109 & -SINWIN*GAIRY(I,J,bi,bj))
110 FORCEY(I,J,bi,bj)=DAIRN(I,J,bi,bj)*(SINWIN*GAIRX(I,J,bi,bj)
111 & +COSWIN*GAIRY(I,J,bi,bj))
112 C IF WIND ON C GRID
113 C FORCEX(I,J,bi,bj)=DAIRN(I,J,bi,bj)*
114 C & (COSWIN*HALF*(GAIRX(I+1,J,bi,bj)+GAIRX(I+1,J+1,bi,bj))
115 C & -SINWIN*HALF*(GAIRY(I,J+1,bi,bj)+GAIRY(I+1,J+1,bi,bj)))
116 C FORCEY(I,J,bi,bj)=DAIRN(I,J,bi,bj)*
117 C & (SINWIN*HALF*(GAIRX(I+1,J,bi,bj)+GAIRX(I+1,J+1,bi,bj))
118 C & +COSWIN*HALF*(GAIRY(I,J+1,bi,bj)+GAIRY(I+1,J+1,bi,bj)))
119
120 C STORE WIND ONLY STRESS
121 WINDX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
122 WINDY(I,J,bi,bj)=FORCEY(I,J,bi,bj)
123
124 C NOW ADD IN TILT
125 FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
126 & -COR_ICE(I,J,bi,bj)*GWATY(I,J,bi,bj)
127 FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)
128 & +COR_ICE(I,J,bi,bj)*GWATX(I,J,bi,bj)
129 C NOW SET UP ICE PRESSURE AND VISCOSITIES
130 PRESS(I,J,bi,bj)=PSTAR*HEFF(I,J,1,bi,bj)
131 & *EXP(-20.0 _d 0*(ONE-AREA(I,J,1,bi,bj)))
132 ZMAX(I,J,bi,bj)=(5.0 _d +12/(2.0 _d +04))*PRESS(I,J,bi,bj)
133 ZMIN(I,J,bi,bj)=4.0 _d +08
134 PRESS(I,J,bi,bj)=PRESS(I,J,bi,bj)*HEFFM(I,J,bi,bj)
135 ENDDO
136 ENDDO
137 ENDDO
138 ENDDO
139
140 #ifdef SEAICE_ALLOW_DYNAMICS
141 IF ( SEAICEuseDYNAMICS ) THEN
142
143 C-- Update overlap regions for PRESS
144 _EXCH_XY_R8(PRESS, myThid)
145
146 DO bj=myByLo(myThid),myByHi(myThid)
147 DO bi=myBxLo(myThid),myBxHi(myThid)
148 DO j=1,sNy
149 DO i=1,sNx
150 C NOW CALCULATE PRESSURE FORCE AND ADD TO EXTERNAL FORCE
151 FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
152 & -(QUART/(DXUICE(I,J,bi,bj)*CSUICE(I,J,bi,bj)))
153 & *(PRESS(I+1,J,bi,bj)+PRESS(I+1,J+1,bi,bj)
154 & -PRESS(I,J,bi,bj)-PRESS(I,J+1,bi,bj))
155 FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)-QUART/DYUICE(I,J,bi,bj)
156 & *(PRESS(I,J+1,bi,bj)+PRESS(I+1,J+1,bi,bj)
157 & -PRESS(I,J,bi,bj)-PRESS(I+1,J,bi,bj))
158 C NOW KEEP FORCEX0
159 FORCEX0(I,J,bi,bj)=FORCEX(I,J,bi,bj)
160 FORCEY0(I,J,bi,bj)=FORCEY(I,J,bi,bj)
161 ENDDO
162 ENDDO
163 ENDDO
164 ENDDO
165
166 C DO PSEUDO-TIMESTEPS TO OBTAIN AN ACCURATE VISCOUS-PLASTIC SOLUTION
167 C 10 PSEUDO-TIMESTEPS OR MORE ARE SUGGESTED FOR HIGH-RESOLUTION (~10KM)
168 C 1 PSEUDO-TIMESTEP CAN BE USED FOR LOW-RESOLUTION GLOBAL MODELING
169 C NPSEUDO is now set in data.seaice input file
170 C TIMESTEP FOR PSEUDO-TIMESTEPPING
171 SEAICE_DT = DELTAT/NPSEUDO
172
173 crg what about DWAIN,DRAGS,DRAGA,ETA,ZETA
174
175 crg later c$taf loop = iteration uice,vice
176
177 DO 5000 KII=1,NPSEUDO
178
179 c$taf store uice,vice = comlev1_seaice_ds,
180 c$taf& key = kii + (ikey_dynamics-1)*NPSEUDO
181
182 C NOW DO PREDICTOR TIME STEP
183 DO bj=myByLo(myThid),myByHi(myThid)
184 DO bi=myBxLo(myThid),myBxHi(myThid)
185 DO j=1-OLy,sNy+OLy
186 DO i=1-OLx,sNx+OLx
187 UICE(I,J,2,bi,bj)=UICE(I,J,1,bi,bj)
188 VICE(I,J,2,bi,bj)=VICE(I,J,1,bi,bj)
189 UICEC(I,J,bi,bj)=UICE(I,J,1,bi,bj)
190 VICEC(I,J,bi,bj)=VICE(I,J,1,bi,bj)
191 ENDDO
192 ENDDO
193 ENDDO
194 ENDDO
195
196 DO bj=myByLo(myThid),myByHi(myThid)
197 DO bi=myBxLo(myThid),myBxHi(myThid)
198 DO j=1,sNy
199 DO i=1,sNx
200 C NOW SET UP NON-LINEAR WATER DRAG, FORCEX, FORCEY
201 DWATN(I,J,bi,bj)=SEAICE_waterDrag*SQRT((UICE(I,J,1,bi,bj)
202 & -GWATX(I,J,bi,bj))**2
203 & +(VICE(I,J,1,bi,bj)-GWATY(I,J,bi,bj))**2)
204 DWATN(I,J,bi,bj)=MAX(DWATN(I,J,bi,bj),QUART)
205 C NOW SET UP SYMMETTRIC DRAG
206 DRAGS(I,J,bi,bj)=DWATN(I,J,bi,bj)*COSWAT
207 C NOW SET UP ANTI SYMMETTRIC DRAG PLUS CORIOLIS
208 DRAGA(I,J,bi,bj)=DWATN(I,J,bi,bj)*SINWAT+COR_ICE(I,J,bi,bj)
209 C NOW ADD IN CURRENT FORCE
210 FORCEX(I,J,bi,bj)=FORCEX0(I,J,bi,bj)+DWATN(I,J,bi,bj)
211 & *(COSWAT*GWATX(I,J,bi,bj)
212 & -SINWAT*GWATY(I,J,bi,bj))
213 FORCEY(I,J,bi,bj)=FORCEY0(I,J,bi,bj)+DWATN(I,J,bi,bj)
214 & *(SINWAT*GWATX(I,J,bi,bj)
215 & +COSWAT*GWATY(I,J,bi,bj))
216 ENDDO
217 ENDDO
218 ENDDO
219 ENDDO
220
221 DO bj=myByLo(myThid),myByHi(myThid)
222 DO bi=myBxLo(myThid),myBxHi(myThid)
223 DO j=1,sNy
224 DO i=1,sNx
225 C NOW EVALUATE STRAIN RATES
226 E11(I,J,bi,bj)=HALF/(DXTICE(I,J,bi,bj)*CSTICE(I,J,bi,bj))
227 & *(UICE(I,J,1,bi,bj)+UICE(I,J-1,1,bi,bj)
228 & -UICE(I-1,J,1,bi,bj)-UICE(I-1,J-1,1,bi,bj))
229 & -QUART*(VICE(I,J,1,bi,bj)+VICE(I-1,J,1,bi,bj)
230 & +VICE(I-1,J-1,1,bi,bj)+VICE(I,J-1,1,bi,bj))
231 & *TNGTICE(I,J,bi,bj)/RADIUS
232 E22(I,J,bi,bj)=HALF/DYTICE(I,J,bi,bj)
233 & *(VICE(I,J,1,bi,bj)+VICE(I-1,J,1,bi,bj)
234 & -VICE(I,J-1,1,bi,bj)-VICE(I-1,J-1,1,bi,bj))
235 E12(I,J,bi,bj)=HALF*(HALF/DYTICE(I,J,bi,bj)
236 & *(UICE(I,J,1,bi,bj)+UICE(I-1,J,1,bi,bj)
237 & -UICE(I,J-1,1,bi,bj)-UICE(I-1,J-1,1,bi,bj))
238 & +HALF/(DXTICE(I,J,bi,bj)*CSTICE(I,J,bi,bj))
239 & *(VICE(I,J,1,bi,bj)+VICE(I,J-1,1,bi,bj)
240 & -VICE(I-1,J,1,bi,bj)-VICE(I-1,J-1,1,bi,bj))
241 & +QUART*(UICE(I,J,1,bi,bj)+UICE(I-1,J,1,bi,bj)
242 & +UICE(I-1,J-1,1,bi,bj)+UICE(I,J-1,1,bi,bj))
243 & *TNGTICE(I,J,bi,bj)/RADIUS)
244 C NOW EVALUATE VISCOSITIES
245 DELT1=(E11(I,J,bi,bj)**2+E22(I,J,bi,bj)**2)*(ONE+ECM2)
246 & +4.0 _d 0*ECM2*E12(I,J,bi,bj)**2
247 1 +TWO*E11(I,J,bi,bj)*E22(I,J,bi,bj)*(ONE-ECM2)
248 DELT2=SQRT(DELT1)
249 DELT2=MAX(GMIN,DELT2)
250 ZETA(I,J,bi,bj)=HALF*PRESS(I,J,bi,bj)/DELT2
251 C NOW PUT MIN AND MAX VISCOSITIES IN
252 ZETA(I,J,bi,bj)=MIN(ZMAX(I,J,bi,bj),ZETA(I,J,bi,bj))
253 ZETA(I,J,bi,bj)=MAX(ZMIN(I,J,bi,bj),ZETA(I,J,bi,bj))
254 C NOW SET VISCOSITIES TO ZERO AT HEFFMFLOW PTS
255 ZETA(I,J,bi,bj)=ZETA(I,J,bi,bj)*HEFFM(I,J,bi,bj)
256 ETA(I,J,bi,bj)=ECM2*ZETA(I,J,bi,bj)
257 ENDDO
258 ENDDO
259 ENDDO
260 ENDDO
261
262 C-- Update overlap regions
263 _EXCH_XY_R8(ETA, myThid)
264 _EXCH_XY_R8(ZETA, myThid)
265
266 C NOW ADI SCHEME (ZHANG-J/ROTHROCK 1999,bi,bj)
267 IF ( SEAICEuseLSR ) THEN
268 CALL LSR( myThid )
269 ELSE
270 CALL ADI( myThid )
271 ENDIF
272
273 C NOW DO MODIFIED EULER STEP
274 DO bj=myByLo(myThid),myByHi(myThid)
275 DO bi=myBxLo(myThid),myBxHi(myThid)
276 DO j=1-OLy,sNy+OLy
277 DO i=1-OLx,sNx+OLx
278 UICE(I,J,1,bi,bj)=HALF*(UICE(I,J,1,bi,bj)+UICE(I,J,2,bi,bj))
279 VICE(I,J,1,bi,bj)=HALF*(VICE(I,J,1,bi,bj)+VICE(I,J,2,bi,bj))
280 UICEC(I,J,bi,bj)=UICE(I,J,1,bi,bj)
281 VICEC(I,J,bi,bj)=VICE(I,J,1,bi,bj)
282 ENDDO
283 ENDDO
284 ENDDO
285 ENDDO
286
287 DO bj=myByLo(myThid),myByHi(myThid)
288 DO bi=myBxLo(myThid),myBxHi(myThid)
289 DO j=1,sNy
290 DO i=1,sNx
291 C NOW SET UP NON-LINEAR WATER DRAG
292 DWATN(I,J,bi,bj)=SEAICE_waterDrag*SQRT((UICE(I,J,1,bi,bj)
293 & -GWATX(I,J,bi,bj))**2
294 & +(VICE(I,J,1,bi,bj)-GWATY(I,J,bi,bj))**2)
295 DWATN(I,J,bi,bj)=MAX(DWATN(I,J,bi,bj),QUART)
296 C NOW SET UP SYMMETTRIC DRAG
297 DRAGS(I,J,bi,bj)=DWATN(I,J,bi,bj)*COSWAT
298 C NOW SET UP ANTI SYMMETTRIC DRAG PLUS CORIOLIS
299 DRAGA(I,J,bi,bj)=DWATN(I,J,bi,bj)*SINWAT+COR_ICE(I,J,bi,bj)
300 C NOW ADD IN CURRENT FORCE
301 FORCEX(I,J,bi,bj)=FORCEX0(I,J,bi,bj)+DWATN(I,J,bi,bj)
302 & *(COSWAT*GWATX(I,J,bi,bj)
303 & -SINWAT*GWATY(I,J,bi,bj))
304 FORCEY(I,J,bi,bj)=FORCEY0(I,J,bi,bj)+DWATN(I,J,bi,bj)
305 & *(SINWAT*GWATX(I,J,bi,bj)
306 & +COSWAT*GWATY(I,J,bi,bj))
307 ENDDO
308 ENDDO
309 ENDDO
310 ENDDO
311
312 DO bj=myByLo(myThid),myByHi(myThid)
313 DO bi=myBxLo(myThid),myBxHi(myThid)
314 DO j=1,sNy
315 DO i=1,sNx
316 C NOW EVALUATE STRAIN RATES
317 E11(I,J,bi,bj)=HALF/(DXTICE(I,J,bi,bj)*CSTICE(I,J,bi,bj))
318 & *(UICE(I,J,1,bi,bj)+UICE(I,J-1,1,bi,bj)
319 & -UICE(I-1,J,1,bi,bj)-UICE(I-1,J-1,1,bi,bj))
320 & -QUART*(VICE(I,J,1,bi,bj)+VICE(I-1,J,1,bi,bj)
321 & +VICE(I-1,J-1,1,bi,bj)+VICE(I,J-1,1,bi,bj))
322 & *TNGTICE(I,J,bi,bj)/RADIUS
323 E22(I,J,bi,bj)=HALF/DYTICE(I,J,bi,bj)
324 & *(VICE(I,J,1,bi,bj)+VICE(I-1,J,1,bi,bj)
325 & -VICE(I,J-1,1,bi,bj)-VICE(I-1,J-1,1,bi,bj))
326 E12(I,J,bi,bj)=HALF*(HALF/DYTICE(I,J,bi,bj)
327 & *(UICE(I,J,1,bi,bj)+UICE(I-1,J,1,bi,bj)
328 & -UICE(I,J-1,1,bi,bj)-UICE(I-1,J-1,1,bi,bj))
329 & +HALF/(DXTICE(I,J,bi,bj)*CSTICE(I,J,bi,bj))
330 & *(VICE(I,J,1,bi,bj)+VICE(I,J-1,1,bi,bj)
331 & -VICE(I-1,J,1,bi,bj)-VICE(I-1,J-1,1,bi,bj))
332 & +QUART*(UICE(I,J,1,bi,bj)+UICE(I-1,J,1,bi,bj)
333 & +UICE(I-1,J-1,1,bi,bj)+UICE(I,J-1,1,bi,bj))
334 & *TNGTICE(I,J,bi,bj)/RADIUS)
335 C NOW EVALUATE VISCOSITIES
336 DELT1=(E11(I,J,bi,bj)**2+E22(I,J,bi,bj)**2)*(ONE+ECM2)
337 & +4. _d 0*ECM2*E12(I,J,bi,bj)**2
338 1 +TWO*E11(I,J,bi,bj)*E22(I,J,bi,bj)*(ONE-ECM2)
339 DELT2=SQRT(DELT1)
340 DELT2=MAX(GMIN,DELT2)
341 ZETA(I,J,bi,bj)=HALF*PRESS(I,J,bi,bj)/DELT2
342 C NOW PUT MIN AND MAX VISCOSITIES IN
343 ZETA(I,J,bi,bj)=MIN(ZMAX(I,J,bi,bj),ZETA(I,J,bi,bj))
344 ZETA(I,J,bi,bj)=MAX(ZMIN(I,J,bi,bj),ZETA(I,J,bi,bj))
345 C NOW SET VISCOSITIES TO ZERO AT HEFFMFLOW PTS
346 ZETA(I,J,bi,bj)=ZETA(I,J,bi,bj)*HEFFM(I,J,bi,bj)
347 ETA(I,J,bi,bj)=ECM2*ZETA(I,J,bi,bj)
348 ENDDO
349 ENDDO
350 ENDDO
351 ENDDO
352
353 C-- Update overlap regions
354 _EXCH_XY_R8(ETA, myThid)
355 _EXCH_XY_R8(ZETA, myThid)
356
357 C GET READY FOR SECOND CALL OF ADI
358 DO bj=myByLo(myThid),myByHi(myThid)
359 DO bi=myBxLo(myThid),myBxHi(myThid)
360 DO j=1-OLy,sNy+OLy
361 DO i=1-OLx,sNx+OLx
362 UICE(I,J,2,bi,bj)=UICEC(I,J,bi,bj)
363 VICE(I,J,2,bi,bj)=VICEC(I,J,bi,bj)
364 ENDDO
365 ENDDO
366 ENDDO
367 ENDDO
368
369 C NOW ADI SCHEME (ZHANG-J/ROTHROCK 1999)
370 IF ( SEAICEuseLSR ) THEN
371 CALL LSR( myThid )
372 ELSE
373 CALL ADI( myThid )
374 ENDIF
375
376 5000 CONTINUE
377
378 c$taf store uice,vice = comlev1, key=ikey_dynamics
379
380 ENDIF
381 #endif SEAICE_ALLOW_DYNAMICS
382
383 C Calculate ocean surface stress
384 CALL OSTRES ( DWATN, COR_ICE, myThid )
385
386 #ifdef SEAICE_ALLOW_DYNAMICS
387 IF ( SEAICEuseDYNAMICS ) THEN
388
389 c Put a cap on ice velocity
390 c limit velocity to 0.40 m s-1 to avoid potential CFL violations
391 c in open water areas (drift of zero thickness ice)
392 DO bj=myByLo(myThid),myByHi(myThid)
393 DO bi=myBxLo(myThid),myBxHi(myThid)
394 DO j=1-OLy,sNy+OLy
395 DO i=1-OLx,sNx+OLx
396 #ifdef SEAICE_DEBUG
397 c write(*,'(2i4,2i2,f7.1,7f12.3)')
398 c & i,j,bi,bj,UVM(I,J,bi,bj),amass(i,j,bi,bj)
399 c & ,gwatx(I,J,bi,bj),gwaty(i,j,bi,bj)
400 c & ,forcex(I,J,bi,bj),forcey(i,j,bi,bj)
401 c & ,uice(i,j,1,bi,bj)
402 c & ,vice(i,j,1,bi,bj)
403 #endif SEAICE_DEBUG
404 UICE(i,j,1,bi,bj)=min(UICE(i,j,1,bi,bj),0.40 _d +00)
405 VICE(i,j,1,bi,bj)=min(VICE(i,j,1,bi,bj),0.40 _d +00)
406 UICE(i,j,1,bi,bj)=max(UICE(i,j,1,bi,bj),-0.40 _d +00)
407 VICE(i,j,1,bi,bj)=max(VICE(i,j,1,bi,bj),-0.40 _d +00)
408 ENDDO
409 ENDDO
410 ENDDO
411 ENDDO
412
413 ENDIF
414 #endif SEAICE_ALLOW_DYNAMICS
415
416 #endif ALLOW_SEAICE
417
418 RETURN
419 END
  ViewVC Help
Powered by ViewVC 1.1.22