/[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.24 - (show annotations) (download)
Wed Dec 5 07:28:29 2007 UTC (16 years, 5 months ago) by dimitri
Branch: MAIN
CVS Tags: checkpoint59m, checkpoint59l
Changes since 1.23: +7 -8 lines
o pkg/seaice: removed SEAICE_FFIELDS.h and seaice_get_forcing.F
  seaice forcing fields can now be read only through pkg/exf

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

  ViewVC Help
Powered by ViewVC 1.1.22