/[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.18 - (show annotations) (download)
Sat Sep 15 01:38:59 2007 UTC (16 years, 8 months ago) by dimitri
Branch: MAIN
CVS Tags: checkpoint59g
Changes since 1.17: +23 -1 lines
adding ice salinity HSALT as a prognostic variable

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

  ViewVC Help
Powered by ViewVC 1.1.22