/[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.11 - (show annotations) (download)
Thu Oct 9 04:19:20 2003 UTC (20 years, 8 months ago) by edhill
Branch: MAIN
CVS Tags: checkpoint51k_post, checkpoint52l_pre, hrcube4, hrcube5, checkpoint52d_pre, checkpoint52j_pre, checkpoint51o_pre, checkpoint51l_post, checkpoint52l_post, checkpoint52k_post, checkpoint52, checkpoint52f_post, checkpoint51t_post, checkpoint51n_post, checkpoint52i_pre, hrcube_1, hrcube_2, hrcube_3, checkpoint51s_post, checkpoint52e_pre, checkpoint52e_post, checkpoint51n_pre, checkpoint52b_pre, checkpoint51l_pre, checkpoint52m_post, checkpoint51q_post, checkpoint52b_post, checkpoint52c_post, checkpoint52f_pre, checkpoint51r_post, checkpoint51i_post, checkpoint52d_post, checkpoint52a_pre, checkpoint52i_post, checkpoint52h_pre, checkpoint52j_post, branch-netcdf, checkpoint51o_post, checkpoint52a_post, ecco_c52_e35, checkpoint51m_post, checkpoint51p_post, checkpoint51u_post
Branch point for: branch-nonh, tg2-branch, netcdf-sm0, checkpoint51n_branch
Changes since 1.10: +2 -1 lines 
 o first check-in for the "branch-genmake2" merge
 o verification suite as run on shelley (gcc 3.2.2):

Wed Oct  8 23:42:29 EDT 2003
                T           S           U           V
G D M    c        m  s        m  s        m  s        m  s
E p a R  g  m  m  e  .  m  m  e  .  m  m  e  .  m  m  e  .
N n k u  2  i  a  a  d  i  a  a  d  i  a  a  d  i  a  a  d
2 d e n  d  n  x  n  .  n  x  n  .  n  x  n  .  n  x  n  .

OPTFILE=NONE

Y Y Y Y 13 16 16 16  0 16 16 16 16 16 16 16 16 13 12  0  0 pass  adjustment.128x64x1
Y Y Y Y 16 16 16 16  0 16 16 16 16 16 16  0  0 16 16  0  0 pass  adjustment.cs-32x32x1
Y Y Y Y 16 16 16 16  0 16 16 16 16 16 16 22  0 16 16 22  0 pass  adjust_nlfs.cs-32x32x1
Y Y Y Y -- 13 13 16 16 13 13 13 13 16 16 16 16 16 16 16 16 N/O   advect_cs
Y Y Y Y -- 22 16 16 16 16 16 16 13 16 16 16 16 16 16 16 16 N/O   advect_xy
Y Y Y Y -- 13 16 13 16 16 16 16 16 16 16 22 16 16 16 16 16 N/O   advect_xz
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 pass  aim.5l_cs
Y Y Y Y 14 16 16 16 16 16 16 16 16 13 16 16 16 16 16 13 16 pass  aim.5l_Equatorial_Channel
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 13 16 16 13 13 16 pass  aim.5l_LatLon
Y Y Y Y 13 16 16 16 16 16 16 16 16 16 13 12 13 13 16 13 16 pass  exp0
Y Y Y Y 14 16 16 16 16 16 16 16 22 16 16 16 13 16 16 22 16 pass  exp1
Y Y Y Y 13 13 16 13 16 16 16 16 16 13 13 16 16 13 13 13 13 pass  exp2
Y Y Y Y 16 16 16 16 16 16 16 16 22 16 16 16 16 16 16 16 16 pass  exp4
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 22 16 16 16 22 16 pass  exp5
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 pass  front_relax
Y Y Y Y 14 16 16 13 13 16 16 13 13 16 13 13 16 12 13 13 16 pass  global_ocean.90x40x15
Y Y Y Y 10 16 16 13 13 16 13 16 16 13 13 13 13 16 16 13 16 FAIL  global_ocean.cs32x15
Y Y Y Y  6 11 12 13 13 12 13 16 13  9  9  9  9 10  9  9 11 FAIL  global_ocean_pressure
Y Y Y Y 14 16 16 13 16 16 16 13 13 13 13 13 16 12 16 13 16 pass  global_with_exf
Y Y Y Y 14 16 16 16 16 16 16 16 16 11 13 22 13 16 16  9 16 pass  hs94.128x64x5
Y Y Y Y 13 16 16 16 16 16 16 16 16 11 16 16 16 13 16 22 13 pass  hs94.1x64x5
Y Y Y Y 14 16 16 16 16 16 16 16 16 13 16 13 13 16 16 22 13 pass  hs94.cs-32x32x5
Y Y Y Y 10 10 16 13 13 16 16 16 22 16 13 13 13 13 13 22 13 FAIL  ideal_2D_oce
Y Y Y Y  8 16 16 16 16 16 16 16 16 13 13  8 16 16 16 16 16 FAIL  internal_wave
Y Y Y Y 14 16 16 16 16 16 16 16 16 13 13 22 13 13 13 22 16 pass  inverted_barometer
Y Y Y Y 12 16 16 16 16 16 16 16 16 16 13 12 13 13 13 13 13 FAIL  lab_sea
Y Y Y Y 11 16 16 16 16 16 16 16 13 13 13 12 13 16 13 12 13 FAIL  natl_box
Y Y Y Y 16 16 16 16 16 16 16 16 22 16 16 16 16 16 16 16 16 pass  plume_on_slope
Y Y Y Y 13 16 16 16 16 13 16 16 16 16 16 16 16 13 16 16 16 pass  solid-body.cs-32x32x1
1 C $Header: /u/u3/gcmpack/MITgcm/pkg/seaice/dynsolver.F,v 1.10.2.1 2003/10/02 18:18:34 adcroft Exp $
2 C $Name: $
3
4 #include "SEAICE_OPTIONS.h"
5
6 CStartOfInterface
7 SUBROUTINE dynsolver( myTime, myIter, myThid )
8 C /==========================================================\
9 C | SUBROUTINE dynsolver |
10 C | o Ice dynamics using either LSR or ADI solvers |
11 C | LSR: Zhang and Hibler, JGR, 102, 8691-8702, 1997 |
12 C | ADI: Zhang and Rothrock, JGR, 105, 3325-3338, 2000 |
13 C | For parallel computing, LSR is preferred |
14 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 "FFIELDS.h"
23 #include "SEAICE.h"
24 #include "SEAICE_GRID.h"
25 #include "SEAICE_PARAMS.h"
26 #include "SEAICE_FFIELDS.h"
27
28 #ifdef ALLOW_AUTODIFF_TAMC
29 # include "tamc.h"
30 #endif
31
32 C === Routine arguments ===
33 C myTime - Simulation time
34 C myIter - Simulation timestep number
35 C myThid - Thread no. that called this routine.
36 _RL myTime
37 INTEGER myIter
38 INTEGER myThid
39 CEndOfInterface
40
41 #ifdef ALLOW_SEAICE
42
43 C === Local variables ===
44 C i,j,bi,bj - Loop counters
45
46 INTEGER i, j, bi, bj, kii
47 _RL DWAT, DAIR, RHOICE, RHOAIR, SINWIN, COSWIN, SINWAT, COSWAT
48 _RL GRAV, ECCEN, ECM2, RADIUS, DELT1, DELT2, PSTAR, AAA
49 _RL TEMPVAR, U1, V1
50
51 _RL PRESS (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
52 _RL PRESS0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
53 _RL DAIRN (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
54 _RL DWATN (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
55 _RL FORCEX0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
56 _RL FORCEY0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
57 _RL E11 (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
58 _RL E22 (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
59 _RL E12 (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
60 _RL COR_ICE (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
61 _RL ZMAX (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
62 _RL ZMIN (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
63
64 C-- FIRST SET UP BASIC CONSTANTS
65 DWAT=0.59 _d 0
66 DAIR=0.01462 _d 0
67 RHOICE=0.91 _d +03
68 RHOAIR=1.3 _d 0
69 GRAV=9.832 _d 0
70 ECCEN=TWO
71 ECM2=ONE/(ECCEN**2)
72 RADIUS=6370. _d 3
73 PSTAR=SEAICE_strength
74
75 C-- 25 DEG GIVES SIN EQUAL TO 0.4226
76 SINWIN=0.4226 _d 0
77 COSWIN=0.9063 _d 0
78 SINWAT=0.4226 _d 0
79 COSWAT=0.9063 _d 0
80
81 C-- Do not introduce turning angle
82 SINWIN=ZERO
83 COSWIN=ONE
84 SINWAT=ZERO
85 COSWAT=ONE
86
87 C-- NOW SET UP MASS PER UNIT AREA AND CORIOLIS TERM
88 DO bj=myByLo(myThid),myByHi(myThid)
89 DO bi=myBxLo(myThid),myBxHi(myThid)
90 DO j=1,sNy
91 DO i=1,sNx
92 AMASS(I,J,bi,bj)=RHOICE*QUART*(HEFF(i,j,1,bi,bj)
93 & +HEFF(i-1,j,1,bi,bj)
94 & +HEFF(i,j-1,1,bi,bj)
95 & +HEFF(i-1,j-1,1,bi,bj))
96 COR_ICE(I,J,bi,bj)=AMASS(I,J,bi,bj)
97 & *TWO*OMEGA*SINEICE(I,J,bi,bj)
98 ENDDO
99 ENDDO
100 ENDDO
101 ENDDO
102
103 C-- NOW SET UP FORCING FIELDS
104
105 C-- Wind stress is computed on South-West B-grid U/V
106 C locations from wind on tracer locations
107 DO bj=myByLo(myThid),myByHi(myThid)
108 DO bi=myBxLo(myThid),myBxHi(myThid)
109 DO j=1,sNy
110 DO i=1,sNx
111 U1=QUART*(UWIND(I-1,J-1,bi,bj)+UWIND(I-1,J,bi,bj)
112 & +UWIND(I ,J-1,bi,bj)+UWIND(I ,J,bi,bj))
113 V1=QUART*(VWIND(I-1,J-1,bi,bj)+VWIND(I-1,J,bi,bj)
114 & +VWIND(I ,J-1,bi,bj)+VWIND(I ,J,bi,bj))
115 AAA=U1**2+V1**2
116 IF ( AAA .LE. SEAICE_EPS_SQ ) THEN
117 AAA=SEAICE_EPS
118 ELSE
119 AAA=SQRT(AAA)
120 ENDIF
121 C first ocean surface stress
122 DAIRN(I,J,bi,bj)=RHOAIR*OCEAN_drag
123 & *(2.70 _d 0+0.142 _d 0*AAA+0.0764 _d 0*AAA*AAA)
124 WINDX(I,J,bi,bj)=DAIRN(I,J,bi,bj)*(COSWIN*U1-SINWIN*V1)
125 WINDY(I,J,bi,bj)=DAIRN(I,J,bi,bj)*(SINWIN*U1+COSWIN*V1)
126
127 C now ice surface stress
128 DAIRN(I,J,bi,bj)=RHOAIR*(SEAICE_drag*AAA*AREA(I,J,1,bi,bj)
129 & +OCEAN_drag*(2.70 _d 0+0.142 _d 0*AAA
130 & +0.0764 _d 0*AAA*AAA)*(ONE-AREA(I,J,1,bi,bj)))
131 FORCEX(I,J,bi,bj)=DAIRN(I,J,bi,bj)*(COSWIN*U1-SINWIN*V1)
132 FORCEY(I,J,bi,bj)=DAIRN(I,J,bi,bj)*(SINWIN*U1+COSWIN*V1)
133 ENDDO
134 ENDDO
135 ENDDO
136 ENDDO
137
138 DO bj=myByLo(myThid),myByHi(myThid)
139 DO bi=myBxLo(myThid),myBxHi(myThid)
140 DO j=1,sNy
141 DO i=1,sNx
142 C-- NOW ADD IN TILT
143 FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
144 & -COR_ICE(I,J,bi,bj)*GWATY(I,J,bi,bj)
145 FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)
146 & +COR_ICE(I,J,bi,bj)*GWATX(I,J,bi,bj)
147 C NOW KEEP FORCEX0
148 FORCEX0(I,J,bi,bj)=FORCEX(I,J,bi,bj)
149 FORCEY0(I,J,bi,bj)=FORCEY(I,J,bi,bj)
150 C-- NOW SET UP ICE PRESSURE AND VISCOSITIES
151 PRESS0(I,J,bi,bj)=PSTAR*HEFF(I,J,1,bi,bj)
152 & *EXP(-20.0 _d 0*(ONE-AREA(I,J,1,bi,bj)))
153 ZMAX(I,J,bi,bj)=(5.0 _d +12/(2.0 _d +04))*PRESS0(I,J,bi,bj)
154 ZMIN(I,J,bi,bj)=4.0 _d +08
155 PRESS0(I,J,bi,bj)=PRESS0(I,J,bi,bj)*HEFFM(I,J,bi,bj)
156 ENDDO
157 ENDDO
158 ENDDO
159 ENDDO
160
161 #ifdef SEAICE_ALLOW_DYNAMICS
162
163 IF ( SEAICEuseDYNAMICS ) THEN
164
165 #ifdef ALLOW_AUTODIFF_TAMC
166 CADJ STORE uice = comlev1, key=ikey_dynamics
167 CADJ STORE vice = comlev1, key=ikey_dynamics
168 #endif /* ALLOW_AUTODIFF_TAMC */
169
170 C DO PSEUDO-TIMESTEPS TO OBTAIN AN ACCURATE VISCOUS-PLASTIC SOLUTION
171 C 1 PSEUDO-TIMESTEP IS SUGGESTED FOR LSR SLOVER
172 C A RANGE OF 5-200 PSEUDO-TIMESTEPS IS SUGGESTED FOR ADI SLOVER
173 C DEPENDING ON ACCURACY REQUIREMENTS OF SPECIFIC APPLICATION
174 C NPSEUDO is set in data.seaice input file
175 C TIMESTEP FOR PSEUDO-TIMESTEPPING
176 SEAICE_DT = DELTAT/NPSEUDO
177
178 crg what about DWAIN,DRAGS,DRAGA,ETA,ZETA
179
180 crg later c$taf loop = iteration uice,vice
181
182 #ifndef ALLOW_AUTODIFF_TAMC
183 DO 5000 KII=1,NPSEUDO
184 #else
185 IF ( NPSEUDO .GT. 1 ) THEN
186 STOP 'S/R DYNSOLVER: NPSEUDO NEEDS TO BE = 1 FOR ADJOINT'
187 ENDIF
188 #endif /* ALLOW_AUTODIFF_TAMC */
189
190 cdm c$taf store uice,vice = comlev1_seaice_ds,
191 cdm c$taf& key = kii + (ikey_dynamics-1)*NPSEUDO
192 C NOW DO PREDICTOR TIME STEP
193 DO bj=myByLo(myThid),myByHi(myThid)
194 DO bi=myBxLo(myThid),myBxHi(myThid)
195 DO j=1-OLy,sNy+OLy
196 DO i=1-OLx,sNx+OLx
197 UICE(I,J,2,bi,bj)=UICE(I,J,1,bi,bj)
198 VICE(I,J,2,bi,bj)=VICE(I,J,1,bi,bj)
199 UICEC(I,J,bi,bj)=UICE(I,J,1,bi,bj)
200 VICEC(I,J,bi,bj)=VICE(I,J,1,bi,bj)
201 ENDDO
202 ENDDO
203 ENDDO
204 ENDDO
205
206 DO bj=myByLo(myThid),myByHi(myThid)
207 DO bi=myBxLo(myThid),myBxHi(myThid)
208 DO j=1,sNy
209 DO i=1,sNx
210 C NOW EVALUATE STRAIN RATES
211 E11(I,J,bi,bj)=HALF/(DXTICE(I,J,bi,bj)*CSTICE(I,J,bi,bj))
212 & *(UICE(I+1,J+1,1,bi,bj)+UICE(I+1,J,1,bi,bj)
213 & -UICE(I,J+1,1,bi,bj)-UICE(I,J,1,bi,bj))
214 & -QUART*(VICE(I+1,J+1,1,bi,bj)+VICE(I,J+1,1,bi,bj)
215 & +VICE(I,J,1,bi,bj)+VICE(I+1,J,1,bi,bj))
216 & *TNGTICE(I,J,bi,bj)/RADIUS
217 E22(I,J,bi,bj)=HALF/DYTICE(I,J,bi,bj)
218 & *(VICE(I+1,J+1,1,bi,bj)+VICE(I,J+1,1,bi,bj)
219 & -VICE(I+1,J,1,bi,bj)-VICE(I,J,1,bi,bj))
220 E12(I,J,bi,bj)=HALF*(HALF/DYTICE(I,J,bi,bj)
221 & *(UICE(I+1,J+1,1,bi,bj)+UICE(I,J+1,1,bi,bj)
222 & -UICE(I+1,J,1,bi,bj)-UICE(I,J,1,bi,bj))
223 & +HALF/(DXTICE(I,J,bi,bj)*CSTICE(I,J,bi,bj))
224 & *(VICE(I+1,J+1,1,bi,bj)+VICE(I+1,J,1,bi,bj)
225 & -VICE(I,J+1,1,bi,bj)-VICE(I,J,1,bi,bj))
226 & +QUART*(UICE(I+1,J+1,1,bi,bj)+UICE(I,J+1,1,bi,bj)
227 & +UICE(I,J,1,bi,bj)+UICE(I+1,J,1,bi,bj))
228 & *TNGTICE(I,J,bi,bj)/RADIUS)
229 C NOW EVALUATE VISCOSITIES
230 DELT1=(E11(I,J,bi,bj)**2+E22(I,J,bi,bj)**2)*(ONE+ECM2)
231 & +4.0 _d 0*ECM2*E12(I,J,bi,bj)**2
232 1 +TWO*E11(I,J,bi,bj)*E22(I,J,bi,bj)*(ONE-ECM2)
233 IF ( DELT1 .LE. SEAICE_EPS_SQ ) THEN
234 DELT2=SEAICE_EPS
235 ELSE
236 DELT2=SQRT(DELT1)
237 ENDIF
238 ZETA(I,J,bi,bj)=HALF*PRESS0(I,J,bi,bj)/DELT2
239 C NOW PUT MIN AND MAX VISCOSITIES IN
240 ZETA(I,J,bi,bj)=MIN(ZMAX(I,J,bi,bj),ZETA(I,J,bi,bj))
241 ZETA(I,J,bi,bj)=MAX(ZMIN(I,J,bi,bj),ZETA(I,J,bi,bj))
242 C NOW SET VISCOSITIES TO ZERO AT HEFFMFLOW PTS
243 ZETA(I,J,bi,bj)=ZETA(I,J,bi,bj)*HEFFM(I,J,bi,bj)
244 ETA(I,J,bi,bj)=ECM2*ZETA(I,J,bi,bj)
245 PRESS(I,J,bi,bj)=TWO*ZETA(I,J,bi,bj)*DELT2
246 ENDDO
247 ENDDO
248 ENDDO
249 ENDDO
250
251 C-- Update overlap regions
252 _EXCH_XY_R8(ETA, myThid)
253 _EXCH_XY_R8(ZETA, myThid)
254 _EXCH_XY_R8(PRESS, myThid)
255
256 DO bj=myByLo(myThid),myByHi(myThid)
257 DO bi=myBxLo(myThid),myBxHi(myThid)
258 DO j=1,sNy
259 DO i=1,sNx
260 C NOW SET UP NON-LINEAR WATER DRAG, FORCEX, FORCEY
261 TEMPVAR=(UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj))**2
262 & +(VICE(I,J,1,bi,bj)-GWATY(I,J,bi,bj))**2
263 IF ( TEMPVAR .LE. (QUART/SEAICE_waterDrag)**2 ) THEN
264 DWATN(I,J,bi,bj)=QUART
265 ELSE
266 DWATN(I,J,bi,bj)=SEAICE_waterDrag*SQRT(TEMPVAR)
267 ENDIF
268 C NOW SET UP SYMMETTRIC DRAG
269 DRAGS(I,J,bi,bj)=DWATN(I,J,bi,bj)*COSWAT
270 C NOW SET UP ANTI SYMMETTRIC DRAG PLUS CORIOLIS
271 DRAGA(I,J,bi,bj)=DWATN(I,J,bi,bj)*SINWAT+COR_ICE(I,J,bi,bj)
272 C NOW ADD IN CURRENT FORCE
273 FORCEX(I,J,bi,bj)=FORCEX0(I,J,bi,bj)+DWATN(I,J,bi,bj)
274 & *(COSWAT*GWATX(I,J,bi,bj)
275 & -SINWAT*GWATY(I,J,bi,bj))
276 FORCEY(I,J,bi,bj)=FORCEY0(I,J,bi,bj)+DWATN(I,J,bi,bj)
277 & *(SINWAT*GWATX(I,J,bi,bj)
278 & +COSWAT*GWATY(I,J,bi,bj))
279 C NOW CALCULATE PRESSURE FORCE AND ADD TO EXTERNAL FORCE
280 FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
281 & -(QUART/(DXUICE(I,J,bi,bj)*CSUICE(I,J,bi,bj)))
282 & *(PRESS(I,J,bi,bj)+PRESS(I,J-1,bi,bj)
283 & -PRESS(I-1,J,bi,bj)-PRESS(I-1,J-1,bi,bj))
284 FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)-QUART/DYUICE(I,J,bi,bj)
285 & *(PRESS(I,J,bi,bj)+PRESS(I-1,J,bi,bj)
286 & -PRESS(I,J-1,bi,bj)-PRESS(I-1,J-1,bi,bj))
287 ENDDO
288 ENDDO
289 ENDDO
290 ENDDO
291
292 C NOW LSR SCHEME (ZHANG-J/HIBLER 1997)
293 C OR ADI SCHEME (ZHANG-J/ROTHROCK 1999)
294 #ifdef ALLOW_AUTODIFF_TAMC
295 CADJ STORE uice = comlev1, key=ikey_dynamics
296 CADJ STORE vice = comlev1, key=ikey_dynamics
297 CALL LSR( 1, myThid )
298 CADJ STORE uice = comlev1, key=ikey_dynamics
299 CADJ STORE vice = comlev1, key=ikey_dynamics
300 #else /* ALLOW_AUTODIFF_TAMC */
301 IF ( SEAICEuseADI ) THEN
302 CALL ADI( myThid )
303 ELSE
304 CALL LSR( 1, myThid )
305 ENDIF
306 #endif /* ALLOW_AUTODIFF_TAMC */
307
308 C NOW DO MODIFIED EULER STEP
309 DO bj=myByLo(myThid),myByHi(myThid)
310 DO bi=myBxLo(myThid),myBxHi(myThid)
311 DO j=1-OLy,sNy+OLy
312 DO i=1-OLx,sNx+OLx
313 UICE(I,J,1,bi,bj)=HALF*(UICE(I,J,1,bi,bj)+UICE(I,J,2,bi,bj))
314 VICE(I,J,1,bi,bj)=HALF*(VICE(I,J,1,bi,bj)+VICE(I,J,2,bi,bj))
315 UICEC(I,J,bi,bj)=UICE(I,J,1,bi,bj)
316 VICEC(I,J,bi,bj)=VICE(I,J,1,bi,bj)
317 ENDDO
318 ENDDO
319 ENDDO
320 ENDDO
321
322 DO bj=myByLo(myThid),myByHi(myThid)
323 DO bi=myBxLo(myThid),myBxHi(myThid)
324 DO j=1,sNy
325 DO i=1,sNx
326 C NOW EVALUATE STRAIN RATES
327 E11(I,J,bi,bj)=HALF/(DXTICE(I,J,bi,bj)*CSTICE(I,J,bi,bj))
328 & *(UICE(I+1,J+1,1,bi,bj)+UICE(I+1,J,1,bi,bj)
329 & -UICE(I,J+1,1,bi,bj)-UICE(I,J,1,bi,bj))
330 & -QUART*(VICE(I+1,J+1,1,bi,bj)+VICE(I,J+1,1,bi,bj)
331 & +VICE(I,J,1,bi,bj)+VICE(I+1,J,1,bi,bj))
332 & *TNGTICE(I,J,bi,bj)/RADIUS
333 E22(I,J,bi,bj)=HALF/DYTICE(I,J,bi,bj)
334 & *(VICE(I+1,J+1,1,bi,bj)+VICE(I,J+1,1,bi,bj)
335 & -VICE(I+1,J,1,bi,bj)-VICE(I,J,1,bi,bj))
336 E12(I,J,bi,bj)=HALF*(HALF/DYTICE(I,J,bi,bj)
337 & *(UICE(I+1,J+1,1,bi,bj)+UICE(I,J+1,1,bi,bj)
338 & -UICE(I+1,J,1,bi,bj)-UICE(I,J,1,bi,bj))
339 & +HALF/(DXTICE(I,J,bi,bj)*CSTICE(I,J,bi,bj))
340 & *(VICE(I+1,J+1,1,bi,bj)+VICE(I+1,J,1,bi,bj)
341 & -VICE(I,J+1,1,bi,bj)-VICE(I,J,1,bi,bj))
342 & +QUART*(UICE(I+1,J+1,1,bi,bj)+UICE(I,J+1,1,bi,bj)
343 & +UICE(I,J,1,bi,bj)+UICE(I+1,J,1,bi,bj))
344 & *TNGTICE(I,J,bi,bj)/RADIUS)
345 C NOW EVALUATE VISCOSITIES
346 DELT1=(E11(I,J,bi,bj)**2+E22(I,J,bi,bj)**2)*(ONE+ECM2)
347 & +4. _d 0*ECM2*E12(I,J,bi,bj)**2
348 1 +TWO*E11(I,J,bi,bj)*E22(I,J,bi,bj)*(ONE-ECM2)
349 IF ( DELT1 .LE. SEAICE_EPS_SQ ) THEN
350 DELT2=SEAICE_EPS
351 ELSE
352 DELT2=SQRT(DELT1)
353 ENDIF
354 ZETA(I,J,bi,bj)=HALF*PRESS0(I,J,bi,bj)/DELT2
355 C NOW PUT MIN AND MAX VISCOSITIES IN
356 ZETA(I,J,bi,bj)=MIN(ZMAX(I,J,bi,bj),ZETA(I,J,bi,bj))
357 ZETA(I,J,bi,bj)=MAX(ZMIN(I,J,bi,bj),ZETA(I,J,bi,bj))
358 C NOW SET VISCOSITIES TO ZERO AT HEFFMFLOW PTS
359 ZETA(I,J,bi,bj)=ZETA(I,J,bi,bj)*HEFFM(I,J,bi,bj)
360 ETA(I,J,bi,bj)=ECM2*ZETA(I,J,bi,bj)
361 PRESS(I,J,bi,bj)=TWO*ZETA(I,J,bi,bj)*DELT2
362 ENDDO
363 ENDDO
364 ENDDO
365 ENDDO
366
367 C-- Update overlap regions
368 _EXCH_XY_R8(ETA, myThid)
369 _EXCH_XY_R8(ZETA, myThid)
370 _EXCH_XY_R8(PRESS, myThid)
371
372 DO bj=myByLo(myThid),myByHi(myThid)
373 DO bi=myBxLo(myThid),myBxHi(myThid)
374 DO j=1,sNy
375 DO i=1,sNx
376 C NOW SET UP NON-LINEAR WATER DRAG, FORCEX, FORCEY
377 TEMPVAR=(UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj))**2
378 & +(VICE(I,J,1,bi,bj)-GWATY(I,J,bi,bj))**2
379 IF ( TEMPVAR .LE. (QUART/SEAICE_waterDrag)**2 ) THEN
380 DWATN(I,J,bi,bj)=QUART
381 ELSE
382 DWATN(I,J,bi,bj)=SEAICE_waterDrag*SQRT(TEMPVAR)
383 ENDIF
384 C NOW SET UP SYMMETTRIC DRAG
385 DRAGS(I,J,bi,bj)=DWATN(I,J,bi,bj)*COSWAT
386 C NOW SET UP ANTI SYMMETTRIC DRAG PLUS CORIOLIS
387 DRAGA(I,J,bi,bj)=DWATN(I,J,bi,bj)*SINWAT+COR_ICE(I,J,bi,bj)
388 C NOW ADD IN CURRENT FORCE
389 FORCEX(I,J,bi,bj)=FORCEX0(I,J,bi,bj)+DWATN(I,J,bi,bj)
390 & *(COSWAT*GWATX(I,J,bi,bj)
391 & -SINWAT*GWATY(I,J,bi,bj))
392 FORCEY(I,J,bi,bj)=FORCEY0(I,J,bi,bj)+DWATN(I,J,bi,bj)
393 & *(SINWAT*GWATX(I,J,bi,bj)
394 & +COSWAT*GWATY(I,J,bi,bj))
395 C NOW CALCULATE PRESSURE FORCE AND ADD TO EXTERNAL FORCE
396 FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
397 & -(QUART/(DXUICE(I,J,bi,bj)*CSUICE(I,J,bi,bj)))
398 & *(PRESS(I,J,bi,bj)+PRESS(I,J-1,bi,bj)
399 & -PRESS(I-1,J,bi,bj)-PRESS(I-1,J-1,bi,bj))
400 FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)-QUART/DYUICE(I,J,bi,bj)
401 & *(PRESS(I,J,bi,bj)+PRESS(I-1,J,bi,bj)
402 & -PRESS(I,J-1,bi,bj)-PRESS(I-1,J-1,bi,bj))
403 ENDDO
404 ENDDO
405 ENDDO
406 ENDDO
407
408 C GET READY FOR SECOND CALL OF ADI
409 IF(SEAICEuseADI ) THEN
410 DO bj=myByLo(myThid),myByHi(myThid)
411 DO bi=myBxLo(myThid),myBxHi(myThid)
412 DO j=1-OLy,sNy+OLy
413 DO i=1-OLx,sNx+OLx
414 UICE(I,J,2,bi,bj)=UICEC(I,J,bi,bj)
415 VICE(I,J,2,bi,bj)=VICEC(I,J,bi,bj)
416 ENDDO
417 ENDDO
418 ENDDO
419 ENDDO
420 ENDIF
421
422 C NOW LSR SCHEME (ZHANG-J/HIBLER 1997)
423 C OR ADI SCHEME (ZHANG-J/ROTHROCK 1999)
424 #ifdef ALLOW_AUTODIFF_TAMC
425 CALL LSR( 2, myThid )
426 #else /* ALLOW_AUTODIFF_TAMC */
427 IF ( SEAICEuseADI ) THEN
428 CALL ADI( myThid )
429 ELSE
430 CALL LSR( 2, myThid )
431 ENDIF
432 #endif /* ALLOW_AUTODIFF_TAMC */
433
434 #ifndef ALLOW_AUTODIFF_TAMC
435 5000 CONTINUE
436 #endif /* ALLOW_AUTODIFF_TAMC */
437
438 cdm c$taf store uice,vice = comlev1, key=ikey_dynamics
439
440 ENDIF
441 #endif /* SEAICE_ALLOW_DYNAMICS */
442
443 C Calculate ocean surface stress
444 CALL OSTRES ( DWATN, COR_ICE, myThid )
445
446 #ifdef SEAICE_ALLOW_DYNAMICS
447
448 IF ( SEAICEuseDYNAMICS ) THEN
449
450 #ifdef ALLOW_AUTODIFF_TAMC
451 CADJ STORE uice = comlev1, key=ikey_dynamics
452 CADJ STORE vice = comlev1, key=ikey_dynamics
453 #endif /* ALLOW_AUTODIFF_TAMC */
454 c Put a cap on ice velocity
455 c limit velocity to 0.40 m s-1 to avoid potential CFL violations
456 c in open water areas (drift of zero thickness ice)
457 DO bj=myByLo(myThid),myByHi(myThid)
458 DO bi=myBxLo(myThid),myBxHi(myThid)
459 DO j=1-OLy,sNy+OLy
460 DO i=1-OLx,sNx+OLx
461 #ifdef SEAICE_DEBUG
462 c write(*,'(2i4,2i2,f7.1,7f12.3)')
463 c & i,j,bi,bj,UVM(I,J,bi,bj),amass(i,j,bi,bj)
464 c & ,gwatx(I,J,bi,bj),gwaty(i,j,bi,bj)
465 c & ,forcex(I,J,bi,bj),forcey(i,j,bi,bj)
466 c & ,uice(i,j,1,bi,bj)
467 c & ,vice(i,j,1,bi,bj)
468 #endif /* SEAICE_DEBUG */
469 UICE(i,j,1,bi,bj)=min(UICE(i,j,1,bi,bj),0.40 _d +00)
470 VICE(i,j,1,bi,bj)=min(VICE(i,j,1,bi,bj),0.40 _d +00)
471 ENDDO
472 ENDDO
473 ENDDO
474 ENDDO
475 #ifdef ALLOW_AUTODIFF_TAMC
476 CADJ STORE uice = comlev1, key=ikey_dynamics
477 CADJ STORE vice = comlev1, key=ikey_dynamics
478 #endif /* ALLOW_AUTODIFF_TAMC */
479 DO bj=myByLo(myThid),myByHi(myThid)
480 DO bi=myBxLo(myThid),myBxHi(myThid)
481 DO j=1-OLy,sNy+OLy
482 DO i=1-OLx,sNx+OLx
483 UICE(i,j,1,bi,bj)=max(UICE(i,j,1,bi,bj),-0.40 _d +00)
484 VICE(i,j,1,bi,bj)=max(VICE(i,j,1,bi,bj),-0.40 _d +00)
485 ENDDO
486 ENDDO
487 ENDDO
488 ENDDO
489
490 ENDIF
491 #endif /* SEAICE_ALLOW_DYNAMICS */
492
493 #endif /* ALLOW_SEAICE */
494
495 RETURN
496 END
  ViewVC Help
Powered by ViewVC 1.1.22