/[MITgcm]/MITgcm/pkg/seaice/seaice_init_varia.F
ViewVC logotype

Contents of /MITgcm/pkg/seaice/seaice_init_varia.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.33 - (show annotations) (download)
Wed Dec 17 03:33:30 2008 UTC (15 years, 5 months ago) by dimitri
Branch: MAIN
CVS Tags: checkpoint61g, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i
Changes since 1.32: +21 -1 lines
added sea ice age tracer;  #define SEAICE_AGE in SEAICE_PARAMS.h

1 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_init_varia.F,v 1.32 2008/10/11 20:37:05 heimbach Exp $
2 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 "DYNVARS.h"
21 #include "SEAICE.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
29 #ifdef ALLOW_OBCS
30 # include "OBCS_OPTIONS.h"
31 # include "OBCS.h"
32 #endif
33
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 _RL PSTAR
44 _RS mask_uice
45 INTEGER myIter, myTile
46
47 #ifdef ALLOW_OBCS
48 INTEGER I_obc, J_obc, kSurface
49 IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
50 kSurface = Nr
51 ELSE
52 kSurface = 1
53 ENDIF
54 #endif /* ALLOW_OBCS */
55
56 C-- Initialise all variables in common blocks:
57 DO bj=myByLo(myThid),myByHi(myThid)
58 DO bi=myBxLo(myThid),myBxHi(myThid)
59 DO k=1,3
60 DO j=1-OLy,sNy+OLy
61 DO i=1-OLx,sNx+OLx
62 HEFF(i,j,k,bi,bj)=0. _d 0
63 AREA(i,j,k,bi,bj)=0. _d 0
64 UICE(i,j,k,bi,bj)=0. _d 0
65 VICE(i,j,k,bi,bj)=0. _d 0
66 ENDDO
67 ENDDO
68 ENDDO
69 #ifdef SEAICE_MULTICATEGORY
70 DO k=1,MULTDIM
71 DO j=1-OLy,sNy+OLy
72 DO i=1-OLx,sNx+OLx
73 TICES(i,j,k,bi,bj)=0. _d 0
74 ENDDO
75 ENDDO
76 ENDDO
77 #endif
78 DO j=1-OLy,sNy+OLy
79 DO i=1-OLx,sNx+OLx
80 ETA (i,j,bi,bj) = 0. _d 0
81 ZETA(i,j,bi,bj) = 0. _d 0
82 DRAGS(i,j,bi,bj) = 0. _d 0
83 DRAGA(i,j,bi,bj) = 0. _d 0
84 FORCEX(i,j,bi,bj) = 0. _d 0
85 FORCEY(i,j,bi,bj) = 0. _d 0
86 UICEC(i,j,bi,bj) = 0. _d 0
87 VICEC(i,j,bi,bj) = 0. _d 0
88 #ifdef SEAICE_CGRID
89 seaiceMassC(i,j,bi,bj)=0. _d 0
90 seaiceMassU(i,j,bi,bj)=0. _d 0
91 seaiceMassV(i,j,bi,bj)=0. _d 0
92 seaiceMaskU(i,j,bi,bj)=0. _d 0
93 seaiceMaskV(i,j,bi,bj)=0. _d 0
94 # ifdef SEAICE_ALLOW_EVP
95 stressDivergenceX(i,j,bi,bj) = 0. _d 0
96 stressDivergenceY(i,j,bi,bj) = 0. _d 0
97 seaice_sigma1 (i,j,bi,bj) = 0. _d 0
98 seaice_sigma2 (i,j,bi,bj) = 0. _d 0
99 seaice_sigma12(i,j,bi,bj) = 0. _d 0
100 # endif /* SEAICE_ALLOW_EVP */
101 #else /* SEAICE_CGRID */
102 AMASS(i,j,bi,bj) = 0. _d 0
103 DAIRN(i,j,bi,bj) = 0. _d 0
104 UVM(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 GWATX(i,j,bi,bj) = 0. _d 0
108 GWATY(i,j,bi,bj) = 0. _d 0
109 KGEO(i,j,bi,bj) = 0
110 #endif /* SEAICE_CGRID */
111 DWATN(i,j,bi,bj) = 0. _d 0
112 PRESS0(i,j,bi,bj) = 0. _d 0
113 FORCEX0(i,j,bi,bj)= 0. _d 0
114 FORCEY0(i,j,bi,bj)= 0. _d 0
115 ZMAX(i,j,bi,bj) = 0. _d 0
116 ZMIN(i,j,bi,bj) = 0. _d 0
117 HSNOW(i,j,bi,bj) = 0. _d 0
118 #ifdef SEAICE_SALINITY
119 HSALT(i,j,bi,bj) = 0. _d 0
120 #endif
121 #ifdef SEAICE_AGE
122 ICEAGE(i,j,bi,bj) = 0. _d 0
123 #endif
124 HEFFM(i,j,bi,bj) = 0. _d 0
125 YNEG (i,j,bi,bj) = 0. _d 0
126 RIVER(i,j,bi,bj) = 0. _d 0
127 TMIX(i,j,bi,bj) = 0. _d 0
128 TICE(i,j,bi,bj) = 0. _d 0
129 TAUX(i,j,bi,bj) = 0. _d 0
130 TAUY(i,j,bi,bj) = 0. _d 0
131 #ifdef ALLOW_SEAICE_COST_EXPORT
132 uHeffExportCell(i,j,bi,bj) = 0. _d 0
133 vHeffExportCell(i,j,bi,bj) = 0. _d 0
134 #endif
135 saltWtrIce(i,j,bi,bj) = 0. _d 0
136 frWtrIce(i,j,bi,bj) = 0. _d 0
137 frWtrAtm(i,j,bi,bj) = 0. _d 0
138 ENDDO
139 ENDDO
140 ENDDO
141 ENDDO
142
143 C-- Initialize grid info
144 DO bj=myByLo(myThid),myByHi(myThid)
145 DO bi=myBxLo(myThid),myBxHi(myThid)
146 DO j=1-OLy,sNy+OLy
147 DO i=1-OLx,sNx+OLx
148 HEFFM(i,j,bi,bj)=ONE
149 IF (_hFacC(i,j,1,bi,bj).eq.0.) HEFFM(i,j,bi,bj)=0. _d 0
150 ENDDO
151 ENDDO
152 DO j=1-OLy+1,sNy+OLy
153 DO i=1-OLx+1,sNx+OLx
154 #ifdef SEAICE_CGRID
155 seaiceMaskU(i,j,bi,bj)= 0.0 _d 0
156 seaiceMaskV(i,j,bi,bj)= 0.0 _d 0
157 mask_uice=HEFFM(i,j,bi,bj)+HEFFM(i-1,j ,bi,bj)
158 IF(mask_uice.GT.1.5) seaiceMaskU(i,j,bi,bj)=ONE
159 mask_uice=HEFFM(i,j,bi,bj)+HEFFM(i ,j-1,bi,bj)
160 IF(mask_uice.GT.1.5) seaiceMaskV(i,j,bi,bj)=ONE
161 #else
162 UVM(i,j,bi,bj)=0. _d 0
163 mask_uice=HEFFM(i,j, bi,bj)+HEFFM(i-1,j-1,bi,bj)
164 & +HEFFM(i,j-1,bi,bj)+HEFFM(i-1,j, bi,bj)
165 IF(mask_uice.GT.3.5) UVM(i,j,bi,bj)=ONE
166 #endif /* SEAICE_CGRID */
167 ENDDO
168 ENDDO
169
170 #ifdef ALLOW_OBCS
171 IF (useOBCS) THEN
172 C-- If OBCS is turned on, close southern and western boundaries
173 #ifdef ALLOW_OBCS_SOUTH
174 DO i=1-Olx,sNx+Olx
175 C Southern boundary
176 J_obc = OB_Js(i,bi,bj)
177 IF (J_obc.NE.0) THEN
178 #ifdef SEAICE_CGRID
179 seaiceMaskU(i,J_obc,bi,bj)= 0.0 _d 0
180 seaiceMaskV(i,J_obc,bi,bj)= 0.0 _d 0
181 #else
182 UVM(i,J_obc,bi,bj)=0. _d 0
183 #endif /* SEAICE_CGRID */
184 ENDIF
185 ENDDO
186 #endif /* ALLOW_OBCS_SOUTH */
187 #ifdef ALLOW_OBCS_WEST
188 DO j=1-Oly,sNy+Oly
189 C Western boundary
190 I_obc=OB_Iw(j,bi,bj)
191 IF (I_obc.NE.0) THEN
192 #ifdef SEAICE_CGRID
193 seaiceMaskU(I_obc,j,bi,bj)= 0.0 _d 0
194 seaiceMaskV(I_obc,j,bi,bj)= 0.0 _d 0
195 #else
196 UVM(I_obc,j,bi,bj)=0. _d 0
197 #endif /* SEAICE_CGRID */
198 ENDIF
199 ENDDO
200 #endif /* ALLOW_OBCS_WEST */
201 ENDIF
202 #endif /* ALLOW_OBCS */
203
204 #ifdef ALLOW_EXCH2
205 #ifndef SEAICE_CGRID
206 C-- Special stuff for cubed sphere: assume grid is rectangular and
207 C set UV mask to zero except for Arctic and Antarctic cube faces.
208 IF (useCubedSphereExchange) THEN
209 myTile = W2_myTileList(bi)
210 IF ( exch2_myFace(myTile) .EQ. 1 .OR.
211 & exch2_myFace(myTile) .EQ. 2 .OR.
212 & exch2_myFace(myTile) .EQ. 4 .OR.
213 & exch2_myFace(myTile) .EQ. 5 ) THEN
214 DO j=1-OLy,sNy+OLy
215 DO i=1-OLx,sNx+OLx
216 UVM(i,j,bi,bj)=0. _d 0
217 ENDDO
218 ENDDO
219 ELSEIF ( exch2_isWedge(myTile) .EQ. 1 ) THEN
220 i=1
221 DO j=1-OLy,sNy+OLy
222 UVM(i,j,bi,bj)=0. _d 0
223 ENDDO
224 ELSEIF ( exch2_isSedge(myTile) .EQ. 1 ) THEN
225 j=1
226 DO i=1-OLx,sNx+OLx
227 UVM(i,j,bi,bj)=0. _d 0
228 ENDDO
229 ENDIF
230 ENDIF
231 #endif /* SEAICE_CGRID */
232 #endif /* ALLOW_EXCH2 */
233
234 DO j=1-OLy,sNy+OLy
235 DO i=1-OLx,sNx+OLx
236 TICE(i,j,bi,bj)=273.0 _d 0
237 #ifdef SEAICE_MULTICATEGORY
238 DO k=1,MULTDIM
239 TICES(i,j,k,bi,bj)=273.0 _d 0
240 ENDDO
241 #endif /* SEAICE_MULTICATEGORY */
242 #ifndef SEAICE_CGRID
243 AMASS (i,j,bi,bj)=1000.0 _d 0
244 #else
245 seaiceMassC(i,j,bi,bj)=1000.0 _d 0
246 seaiceMassU(i,j,bi,bj)=1000.0 _d 0
247 seaiceMassV(i,j,bi,bj)=1000.0 _d 0
248 #endif
249 ENDDO
250 ENDDO
251
252 #ifndef SEAICE_CGRID
253 C-- Choose a proxy level for geostrophic velocity,
254 DO j=1-OLy,sNy+OLy
255 DO i=1-OLx,sNx+OLx
256 #ifdef SEAICE_TEST_ICE_STRESS_1
257 KGEO(i,j,bi,bj) = 1
258 #else /* SEAICE_TEST_ICE_STRESS_1 */
259 IF (klowc(i,j,bi,bj) .LT. 2) THEN
260 KGEO(i,j,bi,bj) = 1
261 ELSE
262 KGEO(i,j,bi,bj) = 2
263 DO WHILE ( abs(rC(KGEO(i,j,bi,bj))) .LT. 50.0 .AND.
264 & KGEO(i,j,bi,bj) .LT. (klowc(i,j,bi,bj)-1) )
265 KGEO(i,j,bi,bj) = KGEO(i,j,bi,bj) + 1
266 ENDDO
267 ENDIF
268 #endif /* SEAICE_TEST_ICE_STRESS_1 */
269 ENDDO
270 ENDDO
271 #endif /* SEAICE_CGRID */
272
273 ENDDO
274 ENDDO
275
276 C-- Update overlap regions
277 #ifdef SEAICE_CGRID
278 CALL EXCH_UV_XY_RL(seaiceMaskU,seaiceMaskV,.FALSE.,myThid)
279 #else
280 _EXCH_XY_R8(UVM, myThid)
281 #endif
282
283 C-- Now lets look at all these beasts
284 IF ( debugLevel .GE. debLevB ) THEN
285 myIter=0
286 CALL PLOT_FIELD_XYRL( HEFFM , 'Current HEFFM ' ,
287 & myIter, myThid )
288 #ifdef SEAICE_CGRID
289 CALL PLOT_FIELD_XYRL( seaiceMaskU, 'Current seaiceMaskU',
290 & myIter, myThid )
291 CALL PLOT_FIELD_XYRL( seaiceMaskV, 'Current seaiceMaskV',
292 & myIter, myThid )
293 #else
294 CALL PLOT_FIELD_XYRL( UVM , 'Current UVM ' ,
295 & myIter, myThid )
296 #endif
297 ENDIF
298
299 C-- Set model variables to initial/restart conditions
300 IF ( .NOT. ( startTime .EQ. baseTime .AND. nIter0 .EQ. 0
301 & .AND. pickupSuff .EQ. ' ') ) THEN
302
303 CALL SEAICE_READ_PICKUP ( myThid )
304
305 ELSE
306
307 DO bj=myByLo(myThid),myByHi(myThid)
308 DO bi=myBxLo(myThid),myBxHi(myThid)
309 DO j=1-OLy,sNy+OLy
310 DO i=1-OLx,sNx+OLx
311 TMIX(i,j,bi,bj)=TICE(i,j,bi,bj)
312 DO k=1,3
313 HEFF(i,j,k,bi,bj)=SEAICE_initialHEFF*HEFFM(i,j,bi,bj)
314 UICE(i,j,k,bi,bj)=ZERO
315 VICE(i,j,k,bi,bj)=ZERO
316 ENDDO
317 ENDDO
318 ENDDO
319 ENDDO
320 ENDDO
321
322 C-- Read initial sea-ice thickness from file if available.
323 IF ( HeffFile .NE. ' ' ) THEN
324 CALL READ_FLD_XY_RL( HeffFile, ' ', ZETA, 0, myThid )
325 _EXCH_XY_R8(ZETA,myThid)
326 DO bj=myByLo(myThid),myByHi(myThid)
327 DO bi=myBxLo(myThid),myBxHi(myThid)
328 DO j=1-OLy,sNy+OLy
329 DO i=1-OLx,sNx+OLx
330 DO k=1,3
331 HEFF(i,j,k,bi,bj) = MAX(ZETA(i,j,bi,bj),ZERO)
332 ENDDO
333 ENDDO
334 ENDDO
335 ENDDO
336 ENDDO
337 ENDIF
338
339 DO bj=myByLo(myThid),myByHi(myThid)
340 DO bi=myBxLo(myThid),myBxHi(myThid)
341 DO j=1-OLy,sNy+OLy
342 DO i=1-OLx,sNx+OLx
343 DO k=1,3
344 IF(HEFF(i,j,k,bi,bj).GT.ZERO)
345 & AREA(i,j,k,bi,bj)=ONE
346 ENDDO
347 ENDDO
348 ENDDO
349 ENDDO
350 ENDDO
351
352 C-- Read initial sea-ice area from file if available.
353 IF ( AreaFile .NE. ' ' ) THEN
354 CALL READ_FLD_XY_RL( AreaFile, ' ', ZETA, 0, 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 HSNOW(i,j,bi,bj)=0.2*AREA(i,j,1,bi,bj)
379 ENDDO
380 ENDDO
381 ENDDO
382 ENDDO
383
384 C-- Read initial snow thickness from file if available.
385 IF ( HsnowFile .NE. ' ' ) THEN
386 CALL READ_FLD_XY_RL( HsnowFile, ' ', ZETA, 0, 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 HSNOW(i,j,bi,bj) = MAX(ZETA(i,j,bi,bj),ZERO)
393 ENDDO
394 ENDDO
395 ENDDO
396 ENDDO
397 ENDIF
398
399 #ifdef SEAICE_SALINITY
400 DO bj=myByLo(myThid),myByHi(myThid)
401 DO bi=myBxLo(myThid),myBxHi(myThid)
402 DO j=1-OLy,sNy+OLy
403 DO i=1-OLx,sNx+OLx
404 HSALT(i,j,bi,bj)=HEFF(i,j,1,bi,bj)*salt(i,j,1,bi,bj)*
405 & ICE2WATR*rhoConstFresh*SEAICE_salinity
406 ENDDO
407 ENDDO
408 ENDDO
409 ENDDO
410
411 C-- Read initial sea ice salinity from file if available.
412 IF ( HsaltFile .NE. ' ' ) THEN
413 CALL READ_FLD_XY_RL( HsaltFile, ' ', ZETA, 0, myThid )
414 _EXCH_XY_R8(ZETA,myThid)
415 DO bj=myByLo(myThid),myByHi(myThid)
416 DO bi=myBxLo(myThid),myBxHi(myThid)
417 DO j=1-OLy,sNy+OLy
418 DO i=1-OLx,sNx+OLx
419 HSALT(i,j,bi,bj) = ZETA(i,j,bi,bj)
420 ENDDO
421 ENDDO
422 ENDDO
423 ENDDO
424 ENDIF
425 #endif /* SEAICE_SALINITY */
426
427 #ifdef SEAICE_AGE
428 C-- Read initial sea ice age from file if available.
429 IF ( IceAgeFile .NE. ' ' ) THEN
430 CALL READ_FLD_XY_RL( IceAgeFile, ' ', ZETA, 0, myThid )
431 _EXCH_XY_R8(ZETA,myThid)
432 DO bj=myByLo(myThid),myByHi(myThid)
433 DO bi=myBxLo(myThid),myBxHi(myThid)
434 DO j=1-OLy,sNy+OLy
435 DO i=1-OLx,sNx+OLx
436 ICEAGE(i,j,bi,bj) = ZETA(i,j,bi,bj)
437 ENDDO
438 ENDDO
439 ENDDO
440 ENDDO
441 ENDIF
442 #endif /* SEAICE_AGE */
443
444 ENDIF
445
446 C--- Complete initialization
447 PSTAR = SEAICE_strength
448 DO bj=myByLo(myThid),myByHi(myThid)
449 DO bi=myBxLo(myThid),myBxHi(myThid)
450 DO j=1-OLy,sNy+OLy
451 DO i=1-OLx,sNx+OLx
452 ZETA(i,j,bi,bj)=HEFF(i,j,1,bi,bj)*(1.0 _d 11)
453 ETA(i,j,bi,bj)=ZETA(i,j,bi,bj)/4.0 _d 0
454 PRESS0(i,j,bi,bj)=PSTAR*HEFF(i,j,1,bi,bj)
455 & *EXP(-20.0 _d 0*(ONE-AREA(i,j,1,bi,bj)))
456 ZMAX(i,j,bi,bj)=(5.0 _d +12/(2.0 _d +04))*PRESS0(i,j,bi,bj)
457 ZMIN(i,j,bi,bj)=SEAICE_zetaMin
458 PRESS0(i,j,bi,bj)=PRESS0(i,j,bi,bj)*HEFFM(i,j,bi,bj)
459 ENDDO
460 ENDDO
461 IF ( useRealFreshWaterFlux .AND. .NOT.useThSIce ) THEN
462 DO j=1-OLy,sNy+OLy
463 DO i=1-OLx,sNx+OLx
464 sIceLoad(i,j,bi,bj) = HEFF(i,j,1,bi,bj)*SEAICE_rhoIce
465 & + HSNOW(i,j,bi,bj)* 330. _d 0
466
467 ENDDO
468 ENDDO
469 ENDIF
470 ENDDO
471 ENDDO
472
473 RETURN
474 END

  ViewVC Help
Powered by ViewVC 1.1.22