/[MITgcm]/MITgcm_contrib/torge/itd/code/seaice_init_varia.F
ViewVC logotype

Annotation of /MITgcm_contrib/torge/itd/code/seaice_init_varia.F

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


Revision 1.1 - (hide annotations) (download)
Fri Apr 27 22:22:17 2012 UTC (13 years, 3 months ago) by dimitri
Branch: MAIN
check-in original code, before itd modifications
seaice_advdiff.F,v 1.60 2012/02/16 01:22:02
seaice_check_pickup.F,v 1.7 2012/03/05 15:21:44
seaice_diagnostics_init.F,v 1.33 2012/02/16 01:22:02
seaice_growth.F,v 1.162 2012/03/15 03:07:31
seaice_init_fixed.F,v 1.19 2012/03/11 13:41:38
seaice_init_varia.F,v 1.72 2012/03/14 22:55:53
seaice_readparms.F,v 1.120 2012/03/14 22:55:53
seaice_write_pickup.F,v 1.14 2012/03/05 15:21:45
seaice_read_pickup.F,v 1.16 2012/03/05 15:21:44
seaice_model.F,v 1.100 2012/03/02 18:56:06
SEAICE.h,v 1.62 2012/03/06 16:51:21
SEAICE_OPTIONS.h,v 1.63 2012/03/08 01:15:02
SEAICE_PARAMS.h,v 1.91 2012/03/11 13:41:38
SEAICE_SIZE.h,v 1.5 2012/03/06 16:51:21
SIZE.h,v 1.28 2009/05/17 21:15:07

1 dimitri 1.1 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_init_varia.F,v 1.72 2012/03/14 22:55:53 heimbach 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