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

Diff of /MITgcm/pkg/seaice/dynsolver.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

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