/[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.21 - (hide annotations) (download)
Fri Sep 21 06:46:30 2007 UTC (16 years, 8 months ago) by dimitri
Branch: MAIN
CVS Tags: checkpoint59h
Changes since 1.20: +3 -2 lines
define SEAICE_salinity as a fraction of the model's surface level salinity

1 dimitri 1.21 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_init_varia.F,v 1.20 2007/09/20 14:32:37 dimitri 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    
15     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     CML#include "SEAICE_GRID.h"
23     #include "SEAICE_DIAGS.h"
24     #include "SEAICE_PARAMS.h"
25     #include "FFIELDS.h"
26     #ifdef ALLOW_EXCH2
27     #include "W2_EXCH2_TOPOLOGY.h"
28     #include "W2_EXCH2_PARAMS.h"
29     #endif /* ALLOW_EXCH2 */
30 dimitri 1.11 #ifdef ALLOW_OBCS
31     #include "OBCS.h"
32 jmc 1.14 #include "OBCS_OPTIONS.h"
33 dimitri 1.11 #endif /* ALLOW_OBCS */
34 heimbach 1.1
35     C === Routine arguments ===
36     C myThid - Thread no. that called this routine.
37     INTEGER myThid
38     CEndOfInterface
39    
40     C === Local variables ===
41     C i,j,k,bi,bj - Loop counters
42    
43     INTEGER i, j, k, bi, bj
44 mlosch 1.8 _RL PSTAR
45 heimbach 1.1 _RS mask_uice
46     INTEGER myIter, myTile
47    
48 dimitri 1.11 #ifdef ALLOW_OBCS
49     INTEGER I_obc, J_obc, kSurface
50     _RL obc_mask
51    
52     if ( buoyancyRelation .eq. 'OCEANICP' ) then
53     kSurface = Nr
54     else
55     kSurface = 1
56     endif
57     #endif /* ALLOW_OBCS */
58    
59 heimbach 1.1 cph(
60     cph make sure TAF sees proper initialisation
61     cph to avoid partial recomputation issues
62     DO bj=myByLo(myThid),myByHi(myThid)
63     DO bi=myBxLo(myThid),myBxHi(myThid)
64     c
65     DO K=1,3
66     DO J=1-OLy,sNy+OLy
67     DO I=1-OLx,sNx+OLx
68 heimbach 1.10 HEFF(I,J,k,bi,bj)=0. _d 0
69     AREA(I,J,k,bi,bj)=0. _d 0
70     UICE(I,J,k,bi,bj)=0. _d 0
71     VICE(I,J,k,bi,bj)=0. _d 0
72 heimbach 1.1 ENDDO
73     ENDDO
74     ENDDO
75     c
76     DO J=1-OLy,sNy+OLy
77     DO I=1-OLx,sNx+OLx
78 heimbach 1.10 ZETA(I,J,bi,bj)=0. _d 0
79     HSNOW(I,J,bi,bj)=0. _d 0
80 dimitri 1.18 #ifdef SEAICE_SALINITY
81     HSALT(I,J,bi,bj)=0. _d 0
82     #endif
83 heimbach 1.1 ENDDO
84     ENDDO
85     c
86     ENDDO
87     ENDDO
88     cph)
89    
90     C-- Initialize grid info
91     DO bj=myByLo(myThid),myByHi(myThid)
92     DO bi=myBxLo(myThid),myBxHi(myThid)
93     DO j=1-OLy,sNy+OLy
94     DO i=1-OLx,sNx+OLx
95     HEFFM(i,j,bi,bj)=ONE
96 heimbach 1.10 IF (_hFacC(i,j,1,bi,bj).eq.0.) HEFFM(i,j,bi,bj)=0. _d 0
97 heimbach 1.1 ENDDO
98     ENDDO
99 dimitri 1.11
100     #ifdef ALLOW_OBCS
101 dimitri 1.16 IF (useOBCS) THEN
102 dimitri 1.11 C Apply OBCS T/S mask (taken from pkg/obcs/obcs_apply_ts.F)
103 dimitri 1.17 C Additionally apply a mask along all open boundaries
104     C Until proper seaice open boundaries are implemented, this
105     C masking is needed to avoid communication across open boundaries.
106 dimitri 1.11 #ifdef ALLOW_OBCS_NORTH
107 dimitri 1.16 DO I=1-Olx,sNx+Olx
108 dimitri 1.11 C Northern boundary
109 dimitri 1.16 J_obc = OB_Jn(I,bi,bj)
110     IF (J_obc.NE.0) THEN
111 dimitri 1.17 HEFFM(I,J_obc,bi,bj)=0. _d 0
112 dimitri 1.16 obc_mask = _maskS(I,J_obc,kSurface,bi,bj)
113     DO J = J_obc, J_obc+Oly
114     IF (obc_mask.eq.0.) HEFFM(i,j,bi,bj)=0. _d 0
115     ENDDO
116     ENDIF
117     ENDDO
118 dimitri 1.11 #endif
119     #ifdef ALLOW_OBCS_SOUTH
120 dimitri 1.16 DO I=1-Olx,sNx+Olx
121 dimitri 1.11 C Southern boundary
122 dimitri 1.16 J_obc = OB_Js(I,bi,bj)
123     IF (J_obc.NE.0) THEN
124 dimitri 1.17 HEFFM(I,J_obc,bi,bj)=0. _d 0
125 dimitri 1.16 obc_mask = _maskS(I,J_obc+1,kSurface,bi,bj)
126     DO J = J_obc-Oly, J_obc
127     IF (obc_mask.eq.0.) HEFFM(i,j,bi,bj)=0. _d 0
128     ENDDO
129     ENDIF
130     ENDDO
131 dimitri 1.11 #endif
132     #ifdef ALLOW_OBCS_EAST
133 dimitri 1.16 DO J=1-Oly,sNy+Oly
134 dimitri 1.11 C Eastern boundary
135 dimitri 1.16 I_obc = OB_Ie(J,bi,bj)
136     IF (I_obc.NE.0) THEN
137 dimitri 1.17 HEFFM(I_obc,J,bi,bj)=0. _d 0
138 dimitri 1.16 obc_mask = _maskW(I_obc,J,kSurface,bi,bj)
139     DO I = I_obc, I_obc+Olx
140     IF (obc_mask.eq.0.) HEFFM(i,j,bi,bj)=0. _d 0
141     ENDDO
142     ENDIF
143     ENDDO
144 dimitri 1.11 #endif
145     #ifdef ALLOW_OBCS_WEST
146 dimitri 1.16 DO J=1-Oly,sNy+Oly
147 dimitri 1.11 C Western boundary
148 dimitri 1.16 I_obc=OB_Iw(J,bi,bj)
149     IF (I_obc.NE.0) THEN
150 dimitri 1.17 HEFFM(I_obc,J,bi,bj)=0. _d 0
151 dimitri 1.16 obc_mask = _maskW(I_obc+1,J,kSurface,bi,bj)
152     DO I = I_obc-Olx, I_obc
153     IF (obc_mask.eq.0.) HEFFM(i,j,bi,bj)=0. _d 0
154     ENDDO
155     ENDIF
156     ENDDO
157 dimitri 1.11 #endif
158 dimitri 1.16 _EXCH_XY_R8(HEFFM, myThid)
159     ENDIF
160 dimitri 1.11 #endif /* ALLOW_OBCS */
161    
162 heimbach 1.1 DO J=1-OLy+1,sNy+OLy
163     DO I=1-OLx+1,sNx+OLx
164 dimitri 1.11 #ifdef SEAICE_CGRID
165 heimbach 1.1 seaiceMaskU(I,J,bi,bj)= 0.0 _d 0
166     seaiceMaskV(I,J,bi,bj)= 0.0 _d 0
167     mask_uice=HEFFM(I,J,bi,bj)+HEFFM(I-1,J ,bi,bj)
168 dimitri 1.12 IF(mask_uice.GT.1.5) seaiceMaskU(I,J,bi,bj)=ONE
169 heimbach 1.1 mask_uice=HEFFM(I,J,bi,bj)+HEFFM(I ,J-1,bi,bj)
170     IF(mask_uice.GT.1.5) seaiceMaskV(I,J,bi,bj)=ONE
171 dimitri 1.11 #else
172     UVM(i,j,bi,bj)=0. _d 0
173     mask_uice=HEFFM(I,J, bi,bj)+HEFFM(I-1,J-1,bi,bj)
174     & +HEFFM(I,J-1,bi,bj)+HEFFM(I-1,J, bi,bj)
175     IF(mask_uice.GT.3.5) UVM(I,J,bi,bj)=ONE
176     #endif /* SEAICE_CGRID */
177 heimbach 1.1 ENDDO
178     ENDDO
179    
180     #ifdef ALLOW_EXCH2
181 dimitri 1.11 #ifndef SEAICE_CGRID
182 heimbach 1.1 C-- Special stuff for cubed sphere: assume grid is rectangular and
183     C set UV mask to zero except for Arctic and Antarctic cube faces.
184     IF (useCubedSphereExchange) THEN
185     myTile = W2_myTileList(bi)
186     IF ( exch2_myFace(myTile) .EQ. 1 .OR.
187     & exch2_myFace(myTile) .EQ. 2 .OR.
188     & exch2_myFace(myTile) .EQ. 4 .OR.
189     & exch2_myFace(myTile) .EQ. 5 ) THEN
190     DO J=1-OLy,sNy+OLy
191     DO I=1-OLx,sNx+OLx
192 heimbach 1.10 UVM(i,j,bi,bj)=0. _d 0
193 heimbach 1.1 ENDDO
194     ENDDO
195     ELSEIF ( exch2_isWedge(myTile) .EQ. 1 ) THEN
196     I=1
197     DO J=1-OLy,sNy+OLy
198 heimbach 1.10 UVM(i,j,bi,bj)=0. _d 0
199 heimbach 1.1 ENDDO
200     ELSEIF ( exch2_isSedge(myTile) .EQ. 1 ) THEN
201     J=1
202     DO I=1-OLx,sNx+OLx
203 heimbach 1.10 UVM(i,j,bi,bj)=0. _d 0
204 heimbach 1.1 ENDDO
205     ENDIF
206     ENDIF
207 dimitri 1.11 #endif /* SEAICE_CGRID */
208 heimbach 1.1 #endif /* ALLOW_EXCH2 */
209    
210     DO j=1-OLy,sNy+OLy
211     DO i=1-OLx,sNx+OLx
212     TICE(I,J,bi,bj)=273.0 _d 0
213 mlosch 1.2 #ifdef SEAICE_MULTICATEGORY
214 heimbach 1.1 DO k=1,MULTDIM
215     TICES(I,J,k,bi,bj)=273.0 _d 0
216     ENDDO
217 mlosch 1.2 #endif /* SEAICE_MULTICATEGORY */
218 heimbach 1.10 UICEC (I,J,bi,bj)=0. _d 0
219     VICEC (I,J,bi,bj)=0. _d 0
220     DWATN (I,J,bi,bj)=0. _d 0
221 heimbach 1.1 #ifndef SEAICE_CGRID
222 heimbach 1.10 DAIRN (I,J,bi,bj)=0. _d 0
223 heimbach 1.1 AMASS (I,J,bi,bj)=1000.0 _d 0
224     #else
225     seaiceMassC(I,J,bi,bj)=1000.0 _d 0
226     seaiceMassU(I,J,bi,bj)=1000.0 _d 0
227     seaiceMassV(I,J,bi,bj)=1000.0 _d 0
228     #endif
229 heimbach 1.10 GWATX (I,J,bi,bj)=0. _d 0
230     GWATY (I,J,bi,bj)=0. _d 0
231 heimbach 1.1 ENDDO
232     ENDDO
233    
234     C-- Choose a proxy level for geostrophic velocity,
235     DO j=1-OLy,sNy+OLy
236     DO i=1-OLx,sNx+OLx
237     #ifdef SEAICE_TEST_ICE_STRESS_1
238     KGEO(I,J,bi,bj) = 1
239     #else /* SEAICE_TEST_ICE_STRESS_1 */
240     IF (klowc(i,j,bi,bj) .LT. 2) THEN
241     KGEO(I,J,bi,bj) = 1
242     ELSE
243     KGEO(I,J,bi,bj) = 2
244     DO WHILE ( abs(rC(KGEO(I,J,bi,bj))) .LT. 50.0 .AND.
245     & KGEO(I,J,bi,bj) .LT. (klowc(i,j,bi,bj)-1) )
246     KGEO(I,J,bi,bj) = KGEO(I,J,bi,bj) + 1
247     ENDDO
248     ENDIF
249     #endif /* SEAICE_TEST_ICE_STRESS_1 */
250     ENDDO
251     ENDDO
252    
253     ENDDO
254     ENDDO
255    
256     C-- Update overlap regions
257     #ifdef SEAICE_CGRID
258     CALL EXCH_UV_XY_RL(seaiceMaskU,seaiceMaskV,.FALSE.,myThid)
259     #else
260     _EXCH_XY_R8(UVM, myThid)
261     #endif
262    
263     C-- Now lets look at all these beasts
264     IF ( debugLevel .GE. debLevB ) THEN
265     myIter=0
266     CALL PLOT_FIELD_XYRL( HEFFM , 'Current HEFFM ' ,
267     & myIter, myThid )
268     #ifdef SEAICE_CGRID
269     CALL PLOT_FIELD_XYRL( seaiceMaskU, 'Current seaiceMaskU',
270     & myIter, myThid )
271     CALL PLOT_FIELD_XYRL( seaiceMaskV, 'Current seaiceMaskV',
272     & myIter, myThid )
273     #else
274     CALL PLOT_FIELD_XYRL( UVM , 'Current UVM ' ,
275     & myIter, myThid )
276     #endif
277     ENDIF
278    
279     #if (defined (SEAICE_CGRID) && defined (SEAICE_ALLOW_EVP))
280     DO bj=myByLo(myThid),myByHi(myThid)
281     DO bi=myBxLo(myThid),myBxHi(myThid)
282     DO j=1-OLy,sNy+OLy
283     DO i=1-OLx,sNx+OLx
284 mlosch 1.5 stressDivergenceX(I,J,bi,bj) = 0. _d 0
285     stressDivergenceY(I,J,bi,bj) = 0. _d 0
286 heimbach 1.1 seaice_sigma1 (I,J,bi,bj) = 0. _d 0
287     seaice_sigma2 (I,J,bi,bj) = 0. _d 0
288     seaice_sigma12(I,J,bi,bj) = 0. _d 0
289     ENDDO
290     ENDDO
291     ENDDO
292     ENDDO
293     #endif /* SEAICE_ALLOW_EVP and SEAICE_CGRID */
294    
295     C-- Set model variables to initial/restart conditions
296     IF ( nIter0 .NE. 0 ) THEN
297    
298     CALL SEAICE_READ_PICKUP ( myThid )
299    
300     ELSE
301    
302     DO bj=myByLo(myThid),myByHi(myThid)
303     DO bi=myBxLo(myThid),myBxHi(myThid)
304     DO j=1-OLy,sNy+OLy
305     DO i=1-OLx,sNx+OLx
306     YNEG(I,J,bi,bj)=ZERO
307     TMIX(I,J,bi,bj)=TICE(I,J,bi,bj)
308     DO k=1,3
309     HEFF(I,J,k,bi,bj)=SEAICE_initialHEFF*HEFFM(i,j,bi,bj)
310     UICE(I,J,k,bi,bj)=ZERO
311     VICE(I,J,k,bi,bj)=ZERO
312     ENDDO
313     ENDDO
314     ENDDO
315     ENDDO
316     ENDDO
317    
318     C-- Read initial sea-ice thickness from file if available.
319     IF ( HeffFile .NE. ' ' ) THEN
320     _BEGIN_MASTER( myThid )
321     CALL READ_FLD_XY_RL( HeffFile, ' ', ZETA, 0, myThid )
322     _END_MASTER(myThid)
323     _EXCH_XY_R8(ZETA,myThid)
324     DO bj=myByLo(myThid),myByHi(myThid)
325     DO bi=myBxLo(myThid),myBxHi(myThid)
326     DO j=1-OLy,sNy+OLy
327     DO i=1-OLx,sNx+OLx
328     DO k=1,3
329     HEFF(I,J,k,bi,bj) = MAX(ZETA(i,j,bi,bj),ZERO)
330     ENDDO
331     ENDDO
332     ENDDO
333     ENDDO
334     ENDDO
335     ENDIF
336    
337     DO bj=myByLo(myThid),myByHi(myThid)
338     DO bi=myBxLo(myThid),myBxHi(myThid)
339     DO j=1-OLy,sNy+OLy
340     DO i=1-OLx,sNx+OLx
341     DO k=1,3
342 mlosch 1.7 IF(HEFF(I,J,k,bi,bj).GT.ZERO)
343     & AREA(I,J,k,bi,bj)=ONE
344 heimbach 1.1 ENDDO
345     ENDDO
346     ENDDO
347     ENDDO
348     ENDDO
349    
350 dimitri 1.9 C-- Read initial area thickness from file if available.
351 mlosch 1.7 IF ( AreaFile .NE. ' ' ) THEN
352     _BEGIN_MASTER( myThid )
353     CALL READ_FLD_XY_RL( AreaFile, ' ', ZETA, 0, myThid )
354     _END_MASTER(myThid)
355     _EXCH_XY_R8(ZETA,myThid)
356     DO bj=myByLo(myThid),myByHi(myThid)
357     DO bi=myBxLo(myThid),myBxHi(myThid)
358     DO j=1-OLy,sNy+OLy
359     DO i=1-OLx,sNx+OLx
360     DO k=1,3
361     AREA(I,J,k,bi,bj) = MAX(ZETA(i,j,bi,bj),ZERO)
362     AREA(I,J,k,bi,bj) = MIN(AREA(i,j,k,bi,bj),ONE)
363     IF ( AREA(I,J,k,bi,bj) .LE. ZERO )
364     & HEFF(I,J,k,bi,bj) = ZERO
365     IF ( HEFF(I,J,k,bi,bj) .LE. ZERO )
366     & AREA(I,J,k,bi,bj) = ZERO
367     ENDDO
368     ENDDO
369     ENDDO
370     ENDDO
371     ENDDO
372     ENDIF
373    
374     DO bj=myByLo(myThid),myByHi(myThid)
375     DO bi=myBxLo(myThid),myBxHi(myThid)
376     DO j=1-OLy,sNy+OLy
377     DO i=1-OLx,sNx+OLx
378 dimitri 1.9 HSNOW(I,J,bi,bj)=0.2*AREA(i,j,1,bi,bj)
379 mlosch 1.7 ENDDO
380     ENDDO
381     ENDDO
382     ENDDO
383 dimitri 1.9
384     C-- Read initial snow thickness from file if available.
385     IF ( HsnowFile .NE. ' ' ) THEN
386     _BEGIN_MASTER( myThid )
387     CALL READ_FLD_XY_RL( HsnowFile, ' ', ZETA, 0, myThid )
388     _END_MASTER(myThid)
389     _EXCH_XY_R8(ZETA,myThid)
390     DO bj=myByLo(myThid),myByHi(myThid)
391     DO bi=myBxLo(myThid),myBxHi(myThid)
392     DO j=1-OLy,sNy+OLy
393     DO i=1-OLx,sNx+OLx
394     HSNOW(I,J,bi,bj) = MAX(ZETA(i,j,bi,bj),ZERO)
395     ENDDO
396     ENDDO
397     ENDDO
398     ENDDO
399     ENDIF
400 dimitri 1.18
401     #ifdef SEAICE_SALINITY
402 dimitri 1.20 DO bj=myByLo(myThid),myByHi(myThid)
403     DO bi=myBxLo(myThid),myBxHi(myThid)
404     DO J=1-OLy,sNy+OLy
405     DO I=1-OLx,sNx+OLx
406 dimitri 1.21 HSALT(I,J,bi,bj)=HEFF(I,J,1,bi,bj)*salt(I,J,1,bi,bj)*
407 dimitri 1.20 & ICE2WATR*rhoConstFresh*SEAICE_salinity
408     ENDDO
409     ENDDO
410     ENDDO
411     ENDDO
412    
413 dimitri 1.18 C-- Read initial sea ice salinity from file if available.
414     IF ( HsaltFile .NE. ' ' ) THEN
415     _BEGIN_MASTER( myThid )
416     CALL READ_FLD_XY_RL( HsaltFile, ' ', ZETA, 0, myThid )
417     _END_MASTER(myThid)
418     _EXCH_XY_R8(ZETA,myThid)
419     DO bj=myByLo(myThid),myByHi(myThid)
420     DO bi=myBxLo(myThid),myBxHi(myThid)
421     DO j=1-OLy,sNy+OLy
422     DO i=1-OLx,sNx+OLx
423 dimitri 1.20 HSALT(I,J,bi,bj) = ZETA(i,j,bi,bj)
424 dimitri 1.18 ENDDO
425     ENDDO
426     ENDDO
427     ENDDO
428     ENDIF
429     #endif /* SEAICE_SALINITY */
430 mlosch 1.7
431 heimbach 1.1 ENDIF
432    
433     C--- Complete initialization
434 mlosch 1.8 PSTAR = SEAICE_strength
435 heimbach 1.1 DO bj=myByLo(myThid),myByHi(myThid)
436     DO bi=myBxLo(myThid),myBxHi(myThid)
437     DO j=1-OLy,sNy+OLy
438     DO i=1-OLx,sNx+OLx
439     ZETA(I,J,bi,bj)=HEFF(I,J,1,bi,bj)*(1.0 _d 11)
440     ETA(I,J,bi,bj)=ZETA(I,J,bi,bj)/4.0 _d 0
441 mlosch 1.8 PRESS0(I,J,bi,bj)=PSTAR*HEFF(I,J,1,bi,bj)
442     & *EXP(-20.0 _d 0*(ONE-AREA(I,J,1,bi,bj)))
443     ZMAX(I,J,bi,bj)=(5.0 _d +12/(2.0 _d +04))*PRESS0(I,J,bi,bj)
444 mlosch 1.19 ZMIN(I,J,bi,bj)=SEAICE_zetaMin
445 mlosch 1.8 PRESS0(I,J,bi,bj)=PRESS0(I,J,bi,bj)*HEFFM(I,J,bi,bj)
446 heimbach 1.1 ENDDO
447     ENDDO
448     IF ( useRealFreshWaterFlux .AND. .NOT.useThSIce ) THEN
449     DO j=1-OLy,sNy+OLy
450     DO i=1-OLx,sNx+OLx
451     sIceLoad(i,j,bi,bj) = HEFF(I,J,1,bi,bj)*SEAICE_rhoIce
452     & + HSNOW(I,J,bi,bj)* 330. _d 0
453    
454     ENDDO
455     ENDDO
456     ENDIF
457     ENDDO
458     ENDDO
459    
460     RETURN
461     END

  ViewVC Help
Powered by ViewVC 1.1.22