/[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.14 - (hide annotations) (download)
Tue Jul 24 23:14:47 2007 UTC (16 years, 10 months ago) by jmc
Branch: MAIN
Changes since 1.13: +2 -2 lines
fix typo in latest check-in

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

  ViewVC Help
Powered by ViewVC 1.1.22