/[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.2 - (hide annotations) (download)
Fri Apr 27 22:25:23 2012 UTC (13 years, 3 months ago) by dimitri
Branch: MAIN
Changes since 1.1: +45 -0 lines
first check in of itd code

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

  ViewVC Help
Powered by ViewVC 1.1.22