/[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.72 - (show annotations) (download)
Wed Mar 14 22:55:53 2012 UTC (12 years, 2 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63k
Changes since 1.71: +15 -2 lines
Merging JPLs
o seaice AREA relaxation code
o modified global mean cost imbalance code

1 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_init_varia.F,v 1.71 2012/03/11 13:43:46 jmc Exp $
2 C $Name: $
3
4 #include "SEAICE_OPTIONS.h"
5 #ifdef ALLOW_OBCS
6 # include "OBCS_OPTIONS.h"
7 #endif
8
9 CStartOfInterface
10 SUBROUTINE SEAICE_INIT_VARIA( myThid )
11 C *==========================================================*
12 C | SUBROUTINE SEAICE_INIT_VARIA |
13 C | o Initialization of sea ice model. |
14 C *==========================================================*
15 C *==========================================================*
16 IMPLICIT NONE
17
18 C === Global variables ===
19 #include "SIZE.h"
20 #include "EEPARAMS.h"
21 #include "PARAMS.h"
22 #include "GRID.h"
23 #include "DYNVARS.h"
24 #include "FFIELDS.h"
25 #include "SEAICE_SIZE.h"
26 #include "SEAICE_PARAMS.h"
27 #include "SEAICE.h"
28 #include "SEAICE_TRACER.h"
29 #include "SEAICE_TAVE.h"
30 #ifdef OBCS_UVICE_OLD
31 # include "OBCS_GRID.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, bi, bj
43 _RL PSTAR
44 INTEGER kSurface
45 #ifdef SEAICE_CGRID
46 _RS mask_uice
47 #endif
48 INTEGER k
49 #ifdef ALLOW_SITRACER
50 INTEGER iTr, jTh
51 #endif
52 #ifdef OBCS_UVICE_OLD
53 INTEGER I_obc, J_obc
54 #endif /* ALLOW_OBCS */
55
56 IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
57 kSurface = Nr
58 ELSE
59 kSurface = 1
60 ENDIF
61
62 C-- Initialise all variables in common blocks:
63 DO bj=myByLo(myThid),myByHi(myThid)
64 DO bi=myBxLo(myThid),myBxHi(myThid)
65 DO j=1-OLy,sNy+OLy
66 DO i=1-OLx,sNx+OLx
67 HEFF(i,j,bi,bj)=0. _d 0
68 AREA(i,j,bi,bj)=0. _d 0
69 UICE(i,j,bi,bj)=0. _d 0
70 VICE(i,j,bi,bj)=0. _d 0
71 #ifdef SEAICE_ALLOW_FREEDRIFT
72 uice_fd(i,j,bi,bj)=0. _d 0
73 vice_fd(i,j,bi,bj)=0. _d 0
74 #endif
75 C
76 uIceNm1(i,j,bi,bj)=0. _d 0
77 vIceNm1(i,j,bi,bj)=0. _d 0
78 #ifdef SEAICE_GROWTH_LEGACY
79 areaNm1(i,j,bi,bj)=0. _d 0
80 hEffNm1(i,j,bi,bj)=0. _d 0
81 #endif
82 ETA (i,j,bi,bj) = 0. _d 0
83 ZETA(i,j,bi,bj) = 0. _d 0
84 DRAGS(i,j,bi,bj) = 0. _d 0
85 DRAGA(i,j,bi,bj) = 0. _d 0
86 FORCEX(i,j,bi,bj) = 0. _d 0
87 FORCEY(i,j,bi,bj) = 0. _d 0
88 uIceC(i,j,bi,bj) = 0. _d 0
89 vIceC(i,j,bi,bj) = 0. _d 0
90 #ifdef SEAICE_CGRID
91 seaiceMassC(i,j,bi,bj)=0. _d 0
92 seaiceMassU(i,j,bi,bj)=0. _d 0
93 seaiceMassV(i,j,bi,bj)=0. _d 0
94 stressDivergenceX(i,j,bi,bj) = 0. _d 0
95 stressDivergenceY(i,j,bi,bj) = 0. _d 0
96 # ifdef SEAICE_ALLOW_EVP
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 WINDX(i,j,bi,bj) = 0. _d 0
105 WINDY(i,j,bi,bj) = 0. _d 0
106 GWATX(i,j,bi,bj) = 0. _d 0
107 GWATY(i,j,bi,bj) = 0. _d 0
108 #endif /* SEAICE_CGRID */
109 DWATN(i,j,bi,bj) = 0. _d 0
110 PRESS0(i,j,bi,bj) = 0. _d 0
111 FORCEX0(i,j,bi,bj)= 0. _d 0
112 FORCEY0(i,j,bi,bj)= 0. _d 0
113 ZMAX(i,j,bi,bj) = 0. _d 0
114 ZMIN(i,j,bi,bj) = 0. _d 0
115 HSNOW(i,j,bi,bj) = 0. _d 0
116 #ifdef SEAICE_VARIABLE_SALINITY
117 HSALT(i,j,bi,bj) = 0. _d 0
118 #endif
119 #ifdef ALLOW_SITRACER
120 DO iTr = 1, SItrMaxNum
121 SItracer(i,j,bi,bj,iTr) = 0. _d 0
122 SItrBucket(i,j,bi,bj,iTr) = 0. _d 0
123 c "ice concentration" tracer that should remain .EQ.1.
124 if (SItrName(iTr).EQ.'one') SItracer(i,j,bi,bj,iTr)=1. _d 0
125 ENDDO
126 DO jTh = 1, 5
127 SItrHEFF (i,j,bi,bj,jTh) = 0. _d 0
128 ENDDO
129 DO jTh = 1, 3
130 SItrAREA (i,j,bi,bj,jTh) = 0. _d 0
131 ENDDO
132 #endif
133 TICE(i,j,bi,bj) = 0. _d 0
134 DO k=1,MULTDIM
135 TICES(i,j,k,bi,bj)=0. _d 0
136 ENDDO
137 TAUX(i,j,bi,bj) = 0. _d 0
138 TAUY(i,j,bi,bj) = 0. _d 0
139 #ifdef ALLOW_SEAICE_COST_EXPORT
140 uHeffExportCell(i,j,bi,bj) = 0. _d 0
141 vHeffExportCell(i,j,bi,bj) = 0. _d 0
142 #endif
143 saltWtrIce(i,j,bi,bj) = 0. _d 0
144 frWtrIce(i,j,bi,bj) = 0. _d 0
145 #if (defined (ALLOW_MEAN_SFLUX_COST_CONTRIBUTION) || defined (ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION))
146 frWtrAtm(i,j,bi,bj) = 0. _d 0
147 AREAforAtmFW(i,j,bi,bj)=0. _d 0
148 #endif
149 ENDDO
150 ENDDO
151 ENDDO
152 ENDDO
153
154 #ifdef ALLOW_TIMEAVE
155 C Initialize averages to zero
156 DO bj = myByLo(myThid), myByHi(myThid)
157 DO bi = myBxLo(myThid), myBxHi(myThid)
158 CALL TIMEAVE_RESET( FUtave , 1, bi, bj, myThid )
159 CALL TIMEAVE_RESET( FVtave , 1, bi, bj, myThid )
160 CALL TIMEAVE_RESET( EmPmRtave, 1, bi, bj, myThid )
161 CALL TIMEAVE_RESET( QNETtave , 1, bi, bj, myThid )
162 CALL TIMEAVE_RESET( QSWtave , 1, bi, bj, myThid )
163 CALL TIMEAVE_RESET( UICEtave , 1, bi, bj, myThid )
164 CALL TIMEAVE_RESET( VICEtave , 1, bi, bj, myThid )
165 CALL TIMEAVE_RESET( HEFFtave , 1, bi, bj, myThid )
166 CALL TIMEAVE_RESET( AREAtave , 1, bi, bj, myThid )
167 SEAICE_timeAve(bi,bj) = ZERO
168 ENDDO
169 ENDDO
170 #endif /* ALLOW_TIMEAVE */
171
172 C-- Initialize (variable) grid info. As long as we allow masking of
173 C-- velocities outside of ice covered areas (in seaice_dynsolver)
174 C-- we need to re-initialize seaiceMaskU/V here for TAF/TAMC
175 #ifdef SEAICE_CGRID
176 DO bj=myByLo(myThid),myByHi(myThid)
177 DO bi=myBxLo(myThid),myBxHi(myThid)
178 DO j=1-OLy+1,sNy+OLy
179 DO i=1-OLx+1,sNx+OLx
180 seaiceMaskU(i,j,bi,bj)= 0.0 _d 0
181 seaiceMaskV(i,j,bi,bj)= 0.0 _d 0
182 mask_uice=HEFFM(i,j,bi,bj)+HEFFM(i-1,j ,bi,bj)
183 IF(mask_uice.GT.1.5 _d 0) seaiceMaskU(i,j,bi,bj)=1.0 _d 0
184 mask_uice=HEFFM(i,j,bi,bj)+HEFFM(i ,j-1,bi,bj)
185 IF(mask_uice.GT.1.5 _d 0) seaiceMaskV(i,j,bi,bj)=1.0 _d 0
186 ENDDO
187 ENDDO
188 ENDDO
189 ENDDO
190 #endif /* SEAICE_CGRID */
191
192 DO bj=myByLo(myThid),myByHi(myThid)
193 DO bi=myBxLo(myThid),myBxHi(myThid)
194 #ifdef OBCS_UVICE_OLD
195 #ifdef SEAICE_CGRID
196 IF (useOBCS) THEN
197 C-- If OBCS is turned on, close southern and western boundaries
198 DO i=1-OLx,sNx+OLx
199 C Southern boundary
200 J_obc = OB_Js(i,bi,bj)
201 IF (J_obc.NE.0) THEN
202 seaiceMaskU(i,J_obc,bi,bj)= 0.0 _d 0
203 seaiceMaskV(i,J_obc,bi,bj)= 0.0 _d 0
204 ENDIF
205 ENDDO
206 DO j=1-OLy,sNy+OLy
207 C Western boundary
208 I_obc=OB_Iw(j,bi,bj)
209 IF (I_obc.NE.0) THEN
210 seaiceMaskU(I_obc,j,bi,bj)= 0.0 _d 0
211 seaiceMaskV(I_obc,j,bi,bj)= 0.0 _d 0
212 ENDIF
213 ENDDO
214 ENDIF
215 #endif /* SEAICE_CGRID */
216 #endif /* OBCS_UVICE_OLD */
217
218 DO j=1-OLy,sNy+OLy
219 DO i=1-OLx,sNx+OLx
220 TICE(i,j,bi,bj)=273.0 _d 0
221 DO k=1,MULTDIM
222 TICES(i,j,k,bi,bj)=273.0 _d 0
223 ENDDO
224 #ifndef SEAICE_CGRID
225 AMASS (i,j,bi,bj)=1000.0 _d 0
226 #else
227 seaiceMassC(i,j,bi,bj)=1000.0 _d 0
228 seaiceMassU(i,j,bi,bj)=1000.0 _d 0
229 seaiceMassV(i,j,bi,bj)=1000.0 _d 0
230 #endif
231 ENDDO
232 ENDDO
233
234 ENDDO
235 ENDDO
236
237 C-- Update overlap regions
238 #ifdef SEAICE_CGRID
239 CALL EXCH_UV_XY_RL(seaiceMaskU,seaiceMaskV,.FALSE.,myThid)
240 #else
241 _EXCH_XY_RS(UVM, myThid)
242 #endif
243
244 C-- Now lets look at all these beasts
245 IF ( debugLevel .GE. debLevC ) THEN
246 CALL PLOT_FIELD_XYRL( HEFFM , 'Current HEFFM ' ,
247 & nIter0, myThid )
248 #ifdef SEAICE_CGRID
249 CALL PLOT_FIELD_XYRL( seaiceMaskU, 'Current seaiceMaskU',
250 & nIter0, myThid )
251 CALL PLOT_FIELD_XYRL( seaiceMaskV, 'Current seaiceMaskV',
252 & nIter0, myThid )
253 #else
254 CALL PLOT_FIELD_XYRS( UVM , 'Current UVM ' ,
255 & nIter0, myThid )
256 #endif
257 ENDIF
258
259 C-- Set model variables to initial/restart conditions
260 IF ( .NOT. ( startTime .EQ. baseTime .AND. nIter0 .EQ. 0
261 & .AND. pickupSuff .EQ. ' ') ) THEN
262
263 CALL SEAICE_READ_PICKUP ( myThid )
264
265 ELSE
266
267 DO bj=myByLo(myThid),myByHi(myThid)
268 DO bi=myBxLo(myThid),myBxHi(myThid)
269 DO j=1-OLy,sNy+OLy
270 DO i=1-OLx,sNx+OLx
271 HEFF(i,j,bi,bj)=SEAICE_initialHEFF*HEFFM(i,j,bi,bj)
272 UICE(i,j,bi,bj)=ZERO
273 VICE(i,j,bi,bj)=ZERO
274 ENDDO
275 ENDDO
276 ENDDO
277 ENDDO
278
279 C-- Read initial sea-ice velocity from file (if available)
280 IF ( uIceFile .NE. ' ' )
281 & CALL READ_FLD_XY_RL( uIceFile, ' ', uIce, 0, myThid )
282 IF ( vIceFile .NE. ' ' )
283 & CALL READ_FLD_XY_RL( vIceFile, ' ', vIce, 0, myThid )
284 IF ( uIceFile .NE. ' ' .OR. vIceFile .NE. ' ' ) THEN
285 #ifdef SEAICE_CGRID
286 DO bj=myByLo(myThid),myByHi(myThid)
287 DO bi=myBxLo(myThid),myBxHi(myThid)
288 DO j=1-OLy,sNy+OLy
289 DO i=1-OLx,sNx+OLx
290 uIce(i,j,bi,bj) = uIce(i,j,bi,bj)*seaiceMaskU(i,j,bi,bj)
291 vIce(i,j,bi,bj) = vIce(i,j,bi,bj)*seaiceMaskV(i,j,bi,bj)
292 ENDDO
293 ENDDO
294 ENDDO
295 ENDDO
296 #endif /* SEAICE_CGRID */
297 CALL EXCH_UV_XY_RL( uIce, vIce, .TRUE., myThid )
298 ENDIF
299
300 C-- Read initial sea-ice thickness from file if available.
301 IF ( HeffFile .NE. ' ' ) THEN
302 CALL READ_FLD_XY_RL( HeffFile, ' ', HEFF, 0, myThid )
303 _EXCH_XY_RL(HEFF,myThid)
304 DO bj=myByLo(myThid),myByHi(myThid)
305 DO bi=myBxLo(myThid),myBxHi(myThid)
306 DO j=1-OLy,sNy+OLy
307 DO i=1-OLx,sNx+OLx
308 HEFF(i,j,bi,bj) = MAX(HEFF(i,j,bi,bj),ZERO)
309 ENDDO
310 ENDDO
311 ENDDO
312 ENDDO
313 ENDIF
314
315 DO bj=myByLo(myThid),myByHi(myThid)
316 DO bi=myBxLo(myThid),myBxHi(myThid)
317 DO j=1-OLy,sNy+OLy
318 DO i=1-OLx,sNx+OLx
319 IF(HEFF(i,j,bi,bj).GT.ZERO) AREA(i,j,bi,bj)=ONE
320 ENDDO
321 ENDDO
322 ENDDO
323 ENDDO
324
325 C-- Read initial sea-ice area from file if available.
326 IF ( AreaFile .NE. ' ' ) THEN
327 CALL READ_FLD_XY_RL( AreaFile, ' ', AREA, 0, myThid )
328 _EXCH_XY_RL(AREA,myThid)
329 DO bj=myByLo(myThid),myByHi(myThid)
330 DO bi=myBxLo(myThid),myBxHi(myThid)
331 DO j=1-OLy,sNy+OLy
332 DO i=1-OLx,sNx+OLx
333 AREA(i,j,bi,bj) = MAX(AREA(i,j,bi,bj),ZERO)
334 AREA(i,j,bi,bj) = MIN(AREA(i,j,bi,bj),ONE)
335 IF ( AREA(i,j,bi,bj) .LE. ZERO ) HEFF(i,j,bi,bj) = ZERO
336 IF ( HEFF(i,j,bi,bj) .LE. ZERO ) AREA(i,j,bi,bj) = ZERO
337 ENDDO
338 ENDDO
339 ENDDO
340 ENDDO
341 ENDIF
342
343 DO bj=myByLo(myThid),myByHi(myThid)
344 DO bi=myBxLo(myThid),myBxHi(myThid)
345 DO j=1-OLy,sNy+OLy
346 DO i=1-OLx,sNx+OLx
347 HSNOW(i,j,bi,bj) = 0.2 _d 0 * AREA(i,j,bi,bj)
348 ENDDO
349 ENDDO
350 ENDDO
351 ENDDO
352
353 C-- Read initial snow thickness from file if available.
354 IF ( HsnowFile .NE. ' ' ) THEN
355 CALL READ_FLD_XY_RL( HsnowFile, ' ', HSNOW, 0, myThid )
356 _EXCH_XY_RL(HSNOW,myThid)
357 DO bj=myByLo(myThid),myByHi(myThid)
358 DO bi=myBxLo(myThid),myBxHi(myThid)
359 DO j=1-OLy,sNy+OLy
360 DO i=1-OLx,sNx+OLx
361 HSNOW(i,j,bi,bj) = MAX(HSNOW(i,j,bi,bj),ZERO)
362 ENDDO
363 ENDDO
364 ENDDO
365 ENDDO
366 ENDIF
367
368 #ifdef SEAICE_VARIABLE_SALINITY
369 DO bj=myByLo(myThid),myByHi(myThid)
370 DO bi=myBxLo(myThid),myBxHi(myThid)
371 DO j=1-OLy,sNy+OLy
372 DO i=1-OLx,sNx+OLx
373 HSALT(i,j,bi,bj)=HEFF(i,j,bi,bj)*salt(i,j,kSurface,bi,bj)*
374 & SEAICE_rhoIce*SEAICE_saltFrac
375 cif & ICE2WATR*rhoConstFresh*SEAICE_saltFrac
376
377 ENDDO
378 ENDDO
379 ENDDO
380 ENDDO
381
382 C-- Read initial sea ice salinity from file if available.
383 IF ( HsaltFile .NE. ' ' ) THEN
384 CALL READ_FLD_XY_RL( HsaltFile, ' ', HSALT, 0, myThid )
385 _EXCH_XY_RL(HSALT,myThid)
386 ENDIF
387 #endif /* SEAICE_VARIABLE_SALINITY */
388
389 #ifdef ALLOW_SITRACER
390 C-- Read initial sea ice age from file if available.
391 DO iTr = 1, SItrMaxNum
392 IF ( SItrFile(iTr) .NE. ' ' ) THEN
393 CALL READ_FLD_XY_RL( siTrFile(iTr), ' ',
394 & SItracer(1-OLx,1-OLy,1,1,iTr), 0, myThid )
395 _EXCH_XY_RL(SItracer(1-OLx,1-OLy,1,1,iTr),myThid)
396 ENDIF
397 ENDDO
398 #endif /* ALLOW_SITRACER */
399
400 ENDIF
401
402 #if (defined (ALLOW_MEAN_SFLUX_COST_CONTRIBUTION) || defined (ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION))
403 DO bj=myByLo(myThid),myByHi(myThid)
404 DO bi=myBxLo(myThid),myBxHi(myThid)
405 DO j=1-OLy,sNy+OLy
406 DO i=1-OLx,sNx+OLx
407 AREAforAtmFW(i,j,bi,bj) = AREA(i,j,bi,bj)
408 ENDDO
409 ENDDO
410 ENDDO
411 ENDDO
412 #endif
413
414 #ifdef ALLOW_OBCS
415 C-- In case we use scheme with a large stencil that extends into overlap:
416 C no longer needed with the right masking in advection & diffusion S/R.
417 c IF ( useOBCS ) THEN
418 c DO bj=myByLo(myThid),myByHi(myThid)
419 c DO bi=myBxLo(myThid),myBxHi(myThid)
420 c CALL OBCS_COPY_TRACER( HEFF(1-OLx,1-OLy,bi,bj),
421 c I 1, bi, bj, myThid )
422 c CALL OBCS_COPY_TRACER( AREA(1-OLx,1-OLy,bi,bj),
423 c I 1, bi, bj, myThid )
424 c CALL OBCS_COPY_TRACER( HSNOW(1-OLx,1-OLy,bi,bj),
425 c I 1, bi, bj, myThid )
426 #ifdef SEAICE_VARIABLE_SALINITY
427 c CALL OBCS_COPY_TRACER( HSALT(1-OLx,1-OLy,bi,bj),
428 c I 1, bi, bj, myThid )
429 #endif
430 c ENDDO
431 c ENDDO
432 c ENDIF
433 #endif /* ALLOW_OBCS */
434
435 C--- Complete initialization
436 PSTAR = SEAICE_strength
437 DO bj=myByLo(myThid),myByHi(myThid)
438 DO bi=myBxLo(myThid),myBxHi(myThid)
439 DO j=1-OLy,sNy+OLy
440 DO i=1-OLx,sNx+OLx
441 ZETA(i,j,bi,bj) = HEFF(i,j,bi,bj)*(1.0 _d 11)
442 ETA(i,j,bi,bj) = ZETA(i,j,bi,bj)/SEAICE_eccen**2
443 PRESS0(i,j,bi,bj) = PSTAR*HEFF(i,j,bi,bj)
444 & *EXP(-20.0 _d 0*(ONE-AREA(i,j,bi,bj)))
445 ZMAX(I,J,bi,bj) = SEAICE_zetaMaxFac*PRESS0(I,J,bi,bj)
446 ZMIN(i,j,bi,bj) = SEAICE_zetaMin
447 PRESS0(i,j,bi,bj) = PRESS0(i,j,bi,bj)*HEFFM(i,j,bi,bj)
448 ENDDO
449 ENDDO
450 IF ( useRealFreshWaterFlux .AND. .NOT.useThSIce ) THEN
451 DO j=1-OLy,sNy+OLy
452 DO i=1-OLx,sNx+OLx
453 sIceLoad(i,j,bi,bj) = HEFF(i,j,bi,bj)*SEAICE_rhoIce
454 & + HSNOW(i,j,bi,bj)*SEAICE_rhoSnow
455
456 ENDDO
457 ENDDO
458 ENDIF
459 ENDDO
460 ENDDO
461
462 RETURN
463 END

  ViewVC Help
Powered by ViewVC 1.1.22