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

Annotation of /MITgcm/pkg/seaice/seaice_init_varia.F

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


Revision 1.36 - (hide annotations) (download)
Wed Mar 18 14:35:38 2009 UTC (15 years, 2 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint61l
Changes since 1.35: +2 -2 lines
forgot a few instances of 330., when replace this number by a runtime
parameter SEAICE_rhoSnow. There is still one instance in
cost_ice_test.F, which I am not touching

1 mlosch 1.36 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_init_varia.F,v 1.35 2009/03/18 13:48:53 mlosch Exp $
2 heimbach 1.1 C $Name: $
3    
4     #include "SEAICE_OPTIONS.h"
5    
6     CStartOfInterface
7     SUBROUTINE SEAICE_INIT_VARIA( myThid )
8     C /==========================================================\
9     C | SUBROUTINE SEAICE_INIT_VARIA |
10     C | o Initialization of sea ice model. |
11     C |==========================================================|
12     C \==========================================================/
13     IMPLICIT NONE
14 jmc 1.23
15 heimbach 1.1 C === Global variables ===
16     #include "SIZE.h"
17     #include "EEPARAMS.h"
18     #include "PARAMS.h"
19     #include "GRID.h"
20 dimitri 1.21 #include "DYNVARS.h"
21 heimbach 1.1 #include "SEAICE.h"
22     #include "SEAICE_DIAGS.h"
23     #include "SEAICE_PARAMS.h"
24     #include "FFIELDS.h"
25     #ifdef ALLOW_EXCH2
26 dimitri 1.24 # include "W2_EXCH2_TOPOLOGY.h"
27     # include "W2_EXCH2_PARAMS.h"
28     #endif
29 dimitri 1.11 #ifdef ALLOW_OBCS
30 dimitri 1.24 # include "OBCS_OPTIONS.h"
31     # include "OBCS.h"
32     #endif
33 heimbach 1.1
34     C === Routine arguments ===
35     C myThid - Thread no. that called this routine.
36     INTEGER myThid
37     CEndOfInterface
38 jmc 1.23
39 heimbach 1.1 C === Local variables ===
40     C i,j,k,bi,bj - Loop counters
41    
42     INTEGER i, j, k, bi, bj
43 mlosch 1.8 _RL PSTAR
44 heimbach 1.1 _RS mask_uice
45     INTEGER myIter, myTile
46    
47 dimitri 1.11 #ifdef ALLOW_OBCS
48     INTEGER I_obc, J_obc, kSurface
49 jmc 1.23 IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
50     kSurface = Nr
51     ELSE
52 dimitri 1.11 kSurface = 1
53 jmc 1.23 ENDIF
54 dimitri 1.11 #endif /* ALLOW_OBCS */
55    
56 jmc 1.23 C-- Initialise all variables in common blocks:
57 heimbach 1.1 DO bj=myByLo(myThid),myByHi(myThid)
58     DO bi=myBxLo(myThid),myBxHi(myThid)
59 jmc 1.23 DO k=1,3
60     DO j=1-OLy,sNy+OLy
61     DO i=1-OLx,sNx+OLx
62     HEFF(i,j,k,bi,bj)=0. _d 0
63     AREA(i,j,k,bi,bj)=0. _d 0
64     UICE(i,j,k,bi,bj)=0. _d 0
65     VICE(i,j,k,bi,bj)=0. _d 0
66     ENDDO
67     ENDDO
68     ENDDO
69     #ifdef SEAICE_MULTICATEGORY
70     DO k=1,MULTDIM
71     DO j=1-OLy,sNy+OLy
72     DO i=1-OLx,sNx+OLx
73     TICES(i,j,k,bi,bj)=0. _d 0
74 heimbach 1.1 ENDDO
75     ENDDO
76     ENDDO
77 jmc 1.23 #endif
78     DO j=1-OLy,sNy+OLy
79     DO i=1-OLx,sNx+OLx
80     ETA (i,j,bi,bj) = 0. _d 0
81     ZETA(i,j,bi,bj) = 0. _d 0
82     DRAGS(i,j,bi,bj) = 0. _d 0
83     DRAGA(i,j,bi,bj) = 0. _d 0
84     FORCEX(i,j,bi,bj) = 0. _d 0
85     FORCEY(i,j,bi,bj) = 0. _d 0
86     UICEC(i,j,bi,bj) = 0. _d 0
87     VICEC(i,j,bi,bj) = 0. _d 0
88     #ifdef SEAICE_CGRID
89     seaiceMassC(i,j,bi,bj)=0. _d 0
90     seaiceMassU(i,j,bi,bj)=0. _d 0
91     seaiceMassV(i,j,bi,bj)=0. _d 0
92     seaiceMaskU(i,j,bi,bj)=0. _d 0
93     seaiceMaskV(i,j,bi,bj)=0. _d 0
94     # ifdef SEAICE_ALLOW_EVP
95     stressDivergenceX(i,j,bi,bj) = 0. _d 0
96     stressDivergenceY(i,j,bi,bj) = 0. _d 0
97     seaice_sigma1 (i,j,bi,bj) = 0. _d 0
98     seaice_sigma2 (i,j,bi,bj) = 0. _d 0
99     seaice_sigma12(i,j,bi,bj) = 0. _d 0
100     # endif /* SEAICE_ALLOW_EVP */
101     #else /* SEAICE_CGRID */
102     AMASS(i,j,bi,bj) = 0. _d 0
103     DAIRN(i,j,bi,bj) = 0. _d 0
104     UVM(i,j,bi,bj) = 0. _d 0
105     WINDX(i,j,bi,bj) = 0. _d 0
106     WINDY(i,j,bi,bj) = 0. _d 0
107 dimitri 1.25 GWATX(i,j,bi,bj) = 0. _d 0
108     GWATY(i,j,bi,bj) = 0. _d 0
109     KGEO(i,j,bi,bj) = 0
110 jmc 1.23 #endif /* SEAICE_CGRID */
111     DWATN(i,j,bi,bj) = 0. _d 0
112     PRESS0(i,j,bi,bj) = 0. _d 0
113     FORCEX0(i,j,bi,bj)= 0. _d 0
114     FORCEY0(i,j,bi,bj)= 0. _d 0
115     ZMAX(i,j,bi,bj) = 0. _d 0
116     ZMIN(i,j,bi,bj) = 0. _d 0
117     HSNOW(i,j,bi,bj) = 0. _d 0
118 dimitri 1.18 #ifdef SEAICE_SALINITY
119 jmc 1.23 HSALT(i,j,bi,bj) = 0. _d 0
120     #endif
121 dimitri 1.33 #ifdef SEAICE_AGE
122     ICEAGE(i,j,bi,bj) = 0. _d 0
123     #endif
124 jmc 1.23 HEFFM(i,j,bi,bj) = 0. _d 0
125     YNEG (i,j,bi,bj) = 0. _d 0
126     RIVER(i,j,bi,bj) = 0. _d 0
127     TMIX(i,j,bi,bj) = 0. _d 0
128     TICE(i,j,bi,bj) = 0. _d 0
129     TAUX(i,j,bi,bj) = 0. _d 0
130     TAUY(i,j,bi,bj) = 0. _d 0
131     #ifdef ALLOW_SEAICE_COST_EXPORT
132     uHeffExportCell(i,j,bi,bj) = 0. _d 0
133     vHeffExportCell(i,j,bi,bj) = 0. _d 0
134 dimitri 1.18 #endif
135 heimbach 1.32 saltWtrIce(i,j,bi,bj) = 0. _d 0
136     frWtrIce(i,j,bi,bj) = 0. _d 0
137     frWtrAtm(i,j,bi,bj) = 0. _d 0
138 heimbach 1.1 ENDDO
139     ENDDO
140     ENDDO
141     ENDDO
142    
143     C-- Initialize grid info
144     DO bj=myByLo(myThid),myByHi(myThid)
145     DO bi=myBxLo(myThid),myBxHi(myThid)
146     DO j=1-OLy,sNy+OLy
147     DO i=1-OLx,sNx+OLx
148     HEFFM(i,j,bi,bj)=ONE
149 heimbach 1.10 IF (_hFacC(i,j,1,bi,bj).eq.0.) HEFFM(i,j,bi,bj)=0. _d 0
150 heimbach 1.1 ENDDO
151     ENDDO
152 jmc 1.23 DO j=1-OLy+1,sNy+OLy
153     DO i=1-OLx+1,sNx+OLx
154 dimitri 1.22 #ifdef SEAICE_CGRID
155 jmc 1.23 seaiceMaskU(i,j,bi,bj)= 0.0 _d 0
156     seaiceMaskV(i,j,bi,bj)= 0.0 _d 0
157     mask_uice=HEFFM(i,j,bi,bj)+HEFFM(i-1,j ,bi,bj)
158     IF(mask_uice.GT.1.5) seaiceMaskU(i,j,bi,bj)=ONE
159     mask_uice=HEFFM(i,j,bi,bj)+HEFFM(i ,j-1,bi,bj)
160     IF(mask_uice.GT.1.5) seaiceMaskV(i,j,bi,bj)=ONE
161 dimitri 1.22 #else
162     UVM(i,j,bi,bj)=0. _d 0
163 jmc 1.23 mask_uice=HEFFM(i,j, bi,bj)+HEFFM(i-1,j-1,bi,bj)
164     & +HEFFM(i,j-1,bi,bj)+HEFFM(i-1,j, bi,bj)
165     IF(mask_uice.GT.3.5) UVM(i,j,bi,bj)=ONE
166 dimitri 1.22 #endif /* SEAICE_CGRID */
167     ENDDO
168     ENDDO
169 mlosch 1.34 #ifdef SEAICE_CGRID
170     C coefficients for metric terms
171     DO j=1-OLy,sNy+OLy
172     DO i=1-OLx,sNx+OLx
173     k1AtC(I,J,bi,bj) = 0.0 _d 0
174     k1AtZ(I,J,bi,bj) = 0.0 _d 0
175     k2AtC(I,J,bi,bj) = 0.0 _d 0
176     k2AtZ(I,J,bi,bj) = 0.0 _d 0
177     ENDDO
178     ENDDO
179     IF ( usingSphericalPolarGrid .AND. SEAICEuseMetricTerms ) THEN
180     C This is the only case where tan(phi) is not zero. In this case
181     C C and U points, and Z and V points have the same phi, so that we
182     C only need a copy here. Do not use tan(YC) and tan(YG), because these
183     C can be the geographical coordinates and not the correct grid
184     C coordinates when the grid is rotated (phi/theta/psiEuler .NE. 0)
185     DO j=1-OLy,sNy+OLy
186     DO i=1-OLx,sNx+OLx
187     k2AtC(I,J,bi,bj) = - _tanPhiAtU(I,J,bi,bj)*recip_rSphere
188     k2AtZ(I,J,bi,bj) = - _tanPhiAtV(I,J,bi,bj)*recip_rSphere
189     ENDDO
190     ENDDO
191     ELSEIF ( usingCurvilinearGrid .AND. SEAICEuseMetricTerms ) THEN
192     C compute metric term coefficients from finite difference approximation
193     DO j=1-OLy,sNy+OLy
194     DO i=1-OLx,sNx+OLx-1
195     k1AtC(I,J,bi,bj) = _recip_dyF(I,J,bi,bj)
196     & * ( _dyG(I+1,J,bi,bj) - _dyG(I,J,bi,bj) )
197     & * _recip_dxF(I,J,bi,bj)
198     ENDDO
199     ENDDO
200     DO j=1-OLy,sNy+OLy
201     DO i=1-OLx+1,sNx+OLx
202     k1AtZ(I,J,bi,bj) = _recip_dyU(I,J,bi,bj)
203     & * ( _dyC(I,J,bi,bj) - _dyC(I-1,J,bi,bj) )
204     & * _recip_dxV(I,J,bi,bj)
205     ENDDO
206     ENDDO
207     DO j=1-OLy,sNy+OLy-1
208     DO i=1-OLx,sNx+OLx
209     k2AtC(I,J,bi,bj) = _recip_dxF(I,J,bi,bj)
210     & * ( _dxG(I,J+1,bi,bj) - _dxG(I,J,bi,bj) )
211     & * _recip_dyF(I,J,bi,bj)
212     ENDDO
213     ENDDO
214     DO j=1-OLy+1,sNy+OLy
215     DO i=1-OLx,sNx+OLx
216     k2AtC(I,J,bi,bj) = _recip_dxV(I,J,bi,bj)
217     & * ( _dxC(I,J,bi,bj) - _dxC(I,J-1,bi,bj) )
218     & * _recip_dyU(I,J,bi,bj)
219     ENDDO
220     ENDDO
221     ENDIF
222     #endif /* SEAICE_CGRID */
223 dimitri 1.11
224     #ifdef ALLOW_OBCS
225 dimitri 1.16 IF (useOBCS) THEN
226 dimitri 1.22 C-- If OBCS is turned on, close southern and western boundaries
227 dimitri 1.11 #ifdef ALLOW_OBCS_SOUTH
228 jmc 1.23 DO i=1-Olx,sNx+Olx
229 dimitri 1.11 C Southern boundary
230 jmc 1.23 J_obc = OB_Js(i,bi,bj)
231 dimitri 1.16 IF (J_obc.NE.0) THEN
232 dimitri 1.22 #ifdef SEAICE_CGRID
233 jmc 1.23 seaiceMaskU(i,J_obc,bi,bj)= 0.0 _d 0
234     seaiceMaskV(i,J_obc,bi,bj)= 0.0 _d 0
235 dimitri 1.22 #else
236 jmc 1.23 UVM(i,J_obc,bi,bj)=0. _d 0
237 dimitri 1.22 #endif /* SEAICE_CGRID */
238 dimitri 1.16 ENDIF
239     ENDDO
240 dimitri 1.22 #endif /* ALLOW_OBCS_SOUTH */
241 dimitri 1.11 #ifdef ALLOW_OBCS_WEST
242 jmc 1.23 DO j=1-Oly,sNy+Oly
243 dimitri 1.11 C Western boundary
244 jmc 1.23 I_obc=OB_Iw(j,bi,bj)
245 dimitri 1.16 IF (I_obc.NE.0) THEN
246 dimitri 1.22 #ifdef SEAICE_CGRID
247 jmc 1.23 seaiceMaskU(I_obc,j,bi,bj)= 0.0 _d 0
248     seaiceMaskV(I_obc,j,bi,bj)= 0.0 _d 0
249 dimitri 1.22 #else
250 jmc 1.23 UVM(I_obc,j,bi,bj)=0. _d 0
251 dimitri 1.22 #endif /* SEAICE_CGRID */
252 dimitri 1.16 ENDIF
253     ENDDO
254 dimitri 1.22 #endif /* ALLOW_OBCS_WEST */
255 dimitri 1.16 ENDIF
256 dimitri 1.11 #endif /* ALLOW_OBCS */
257    
258 heimbach 1.1 #ifdef ALLOW_EXCH2
259 dimitri 1.11 #ifndef SEAICE_CGRID
260 heimbach 1.1 C-- Special stuff for cubed sphere: assume grid is rectangular and
261     C set UV mask to zero except for Arctic and Antarctic cube faces.
262     IF (useCubedSphereExchange) THEN
263     myTile = W2_myTileList(bi)
264     IF ( exch2_myFace(myTile) .EQ. 1 .OR.
265     & exch2_myFace(myTile) .EQ. 2 .OR.
266     & exch2_myFace(myTile) .EQ. 4 .OR.
267     & exch2_myFace(myTile) .EQ. 5 ) THEN
268 jmc 1.23 DO j=1-OLy,sNy+OLy
269     DO i=1-OLx,sNx+OLx
270 heimbach 1.10 UVM(i,j,bi,bj)=0. _d 0
271 heimbach 1.1 ENDDO
272     ENDDO
273     ELSEIF ( exch2_isWedge(myTile) .EQ. 1 ) THEN
274 jmc 1.23 i=1
275     DO j=1-OLy,sNy+OLy
276 heimbach 1.10 UVM(i,j,bi,bj)=0. _d 0
277 heimbach 1.1 ENDDO
278     ELSEIF ( exch2_isSedge(myTile) .EQ. 1 ) THEN
279 jmc 1.23 j=1
280     DO i=1-OLx,sNx+OLx
281 heimbach 1.10 UVM(i,j,bi,bj)=0. _d 0
282 heimbach 1.1 ENDDO
283     ENDIF
284     ENDIF
285 dimitri 1.11 #endif /* SEAICE_CGRID */
286 heimbach 1.1 #endif /* ALLOW_EXCH2 */
287    
288     DO j=1-OLy,sNy+OLy
289     DO i=1-OLx,sNx+OLx
290 jmc 1.23 TICE(i,j,bi,bj)=273.0 _d 0
291 mlosch 1.2 #ifdef SEAICE_MULTICATEGORY
292 heimbach 1.1 DO k=1,MULTDIM
293 jmc 1.23 TICES(i,j,k,bi,bj)=273.0 _d 0
294 heimbach 1.1 ENDDO
295 mlosch 1.2 #endif /* SEAICE_MULTICATEGORY */
296 heimbach 1.1 #ifndef SEAICE_CGRID
297 jmc 1.23 AMASS (i,j,bi,bj)=1000.0 _d 0
298 heimbach 1.1 #else
299 jmc 1.23 seaiceMassC(i,j,bi,bj)=1000.0 _d 0
300     seaiceMassU(i,j,bi,bj)=1000.0 _d 0
301     seaiceMassV(i,j,bi,bj)=1000.0 _d 0
302 heimbach 1.1 #endif
303     ENDDO
304     ENDDO
305    
306 dimitri 1.25 #ifndef SEAICE_CGRID
307 heimbach 1.1 C-- Choose a proxy level for geostrophic velocity,
308     DO j=1-OLy,sNy+OLy
309     DO i=1-OLx,sNx+OLx
310 mlosch 1.35 #ifdef SEAICE_BICE_STRESS
311 jmc 1.23 KGEO(i,j,bi,bj) = 1
312 mlosch 1.35 #else /* SEAICE_BICE_STRESS */
313 heimbach 1.1 IF (klowc(i,j,bi,bj) .LT. 2) THEN
314 jmc 1.23 KGEO(i,j,bi,bj) = 1
315 heimbach 1.1 ELSE
316 jmc 1.23 KGEO(i,j,bi,bj) = 2
317     DO WHILE ( abs(rC(KGEO(i,j,bi,bj))) .LT. 50.0 .AND.
318     & KGEO(i,j,bi,bj) .LT. (klowc(i,j,bi,bj)-1) )
319     KGEO(i,j,bi,bj) = KGEO(i,j,bi,bj) + 1
320 heimbach 1.1 ENDDO
321     ENDIF
322 mlosch 1.35 #endif /* SEAICE_BICE_STRESS */
323 heimbach 1.1 ENDDO
324     ENDDO
325 dimitri 1.25 #endif /* SEAICE_CGRID */
326 heimbach 1.1
327     ENDDO
328     ENDDO
329    
330     C-- Update overlap regions
331     #ifdef SEAICE_CGRID
332     CALL EXCH_UV_XY_RL(seaiceMaskU,seaiceMaskV,.FALSE.,myThid)
333     #else
334     _EXCH_XY_R8(UVM, myThid)
335     #endif
336    
337     C-- Now lets look at all these beasts
338     IF ( debugLevel .GE. debLevB ) THEN
339     myIter=0
340     CALL PLOT_FIELD_XYRL( HEFFM , 'Current HEFFM ' ,
341     & myIter, myThid )
342     #ifdef SEAICE_CGRID
343     CALL PLOT_FIELD_XYRL( seaiceMaskU, 'Current seaiceMaskU',
344     & myIter, myThid )
345     CALL PLOT_FIELD_XYRL( seaiceMaskV, 'Current seaiceMaskV',
346     & myIter, myThid )
347     #else
348     CALL PLOT_FIELD_XYRL( UVM , 'Current UVM ' ,
349     & myIter, myThid )
350     #endif
351     ENDIF
352    
353     C-- Set model variables to initial/restart conditions
354 mlosch 1.26 IF ( .NOT. ( startTime .EQ. baseTime .AND. nIter0 .EQ. 0
355     & .AND. pickupSuff .EQ. ' ') ) THEN
356 heimbach 1.1
357     CALL SEAICE_READ_PICKUP ( myThid )
358    
359     ELSE
360    
361     DO bj=myByLo(myThid),myByHi(myThid)
362     DO bi=myBxLo(myThid),myBxHi(myThid)
363     DO j=1-OLy,sNy+OLy
364     DO i=1-OLx,sNx+OLx
365 jmc 1.23 TMIX(i,j,bi,bj)=TICE(i,j,bi,bj)
366 heimbach 1.1 DO k=1,3
367 jmc 1.23 HEFF(i,j,k,bi,bj)=SEAICE_initialHEFF*HEFFM(i,j,bi,bj)
368     UICE(i,j,k,bi,bj)=ZERO
369     VICE(i,j,k,bi,bj)=ZERO
370 heimbach 1.1 ENDDO
371     ENDDO
372     ENDDO
373     ENDDO
374     ENDDO
375    
376     C-- Read initial sea-ice thickness from file if available.
377     IF ( HeffFile .NE. ' ' ) THEN
378     CALL READ_FLD_XY_RL( HeffFile, ' ', ZETA, 0, myThid )
379     _EXCH_XY_R8(ZETA,myThid)
380     DO bj=myByLo(myThid),myByHi(myThid)
381     DO bi=myBxLo(myThid),myBxHi(myThid)
382     DO j=1-OLy,sNy+OLy
383     DO i=1-OLx,sNx+OLx
384     DO k=1,3
385 jmc 1.23 HEFF(i,j,k,bi,bj) = MAX(ZETA(i,j,bi,bj),ZERO)
386 heimbach 1.1 ENDDO
387     ENDDO
388     ENDDO
389     ENDDO
390     ENDDO
391     ENDIF
392    
393     DO bj=myByLo(myThid),myByHi(myThid)
394     DO bi=myBxLo(myThid),myBxHi(myThid)
395     DO j=1-OLy,sNy+OLy
396     DO i=1-OLx,sNx+OLx
397     DO k=1,3
398 jmc 1.23 IF(HEFF(i,j,k,bi,bj).GT.ZERO)
399     & AREA(i,j,k,bi,bj)=ONE
400 heimbach 1.1 ENDDO
401     ENDDO
402     ENDDO
403     ENDDO
404     ENDDO
405 jmc 1.23
406 dimitri 1.27 C-- Read initial sea-ice area from file if available.
407 mlosch 1.7 IF ( AreaFile .NE. ' ' ) THEN
408     CALL READ_FLD_XY_RL( AreaFile, ' ', ZETA, 0, myThid )
409     _EXCH_XY_R8(ZETA,myThid)
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     DO k=1,3
415 jmc 1.23 AREA(i,j,k,bi,bj) = MAX(ZETA(i,j,bi,bj),ZERO)
416     AREA(i,j,k,bi,bj) = MIN(AREA(i,j,k,bi,bj),ONE)
417     IF ( AREA(i,j,k,bi,bj) .LE. ZERO )
418     & HEFF(i,j,k,bi,bj) = ZERO
419     IF ( HEFF(i,j,k,bi,bj) .LE. ZERO )
420     & AREA(i,j,k,bi,bj) = ZERO
421 mlosch 1.7 ENDDO
422     ENDDO
423     ENDDO
424     ENDDO
425     ENDDO
426     ENDIF
427    
428     DO bj=myByLo(myThid),myByHi(myThid)
429     DO bi=myBxLo(myThid),myBxHi(myThid)
430     DO j=1-OLy,sNy+OLy
431     DO i=1-OLx,sNx+OLx
432 jmc 1.23 HSNOW(i,j,bi,bj)=0.2*AREA(i,j,1,bi,bj)
433 mlosch 1.7 ENDDO
434     ENDDO
435     ENDDO
436     ENDDO
437 dimitri 1.9
438     C-- Read initial snow thickness from file if available.
439     IF ( HsnowFile .NE. ' ' ) THEN
440     CALL READ_FLD_XY_RL( HsnowFile, ' ', ZETA, 0, myThid )
441     _EXCH_XY_R8(ZETA,myThid)
442     DO bj=myByLo(myThid),myByHi(myThid)
443     DO bi=myBxLo(myThid),myBxHi(myThid)
444     DO j=1-OLy,sNy+OLy
445     DO i=1-OLx,sNx+OLx
446 jmc 1.23 HSNOW(i,j,bi,bj) = MAX(ZETA(i,j,bi,bj),ZERO)
447 dimitri 1.9 ENDDO
448     ENDDO
449     ENDDO
450     ENDDO
451     ENDIF
452 dimitri 1.18
453     #ifdef SEAICE_SALINITY
454 dimitri 1.20 DO bj=myByLo(myThid),myByHi(myThid)
455     DO bi=myBxLo(myThid),myBxHi(myThid)
456 jmc 1.23 DO j=1-OLy,sNy+OLy
457     DO i=1-OLx,sNx+OLx
458     HSALT(i,j,bi,bj)=HEFF(i,j,1,bi,bj)*salt(i,j,1,bi,bj)*
459 dimitri 1.20 & ICE2WATR*rhoConstFresh*SEAICE_salinity
460     ENDDO
461     ENDDO
462     ENDDO
463     ENDDO
464    
465 dimitri 1.18 C-- Read initial sea ice salinity from file if available.
466     IF ( HsaltFile .NE. ' ' ) THEN
467     CALL READ_FLD_XY_RL( HsaltFile, ' ', ZETA, 0, myThid )
468     _EXCH_XY_R8(ZETA,myThid)
469     DO bj=myByLo(myThid),myByHi(myThid)
470     DO bi=myBxLo(myThid),myBxHi(myThid)
471     DO j=1-OLy,sNy+OLy
472     DO i=1-OLx,sNx+OLx
473 jmc 1.23 HSALT(i,j,bi,bj) = ZETA(i,j,bi,bj)
474 dimitri 1.18 ENDDO
475     ENDDO
476     ENDDO
477     ENDDO
478     ENDIF
479     #endif /* SEAICE_SALINITY */
480 jmc 1.23
481 dimitri 1.33 #ifdef SEAICE_AGE
482     C-- Read initial sea ice age from file if available.
483     IF ( IceAgeFile .NE. ' ' ) THEN
484     CALL READ_FLD_XY_RL( IceAgeFile, ' ', ZETA, 0, myThid )
485     _EXCH_XY_R8(ZETA,myThid)
486     DO bj=myByLo(myThid),myByHi(myThid)
487     DO bi=myBxLo(myThid),myBxHi(myThid)
488     DO j=1-OLy,sNy+OLy
489     DO i=1-OLx,sNx+OLx
490     ICEAGE(i,j,bi,bj) = ZETA(i,j,bi,bj)
491     ENDDO
492     ENDDO
493     ENDDO
494     ENDDO
495     ENDIF
496     #endif /* SEAICE_AGE */
497    
498 heimbach 1.1 ENDIF
499    
500     C--- Complete initialization
501 mlosch 1.8 PSTAR = SEAICE_strength
502 heimbach 1.1 DO bj=myByLo(myThid),myByHi(myThid)
503     DO bi=myBxLo(myThid),myBxHi(myThid)
504     DO j=1-OLy,sNy+OLy
505     DO i=1-OLx,sNx+OLx
506 jmc 1.23 ZETA(i,j,bi,bj)=HEFF(i,j,1,bi,bj)*(1.0 _d 11)
507     ETA(i,j,bi,bj)=ZETA(i,j,bi,bj)/4.0 _d 0
508     PRESS0(i,j,bi,bj)=PSTAR*HEFF(i,j,1,bi,bj)
509     & *EXP(-20.0 _d 0*(ONE-AREA(i,j,1,bi,bj)))
510     ZMAX(i,j,bi,bj)=(5.0 _d +12/(2.0 _d +04))*PRESS0(i,j,bi,bj)
511     ZMIN(i,j,bi,bj)=SEAICE_zetaMin
512     PRESS0(i,j,bi,bj)=PRESS0(i,j,bi,bj)*HEFFM(i,j,bi,bj)
513 heimbach 1.1 ENDDO
514     ENDDO
515     IF ( useRealFreshWaterFlux .AND. .NOT.useThSIce ) THEN
516     DO j=1-OLy,sNy+OLy
517     DO i=1-OLx,sNx+OLx
518 jmc 1.23 sIceLoad(i,j,bi,bj) = HEFF(i,j,1,bi,bj)*SEAICE_rhoIce
519 mlosch 1.36 & + HSNOW(i,j,bi,bj)*SEAICE_rhoSnow
520 jmc 1.23
521 heimbach 1.1 ENDDO
522     ENDDO
523     ENDIF
524     ENDDO
525     ENDDO
526    
527     RETURN
528     END

  ViewVC Help
Powered by ViewVC 1.1.22