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

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

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


Revision 1.3 - (hide annotations) (download)
Wed Sep 26 17:03:29 2012 UTC (12 years, 10 months ago) by torge
Branch: MAIN
Changes since 1.2: +3 -2 lines
change documentation of Hlimit calculation

1 dimitri 1.1 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_init_fixed.F,v 1.19 2012/03/11 13:41:38 jmc Exp $
2     C $Name: $
3    
4     #include "SEAICE_OPTIONS.h"
5    
6     CStartOfInterface
7     SUBROUTINE SEAICE_INIT_FIXED( myThid )
8     C *==========================================================*
9     C | SUBROUTINE SEAICE_INIT_FIXED
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 "FFIELDS.h"
21     #include "SEAICE_SIZE.h"
22     #include "SEAICE_PARAMS.h"
23     #include "SEAICE.h"
24     #include "SEAICE_TRACER.h"
25    
26     C === Routine arguments ===
27     C myThid - Thread no. that called this routine.
28     INTEGER myThid
29     CEndOfInterface
30    
31     C === Local variables ===
32 dimitri 1.2 CToM<<<
33     #ifdef SEAICE_ITD
34     C msgBuf :: Informational/error message buffer
35     CHARACTER*(MAX_LEN_MBUF) msgBuf
36     CHARACTER*15 HlimitMsgFormat
37     C k - loop counter for ITD categories
38     INTEGER k
39     #endif
40     C>>>ToM
41 dimitri 1.1 C i,j,k,bi,bj - Loop counters
42    
43     INTEGER i, j, bi, bj
44     INTEGER kSurface
45     #ifdef ALLOW_SITRACER
46     INTEGER iTracer
47     #endif
48     #ifndef SEAICE_CGRID
49     _RS mask_uice
50     #endif
51     #ifdef SHORTWAVE_HEATING
52     cif Helper variable for determining the fraction of sw radiation
53     cif penetrating the model shallowest layer
54     _RL dummyTime
55     _RL swfracba(2)
56     _RL tmpFac
57     #endif /* SHORTWAVE_HEATING */
58    
59     IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
60     kSurface = Nr
61     ELSE
62     kSurface = 1
63     ENDIF
64    
65     C Initialize MNC variable information for SEAICE
66     IF ( useMNC .AND.
67     & (seaice_tave_mnc.OR.seaice_dump_mnc.OR.SEAICE_mon_mnc)
68     & ) THEN
69     CALL SEAICE_MNC_INIT( myThid )
70     ENDIF
71    
72     _BEGIN_MASTER(myThid)
73     #ifdef SHORTWAVE_HEATING
74     tmpFac = -1.0
75     dummyTime = 1.0
76     swfracba(1) = ABS(rF(1))
77     swfracba(2) = ABS(rF(2))
78     CALL SWFRAC(
79     I 2, tmpFac,
80     U swfracba,
81     I dummyTime, 0, myThid )
82     SWFracB = swfracba(2)
83     #else /* SHORTWAVE_HEATING */
84     SWFracB = 0. _d 0
85     #endif /* SHORTWAVE_HEATING */
86     _END_MASTER(myThid)
87    
88     C-- Set mcPheePiston coeff (if still unset)
89     _BEGIN_MASTER(myThid)
90     IF ( SEAICE_mcPheePiston.EQ.UNSET_RL ) THEN
91     IF ( SEAICE_availHeatFrac.NE.UNSET_RL ) THEN
92     SEAICE_mcPheePiston = SEAICE_availHeatFrac
93     & * drF(kSurface)/SEAICE_deltaTtherm
94     ELSE
95     SEAICE_mcPheePiston = MCPHEE_TAPER_FAC
96     & * STANTON_NUMBER * USTAR_BASE
97     SEAICE_mcPheePiston = MIN( SEAICE_mcPheePiston,
98     & drF(kSurface)/SEAICE_deltaTtherm )
99     ENDIF
100     ENDIF
101     _END_MASTER(myThid)
102    
103     C-- SItracer specifications for basic tracers
104     #ifdef ALLOW_SITRACER
105     _BEGIN_MASTER(myThid)
106     DO iTracer = 1, SItrNumInUse
107     C "ice concentration" tracer that should remain .EQ.1.
108     IF (SItrName(iTracer).EQ.'one') THEN
109     SItrFromOcean0(iTracer) =ONE
110     SItrFromFlood0(iTracer) =ONE
111     SItrExpand0(iTracer) =ONE
112     SItrFromOceanFrac(iTracer) =ZERO
113     SItrFromFloodFrac(iTracer) =ZERO
114     ENDIF
115     C age tracer: no age in ocean, or effect from ice cover changes
116     IF (SItrName(iTracer).EQ.'age') THEN
117     SItrFromOcean0(iTracer) =ZERO
118     SItrFromFlood0(iTracer) =ZERO
119     SItrExpand0(iTracer) =ZERO
120     SItrFromOceanFrac(iTracer) =ZERO
121     SItrFromFloodFrac(iTracer) =ZERO
122     ENDIf
123     C salinity tracer:
124     IF (SItrName(iTracer).EQ.'salinity') THEN
125     SItrMate(iTracer) ='HEFF'
126     SItrExpand0(iTracer) =ZERO
127     IF ( SEAICE_salinityTracer ) THEN
128     SEAICE_salt0 = ZERO
129     SEAICE_saltFrac = ZERO
130     ENDIF
131     ENDIF
132     C simple, made up, ice surface roughness index prototype
133     IF (SItrName(iTracer).EQ.'ridge') THEN
134     SItrMate(iTracer) ='AREA'
135     SItrFromOcean0(iTracer) =ZERO
136     SItrFromFlood0(iTracer) =ZERO
137     SItrExpand0(iTracer) =ZERO
138     SItrFromOceanFrac(iTracer) =ZERO
139     SItrFromFloodFrac(iTracer) =ZERO
140     ENDIF
141     ENDDO
142     _END_MASTER(myThid)
143     #endif
144    
145     C-- all threads wait for master to finish initialisation of shared params
146     _BARRIER
147    
148     C-- Initialize grid info
149     DO bj=myByLo(myThid),myByHi(myThid)
150     DO bi=myBxLo(myThid),myBxHi(myThid)
151     DO j=1-OLy,sNy+OLy
152     DO i=1-OLx,sNx+OLx
153     HEFFM(i,j,bi,bj) = 0. _d 0
154     ENDDO
155     ENDDO
156     DO j=1-OLy,sNy+OLy
157     DO i=1-OLx,sNx+OLx
158     HEFFM(i,j,bi,bj)= 1. _d 0
159     IF (_hFacC(i,j,kSurface,bi,bj).eq.0.)
160     & HEFFM(i,j,bi,bj)= 0. _d 0
161     ENDDO
162     ENDDO
163     #ifndef SEAICE_CGRID
164     DO j=1-OLy+1,sNy+OLy
165     DO i=1-OLx+1,sNx+OLx
166     UVM(i,j,bi,bj)=0. _d 0
167     mask_uice=HEFFM(i,j, bi,bj)+HEFFM(i-1,j-1,bi,bj)
168     & +HEFFM(i,j-1,bi,bj)+HEFFM(i-1,j, bi,bj)
169     IF(mask_uice.GT.3.5 _d 0) UVM(i,j,bi,bj)=1. _d 0
170     ENDDO
171     ENDDO
172     #endif /* SEAICE_CGRID */
173     ENDDO
174     ENDDO
175    
176 dimitri 1.2 CToM<<<
177     #ifdef SEAICE_ITD
178     C use Equ. 22 of Lipscomb et al. (2001, JGR) to generate ice
179     C thickness category limits:
180 torge 1.3 C - dependends on given number of categories nITD
181 dimitri 1.2 C - choose between original parameters of Lipscomb et al. (2001):
182     C Lipscomb: c1=3.0/N, c2=15*c1, c3=3.0
183 torge 1.3 C or emphasize thin end of ITD in order to enhance ice growth:
184 dimitri 1.2 C ToM: c1=1.5/N, c2=42*c1, c3=3.3
185 torge 1.3 C (parameters c1 to c3 are set in seaice_readparms.F)
186 dimitri 1.2 C---+-|--1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
187     Hlimit(0) = 0.0
188     Hlimit_c1=Hlimit_c1/float(nITD)
189     Hlimit_c2=Hlimit_c2*Hlimit_c1
190     DO k=1,nITD-1
191     Hlimit(k) = Hlimit(k-1)
192     & + Hlimit_c1
193     & + Hlimit_c2
194     & *( 1.+tanh( Hlimit_c3 *(float(k-1)/float(nITD) -1. ) ) )
195     ENDDO
196     C thickest category is unlimited
197     Hlimit(nITD) = 999.9
198    
199     WRITE(msgBuf,'(A,I2,A)')
200     & ' SEAICE_INIT_FIXED: ', nITD,
201     & ' sea ice thickness categories'
202     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
203     & SQUEEZE_RIGHT , myThid)
204     WRITE(HlimitMsgFormat,'(A,I2,A)') '(A,',nITD,'F6.2,F6.1)'
205     WRITE(msgBuf,HlimitMsgFormat)
206     & ' SEAICE_INIT_FIXED: Hlimit = ',
207     & Hlimit(0:nITD-1),Hlimit(nITD)
208     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
209     & SQUEEZE_RIGHT , myThid)
210    
211     #endif
212     C>>>ToM
213 dimitri 1.1 C coefficients for metric terms
214     DO bj=myByLo(myThid),myByHi(myThid)
215     DO bi=myBxLo(myThid),myBxHi(myThid)
216     #ifdef SEAICE_CGRID
217     DO j=1-OLy,sNy+OLy
218     DO i=1-OLx,sNx+OLx
219     k1AtC(I,J,bi,bj) = 0.0 _d 0
220     k1AtZ(I,J,bi,bj) = 0.0 _d 0
221     k2AtC(I,J,bi,bj) = 0.0 _d 0
222     k2AtZ(I,J,bi,bj) = 0.0 _d 0
223     ENDDO
224     ENDDO
225     IF ( usingSphericalPolarGrid .AND. SEAICEuseMetricTerms ) THEN
226     C This is the only case where tan(phi) is not zero. In this case
227     C C and U points, and Z and V points have the same phi, so that we
228     C only need a copy here. Do not use tan(YC) and tan(YG), because these
229     C can be the geographical coordinates and not the correct grid
230     C coordinates when the grid is rotated (phi/theta/psiEuler .NE. 0)
231     DO j=1-OLy,sNy+OLy
232     DO i=1-OLx,sNx+OLx
233     k2AtC(I,J,bi,bj) = - _tanPhiAtU(I,J,bi,bj)*recip_rSphere
234     k2AtZ(I,J,bi,bj) = - _tanPhiAtV(I,J,bi,bj)*recip_rSphere
235     ENDDO
236     ENDDO
237     ELSEIF ( usingCurvilinearGrid .AND. SEAICEuseMetricTerms ) THEN
238     C compute metric term coefficients from finite difference approximation
239     DO j=1-OLy,sNy+OLy
240     DO i=1-OLx,sNx+OLx-1
241     k1AtC(I,J,bi,bj) = _recip_dyF(I,J,bi,bj)
242     & * ( _dyG(I+1,J,bi,bj) - _dyG(I,J,bi,bj) )
243     & * _recip_dxF(I,J,bi,bj)
244     ENDDO
245     ENDDO
246     DO j=1-OLy,sNy+OLy
247     DO i=1-OLx+1,sNx+OLx
248     k1AtZ(I,J,bi,bj) = _recip_dyU(I,J,bi,bj)
249     & * ( _dyC(I,J,bi,bj) - _dyC(I-1,J,bi,bj) )
250     & * _recip_dxV(I,J,bi,bj)
251     ENDDO
252     ENDDO
253     DO j=1-OLy,sNy+OLy-1
254     DO i=1-OLx,sNx+OLx
255     k2AtC(I,J,bi,bj) = _recip_dxF(I,J,bi,bj)
256     & * ( _dxG(I,J+1,bi,bj) - _dxG(I,J,bi,bj) )
257     & * _recip_dyF(I,J,bi,bj)
258     ENDDO
259     ENDDO
260     DO j=1-OLy+1,sNy+OLy
261     DO i=1-OLx,sNx+OLx
262     k2AtZ(I,J,bi,bj) = _recip_dxV(I,J,bi,bj)
263     & * ( _dxC(I,J,bi,bj) - _dxC(I,J-1,bi,bj) )
264     & * _recip_dyU(I,J,bi,bj)
265     ENDDO
266     ENDDO
267     ENDIF
268     #else /* not SEAICE_CGRID */
269     DO j=1-OLy,sNy+OLy
270     DO i=1-OLx,sNx+OLx
271     k1AtC(I,J,bi,bj) = 0.0 _d 0
272     k1AtU(I,J,bi,bj) = 0.0 _d 0
273     k1AtV(I,J,bi,bj) = 0.0 _d 0
274     k2AtC(I,J,bi,bj) = 0.0 _d 0
275     k2AtU(I,J,bi,bj) = 0.0 _d 0
276     k2AtV(I,J,bi,bj) = 0.0 _d 0
277     ENDDO
278     ENDDO
279     IF ( usingSphericalPolarGrid .AND. SEAICEuseMetricTerms ) THEN
280     C This is the only case where tan(phi) is not zero. In this case
281     C C and U points, and Z and V points have the same phi, so that we
282     C only need a copy here. Do not use tan(YC) and tan(YG), because these
283     C can be the geographical coordinates and not the correct grid
284     C coordinates when the grid is rotated (phi/theta/psiEuler .NE. 0)
285     DO j=1-OLy,sNy+OLy
286     DO i=1-OLx,sNx+OLx
287     k2AtC(I,J,bi,bj) = - _tanPhiAtU(I,J,bi,bj)*recip_rSphere
288     k2AtU(I,J,bi,bj) = - _tanPhiAtU(I,J,bi,bj)*recip_rSphere
289     k2AtV(I,J,bi,bj) = - _tanPhiAtV(I,J,bi,bj)*recip_rSphere
290     ENDDO
291     ENDDO
292     ELSEIF ( usingCurvilinearGrid .AND. SEAICEuseMetricTerms ) THEN
293     C compute metric term coefficients from finite difference approximation
294     DO j=1-OLy,sNy+OLy
295     DO i=1-OLx,sNx+OLx-1
296     k1AtC(I,J,bi,bj) = _recip_dyF(I,J,bi,bj)
297     & * ( _dyG(I+1,J,bi,bj) - _dyG(I,J,bi,bj) )
298     & * _recip_dxF(I,J,bi,bj)
299     ENDDO
300     ENDDO
301     DO j=1-OLy,sNy+OLy
302     DO i=1-OLx+1,sNx+OLx
303     k1AtU(I,J,bi,bj) = _recip_dyG(I,J,bi,bj)
304     & * ( _dyF(I,J,bi,bj) - _dyF(I-1,J,bi,bj) )
305     & * _recip_dxC(I,J,bi,bj)
306     ENDDO
307     ENDDO
308     DO j=1-OLy,sNy+OLy
309     DO i=1-OLx,sNx+OLx-1
310     k1AtV(I,J,bi,bj) = _recip_dyC(I,J,bi,bj)
311     & * ( _dyU(I+1,J,bi,bj) - _dyU(I,J,bi,bj) )
312     & * _recip_dxG(I,J,bi,bj)
313     ENDDO
314     ENDDO
315     DO j=1-OLy,sNy+OLy-1
316     DO i=1-OLx,sNx+OLx
317     k2AtC(I,J,bi,bj) = _recip_dxF(I,J,bi,bj)
318     & * ( _dxG(I,J+1,bi,bj) - _dxG(I,J,bi,bj) )
319     & * _recip_dyF(I,J,bi,bj)
320     ENDDO
321     ENDDO
322     DO j=1-OLy,sNy+OLy-1
323     DO i=1-OLx,sNx+OLx
324     k2AtU(I,J,bi,bj) = _recip_dxC(I,J,bi,bj)
325     & * ( _dxV(I,J+1,bi,bj) - _dxV(I,J,bi,bj) )
326     & * _recip_dyG(I,J,bi,bj)
327     ENDDO
328     ENDDO
329     DO j=1-OLy+1,sNy+OLy
330     DO i=1-OLx,sNx+OLx
331     k2AtV(I,J,bi,bj) = _recip_dxG(I,J,bi,bj)
332     & * ( _dxF(I,J,bi,bj) - _dxF(I,J-1,bi,bj) )
333     & * _recip_dyC(I,J,bi,bj)
334     ENDDO
335     ENDDO
336     ENDIF
337     #endif /* not SEAICE_CGRID */
338     ENDDO
339     ENDDO
340    
341     #ifndef SEAICE_CGRID
342     C-- Choose a proxy level for geostrophic velocity,
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     KGEO(i,j,bi,bj) = 0
348     ENDDO
349     ENDDO
350     DO j=1-OLy,sNy+OLy
351     DO i=1-OLx,sNx+OLx
352     #ifdef SEAICE_BICE_STRESS
353     KGEO(i,j,bi,bj) = 1
354     #else /* SEAICE_BICE_STRESS */
355     IF (klowc(i,j,bi,bj) .LT. 2) THEN
356     KGEO(i,j,bi,bj) = 1
357     ELSE
358     KGEO(i,j,bi,bj) = 2
359     DO WHILE ( abs(rC(KGEO(i,j,bi,bj))) .LT. 50.0 _d 0 .AND.
360     & KGEO(i,j,bi,bj) .LT. (klowc(i,j,bi,bj)-1) )
361     KGEO(i,j,bi,bj) = KGEO(i,j,bi,bj) + 1
362     ENDDO
363     ENDIF
364     #endif /* SEAICE_BICE_STRESS */
365     ENDDO
366     ENDDO
367     ENDDO
368     ENDDO
369     #endif /* SEAICE_CGRID */
370    
371     #ifdef ALLOW_DIAGNOSTICS
372     IF ( useDiagnostics ) THEN
373     CALL SEAICE_DIAGNOSTICS_INIT( myThid )
374     ENDIF
375     #endif
376    
377     C-- Summarise pkg/seaice configuration
378     CALL SEAICE_SUMMARY( myThid )
379    
380     RETURN
381     END

  ViewVC Help
Powered by ViewVC 1.1.22