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

Contents of /MITgcm_contrib/torge/grease/code/seaice_init_fixed.F

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


Revision 1.1 - (show annotations) (download)
Thu Jan 16 23:37:42 2014 UTC (11 years, 6 months ago) by torge
Branch: MAIN
CVS Tags: HEAD
introducing a grease ice category
by making use of the seaice_tracer infrastructure

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

  ViewVC Help
Powered by ViewVC 1.1.22