/[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.46 - (hide annotations) (download)
Sun Jun 28 01:05:41 2009 UTC (14 years, 11 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint61t, checkpoint61s
Changes since 1.45: +2 -2 lines
add bj in exch2 arrays and S/R

1 jmc 1.46 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_init_varia.F,v 1.45 2009/06/25 14:36:15 heimbach 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 jmc 1.38 # include "W2_EXCH2_SIZE.h"
27 dimitri 1.24 # include "W2_EXCH2_TOPOLOGY.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 mlosch 1.44 INTEGER myIter, myTile, kSurface
46     #ifdef ALLOW_OBCS
47     INTEGER I_obc, J_obc
48     #endif /* ALLOW_OBCS */
49 heimbach 1.1
50 jmc 1.23 IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
51 mlosch 1.44 kSurface = Nr
52 jmc 1.23 ELSE
53 mlosch 1.44 kSurface = 1
54 jmc 1.23 ENDIF
55 dimitri 1.11
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 mlosch 1.42 DO j=1-OLy,sNy+OLy
60     DO i=1-OLx,sNx+OLx
61 mlosch 1.43 HEFF(i,j,bi,bj)=0. _d 0
62     AREA(i,j,bi,bj)=0. _d 0
63     UICE(i,j,bi,bj)=0. _d 0
64     VICE(i,j,bi,bj)=0. _d 0
65     C
66 mlosch 1.42 uIceNm1(i,j,bi,bj)=0. _d 0
67     vIceNm1(i,j,bi,bj)=0. _d 0
68     areaNm1(i,j,bi,bj)=0. _d 0
69     hEffNm1(i,j,bi,bj)=0. _d 0
70     ENDDO
71     ENDDO
72 jmc 1.23 #ifdef SEAICE_MULTICATEGORY
73     DO k=1,MULTDIM
74     DO j=1-OLy,sNy+OLy
75     DO i=1-OLx,sNx+OLx
76     TICES(i,j,k,bi,bj)=0. _d 0
77 heimbach 1.1 ENDDO
78     ENDDO
79     ENDDO
80 jmc 1.23 #endif
81     DO j=1-OLy,sNy+OLy
82     DO i=1-OLx,sNx+OLx
83     ETA (i,j,bi,bj) = 0. _d 0
84     ZETA(i,j,bi,bj) = 0. _d 0
85     DRAGS(i,j,bi,bj) = 0. _d 0
86     DRAGA(i,j,bi,bj) = 0. _d 0
87     FORCEX(i,j,bi,bj) = 0. _d 0
88     FORCEY(i,j,bi,bj) = 0. _d 0
89     UICEC(i,j,bi,bj) = 0. _d 0
90     VICEC(i,j,bi,bj) = 0. _d 0
91     #ifdef SEAICE_CGRID
92     seaiceMassC(i,j,bi,bj)=0. _d 0
93     seaiceMassU(i,j,bi,bj)=0. _d 0
94     seaiceMassV(i,j,bi,bj)=0. _d 0
95     stressDivergenceX(i,j,bi,bj) = 0. _d 0
96     stressDivergenceY(i,j,bi,bj) = 0. _d 0
97 heimbach 1.39 # ifdef SEAICE_ALLOW_EVP
98 jmc 1.23 seaice_sigma1 (i,j,bi,bj) = 0. _d 0
99     seaice_sigma2 (i,j,bi,bj) = 0. _d 0
100     seaice_sigma12(i,j,bi,bj) = 0. _d 0
101     # endif /* SEAICE_ALLOW_EVP */
102     #else /* SEAICE_CGRID */
103     AMASS(i,j,bi,bj) = 0. _d 0
104     DAIRN(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 jmc 1.23 #endif /* SEAICE_CGRID */
110     DWATN(i,j,bi,bj) = 0. _d 0
111     PRESS0(i,j,bi,bj) = 0. _d 0
112     FORCEX0(i,j,bi,bj)= 0. _d 0
113     FORCEY0(i,j,bi,bj)= 0. _d 0
114     ZMAX(i,j,bi,bj) = 0. _d 0
115     ZMIN(i,j,bi,bj) = 0. _d 0
116     HSNOW(i,j,bi,bj) = 0. _d 0
117 dimitri 1.18 #ifdef SEAICE_SALINITY
118 jmc 1.23 HSALT(i,j,bi,bj) = 0. _d 0
119     #endif
120 dimitri 1.33 #ifdef SEAICE_AGE
121     ICEAGE(i,j,bi,bj) = 0. _d 0
122     #endif
123 jmc 1.23 YNEG (i,j,bi,bj) = 0. _d 0
124     RIVER(i,j,bi,bj) = 0. _d 0
125     TMIX(i,j,bi,bj) = 0. _d 0
126     TICE(i,j,bi,bj) = 0. _d 0
127     TAUX(i,j,bi,bj) = 0. _d 0
128     TAUY(i,j,bi,bj) = 0. _d 0
129     #ifdef ALLOW_SEAICE_COST_EXPORT
130     uHeffExportCell(i,j,bi,bj) = 0. _d 0
131     vHeffExportCell(i,j,bi,bj) = 0. _d 0
132 dimitri 1.18 #endif
133 heimbach 1.32 saltWtrIce(i,j,bi,bj) = 0. _d 0
134     frWtrIce(i,j,bi,bj) = 0. _d 0
135 heimbach 1.1 ENDDO
136     ENDDO
137     ENDDO
138     ENDDO
139    
140 mlosch 1.40 C-- Initialize (variable) grid info. As long as we allow masking of
141     C-- velocities outside of ice covered areas (in seaice_dynsolver)
142     C-- we need to re-initialize seaiceMaskU/V here for TAF/TAMC
143     #ifdef SEAICE_CGRID
144 heimbach 1.1 DO bj=myByLo(myThid),myByHi(myThid)
145     DO bi=myBxLo(myThid),myBxHi(myThid)
146 jmc 1.23 DO j=1-OLy+1,sNy+OLy
147     DO i=1-OLx+1,sNx+OLx
148     seaiceMaskU(i,j,bi,bj)= 0.0 _d 0
149     seaiceMaskV(i,j,bi,bj)= 0.0 _d 0
150     mask_uice=HEFFM(i,j,bi,bj)+HEFFM(i-1,j ,bi,bj)
151 heimbach 1.45 IF(mask_uice.GT.1.5 _d 0) seaiceMaskU(i,j,bi,bj)=1.0 _d 0
152 jmc 1.23 mask_uice=HEFFM(i,j,bi,bj)+HEFFM(i ,j-1,bi,bj)
153 heimbach 1.45 IF(mask_uice.GT.1.5 _d 0) seaiceMaskV(i,j,bi,bj)=1.0 _d 0
154 dimitri 1.22 ENDDO
155     ENDDO
156 mlosch 1.40 ENDDO
157     ENDDO
158 mlosch 1.34 #endif /* SEAICE_CGRID */
159 dimitri 1.11
160 mlosch 1.40 DO bj=myByLo(myThid),myByHi(myThid)
161     DO bi=myBxLo(myThid),myBxHi(myThid)
162 dimitri 1.11 #ifdef ALLOW_OBCS
163 dimitri 1.16 IF (useOBCS) THEN
164 dimitri 1.22 C-- If OBCS is turned on, close southern and western boundaries
165 dimitri 1.11 #ifdef ALLOW_OBCS_SOUTH
166 jmc 1.23 DO i=1-Olx,sNx+Olx
167 dimitri 1.11 C Southern boundary
168 jmc 1.23 J_obc = OB_Js(i,bi,bj)
169 dimitri 1.16 IF (J_obc.NE.0) THEN
170 dimitri 1.22 #ifdef SEAICE_CGRID
171 jmc 1.23 seaiceMaskU(i,J_obc,bi,bj)= 0.0 _d 0
172     seaiceMaskV(i,J_obc,bi,bj)= 0.0 _d 0
173 dimitri 1.22 #else
174 jmc 1.23 UVM(i,J_obc,bi,bj)=0. _d 0
175 dimitri 1.22 #endif /* SEAICE_CGRID */
176 dimitri 1.16 ENDIF
177     ENDDO
178 dimitri 1.22 #endif /* ALLOW_OBCS_SOUTH */
179 dimitri 1.11 #ifdef ALLOW_OBCS_WEST
180 jmc 1.23 DO j=1-Oly,sNy+Oly
181 dimitri 1.11 C Western boundary
182 jmc 1.23 I_obc=OB_Iw(j,bi,bj)
183 dimitri 1.16 IF (I_obc.NE.0) THEN
184 dimitri 1.22 #ifdef SEAICE_CGRID
185 jmc 1.23 seaiceMaskU(I_obc,j,bi,bj)= 0.0 _d 0
186     seaiceMaskV(I_obc,j,bi,bj)= 0.0 _d 0
187 dimitri 1.22 #else
188 jmc 1.23 UVM(I_obc,j,bi,bj)=0. _d 0
189 dimitri 1.22 #endif /* SEAICE_CGRID */
190 dimitri 1.16 ENDIF
191     ENDDO
192 dimitri 1.22 #endif /* ALLOW_OBCS_WEST */
193 dimitri 1.16 ENDIF
194 dimitri 1.11 #endif /* ALLOW_OBCS */
195    
196 heimbach 1.1 #ifdef ALLOW_EXCH2
197 dimitri 1.11 #ifndef SEAICE_CGRID
198 heimbach 1.1 C-- Special stuff for cubed sphere: assume grid is rectangular and
199     C set UV mask to zero except for Arctic and Antarctic cube faces.
200     IF (useCubedSphereExchange) THEN
201 jmc 1.46 myTile = W2_myTileList(bi,bj)
202 heimbach 1.1 IF ( exch2_myFace(myTile) .EQ. 1 .OR.
203     & exch2_myFace(myTile) .EQ. 2 .OR.
204     & exch2_myFace(myTile) .EQ. 4 .OR.
205     & exch2_myFace(myTile) .EQ. 5 ) THEN
206 jmc 1.23 DO j=1-OLy,sNy+OLy
207     DO i=1-OLx,sNx+OLx
208 heimbach 1.10 UVM(i,j,bi,bj)=0. _d 0
209 heimbach 1.1 ENDDO
210     ENDDO
211     ELSEIF ( exch2_isWedge(myTile) .EQ. 1 ) THEN
212 jmc 1.23 i=1
213     DO j=1-OLy,sNy+OLy
214 heimbach 1.10 UVM(i,j,bi,bj)=0. _d 0
215 heimbach 1.1 ENDDO
216     ELSEIF ( exch2_isSedge(myTile) .EQ. 1 ) THEN
217 jmc 1.23 j=1
218     DO i=1-OLx,sNx+OLx
219 heimbach 1.10 UVM(i,j,bi,bj)=0. _d 0
220 heimbach 1.1 ENDDO
221     ENDIF
222     ENDIF
223 dimitri 1.11 #endif /* SEAICE_CGRID */
224 heimbach 1.1 #endif /* ALLOW_EXCH2 */
225    
226     DO j=1-OLy,sNy+OLy
227     DO i=1-OLx,sNx+OLx
228 jmc 1.23 TICE(i,j,bi,bj)=273.0 _d 0
229 mlosch 1.2 #ifdef SEAICE_MULTICATEGORY
230 heimbach 1.1 DO k=1,MULTDIM
231 jmc 1.23 TICES(i,j,k,bi,bj)=273.0 _d 0
232 heimbach 1.1 ENDDO
233 mlosch 1.2 #endif /* SEAICE_MULTICATEGORY */
234 heimbach 1.1 #ifndef SEAICE_CGRID
235 jmc 1.23 AMASS (i,j,bi,bj)=1000.0 _d 0
236 heimbach 1.1 #else
237 jmc 1.23 seaiceMassC(i,j,bi,bj)=1000.0 _d 0
238     seaiceMassU(i,j,bi,bj)=1000.0 _d 0
239     seaiceMassV(i,j,bi,bj)=1000.0 _d 0
240 heimbach 1.1 #endif
241     ENDDO
242     ENDDO
243    
244     ENDDO
245     ENDDO
246    
247     C-- Update overlap regions
248     #ifdef SEAICE_CGRID
249     CALL EXCH_UV_XY_RL(seaiceMaskU,seaiceMaskV,.FALSE.,myThid)
250     #else
251 jmc 1.37 _EXCH_XY_RL(UVM, myThid)
252 heimbach 1.1 #endif
253    
254     C-- Now lets look at all these beasts
255     IF ( debugLevel .GE. debLevB ) THEN
256     myIter=0
257     CALL PLOT_FIELD_XYRL( HEFFM , 'Current HEFFM ' ,
258     & myIter, myThid )
259     #ifdef SEAICE_CGRID
260     CALL PLOT_FIELD_XYRL( seaiceMaskU, 'Current seaiceMaskU',
261     & myIter, myThid )
262     CALL PLOT_FIELD_XYRL( seaiceMaskV, 'Current seaiceMaskV',
263     & myIter, myThid )
264     #else
265     CALL PLOT_FIELD_XYRL( UVM , 'Current UVM ' ,
266     & myIter, myThid )
267     #endif
268     ENDIF
269    
270     C-- Set model variables to initial/restart conditions
271 mlosch 1.26 IF ( .NOT. ( startTime .EQ. baseTime .AND. nIter0 .EQ. 0
272     & .AND. pickupSuff .EQ. ' ') ) THEN
273 heimbach 1.1
274     CALL SEAICE_READ_PICKUP ( myThid )
275    
276     ELSE
277    
278     DO bj=myByLo(myThid),myByHi(myThid)
279     DO bi=myBxLo(myThid),myBxHi(myThid)
280     DO j=1-OLy,sNy+OLy
281     DO i=1-OLx,sNx+OLx
282 jmc 1.23 TMIX(i,j,bi,bj)=TICE(i,j,bi,bj)
283 mlosch 1.43 HEFF(i,j,bi,bj)=SEAICE_initialHEFF*HEFFM(i,j,bi,bj)
284     UICE(i,j,bi,bj)=ZERO
285     VICE(i,j,bi,bj)=ZERO
286 heimbach 1.1 ENDDO
287     ENDDO
288     ENDDO
289     ENDDO
290    
291     C-- Read initial sea-ice thickness from file if available.
292     IF ( HeffFile .NE. ' ' ) THEN
293 mlosch 1.43 CALL READ_FLD_XY_RL( HeffFile, ' ', HEFF, 0, myThid )
294     _EXCH_XY_RL(HEFF,myThid)
295 heimbach 1.1 DO bj=myByLo(myThid),myByHi(myThid)
296     DO bi=myBxLo(myThid),myBxHi(myThid)
297     DO j=1-OLy,sNy+OLy
298     DO i=1-OLx,sNx+OLx
299 mlosch 1.43 HEFF(i,j,bi,bj) = MAX(HEFF(i,j,bi,bj),ZERO)
300 heimbach 1.1 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 mlosch 1.43 IF(HEFF(i,j,bi,bj).GT.ZERO) AREA(i,j,bi,bj)=ONE
311 heimbach 1.1 ENDDO
312     ENDDO
313     ENDDO
314     ENDDO
315 jmc 1.23
316 dimitri 1.27 C-- Read initial sea-ice area from file if available.
317 mlosch 1.7 IF ( AreaFile .NE. ' ' ) THEN
318 mlosch 1.43 CALL READ_FLD_XY_RL( AreaFile, ' ', AREA, 0, myThid )
319     _EXCH_XY_RL(AREA,myThid)
320 mlosch 1.7 DO bj=myByLo(myThid),myByHi(myThid)
321     DO bi=myBxLo(myThid),myBxHi(myThid)
322     DO j=1-OLy,sNy+OLy
323     DO i=1-OLx,sNx+OLx
324 mlosch 1.43 AREA(i,j,bi,bj) = MAX(AREA(i,j,bi,bj),ZERO)
325     AREA(i,j,bi,bj) = MIN(AREA(i,j,bi,bj),ONE)
326     IF ( AREA(i,j,bi,bj) .LE. ZERO ) HEFF(i,j,bi,bj) = ZERO
327     IF ( HEFF(i,j,bi,bj) .LE. ZERO ) AREA(i,j,bi,bj) = ZERO
328 mlosch 1.7 ENDDO
329     ENDDO
330     ENDDO
331     ENDDO
332     ENDIF
333    
334     DO bj=myByLo(myThid),myByHi(myThid)
335     DO bi=myBxLo(myThid),myBxHi(myThid)
336     DO j=1-OLy,sNy+OLy
337     DO i=1-OLx,sNx+OLx
338 mlosch 1.43 HSNOW(i,j,bi,bj)=0.2*AREA(i,j,bi,bj)
339 mlosch 1.7 ENDDO
340     ENDDO
341     ENDDO
342     ENDDO
343 dimitri 1.9
344     C-- Read initial snow thickness from file if available.
345     IF ( HsnowFile .NE. ' ' ) THEN
346 mlosch 1.43 CALL READ_FLD_XY_RL( HsnowFile, ' ', HSNOW, 0, myThid )
347     _EXCH_XY_RL(HSNOW,myThid)
348 dimitri 1.9 DO bj=myByLo(myThid),myByHi(myThid)
349     DO bi=myBxLo(myThid),myBxHi(myThid)
350     DO j=1-OLy,sNy+OLy
351     DO i=1-OLx,sNx+OLx
352 mlosch 1.43 HSNOW(i,j,bi,bj) = MAX(HSNOW(i,j,bi,bj),ZERO)
353 dimitri 1.9 ENDDO
354     ENDDO
355     ENDDO
356     ENDDO
357     ENDIF
358 dimitri 1.18
359     #ifdef SEAICE_SALINITY
360 dimitri 1.20 DO bj=myByLo(myThid),myByHi(myThid)
361     DO bi=myBxLo(myThid),myBxHi(myThid)
362 jmc 1.23 DO j=1-OLy,sNy+OLy
363     DO i=1-OLx,sNx+OLx
364 mlosch 1.44 HSALT(i,j,bi,bj)=HEFF(i,j,bi,bj)*salt(i,j,kSurface,bi,bj)*
365 dimitri 1.20 & ICE2WATR*rhoConstFresh*SEAICE_salinity
366     ENDDO
367     ENDDO
368     ENDDO
369     ENDDO
370    
371 dimitri 1.18 C-- Read initial sea ice salinity from file if available.
372     IF ( HsaltFile .NE. ' ' ) THEN
373 mlosch 1.43 CALL READ_FLD_XY_RL( HsaltFile, ' ', HSALT, 0, myThid )
374     _EXCH_XY_RL(HSALT,myThid)
375 dimitri 1.18 ENDIF
376     #endif /* SEAICE_SALINITY */
377 jmc 1.23
378 dimitri 1.33 #ifdef SEAICE_AGE
379     C-- Read initial sea ice age from file if available.
380     IF ( IceAgeFile .NE. ' ' ) THEN
381 mlosch 1.43 CALL READ_FLD_XY_RL( IceAgeFile, ' ', ICEAGE, 0, myThid )
382     _EXCH_XY_RL(ICEAGE,myThid)
383 dimitri 1.33 ENDIF
384     #endif /* SEAICE_AGE */
385    
386 heimbach 1.1 ENDIF
387    
388     C--- Complete initialization
389 mlosch 1.8 PSTAR = SEAICE_strength
390 heimbach 1.1 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 mlosch 1.43 ZETA(i,j,bi,bj)=HEFF(i,j,bi,bj)*(1.0 _d 11)
395 jmc 1.23 ETA(i,j,bi,bj)=ZETA(i,j,bi,bj)/4.0 _d 0
396 mlosch 1.43 PRESS0(i,j,bi,bj)=PSTAR*HEFF(i,j,bi,bj)
397     & *EXP(-20.0 _d 0*(ONE-AREA(i,j,bi,bj)))
398 jmc 1.23 ZMAX(i,j,bi,bj)=(5.0 _d +12/(2.0 _d +04))*PRESS0(i,j,bi,bj)
399     ZMIN(i,j,bi,bj)=SEAICE_zetaMin
400     PRESS0(i,j,bi,bj)=PRESS0(i,j,bi,bj)*HEFFM(i,j,bi,bj)
401 heimbach 1.1 ENDDO
402     ENDDO
403     IF ( useRealFreshWaterFlux .AND. .NOT.useThSIce ) THEN
404     DO j=1-OLy,sNy+OLy
405     DO i=1-OLx,sNx+OLx
406 mlosch 1.43 sIceLoad(i,j,bi,bj) = HEFF(i,j,bi,bj)*SEAICE_rhoIce
407 mlosch 1.36 & + HSNOW(i,j,bi,bj)*SEAICE_rhoSnow
408 jmc 1.23
409 heimbach 1.1 ENDDO
410     ENDDO
411     ENDIF
412     ENDDO
413     ENDDO
414    
415     RETURN
416     END

  ViewVC Help
Powered by ViewVC 1.1.22