/[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.5 - (hide annotations) (download)
Mon Dec 10 22:19:49 2012 UTC (12 years, 7 months ago) by torge
Branch: MAIN
Changes since 1.4: +10 -1 lines
include updates from main branch

1 torge 1.5 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_init_fixed.F,v 1.21 2012/11/07 09:46:18 mlosch Exp $
2 dimitri 1.1 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 torge 1.4 IF (nITD .gt. 1) THEN
191     DO k=1,nITD-1
192     Hlimit(k) = Hlimit(k-1)
193     & + Hlimit_c1
194     & + Hlimit_c2
195     & *( 1.+tanh( Hlimit_c3 *(float(k-1)/float(nITD) -1. ) ) )
196     ENDDO
197     ENDIF
198 dimitri 1.2 C thickest category is unlimited
199     Hlimit(nITD) = 999.9
200    
201     WRITE(msgBuf,'(A,I2,A)')
202     & ' SEAICE_INIT_FIXED: ', nITD,
203     & ' sea ice thickness categories'
204     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
205     & SQUEEZE_RIGHT , myThid)
206     WRITE(HlimitMsgFormat,'(A,I2,A)') '(A,',nITD,'F6.2,F6.1)'
207     WRITE(msgBuf,HlimitMsgFormat)
208     & ' SEAICE_INIT_FIXED: Hlimit = ',
209     & Hlimit(0:nITD-1),Hlimit(nITD)
210     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
211     & SQUEEZE_RIGHT , myThid)
212    
213     #endif
214     C>>>ToM
215 dimitri 1.1 C coefficients for metric terms
216     DO bj=myByLo(myThid),myByHi(myThid)
217     DO bi=myBxLo(myThid),myBxHi(myThid)
218     #ifdef SEAICE_CGRID
219     DO j=1-OLy,sNy+OLy
220     DO i=1-OLx,sNx+OLx
221     k1AtC(I,J,bi,bj) = 0.0 _d 0
222     k1AtZ(I,J,bi,bj) = 0.0 _d 0
223     k2AtC(I,J,bi,bj) = 0.0 _d 0
224     k2AtZ(I,J,bi,bj) = 0.0 _d 0
225     ENDDO
226     ENDDO
227     IF ( usingSphericalPolarGrid .AND. SEAICEuseMetricTerms ) THEN
228     C This is the only case where tan(phi) is not zero. In this case
229     C C and U points, and Z and V points have the same phi, so that we
230     C only need a copy here. Do not use tan(YC) and tan(YG), because these
231     C can be the geographical coordinates and not the correct grid
232     C coordinates when the grid is rotated (phi/theta/psiEuler .NE. 0)
233     DO j=1-OLy,sNy+OLy
234     DO i=1-OLx,sNx+OLx
235     k2AtC(I,J,bi,bj) = - _tanPhiAtU(I,J,bi,bj)*recip_rSphere
236     k2AtZ(I,J,bi,bj) = - _tanPhiAtV(I,J,bi,bj)*recip_rSphere
237     ENDDO
238     ENDDO
239     ELSEIF ( usingCurvilinearGrid .AND. SEAICEuseMetricTerms ) THEN
240     C compute metric term coefficients from finite difference approximation
241     DO j=1-OLy,sNy+OLy
242     DO i=1-OLx,sNx+OLx-1
243     k1AtC(I,J,bi,bj) = _recip_dyF(I,J,bi,bj)
244     & * ( _dyG(I+1,J,bi,bj) - _dyG(I,J,bi,bj) )
245     & * _recip_dxF(I,J,bi,bj)
246     ENDDO
247     ENDDO
248     DO j=1-OLy,sNy+OLy
249     DO i=1-OLx+1,sNx+OLx
250     k1AtZ(I,J,bi,bj) = _recip_dyU(I,J,bi,bj)
251     & * ( _dyC(I,J,bi,bj) - _dyC(I-1,J,bi,bj) )
252     & * _recip_dxV(I,J,bi,bj)
253     ENDDO
254     ENDDO
255     DO j=1-OLy,sNy+OLy-1
256     DO i=1-OLx,sNx+OLx
257     k2AtC(I,J,bi,bj) = _recip_dxF(I,J,bi,bj)
258     & * ( _dxG(I,J+1,bi,bj) - _dxG(I,J,bi,bj) )
259     & * _recip_dyF(I,J,bi,bj)
260     ENDDO
261     ENDDO
262     DO j=1-OLy+1,sNy+OLy
263     DO i=1-OLx,sNx+OLx
264     k2AtZ(I,J,bi,bj) = _recip_dxV(I,J,bi,bj)
265     & * ( _dxC(I,J,bi,bj) - _dxC(I,J-1,bi,bj) )
266     & * _recip_dyU(I,J,bi,bj)
267     ENDDO
268     ENDDO
269     ENDIF
270     #else /* not SEAICE_CGRID */
271     DO j=1-OLy,sNy+OLy
272     DO i=1-OLx,sNx+OLx
273     k1AtC(I,J,bi,bj) = 0.0 _d 0
274     k1AtU(I,J,bi,bj) = 0.0 _d 0
275     k1AtV(I,J,bi,bj) = 0.0 _d 0
276     k2AtC(I,J,bi,bj) = 0.0 _d 0
277     k2AtU(I,J,bi,bj) = 0.0 _d 0
278     k2AtV(I,J,bi,bj) = 0.0 _d 0
279     ENDDO
280     ENDDO
281     IF ( usingSphericalPolarGrid .AND. SEAICEuseMetricTerms ) THEN
282     C This is the only case where tan(phi) is not zero. In this case
283     C C and U points, and Z and V points have the same phi, so that we
284     C only need a copy here. Do not use tan(YC) and tan(YG), because these
285     C can be the geographical coordinates and not the correct grid
286     C coordinates when the grid is rotated (phi/theta/psiEuler .NE. 0)
287     DO j=1-OLy,sNy+OLy
288     DO i=1-OLx,sNx+OLx
289     k2AtC(I,J,bi,bj) = - _tanPhiAtU(I,J,bi,bj)*recip_rSphere
290     k2AtU(I,J,bi,bj) = - _tanPhiAtU(I,J,bi,bj)*recip_rSphere
291     k2AtV(I,J,bi,bj) = - _tanPhiAtV(I,J,bi,bj)*recip_rSphere
292     ENDDO
293     ENDDO
294     ELSEIF ( usingCurvilinearGrid .AND. SEAICEuseMetricTerms ) THEN
295     C compute metric term coefficients from finite difference approximation
296     DO j=1-OLy,sNy+OLy
297     DO i=1-OLx,sNx+OLx-1
298     k1AtC(I,J,bi,bj) = _recip_dyF(I,J,bi,bj)
299     & * ( _dyG(I+1,J,bi,bj) - _dyG(I,J,bi,bj) )
300     & * _recip_dxF(I,J,bi,bj)
301     ENDDO
302     ENDDO
303     DO j=1-OLy,sNy+OLy
304     DO i=1-OLx+1,sNx+OLx
305     k1AtU(I,J,bi,bj) = _recip_dyG(I,J,bi,bj)
306     & * ( _dyF(I,J,bi,bj) - _dyF(I-1,J,bi,bj) )
307     & * _recip_dxC(I,J,bi,bj)
308     ENDDO
309     ENDDO
310     DO j=1-OLy,sNy+OLy
311     DO i=1-OLx,sNx+OLx-1
312     k1AtV(I,J,bi,bj) = _recip_dyC(I,J,bi,bj)
313     & * ( _dyU(I+1,J,bi,bj) - _dyU(I,J,bi,bj) )
314     & * _recip_dxG(I,J,bi,bj)
315     ENDDO
316     ENDDO
317     DO j=1-OLy,sNy+OLy-1
318     DO i=1-OLx,sNx+OLx
319     k2AtC(I,J,bi,bj) = _recip_dxF(I,J,bi,bj)
320     & * ( _dxG(I,J+1,bi,bj) - _dxG(I,J,bi,bj) )
321     & * _recip_dyF(I,J,bi,bj)
322     ENDDO
323     ENDDO
324     DO j=1-OLy,sNy+OLy-1
325     DO i=1-OLx,sNx+OLx
326     k2AtU(I,J,bi,bj) = _recip_dxC(I,J,bi,bj)
327     & * ( _dxV(I,J+1,bi,bj) - _dxV(I,J,bi,bj) )
328     & * _recip_dyG(I,J,bi,bj)
329     ENDDO
330     ENDDO
331     DO j=1-OLy+1,sNy+OLy
332     DO i=1-OLx,sNx+OLx
333     k2AtV(I,J,bi,bj) = _recip_dxG(I,J,bi,bj)
334     & * ( _dxF(I,J,bi,bj) - _dxF(I,J-1,bi,bj) )
335     & * _recip_dyC(I,J,bi,bj)
336     ENDDO
337     ENDDO
338     ENDIF
339     #endif /* not SEAICE_CGRID */
340     ENDDO
341     ENDDO
342    
343     #ifndef SEAICE_CGRID
344     C-- Choose a proxy level for geostrophic velocity,
345     DO bj=myByLo(myThid),myByHi(myThid)
346     DO bi=myBxLo(myThid),myBxHi(myThid)
347     DO j=1-OLy,sNy+OLy
348     DO i=1-OLx,sNx+OLx
349     KGEO(i,j,bi,bj) = 0
350     ENDDO
351     ENDDO
352     DO j=1-OLy,sNy+OLy
353     DO i=1-OLx,sNx+OLx
354     #ifdef SEAICE_BICE_STRESS
355     KGEO(i,j,bi,bj) = 1
356     #else /* SEAICE_BICE_STRESS */
357     IF (klowc(i,j,bi,bj) .LT. 2) THEN
358     KGEO(i,j,bi,bj) = 1
359     ELSE
360     KGEO(i,j,bi,bj) = 2
361     DO WHILE ( abs(rC(KGEO(i,j,bi,bj))) .LT. 50.0 _d 0 .AND.
362     & KGEO(i,j,bi,bj) .LT. (klowc(i,j,bi,bj)-1) )
363     KGEO(i,j,bi,bj) = KGEO(i,j,bi,bj) + 1
364     ENDDO
365     ENDIF
366     #endif /* SEAICE_BICE_STRESS */
367     ENDDO
368     ENDDO
369     ENDDO
370     ENDDO
371     #endif /* SEAICE_CGRID */
372    
373 torge 1.5 #ifdef SEAICE_ALLOW_JFNK
374     C initialise some diagnostic counters for the JFNK solver
375     totalNewtonIters = 0
376     totalNewtonFails = 0
377     totalKrylovIters = 0
378     totalKrylovFails = 0
379     totalJFNKtimeSteps = 0
380     #endif /* SEAICE_ALLOW_JFNK */
381    
382 dimitri 1.1 #ifdef ALLOW_DIAGNOSTICS
383     IF ( useDiagnostics ) THEN
384     CALL SEAICE_DIAGNOSTICS_INIT( myThid )
385     ENDIF
386     #endif
387    
388     C-- Summarise pkg/seaice configuration
389     CALL SEAICE_SUMMARY( myThid )
390    
391     RETURN
392     END

  ViewVC Help
Powered by ViewVC 1.1.22