/[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.22 - (hide annotations) (download)
Thu Oct 11 10:30:34 2007 UTC (16 years, 7 months ago) by dimitri
Branch: MAIN
CVS Tags: checkpoint59i, checkpoint59j
Changes since 1.21: +33 -64 lines
o Added open boundary conditions capability for seaice HEFF and AREA
For time being this capability requires pkg/exf, pkg/obcs, and pkg/seaice
UICE and VICE are masked at the southern and western open boundaries, if
any, of the domain.  Eventually needs to be extended to SALT and SNOW.

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

  ViewVC Help
Powered by ViewVC 1.1.22