/[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.16 - (hide annotations) (download)
Mon Aug 6 19:36:53 2007 UTC (16 years, 10 months ago) by dimitri
Branch: MAIN
CVS Tags: checkpoint59f
Changes since 1.15: +74 -38 lines
Close southern and western boundaries for UICE/VICE when useOBCS=.TRUE.

1 dimitri 1.16 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_init_varia.F,v 1.15 2007/07/25 08:51:29 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 dimitri 1.16 IF (useOBCS) THEN
98 dimitri 1.11 C Apply OBCS T/S mask (taken from pkg/obcs/obcs_apply_ts.F)
99     #ifdef ALLOW_OBCS_NORTH
100 dimitri 1.16 DO I=1-Olx,sNx+Olx
101 dimitri 1.11 C Northern boundary
102 dimitri 1.16 J_obc = OB_Jn(I,bi,bj)
103     IF (J_obc.NE.0) THEN
104     obc_mask = _maskS(I,J_obc,kSurface,bi,bj)
105     DO J = J_obc, J_obc+Oly
106     IF (obc_mask.eq.0.) HEFFM(i,j,bi,bj)=0. _d 0
107     ENDDO
108     ENDIF
109     ENDDO
110 dimitri 1.11 #endif
111     #ifdef ALLOW_OBCS_SOUTH
112 dimitri 1.16 DO I=1-Olx,sNx+Olx
113 dimitri 1.11 C Southern boundary
114 dimitri 1.16 J_obc = OB_Js(I,bi,bj)
115     IF (J_obc.NE.0) THEN
116     obc_mask = _maskS(I,J_obc+1,kSurface,bi,bj)
117     DO J = J_obc-Oly, J_obc
118     IF (obc_mask.eq.0.) HEFFM(i,j,bi,bj)=0. _d 0
119     ENDDO
120     ENDIF
121     ENDDO
122 dimitri 1.11 #endif
123     #ifdef ALLOW_OBCS_EAST
124 dimitri 1.16 DO J=1-Oly,sNy+Oly
125 dimitri 1.11 C Eastern boundary
126 dimitri 1.16 I_obc = OB_Ie(J,bi,bj)
127     IF (I_obc.NE.0) THEN
128     obc_mask = _maskW(I_obc,J,kSurface,bi,bj)
129     DO I = I_obc, I_obc+Olx
130     IF (obc_mask.eq.0.) HEFFM(i,j,bi,bj)=0. _d 0
131     ENDDO
132     ENDIF
133     ENDDO
134 dimitri 1.11 #endif
135     #ifdef ALLOW_OBCS_WEST
136 dimitri 1.16 DO J=1-Oly,sNy+Oly
137 dimitri 1.11 C Western boundary
138 dimitri 1.16 I_obc=OB_Iw(J,bi,bj)
139     IF (I_obc.NE.0) THEN
140     obc_mask = _maskW(I_obc+1,J,kSurface,bi,bj)
141     DO I = I_obc-Olx, I_obc
142     IF (obc_mask.eq.0.) HEFFM(i,j,bi,bj)=0. _d 0
143     ENDDO
144     ENDIF
145     ENDDO
146 dimitri 1.11 #endif
147 dimitri 1.16 _EXCH_XY_R8(HEFFM, myThid)
148     ENDIF
149 dimitri 1.11 #endif /* ALLOW_OBCS */
150    
151 heimbach 1.1 DO J=1-OLy+1,sNy+OLy
152     DO I=1-OLx+1,sNx+OLx
153 dimitri 1.11 #ifdef SEAICE_CGRID
154 heimbach 1.1 seaiceMaskU(I,J,bi,bj)= 0.0 _d 0
155     seaiceMaskV(I,J,bi,bj)= 0.0 _d 0
156     mask_uice=HEFFM(I,J,bi,bj)+HEFFM(I-1,J ,bi,bj)
157 dimitri 1.12 IF(mask_uice.GT.1.5) seaiceMaskU(I,J,bi,bj)=ONE
158 heimbach 1.1 mask_uice=HEFFM(I,J,bi,bj)+HEFFM(I ,J-1,bi,bj)
159     IF(mask_uice.GT.1.5) seaiceMaskV(I,J,bi,bj)=ONE
160 dimitri 1.11 #else
161     UVM(i,j,bi,bj)=0. _d 0
162     mask_uice=HEFFM(I,J, bi,bj)+HEFFM(I-1,J-1,bi,bj)
163     & +HEFFM(I,J-1,bi,bj)+HEFFM(I-1,J, bi,bj)
164     IF(mask_uice.GT.3.5) UVM(I,J,bi,bj)=ONE
165     #endif /* SEAICE_CGRID */
166 heimbach 1.1 ENDDO
167     ENDDO
168    
169 dimitri 1.16 #ifdef ALLOW_OBCS
170     IF (useOBCS) THEN
171     C-- If OBCS is turned on, close southern and western boundaries
172     #ifdef ALLOW_OBCS_SOUTH
173     DO I=1-Olx,sNx+Olx
174     C Southern boundary
175     J_obc = OB_Js(I,bi,bj)
176     IF (J_obc.NE.0) THEN
177     #ifdef SEAICE_CGRID
178     seaiceMaskU(I,J_obc,bi,bj)= 0.0 _d 0
179     seaiceMaskV(I,J_obc,bi,bj)= 0.0 _d 0
180     #else
181     UVM(I,J_obc,bi,bj)=0. _d 0
182     #endif
183     ENDIF
184     ENDDO
185     #endif /* ALLOW_OBCS_SOUTH */
186     #ifdef ALLOW_OBCS_WEST
187     DO J=1-Oly,sNy+Oly
188     C Western boundary
189     I_obc=OB_Iw(J,bi,bj)
190     IF (I_obc.NE.0) THEN
191     #ifdef SEAICE_CGRID
192     seaiceMaskU(I_obc,J,bi,bj)= 0.0 _d 0
193     seaiceMaskV(I_obc,J,bi,bj)= 0.0 _d 0
194     #else
195     UVM(I_obc,J,bi,bj)=0. _d 0
196     #endif
197     ENDIF
198     ENDDO
199     #endif /* ALLOW_OBCS_WEST */
200     ENDIF
201     #endif /* ALLOW_OBCS */
202    
203 heimbach 1.1 #ifdef ALLOW_EXCH2
204 dimitri 1.11 #ifndef SEAICE_CGRID
205 heimbach 1.1 C-- Special stuff for cubed sphere: assume grid is rectangular and
206     C set UV mask to zero except for Arctic and Antarctic cube faces.
207     IF (useCubedSphereExchange) THEN
208     myTile = W2_myTileList(bi)
209     IF ( exch2_myFace(myTile) .EQ. 1 .OR.
210     & exch2_myFace(myTile) .EQ. 2 .OR.
211     & exch2_myFace(myTile) .EQ. 4 .OR.
212     & exch2_myFace(myTile) .EQ. 5 ) THEN
213     DO J=1-OLy,sNy+OLy
214     DO I=1-OLx,sNx+OLx
215 heimbach 1.10 UVM(i,j,bi,bj)=0. _d 0
216 heimbach 1.1 ENDDO
217     ENDDO
218     ELSEIF ( exch2_isWedge(myTile) .EQ. 1 ) THEN
219     I=1
220     DO J=1-OLy,sNy+OLy
221 heimbach 1.10 UVM(i,j,bi,bj)=0. _d 0
222 heimbach 1.1 ENDDO
223     ELSEIF ( exch2_isSedge(myTile) .EQ. 1 ) THEN
224     J=1
225     DO I=1-OLx,sNx+OLx
226 heimbach 1.10 UVM(i,j,bi,bj)=0. _d 0
227 heimbach 1.1 ENDDO
228     ENDIF
229     ENDIF
230 dimitri 1.11 #endif /* SEAICE_CGRID */
231 heimbach 1.1 #endif /* ALLOW_EXCH2 */
232    
233     DO j=1-OLy,sNy+OLy
234     DO i=1-OLx,sNx+OLx
235     TICE(I,J,bi,bj)=273.0 _d 0
236 mlosch 1.2 #ifdef SEAICE_MULTICATEGORY
237 heimbach 1.1 DO k=1,MULTDIM
238     TICES(I,J,k,bi,bj)=273.0 _d 0
239     ENDDO
240 mlosch 1.2 #endif /* SEAICE_MULTICATEGORY */
241 heimbach 1.10 UICEC (I,J,bi,bj)=0. _d 0
242     VICEC (I,J,bi,bj)=0. _d 0
243     DWATN (I,J,bi,bj)=0. _d 0
244 heimbach 1.1 #ifndef SEAICE_CGRID
245 heimbach 1.10 DAIRN (I,J,bi,bj)=0. _d 0
246 heimbach 1.1 AMASS (I,J,bi,bj)=1000.0 _d 0
247     #else
248     seaiceMassC(I,J,bi,bj)=1000.0 _d 0
249     seaiceMassU(I,J,bi,bj)=1000.0 _d 0
250     seaiceMassV(I,J,bi,bj)=1000.0 _d 0
251     #endif
252 heimbach 1.10 GWATX (I,J,bi,bj)=0. _d 0
253     GWATY (I,J,bi,bj)=0. _d 0
254 heimbach 1.1 ENDDO
255     ENDDO
256    
257     C-- Choose a proxy level for geostrophic velocity,
258     DO j=1-OLy,sNy+OLy
259     DO i=1-OLx,sNx+OLx
260     #ifdef SEAICE_TEST_ICE_STRESS_1
261     KGEO(I,J,bi,bj) = 1
262     #else /* SEAICE_TEST_ICE_STRESS_1 */
263     IF (klowc(i,j,bi,bj) .LT. 2) THEN
264     KGEO(I,J,bi,bj) = 1
265     ELSE
266     KGEO(I,J,bi,bj) = 2
267     DO WHILE ( abs(rC(KGEO(I,J,bi,bj))) .LT. 50.0 .AND.
268     & KGEO(I,J,bi,bj) .LT. (klowc(i,j,bi,bj)-1) )
269     KGEO(I,J,bi,bj) = KGEO(I,J,bi,bj) + 1
270     ENDDO
271     ENDIF
272     #endif /* SEAICE_TEST_ICE_STRESS_1 */
273     ENDDO
274     ENDDO
275    
276     ENDDO
277     ENDDO
278    
279     C-- Update overlap regions
280     #ifdef SEAICE_CGRID
281     CALL EXCH_UV_XY_RL(seaiceMaskU,seaiceMaskV,.FALSE.,myThid)
282     #else
283     _EXCH_XY_R8(UVM, myThid)
284     #endif
285    
286     C-- Now lets look at all these beasts
287     IF ( debugLevel .GE. debLevB ) THEN
288     myIter=0
289     CALL PLOT_FIELD_XYRL( HEFFM , 'Current HEFFM ' ,
290     & myIter, myThid )
291     #ifdef SEAICE_CGRID
292     CALL PLOT_FIELD_XYRL( seaiceMaskU, 'Current seaiceMaskU',
293     & myIter, myThid )
294     CALL PLOT_FIELD_XYRL( seaiceMaskV, 'Current seaiceMaskV',
295     & myIter, myThid )
296     #else
297     CALL PLOT_FIELD_XYRL( UVM , 'Current UVM ' ,
298     & myIter, myThid )
299     #endif
300     ENDIF
301    
302     #if (defined (SEAICE_CGRID) && defined (SEAICE_ALLOW_EVP))
303     DO bj=myByLo(myThid),myByHi(myThid)
304     DO bi=myBxLo(myThid),myBxHi(myThid)
305     DO j=1-OLy,sNy+OLy
306     DO i=1-OLx,sNx+OLx
307 mlosch 1.5 stressDivergenceX(I,J,bi,bj) = 0. _d 0
308     stressDivergenceY(I,J,bi,bj) = 0. _d 0
309 heimbach 1.1 seaice_sigma1 (I,J,bi,bj) = 0. _d 0
310     seaice_sigma2 (I,J,bi,bj) = 0. _d 0
311     seaice_sigma12(I,J,bi,bj) = 0. _d 0
312     ENDDO
313     ENDDO
314     ENDDO
315     ENDDO
316     #endif /* SEAICE_ALLOW_EVP and SEAICE_CGRID */
317    
318     C-- Set model variables to initial/restart conditions
319     IF ( nIter0 .NE. 0 ) THEN
320    
321     CALL SEAICE_READ_PICKUP ( myThid )
322    
323     ELSE
324    
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     YNEG(I,J,bi,bj)=ZERO
330     TMIX(I,J,bi,bj)=TICE(I,J,bi,bj)
331     DO k=1,3
332     HEFF(I,J,k,bi,bj)=SEAICE_initialHEFF*HEFFM(i,j,bi,bj)
333     UICE(I,J,k,bi,bj)=ZERO
334     VICE(I,J,k,bi,bj)=ZERO
335     ENDDO
336     ENDDO
337     ENDDO
338     ENDDO
339     ENDDO
340    
341     C-- Read initial sea-ice thickness from file if available.
342     IF ( HeffFile .NE. ' ' ) THEN
343     _BEGIN_MASTER( myThid )
344     CALL READ_FLD_XY_RL( HeffFile, ' ', ZETA, 0, myThid )
345     _END_MASTER(myThid)
346     _EXCH_XY_R8(ZETA,myThid)
347     DO bj=myByLo(myThid),myByHi(myThid)
348     DO bi=myBxLo(myThid),myBxHi(myThid)
349     DO j=1-OLy,sNy+OLy
350     DO i=1-OLx,sNx+OLx
351     DO k=1,3
352     HEFF(I,J,k,bi,bj) = MAX(ZETA(i,j,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     DO k=1,3
365 mlosch 1.7 IF(HEFF(I,J,k,bi,bj).GT.ZERO)
366     & AREA(I,J,k,bi,bj)=ONE
367 heimbach 1.1 ENDDO
368     ENDDO
369     ENDDO
370     ENDDO
371     ENDDO
372    
373 dimitri 1.9 C-- Read initial area thickness from file if available.
374 mlosch 1.7 IF ( AreaFile .NE. ' ' ) THEN
375     _BEGIN_MASTER( myThid )
376     CALL READ_FLD_XY_RL( AreaFile, ' ', ZETA, 0, myThid )
377     _END_MASTER(myThid)
378     _EXCH_XY_R8(ZETA,myThid)
379     DO bj=myByLo(myThid),myByHi(myThid)
380     DO bi=myBxLo(myThid),myBxHi(myThid)
381     DO j=1-OLy,sNy+OLy
382     DO i=1-OLx,sNx+OLx
383     DO k=1,3
384     AREA(I,J,k,bi,bj) = MAX(ZETA(i,j,bi,bj),ZERO)
385     AREA(I,J,k,bi,bj) = MIN(AREA(i,j,k,bi,bj),ONE)
386     IF ( AREA(I,J,k,bi,bj) .LE. ZERO )
387     & HEFF(I,J,k,bi,bj) = ZERO
388     IF ( HEFF(I,J,k,bi,bj) .LE. ZERO )
389     & AREA(I,J,k,bi,bj) = ZERO
390     ENDDO
391     ENDDO
392     ENDDO
393     ENDDO
394     ENDDO
395     ENDIF
396    
397     DO bj=myByLo(myThid),myByHi(myThid)
398     DO bi=myBxLo(myThid),myBxHi(myThid)
399     DO j=1-OLy,sNy+OLy
400     DO i=1-OLx,sNx+OLx
401 dimitri 1.9 HSNOW(I,J,bi,bj)=0.2*AREA(i,j,1,bi,bj)
402 mlosch 1.7 ENDDO
403     ENDDO
404     ENDDO
405     ENDDO
406 dimitri 1.9
407     C-- Read initial snow thickness from file if available.
408     IF ( HsnowFile .NE. ' ' ) THEN
409     _BEGIN_MASTER( myThid )
410     CALL READ_FLD_XY_RL( HsnowFile, ' ', ZETA, 0, myThid )
411     _END_MASTER(myThid)
412     _EXCH_XY_R8(ZETA,myThid)
413     DO bj=myByLo(myThid),myByHi(myThid)
414     DO bi=myBxLo(myThid),myBxHi(myThid)
415     DO j=1-OLy,sNy+OLy
416     DO i=1-OLx,sNx+OLx
417     HSNOW(I,J,bi,bj) = MAX(ZETA(i,j,bi,bj),ZERO)
418     ENDDO
419     ENDDO
420     ENDDO
421     ENDDO
422     ENDIF
423 mlosch 1.7
424 heimbach 1.1 ENDIF
425    
426     C--- Complete initialization
427 mlosch 1.8 PSTAR = SEAICE_strength
428 heimbach 1.1 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     ZETA(I,J,bi,bj)=HEFF(I,J,1,bi,bj)*(1.0 _d 11)
433     ETA(I,J,bi,bj)=ZETA(I,J,bi,bj)/4.0 _d 0
434 mlosch 1.8 PRESS0(I,J,bi,bj)=PSTAR*HEFF(I,J,1,bi,bj)
435     & *EXP(-20.0 _d 0*(ONE-AREA(I,J,1,bi,bj)))
436     ZMAX(I,J,bi,bj)=(5.0 _d +12/(2.0 _d +04))*PRESS0(I,J,bi,bj)
437     ZMIN(I,J,bi,bj)=0.0 _d 0
438     PRESS0(I,J,bi,bj)=PRESS0(I,J,bi,bj)*HEFFM(I,J,bi,bj)
439 heimbach 1.1 ENDDO
440     ENDDO
441     IF ( useRealFreshWaterFlux .AND. .NOT.useThSIce ) THEN
442     DO j=1-OLy,sNy+OLy
443     DO i=1-OLx,sNx+OLx
444     sIceLoad(i,j,bi,bj) = HEFF(I,J,1,bi,bj)*SEAICE_rhoIce
445     & + HSNOW(I,J,bi,bj)* 330. _d 0
446    
447     ENDDO
448     ENDDO
449     ENDIF
450     ENDDO
451     ENDDO
452    
453     RETURN
454     END

  ViewVC Help
Powered by ViewVC 1.1.22