/[MITgcm]/MITgcm/pkg/seaice/seaice_growth.F
ViewVC logotype

Contents of /MITgcm/pkg/seaice/seaice_growth.F

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

Revision 1.124 - (show annotations) (download)
Tue Jun 7 03:58:23 2011 UTC (13 years ago) by gforget
Branch: MAIN
Changes since 1.123: +29 -6 lines 
- introducing ALLOW_SITRACER and seaice_tracer_phys.F to handle generic seaice tracer.
  For now it covers, and was tested for, salinity and age (work in progress).
- introducing siEps (1e-5, parameter, defined in SEAICE_PARAMS.h).
1 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_growth.F,v 1.123 2011/05/28 22:39:20 gforget Exp $
2 C $Name: $
3
4 #include "SEAICE_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: SEAICE_GROWTH
8 C !INTERFACE:
9 SUBROUTINE SEAICE_GROWTH( myTime, myIter, myThid )
10 C !DESCRIPTION: \bv
11 C *==========================================================*
12 C | SUBROUTINE seaice_growth
13 C | o Updata ice thickness and snow depth
14 C *==========================================================*
15 C \ev
16
17 C !USES:
18 IMPLICIT NONE
19 C === Global variables ===
20 #include "SIZE.h"
21 #include "EEPARAMS.h"
22 #include "PARAMS.h"
23 #include "DYNVARS.h"
24 #include "GRID.h"
25 #include "FFIELDS.h"
26 #include "SEAICE_SIZE.h"
27 #include "SEAICE_PARAMS.h"
28 #include "SEAICE.h"
29 #include "SEAICE_TRACER.h"
30 #ifdef ALLOW_EXF
31 # include "EXF_OPTIONS.h"
32 # include "EXF_FIELDS.h"
33 # include "EXF_PARAM.h"
34 #endif
35 #ifdef ALLOW_SALT_PLUME
36 # include "SALT_PLUME.h"
37 #endif
38 #ifdef ALLOW_AUTODIFF_TAMC
39 # include "tamc.h"
40 #endif
41
42 C !INPUT/OUTPUT PARAMETERS:
43 C === Routine arguments ===
44 C myTime :: Simulation time
45 C myIter :: Simulation timestep number
46 C myThid :: Thread no. that called this routine.
47 _RL myTime
48 INTEGER myIter, myThid
49
50 C !FUNCTIONS:
51 #ifdef ALLOW_DIAGNOSTICS
52 LOGICAL DIAGNOSTICS_IS_ON
53 EXTERNAL DIAGNOSTICS_IS_ON
54 #endif
55
56 C !LOCAL VARIABLES:
57 C === Local variables ===
58 C
59 C unit/sign convention:
60 C Within the thermodynamic computation all stocks, except HSNOW,
61 C are in 'effective ice meters' units, and >0 implies more ice.
62 C This holds for stocks due to ocean and atmosphere heat,
63 C at the outset of 'PART 2: determine heat fluxes/stocks'
64 C and until 'PART 7: determine ocean model forcing'
65 C This strategy minimizes the need for multiplications/divisions
66 C by ice fraction, heat capacity, etc. The only conversions that
67 C occurs are for the HSNOW (in effective snow meters) and
68 C PRECIP (fresh water m/s).
69 C
70 C HEFF is effective Hice thickness (m3/m2)
71 C HSNOW is Heffective snow thickness (m3/m2)
72 C HSALT is Heffective salt content (g/m2)
73 C AREA is the seaice cover fraction (0<=AREA<=1)
74 C Q denotes heat stocks -- converted to ice stocks (m3/m2) early on
75 C
76 C For all other stocks/increments, such as d_HEFFbyATMonOCN
77 C or a_QbyATM_cover, the naming convention is as follows:
78 C The prefix 'a_' means available, the prefix 'd_' means delta
79 C (i.e. increment), and the prefix 'r_' means residual.
80 C The suffix '_cover' denotes a value for the ice covered fraction
81 C of the grid cell, whereas '_open' is for the open water fraction.
82 C The main part of the name states what ice/snow stock is concerned
83 C (e.g. QbyATM or HEFF), and how it is affected (e.g. d_HEFFbyATMonOCN
84 C is the increment of HEFF due to the ATMosphere extracting heat from the
85 C OCeaN surface, or providing heat to the OCeaN surface).
86
87 CEOP
88 C i,j,bi,bj :: Loop counters
89 INTEGER i, j, bi, bj
90 C number of surface interface layer
91 INTEGER kSurface
92 C constants
93 _RL TBC, ICE2SNOW
94 _RL QI, QS
95
96 C a_QbyATM_cover :: available heat (in W/m^2) due to the interaction of
97 C the atmosphere and the ocean surface - for ice covered water
98 C a_QbyATM_open :: same but for open water
99 C r_QbyATM_cover :: residual of a_QbyATM_cover after freezing/melting processes
100 C r_QbyATM_open :: same but for open water
101 _RL a_QbyATM_cover (1:sNx,1:sNy)
102 _RL a_QbyATM_open (1:sNx,1:sNy)
103 _RL r_QbyATM_cover (1:sNx,1:sNy)
104 _RL r_QbyATM_open (1:sNx,1:sNy)
105 C a_QSWbyATM_open - short wave heat flux over ocean in W/m^2
106 C a_QSWbyATM_cover - short wave heat flux under ice in W/m^2
107 _RL a_QSWbyATM_open (1:sNx,1:sNy)
108 _RL a_QSWbyATM_cover (1:sNx,1:sNy)
109 C a_QbyOCN :: available heat (in in W/m^2) due to the
110 C interaction of the ice pack and the ocean surface
111 C r_QbyOCN :: residual of a_QbyOCN after freezing/melting
112 C processes have been accounted for
113 _RL a_QbyOCN (1:sNx,1:sNy)
114 _RL r_QbyOCN (1:sNx,1:sNy)
115
116 C conversion factors to go from Q (W/m2) to HEFF (ice meters)
117 _RL convertQ2HI, convertHI2Q
118 C conversion factors to go from precip (m/s) unit to HEFF (ice meters)
119 _RL convertPRECIP2HI, convertHI2PRECIP
120
121 C ICE/SNOW stocks tendencies associated with the various melt/freeze processes
122 _RL d_AREAbyATM (1:sNx,1:sNy)
123 _RL d_AREAbyOCN (1:sNx,1:sNy)
124 _RL d_AREAbyICE (1:sNx,1:sNy)
125
126 c The change of mean ice thickness due to out-of-bounds values following
127 c sea ice dyhnamics
128 _RL d_HEFFbyNEG (1:sNx,1:sNy)
129
130 c The change of mean ice thickness due to turbulent ocean-sea ice heat fluxes
131 _RL d_HEFFbyOCNonICE (1:sNx,1:sNy)
132
133 c The sum of mean ice thickness increments due to atmospheric fluxes over the open water
134 c fraction and ice-covered fractions of the grid cell
135 _RL d_HEFFbyATMonOCN (1:sNx,1:sNy)
136 c The change of mean ice thickness due to flooding by snow
137 _RL d_HEFFbyFLOODING (1:sNx,1:sNy)
138
139 c The mean ice thickness increments due to atmospheric fluxes over the open water
140 c fraction and ice-covered fractions of the grid cell, respectively
141 _RL d_HEFFbyATMonOCN_open(1:sNx,1:sNy)
142 _RL d_HEFFbyATMonOCN_cover(1:sNx,1:sNy)
143
144 _RL d_HSNWbyNEG (1:sNx,1:sNy)
145 _RL d_HSNWbyATMonSNW (1:sNx,1:sNy)
146 _RL d_HSNWbyOCNonSNW (1:sNx,1:sNy)
147 _RL d_HSNWbyRAIN (1:sNx,1:sNy)
148
149 _RL d_HFRWbyRAIN (1:sNx,1:sNy)
150 C
151 C a_FWbySublim :: fresh water flux implied by latent heat of
152 C sublimation to atmosphere, same sign convention
153 C as EVAP (positive upward)
154 _RL a_FWbySublim (1:sNx,1:sNy)
155 #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
156 _RL d_HEFFbySublim (1:sNx,1:sNy)
157 _RL d_HSNWbySublim (1:sNx,1:sNy)
158 _RL rodt, rrodt
159 #endif /* SEAICE_ADD_SUBLIMATION_TO_FWBUDGET */
160
161 C actual ice thickness (with upper and lower limit)
162 _RL heffActual (1:sNx,1:sNy)
163 C actual snow thickness
164 _RL hsnowActual (1:sNx,1:sNy)
165 C actual ice thickness (with lower limit only) Reciprocal
166 _RL recip_heffActual (1:sNx,1:sNy)
167 C local value (=HO or HO_south)
168 _RL HO_loc
169
170 C AREA_PRE :: hold sea-ice fraction field before any seaice-thermo update
171 _RL AREApreTH (1:sNx,1:sNy)
172 _RL HEFFpreTH (1:sNx,1:sNy)
173 _RL HSNWpreTH (1:sNx,1:sNy)
174
175 C wind speed
176 _RL UG (1:sNx,1:sNy)
177 #ifdef ALLOW_ATM_WIND
178 _RL SPEED_SQ
179 #endif
180
181 C pathological cases thresholds
182 _RL heffTooHeavy
183
184 C temporary variables available for the various computations
185 #ifdef SEAICE_GROWTH_LEGACY
186 _RL tmpscal0
187 #endif
188 _RL tmpscal1, tmpscal2, tmpscal3, tmpscal4
189 _RL tmparr1 (1:sNx,1:sNy)
190
191 #ifdef ALLOW_SEAICE_FLOODING
192 _RL hDraft
193 #endif /* ALLOW_SEAICE_FLOODING */
194
195 #ifdef SEAICE_VARIABLE_SALINITY
196 _RL saltFluxAdjust (1:sNx,1:sNy)
197 #endif
198
199 INTEGER nDim
200 #ifdef SEAICE_MULTICATEGORY
201 INTEGER ilockey
202 PARAMETER ( nDim = MULTDIM )
203 INTEGER it
204 _RL pFac
205 _RL heffActualP (1:sNx,1:sNy)
206 _RL a_QbyATMmult_cover (1:sNx,1:sNy)
207 _RL a_QSWbyATMmult_cover(1:sNx,1:sNy)
208 _RL a_FWbySublimMult (1:sNx,1:sNy)
209 #else
210 PARAMETER ( nDim = 1 )
211 #endif /* SEAICE_MULTICATEGORY */
212
213 #ifdef MCPHEE_OCEAN_ICE_HEAT_FLUX
214 C Factor by which we increase the upper ocean friction velocity (u*) when
215 C ice is absent in a grid cell (dimensionless)
216 _RL MixedLayerTurbulenceFactor
217
218 c The Stanton number for the McPhee
219 c ocean-ice heat flux parameterization (dimensionless)
220 _RL STANTON_NUMBER
221
222 c A typical friction velocity beneath sea ice for the
223 c McPhee heat flux parameterization (m/s)
224 _RL USTAR_BASE
225
226 _RL surf_theta
227 #endif
228
229 #ifdef ALLOW_DIAGNOSTICS
230 c Helper variables for diagnostics
231 _RL DIAGarray (1:sNx,1:sNy)
232 _RL DIAGarrayA (1:sNx,1:sNy)
233 _RL DIAGarrayB (1:sNx,1:sNy)
234 _RL DIAGarrayC (1:sNx,1:sNy)
235 _RL DIAGarrayD (1:sNx,1:sNy)
236 #endif
237
238 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
239
240 C ===================================================================
241 C =================PART 0: constants and initializations=============
242 C ===================================================================
243
244 IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
245 kSurface = Nr
246 ELSE
247 kSurface = 1
248 ENDIF
249
250 C Cutoff for iceload
251 heffTooHeavy=dRf(kSurface) / 5. _d 0
252
253 C FREEZING TEMP. OF SEA WATER (deg C)
254 TBC = SEAICE_freeze
255
256 C RATIO OF SEA ICE DENSITY to SNOW DENSITY
257 ICE2SNOW = SEAICE_rhoIce/SEAICE_rhoSnow
258
259 C HEAT OF FUSION OF ICE (J/m^3)
260 QI = SEAICE_rhoIce*SEAICE_lhFusion
261 C HEAT OF FUSION OF SNOW (J/m^3)
262 QS = SEAICE_rhoSnow*SEAICE_lhFusion
263 C
264 C note that QI/QS=ICE2SNOW -- except it wasnt in old code
265
266 #ifdef MCPHEE_OCEAN_ICE_HEAT_FLUX
267 STANTON_NUMBER = 0.0056 _d 0
268 USTAR_BASE = 0.0125 _d 0
269 #endif
270
271 C conversion factors to go from Q (W/m2) to HEFF (ice meters)
272 convertQ2HI=SEAICE_deltaTtherm/QI
273 convertHI2Q=1/convertQ2HI
274 C conversion factors to go from precip (m/s) unit to HEFF (ice meters)
275 convertPRECIP2HI=SEAICE_deltaTtherm*rhoConstFresh/SEAICE_rhoIce
276 convertHI2PRECIP=1./convertPRECIP2HI
277
278 DO bj=myByLo(myThid),myByHi(myThid)
279 DO bi=myBxLo(myThid),myBxHi(myThid)
280
281 #ifdef ALLOW_AUTODIFF_TAMC
282 act1 = bi - myBxLo(myThid)
283 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
284 act2 = bj - myByLo(myThid)
285 max2 = myByHi(myThid) - myByLo(myThid) + 1
286 act3 = myThid - 1
287 max3 = nTx*nTy
288 act4 = ikey_dynamics - 1
289 iicekey = (act1 + 1) + act2*max1
290 & + act3*max1*max2
291 & + act4*max1*max2*max3
292 #endif /* ALLOW_AUTODIFF_TAMC */
293
294
295 C array initializations
296 C =====================
297
298 DO J=1,sNy
299 DO I=1,sNx
300 a_QbyATM_cover (I,J) = 0.0 _d 0
301 a_QbyATM_open(I,J) = 0.0 _d 0
302 r_QbyATM_cover (I,J) = 0.0 _d 0
303 r_QbyATM_open (I,J) = 0.0 _d 0
304
305 a_QSWbyATM_open (I,J) = 0.0 _d 0
306 a_QSWbyATM_cover (I,J) = 0.0 _d 0
307
308 a_QbyOCN (I,J) = 0.0 _d 0
309 r_QbyOCN (I,J) = 0.0 _d 0
310
311 d_AREAbyATM(I,J) = 0.0 _d 0
312 d_AREAbyICE(I,J) = 0.0 _d 0
313 d_AREAbyOCN(I,J) = 0.0 _d 0
314
315 d_HEFFbyNEG(I,J) = 0.0 _d 0
316 d_HEFFbyOCNonICE(I,J) = 0.0 _d 0
317 d_HEFFbyATMonOCN(I,J) = 0.0 _d 0
318 d_HEFFbyFLOODING(I,J) = 0.0 _d 0
319
320 d_HEFFbyATMonOCN_open(I,J) = 0.0 _d 0
321 d_HEFFbyATMonOCN_cover(I,J) = 0.0 _d 0
322
323 d_HSNWbyNEG(I,J) = 0.0 _d 0
324 d_HSNWbyATMonSNW(I,J) = 0.0 _d 0
325 d_HSNWbyOCNonSNW(I,J) = 0.0 _d 0
326 d_HSNWbyRAIN(I,J) = 0.0 _d 0
327 a_FWbySublim(I,J) = 0.0 _d 0
328 #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
329 d_HEFFbySublim(I,J) = 0.0 _d 0
330 d_HSNWbySublim(I,J) = 0.0 _d 0
331 #endif /* SEAICE_ADD_SUBLIMATION_TO_FWBUDGET */
332 c
333 d_HFRWbyRAIN(I,J) = 0.0 _d 0
334
335 tmparr1(I,J) = 0.0 _d 0
336
337 #ifdef SEAICE_VARIABLE_SALINITY
338 saltFluxAdjust(I,J) = 0.0 _d 0
339 #endif
340 #ifdef SEAICE_MULTICATEGORY
341 a_QbyATMmult_cover(I,J) = 0.0 _d 0
342 a_QSWbyATMmult_cover(I,J) = 0.0 _d 0
343 a_FWbySublimMult(I,J) = 0.0 _d 0
344 #endif /* SEAICE_MULTICATEGORY */
345 ENDDO
346 ENDDO
347 #ifdef ALLOW_MEAN_SFLUX_COST_CONTRIBUTION
348 DO J=1-oLy,sNy+oLy
349 DO I=1-oLx,sNx+oLx
350 frWtrAtm(I,J,bi,bj) = 0.0 _d 0
351 ENDDO
352 ENDDO
353 #endif
354
355
356 C =====================================================================
357 C ===========PART 1: treat pathological cases (post advdiff)===========
358 C =====================================================================
359
360 #ifdef SEAICE_GROWTH_LEGACY
361
362 DO J=1,sNy
363 DO I=1,sNx
364 HEFFpreTH(I,J)=HEFFNM1(I,J,bi,bj)
365 HSNWpreTH(I,J)=HSNOW(I,J,bi,bj)
366 AREApreTH(I,J)=AREANM1(I,J,bi,bj)
367 d_HEFFbyNEG(I,J)=0. _d 0
368 d_HSNWbyNEG(I,J)=0. _d 0
369 ENDDO
370 ENDDO
371
372 #else /* SEAICE_GROWTH_LEGACY */
373
374 #ifdef ALLOW_AUTODIFF_TAMC
375 #ifdef SEAICE_MODIFY_GROWTH_ADJ
376 Cgf no dependency through pathological cases treatment
377 if ( SEAICEadjMODE.EQ.0 ) then
378 #endif
379 #endif
380
381 C 1) treat the case of negative values:
382
383 #ifdef ALLOW_AUTODIFF_TAMC
384 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
385 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
386 CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
387 #endif /* ALLOW_AUTODIFF_TAMC */
388 DO J=1,sNy
389 DO I=1,sNx
390 d_HEFFbyNEG(I,J)=MAX(-HEFF(I,J,bi,bj),0. _d 0)
391 HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J)
392 d_HSNWbyNEG(I,J)=MAX(-HSNOW(I,J,bi,bj),0. _d 0)
393 HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J)
394 AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),0. _d 0)
395 ENDDO
396 ENDDO
397
398 C 1.25) treat the case of very thin ice:
399
400 #ifdef ALLOW_AUTODIFF_TAMC
401 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
402 #endif /* ALLOW_AUTODIFF_TAMC */
403 DO J=1,sNy
404 DO I=1,sNx
405 if (HEFF(I,J,bi,bj).LE.siEps) then
406 tmpscal2=-HEFF(I,J,bi,bj)
407 tmpscal3=-HSNOW(I,J,bi,bj)
408 TICE(I,J,bi,bj)=celsius2K
409 #ifdef SEAICE_AGE
410 IceAgeTr(i,j,bi,bj,2)=0. _d 0
411 #endif /* SEAICE_AGE */
412 else
413 tmpscal2=0. _d 0
414 tmpscal3=0. _d 0
415 endif
416 HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+tmpscal2
417 d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
418 HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3
419 d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
420 ENDDO
421 ENDDO
422
423 C 1.5) treat the case of area but no ice/snow:
424
425 #ifdef ALLOW_AUTODIFF_TAMC
426 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
427 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
428 #endif /* ALLOW_AUTODIFF_TAMC */
429 DO J=1,sNy
430 DO I=1,sNx
431 IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.
432 & (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0
433 ENDDO
434 ENDDO
435
436 C 2) treat the case of very small area:
437
438 #ifdef ALLOW_PRECLUDE_INFINITESIMAL_AREA
439 #ifdef ALLOW_AUTODIFF_TAMC
440 CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
441 #endif /* ALLOW_AUTODIFF_TAMC */
442 DO J=1,sNy
443 DO I=1,sNx
444 IF ((HEFF(i,j,bi,bj).GT.0).OR.(HSNOW(i,j,bi,bj).GT.0)) THEN
445 #ifdef SEAICE_AGE
446 IF (AREA(I,J,bi,bj).LT.areaMin) THEN
447 IceAgeTr(i,j,bi,bj,1)=IceAgeTr(i,j,bi,bj,1)
448 & + IceAgeTr(i,j,bi,bj,1)
449 & *(areaMin-AREA(I,J,bi,bj)) / areaMin
450 ENDIF
451 #endif /* SEAICE_AGE */
452 AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),areaMin)
453 ENDIF
454 ENDDO
455 ENDDO
456 #endif
457
458 C 2.5) treat case of excessive ice cover, e.g., due to ridging:
459
460 #ifdef ALLOW_AUTODIFF_TAMC
461 CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
462 #endif /* ALLOW_AUTODIFF_TAMC */
463 DO J=1,sNy
464 DO I=1,sNx
465 #ifdef SEAICE_AGE
466 IF (AREA(I,J,bi,bj).GT.areaMax) THEN
467 IceAgeTr(i,j,bi,bj,1)=IceAgeTr(i,j,bi,bj,1)
468 & - IceAgeTr(i,j,bi,bj,1)
469 & *(AREA(I,J,bi,bj)-areaMax) / areaMax
470 ENDIF
471 #endif /* SEAICE_AGE */
472 AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),areaMax)
473 ENDDO
474 ENDDO
475
476 #ifdef ALLOW_AUTODIFF_TAMC
477 #ifdef SEAICE_MODIFY_GROWTH_ADJ
478 endif
479 #endif
480 #endif
481
482 C 3) store regularized values of heff, hsnow, area at the onset of thermo.
483 DO J=1,sNy
484 DO I=1,sNx
485 HEFFpreTH(I,J)=HEFF(I,J,bi,bj)
486 HSNWpreTH(I,J)=HSNOW(I,J,bi,bj)
487 AREApreTH(I,J)=AREA(I,J,bi,bj)
488 #ifdef ALLOW_SITRACER
489 SItrHEFF(I,J,bi,bj,1)=HEFF(I,J,bi,bj)
490 #endif
491 ENDDO
492 ENDDO
493
494 C 4) treat sea ice salinity pathological cases
495 #ifdef SEAICE_VARIABLE_SALINITY
496 #ifdef ALLOW_AUTODIFF_TAMC
497 CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
498 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
499 #endif /* ALLOW_AUTODIFF_TAMC */
500 DO J=1,sNy
501 DO I=1,sNx
502 IF ( (HSALT(I,J,bi,bj) .LT. 0.0).OR.
503 & (HEFF(I,J,bi,bj) .EQ. 0.0) ) THEN
504 saltFluxAdjust(I,J) = - HEFFM(I,J,bi,bj) *
505 & HSALT(I,J,bi,bj) / SEAICE_deltaTtherm
506 HSALT(I,J,bi,bj) = 0.0 _d 0
507 ENDIF
508 ENDDO
509 ENDDO
510 #endif /* SEAICE_VARIABLE_SALINITY */
511
512 C 5) treat sea ice age pathological cases
513 C ...
514
515 #ifdef SEAICE_AGE
516
517 #ifdef ALLOW_AUTODIFF_TAMC
518 CADJ STORE IceAgeTr(:,:,bi,bj,1) = comlev1_bibj,
519 CADJ & key = iicekey,byte=isbyte
520 CADJ STORE IceAgeTr(:,:,bi,bj,2) = comlev1_bibj,
521 CADJ & key = iicekey,byte=isbyte
522 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
523 CADJ & key = iicekey,byte=isbyte
524 #endif /* ALLOW_AUTODIFF_TAMC */
525 DO J=1,sNy
526 DO I=1,sNx
527 IF (HEFF(i,j,bi,bj).EQ.0.) THEN
528 IceAgeTr(i,j,bi,bj,1)=0. _d 0
529 IceAgeTr(i,j,bi,bj,2)=0. _d 0
530 ENDIF
531 IF (IceAgeTr(i,j,bi,bj,1).LT.0.) IceAgeTr(i,j,bi,bj,1)=0. _d 0
532 IF (IceAgeTr(i,j,bi,bj,2).LT.0.) IceAgeTr(i,j,bi,bj,2)=0. _d 0
533 ENDDO
534 ENDDO
535 #endif /* SEAICE_AGE */
536
537 #endif /* SEAICE_GROWTH_LEGACY */
538
539 #ifdef ALLOW_AUTODIFF_TAMC
540 #ifdef SEAICE_MODIFY_GROWTH_ADJ
541 Cgf no additional dependency of air-sea fluxes to ice
542 if ( SEAICEadjMODE.GE.1 ) then
543 DO J=1,sNy
544 DO I=1,sNx
545 HEFFpreTH(I,J) = 0. _d 0
546 HSNWpreTH(I,J) = 0. _d 0
547 AREApreTH(I,J) = 0. _d 0
548 ENDDO
549 ENDDO
550 endif
551 #endif
552 #endif
553
554 C 4) COMPUTE ACTUAL ICE/SNOW THICKNESS; USE MIN/MAX VALUES
555 C TO REGULARIZE SEAICE_SOLVE4TEMP/d_AREA COMPUTATIONS
556
557 #ifdef ALLOW_AUTODIFF_TAMC
558 CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
559 CADJ STORE HEFFpreTH = comlev1_bibj, key = iicekey, byte = isbyte
560 CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte
561 #endif /* ALLOW_AUTODIFF_TAMC */
562 DO J=1,sNy
563 DO I=1,sNx
564 tmpscal1 = MAX(areaMin,AREApreTH(I,J))
565 hsnowActual(I,J) = HSNWpreTH(I,J)/tmpscal1
566 tmpscal2 = HEFFpreTH(I,J)/tmpscal1
567 #ifdef SEAICE_GROWTH_LEGACY
568 heffActual(I,J) = MAX(tmpscal2,hiceMin)
569 #else
570 heffActual(I,J) = sqrt(tmpscal2*tmpscal2 + hiceMin*hiceMin)
571 c here we dont impose the area minimum
572 recip_heffActual(I,J) = AREApreTH(I,J) /
573 & sqrt(HEFFpreTH(I,J)*HEFFpreTH(I,J) + hiceMin*hiceMin)
574 #endif
575 ENDDO
576 ENDDO
577
578 #ifdef ALLOW_AUTODIFF_TAMC
579 #ifdef SEAICE_SIMPLIFY_GROWTH_ADJ
580 CALL ZERO_ADJ_1D( sNx*sNy, heffActual, myThid)
581 CALL ZERO_ADJ_1D( sNx*sNy, hsnowActual, myThid)
582 CALL ZERO_ADJ_1D( sNx*sNy, recip_heffActual, myThid)
583 #endif
584 #endif
585
586
587
588 C ===================================================================
589 C ================PART 2: determine heat fluxes/stocks===============
590 C ===================================================================
591
592 C determine available heat due to the atmosphere -- for open water
593 C ================================================================
594
595 C ocean surface/mixed layer temperature
596 DO J=1,sNy
597 DO I=1,sNx
598 TMIX(I,J,bi,bj)=theta(I,J,kSurface,bi,bj)+celsius2K
599 ENDDO
600 ENDDO
601
602 C wind speed from exf
603 DO J=1,sNy
604 DO I=1,sNx
605 UG(I,J) = MAX(SEAICE_EPS,wspeed(I,J,bi,bj))
606 ENDDO
607 ENDDO
608
609 #ifdef ALLOW_AUTODIFF_TAMC
610 CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
611 CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
612 cCADJ STORE UG = comlev1_bibj, key = iicekey,byte=isbyte
613 cCADJ STORE TMIX(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
614 #endif /* ALLOW_AUTODIFF_TAMC */
615
616 CALL SEAICE_BUDGET_OCEAN(
617 I UG,
618 U TMIX,
619 O a_QbyATM_open, a_QSWbyATM_open,
620 I bi, bj, myTime, myIter, myThid )
621
622 C determine available heat due to the atmosphere -- for ice covered water
623 C =======================================================================
624
625 #ifdef ALLOW_ATM_WIND
626 IF (useRelativeWind) THEN
627 C Compute relative wind speed over sea ice.
628 DO J=1,sNy
629 DO I=1,sNx
630 SPEED_SQ =
631 & (uWind(I,J,bi,bj)
632 & +0.5 _d 0*(uVel(i,j,kSurface,bi,bj)
633 & +uVel(i+1,j,kSurface,bi,bj))
634 & -0.5 _d 0*(uice(i,j,bi,bj)+uice(i+1,j,bi,bj)))**2
635 & +(vWind(I,J,bi,bj)
636 & +0.5 _d 0*(vVel(i,j,kSurface,bi,bj)
637 & +vVel(i,j+1,kSurface,bi,bj))
638 & -0.5 _d 0*(vice(i,j,bi,bj)+vice(i,j+1,bi,bj)))**2
639 IF ( SPEED_SQ .LE. SEAICE_EPS_SQ ) THEN
640 UG(I,J)=SEAICE_EPS
641 ELSE
642 UG(I,J)=SQRT(SPEED_SQ)
643 ENDIF
644 ENDDO
645 ENDDO
646 ENDIF
647 #endif
648
649 #ifdef ALLOW_AUTODIFF_TAMC
650 CADJ STORE tice = comlev1, key = ikey_dynamics, byte = isbyte
651 CADJ STORE hsnowActual = comlev1_bibj, key = iicekey, byte = isbyte
652 CADJ STORE heffActual = comlev1_bibj, key = iicekey, byte = isbyte
653 CADJ STORE UG = comlev1_bibj, key = iicekey, byte = isbyte
654 # ifdef SEAICE_MULTICATEGORY
655 CADJ STORE tices = comlev1, key = ikey_dynamics, byte = isbyte
656 # endif /* SEAICE_MULTICATEGORY */
657 #endif /* ALLOW_AUTODIFF_TAMC */
658
659 C-- Start loop over multi-categories, if SEAICE_MULTICATEGORY is undefined
660 C nDim = 1, and there is only one loop iteration
661 #ifdef SEAICE_MULTICATEGORY
662 DO IT=1,nDim
663 #ifdef ALLOW_AUTODIFF_TAMC
664 C Why do we need this store directive when we have just stored
665 C TICES before the loop?
666 ilockey = (iicekey-1)*nDim + IT
667 CADJ STORE tices(:,:,it,bi,bj) = comlev1_multdim,
668 CADJ & key = ilockey, byte = isbyte
669 #endif /* ALLOW_AUTODIFF_TAMC */
670 pFac = (2.0 _d 0*real(IT)-1.0 _d 0)/nDim
671 DO J=1,sNy
672 DO I=1,sNx
673 heffActualP(I,J)= heffActual(I,J)*pFac
674 TICE(I,J,bi,bj)=TICES(I,J,IT,bi,bj)
675 ENDDO
676 ENDDO
677 CALL SEAICE_SOLVE4TEMP(
678 I UG, heffActualP, hsnowActual,
679 U TICE,
680 O a_QbyATMmult_cover, a_QSWbyATMmult_cover,
681 O a_FWbySublimMult,
682 I bi, bj, myTime, myIter, myThid )
683 DO J=1,sNy
684 DO I=1,sNx
685 C average over categories
686 a_QbyATM_cover (I,J) = a_QbyATM_cover(I,J)
687 & + a_QbyATMmult_cover(I,J)/nDim
688 a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
689 & + a_QSWbyATMmult_cover(I,J)/nDim
690 a_FWbySublim (I,J) = a_FWbySublim(I,J)
691 & + a_FWbySublimMult(I,J)/nDim
692 TICES(I,J,IT,bi,bj) = TICE(I,J,bi,bj)
693 ENDDO
694 ENDDO
695 ENDDO
696 #else
697 CALL SEAICE_SOLVE4TEMP(
698 I UG, heffActual, hsnowActual,
699 U TICE,
700 O a_QbyATM_cover, a_QSWbyATM_cover, a_FWbySublim,
701 I bi, bj, myTime, myIter, myThid )
702 #endif /* SEAICE_MULTICATEGORY */
703 C-- End loop over multi-categories
704
705 #ifdef ALLOW_DIAGNOSTICS
706 IF ( useDiagnostics ) THEN
707 IF ( DIAGNOSTICS_IS_ON('SIatmQnt',myThid) ) THEN
708 DO J=1,sNy
709 DO I=1,sNx
710 CML If I consider the atmosphere above the ice, the surface flux
711 CML which is relevant for the air temperature dT/dt Eq
712 CML accounts for sensible and radiation (with different treatment
713 CML according to wave-length) fluxes but not for "latent heat flux",
714 CML since it does not contribute to heating the air.
715 CML So this diagnostic is only good for heat budget calculations within
716 CML the ice-ocean system.
717 DIAGarray(I,J) = maskC(I,J,kSurface,bi,bj) * (
718 & a_QbyATM_cover(I,J) * AREApreTH(I,J)
719 #ifndef SEAICE_GROWTH_LEGACY
720 & + a_QSWbyATM_cover(I,J) * AREApreTH(I,J)
721 #endif /* SEAICE_GROWTH_LEGACY */
722 & + a_QbyATM_open(I,J) * ( ONE - AREApreTH(I,J) ) )
723 ENDDO
724 ENDDO
725 CALL DIAGNOSTICS_FILL(DIAGarray,'SIatmQnt',0,1,3,bi,bj,myThid)
726 ENDIF
727 IF ( DIAGNOSTICS_IS_ON('SIfwSubl',myThid) ) THEN
728 DO J=1,sNy
729 DO I=1,sNx
730 DIAGarray(I,J) = maskC(I,J,kSurface,bi,bj) *
731 & a_FWbySublim(I,J) * AREApreTH(I,J)
732 ENDDO
733 ENDDO
734 CALL DIAGNOSTICS_FILL(DIAGarray,'SIfwSubl',0,1,3,bi,bj,myThid)
735 ENDIF
736 IF ( DIAGNOSTICS_IS_ON('SIatmFW ',myThid) ) THEN
737 DO J=1,sNy
738 DO I=1,sNx
739 DIAGarray(I,J) = maskC(I,J,kSurface,bi,bj)*(
740 & PRECIP(I,J,bi,bj)
741 & - EVAP(I,J,bi,bj)
742 & *( ONE - AREApreTH(I,J) )
743 #ifdef ALLOW_RUNOFF
744 & + RUNOFF(I,J,bi,bj)
745 #endif /* ALLOW_RUNOFF */
746 & )*rhoConstFresh
747 #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
748 & - a_FWbySublim(I,J)*AREApreTH(I,J)
749 #endif /* SEAICE_ADD_SUBLIMATION_TO_FWBUDGET */
750 ENDDO
751 ENDDO
752 CALL DIAGNOSTICS_FILL(DIAGarray,'SIatmFW ',0,1,3,bi,bj,myThid)
753 ENDIF
754 ENDIF
755 #endif /* ALLOW_DIAGNOSTICS */
756
757 C switch heat fluxes from W/m2 to 'effective' ice meters
758 DO J=1,sNy
759 DO I=1,sNx
760 a_QbyATM_cover(I,J) = a_QbyATM_cover(I,J)
761 & * convertQ2HI * AREApreTH(I,J)
762 a_QSWbyATM_cover(I,J) = a_QSWbyATM_cover(I,J)
763 & * convertQ2HI * AREApreTH(I,J)
764 a_QbyATM_open(I,J) = a_QbyATM_open(I,J)
765 & * convertQ2HI * ( ONE - AREApreTH(I,J) )
766 a_QSWbyATM_open(I,J) = a_QSWbyATM_open(I,J)
767 & * convertQ2HI * ( ONE - AREApreTH(I,J) )
768 C and initialize r_QbyATM_cover/r_QbyATM_open
769 r_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)
770 r_QbyATM_open(I,J)=a_QbyATM_open(I,J)
771 #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
772 C Convert fresh water flux by sublimation to 'effective' ice meters.
773 C Negative sublimation is resublimation and will be added as snow.
774 a_FWbySublim(I,J) = SEAICE_deltaTtherm/SEAICE_rhoIce
775 & * a_FWbySublim(I,J)*AREApreTH(I,J)
776 #endif /* SEAICE_ADD_SUBLIMATION_TO_FWBUDGET */
777 ENDDO
778 ENDDO
779
780 #ifdef ALLOW_DIAGNOSTICS
781 IF ( useDiagnostics ) THEN
782 IF ( DIAGNOSTICS_IS_ON('SIaQbATC',myThid) ) THEN
783 CALL DIAGNOSTICS_FILL(a_QbyATM_cover,
784 & 'SIaQbATC',0,1,3,bi,bj,myThid)
785 ENDIF
786 IF ( DIAGNOSTICS_IS_ON('SIaQbATO',myThid) ) THEN
787 CALL DIAGNOSTICS_FILL(a_QbyATM_open,
788 & 'SIaQbATO',0,1,3,bi,bj,myThid)
789 ENDIF
790 ENDIF
791 #endif
792 #ifdef ALLOW_AUTODIFF_TAMC
793 #ifdef SEAICE_MODIFY_GROWTH_ADJ
794 Cgf no additional dependency through ice cover
795 if ( SEAICEadjMODE.GE.3 ) then
796 DO J=1,sNy
797 DO I=1,sNx
798 a_QbyATM_cover(I,J) = 0. _d 0
799 r_QbyATM_cover(I,J) = 0. _d 0
800 a_QSWbyATM_cover(I,J) = 0. _d 0
801 ENDDO
802 ENDDO
803 endif
804 #endif
805 #endif
806
807 C determine available heat due to the ice pack tying the
808 C underlying surface water temperature to freezing point
809 C ======================================================
810
811 #ifdef ALLOW_AUTODIFF_TAMC
812 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
813 #endif
814
815 DO J=1,sNy
816 DO I=1,sNx
817
818 #ifdef SEAICE_VARIABLE_FREEZING_POINT
819 TBC = -0.0575 _d 0*salt(I,J,kSurface,bi,bj) + 0.0901 _d 0
820 #endif /* SEAICE_VARIABLE_FREEZING_POINT */
821
822 #ifdef MCPHEE_OCEAN_ICE_HEAT_FLUX
823
824 c Bound the ocean temperature to be at or above the freezing point.
825 surf_theta = max(theta(I,J,kSurface,bi,bj), TBC)
826 #ifdef GRADIENT_MIXED_LAYER_TURBULENCE_FACTOR
827 MixedLayerTurbulenceFactor = 12.5 _d 0 -
828 & 11.5 _d 0 *AREApreTH(I,J)
829 #else
830 c If ice is present, MixedLayerTurbulenceFactor = 1.0, else 12.50
831 IF (AREApreTH(I,J) .GT. 0. _d 0) THEN
832 MixedLayerTurbulenceFactor = ONE
833 ELSE
834 MixedLayerTurbulenceFactor = 12.5 _d 0
835 ENDIF
836 #endif /* GRADIENT_MIXED_LAYER_TURBULENCE_FACTOR */
837
838 c The turbulent ocean-ice heat flux
839 a_QbyOCN(I,J) = -STANTON_NUMBER * USTAR_BASE * rhoConst *
840 & HeatCapacity_Cp *(surf_theta - TBC)*
841 & MixedLayerTurbulenceFactor*maskC(i,j,kSurface,bi,bj)
842
843 c The turbulent ocean-ice heat flux converted to meters
844 c of potential ice melt
845 a_QbyOCN(I,J) = a_QbyOCN(I,J) * convertQ2HI
846
847 c by design a_QbyOCN .LE. 0. so that initial ice growth cannot
848 c be triggered by this term, which Ian says is better for adjoint
849 #else
850
851 IF ( theta(I,J,kSurface,bi,bj) .GE. TBC ) THEN
852 tmpscal1 = SEAICE_availHeatFrac
853 ELSE
854 tmpscal1 = SEAICE_availHeatFracFrz
855 ENDIF
856 tmpscal1 = tmpscal1
857 & * dRf(kSurface) * maskC(i,j,kSurface,bi,bj)
858
859 a_QbyOCN(i,j) = -tmpscal1 * (HeatCapacity_Cp*rhoConst/QI)
860 & * (theta(I,J,kSurface,bi,bj)-TBC)
861
862 #endif /* MCPHEE_OCEAN_ICE_HEAT_FLUX */
863
864 #ifdef ALLOW_DIAGNOSTICS
865 DIAGarrayA (I,J) = a_QbyOCN(I,J)
866 #endif
867 r_QbyOCN(i,j) = a_QbyOCN(i,j)
868
869 ENDDO
870 ENDDO
871
872 #ifdef ALLOW_DIAGNOSTICS
873 IF ( useDiagnostics ) THEN
874 IF ( DIAGNOSTICS_IS_ON('SIaQbOCN',myThid) ) THEN
875 CALL DIAGNOSTICS_FILL(DIAGarrayA,
876 & 'SIaQbOCN',0,1,3,bi,bj,myThid)
877 ENDIF
878 ENDIF
879 #endif
880
881
882
883 #ifdef ALLOW_AUTODIFF_TAMC
884 #ifdef SEAICE_SIMPLIFY_GROWTH_ADJ
885 CALL ZERO_ADJ_1D( sNx*sNy, r_QbyOCN, myThid)
886 #endif
887 #endif
888
889
890 C ===================================================================
891 C =========PART 3: determine effective thicknesses increments========
892 C ===================================================================
893
894 C compute ice thickness tendency due to ice-ocean interaction
895 C ===========================================================
896
897 #ifdef ALLOW_AUTODIFF_TAMC
898 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
899 CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
900 #endif /* ALLOW_AUTODIFF_TAMC */
901
902 DO J=1,sNy
903 DO I=1,sNx
904 d_HEFFbyOCNonICE(I,J)=MAX(r_QbyOCN(i,j), -HEFF(I,J,bi,bj))
905 r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)
906 HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj) + d_HEFFbyOCNonICE(I,J)
907 #ifdef ALLOW_DIAGNOSTICS
908 DIAGarrayA(I,J) = d_HEFFbyOCNonICE(i,j)
909 #endif
910 #ifdef ALLOW_SITRACER
911 SItrHEFF(I,J,bi,bj,2)=HEFF(I,J,bi,bj)
912 #endif
913 ENDDO
914 ENDDO
915
916 #ifdef ALLOW_DIAGNOSTICS
917 IF ( useDiagnostics ) THEN
918 IF ( DIAGNOSTICS_IS_ON('SIdHbOCN',myThid) ) THEN
919 CALL DIAGNOSTICS_FILL(DIAGarrayA,
920 & 'SIdHbOCN',0,1,3,bi,bj,myThid)
921 ENDIF
922 ENDIF
923 #endif
924
925 #ifdef SEAICE_GROWTH_LEGACY
926 #ifdef ALLOW_DIAGNOSTICS
927 IF ( useDiagnostics ) THEN
928 IF ( DIAGNOSTICS_IS_ON('SIyneg ',myThid) ) THEN
929 CALL DIAGNOSTICS_FILL(d_HEFFbyOCNonICE,
930 & 'SIyneg ',0,1,1,bi,bj,myThid)
931 ENDIF
932 ENDIF
933 #endif
934 #endif
935
936 C compute snow melt tendency due to snow-atmosphere interaction
937 C ==================================================================
938
939 #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
940 #ifdef ALLOW_AUTODIFF_TAMC
941 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
942 CADJ STORE a_FWbySublim = comlev1_bibj,key=iicekey,byte=isbyte
943 #endif /* ALLOW_AUTODIFF_TAMC */
944 C First apply sublimation to snow
945 rodt = ICE2SNOW
946 rrodt = 1./rodt
947 DO J=1,sNy
948 DO I=1,sNx
949 IF ( a_FWbySublim(I,J) .LT. 0. _d 0 ) THEN
950 C resublimate as snow
951 d_HSNWbySublim(I,J) = -a_FWbySublim(I,J)*rodt
952 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbySublim(I,J)
953 a_FWbySublim(I,J) = 0. _d 0
954 ENDIF
955 C sublimate snow first
956 tmpscal1 = MIN(a_FWbySublim(I,J)*rodt,HSNOW(I,J,bi,bj))
957 tmpscal2 = MAX(tmpscal1,0. _d 0)
958 d_HSNWbySublim(I,J) = - tmpscal2
959 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) - tmpscal2
960 a_FWbySublim(I,J) = a_FWbySublim(I,J) - tmpscal2*rrodt
961 ENDDO
962 ENDDO
963 #endif /* SEAICE_ADD_SUBLIMATION_TO_FWBUDGET */
964
965 #ifdef ALLOW_AUTODIFF_TAMC
966 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
967 CADJ STORE r_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
968 #endif /* ALLOW_AUTODIFF_TAMC */
969
970 DO J=1,sNy
971 DO I=1,sNx
972 tmpscal1=MAX(r_QbyATM_cover(I,J)*ICE2SNOW,-HSNOW(I,J,bi,bj))
973 tmpscal2=MIN(tmpscal1,0. _d 0)
974 #ifdef SEAICE_MODIFY_GROWTH_ADJ
975 Cgf no additional dependency through snow
976 if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
977 #endif
978 d_HSNWbyATMonSNW(I,J)= tmpscal2
979 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + tmpscal2
980 r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2/ICE2SNOW
981 ENDDO
982 ENDDO
983
984 C compute ice thickness tendency due to the atmosphere
985 C ====================================================
986
987 #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
988 #ifdef ALLOW_AUTODIFF_TAMC
989 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
990 CADJ STORE a_FWbySublim = comlev1_bibj,key=iicekey,byte=isbyte
991 #endif /* ALLOW_AUTODIFF_TAMC */
992 C Apply sublimation to ice
993 rodt = 1. _d 0
994 rrodt = 1./rodt
995 DO J=1,sNy
996 DO I=1,sNx
997 C If anything is left, sublimate ice
998 tmpscal1 = MIN(a_FWbySublim(I,J)*rodt,HEFF(I,J,bi,bj))
999 tmpscal2 = MAX(tmpscal1,0. _d 0)
1000 d_HEFFbySublim(I,J) = - tmpscal2
1001 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) - tmpscal2
1002 a_FWbySublim(I,J) = a_FWbySublim(I,J) - tmpscal2*rrodt
1003 ENDDO
1004 ENDDO
1005 #endif /* SEAICE_ADD_SUBLIMATION_TO_FWBUDGET */
1006
1007 #ifdef ALLOW_AUTODIFF_TAMC
1008 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1009 CADJ STORE r_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1010 #endif /* ALLOW_AUTODIFF_TAMC */
1011
1012 Cgf note: this block is not actually tested by lab_sea
1013 Cgf where all experiments start in January. So even though
1014 Cgf the v1.81=>v1.82 revision would change results in
1015 Cgf warming conditions, the lab_sea results were not changed.
1016
1017 DO J=1,sNy
1018 DO I=1,sNx
1019
1020 #ifdef SEAICE_GROWTH_LEGACY
1021 tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J))
1022 #else
1023 tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J)+
1024 c Limit ice growth by potential melt by ocean
1025 & AREApreTH(I,J) * r_QbyOCN(I,J))
1026 #endif /* SEAICE_GROWTH_LEGACY */
1027
1028 d_HEFFbyATMonOCN_cover(I,J)=tmpscal2
1029 d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal2
1030 r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)-tmpscal2
1031 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal2
1032
1033 #ifdef ALLOW_DIAGNOSTICS
1034 DIAGarrayA(I,J) = d_HEFFbyATMonOCN_cover(I,J)
1035 #endif
1036 #ifdef ALLOW_SITRACER
1037 SItrHEFF(I,J,bi,bj,3)=HEFF(I,J,bi,bj)
1038 #endif
1039 ENDDO
1040 ENDDO
1041
1042 #ifdef ALLOW_DIAGNOSTICS
1043 IF ( useDiagnostics ) THEN
1044 IF ( DIAGNOSTICS_IS_ON('SIdHbATC',myThid) ) THEN
1045 CALL DIAGNOSTICS_FILL(DIAGarrayA,
1046 & 'SIdHbATC',0,1,3,bi,bj,myThid)
1047 ENDIF
1048 ENDIF
1049 #endif /* ALLOW DIAGNOSTICS */
1050
1051
1052 C attribute precip to fresh water or snow stock,
1053 C depending on atmospheric conditions.
1054 C =================================================
1055 #ifdef ALLOW_ATM_TEMP
1056 #ifdef ALLOW_AUTODIFF_TAMC
1057 CADJ STORE a_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1058 CADJ STORE PRECIP(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1059 CADJ STORE AREApreTH = comlev1_bibj,key=iicekey,byte=isbyte
1060 #endif /* ALLOW_AUTODIFF_TAMC */
1061 DO J=1,sNy
1062 DO I=1,sNx
1063 C possible alternatives to the a_QbyATM_cover criterium
1064 c IF (TICE(I,J,bi,bj) .LT. TMIX) THEN
1065 c IF (atemp(I,J,bi,bj) .LT. celsius2K) THEN
1066 IF ( a_QbyATM_cover(I,J).GE. 0. _d 0 ) THEN
1067 C add precip as snow
1068 d_HFRWbyRAIN(I,J)=0. _d 0
1069 d_HSNWbyRAIN(I,J)=convertPRECIP2HI*ICE2SNOW*
1070 & PRECIP(I,J,bi,bj)*AREApreTH(I,J)
1071 ELSE
1072 C add precip to the fresh water bucket
1073 d_HFRWbyRAIN(I,J)=-convertPRECIP2HI*
1074 & PRECIP(I,J,bi,bj)*AREApreTH(I,J)
1075 d_HSNWbyRAIN(I,J)=0. _d 0
1076 ENDIF
1077 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyRAIN(I,J)
1078 ENDDO
1079 ENDDO
1080 Cgf note: this does not affect air-sea heat flux,
1081 Cgf since the implied air heat gain to turn
1082 Cgf rain to snow is not a surface process.
1083 #ifdef ALLOW_DIAGNOSTICS
1084 IF ( useDiagnostics ) THEN
1085 IF ( DIAGNOSTICS_IS_ON('SIsnPrcp',myThid) ) THEN
1086 DO J=1,sNy
1087 DO I=1,sNx
1088 DIAGarray(I,J) = maskC(I,J,kSurface,bi,bj)
1089 & * d_HSNWbyRAIN(I,J)*SEAICE_rhoSnow/SEAICE_deltaTtherm
1090 ENDDO
1091 ENDDO
1092 CALL DIAGNOSTICS_FILL(DIAGarray,'SIsnPrcp',0,1,3,bi,bj,myThid)
1093 ENDIF
1094 ENDIF
1095 #endif /* ALLOW_DIAGNOSTICS */
1096 #endif /* ALLOW_ATM_TEMP */
1097
1098 C compute snow melt due to heat available from ocean.
1099 C =================================================================
1100
1101 Cgf do we need to keep this comment and cpp bracket?
1102 Cph( very sensitive bit here by JZ
1103 #ifndef SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING
1104 #ifdef ALLOW_AUTODIFF_TAMC
1105 CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1106 CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1107 #endif /* ALLOW_AUTODIFF_TAMC */
1108 DO J=1,sNy
1109 DO I=1,sNx
1110 tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW, -HSNOW(I,J,bi,bj))
1111 tmpscal2=MIN(tmpscal1,0. _d 0)
1112 #ifdef SEAICE_MODIFY_GROWTH_ADJ
1113 Cgf no additional dependency through snow
1114 if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1115 #endif
1116 d_HSNWbyOCNonSNW(I,J) = tmpscal2
1117 r_QbyOCN(I,J)=r_QbyOCN(I,J)
1118 & -d_HSNWbyOCNonSNW(I,J)/ICE2SNOW
1119 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+d_HSNWbyOCNonSNW(I,J)
1120 ENDDO
1121 ENDDO
1122 #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */
1123 Cph)
1124
1125 C gain of new ice over open water
1126 C ===============================
1127 #ifndef SEAICE_GROWTH_LEGACY
1128 #ifdef ALLOW_AUTODIFF_TAMC
1129 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1130 CADJ STORE r_QbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte
1131 CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1132 CADJ STORE a_QSWbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1133 CADJ STORE a_QSWbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte
1134 #endif /* ALLOW_AUTODIFF_TAMC */
1135 DO J=1,sNy
1136 DO I=1,sNx
1137 #ifdef ALLOW_DIAGNOSTICS
1138 DIAGarrayA(I,J) = ZERO
1139 #endif
1140
1141 #ifdef SEAICE_DO_OPEN_WATER_GROWTH
1142
1143 c Initial ice growth is triggered by open water
1144 c heat flux overcoming potential melt by ocean
1145 tmpscal1=r_QbyATM_open(I,J)+r_QbyOCN(i,j) *
1146 & (1.0 _d 0 - AREApreTH(i,J))
1147 c Penetrative shortwave flux beyond first layer
1148 c that is therefore not available to ice growth/melt
1149 tmpscal2=SWFRACB * a_QSWbyATM_open(I,J)
1150 #ifdef SEAICE_DO_OPEN_WATER_MELT
1151 C allow not only growth but also melt by open ocean heat flux
1152 tmpscal3=MAX(tmpscal1-tmpscal2, -HEFF(I,J,bi,bj))
1153 #else
1154 tmpscal3=MAX(tmpscal1-tmpscal2, 0. _d 0)
1155 #endif /* SEAICE_DO_OPEN_WATER_MELT */
1156 d_HEFFbyATMonOCN_open(I,J)=tmpscal3
1157 d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3
1158 r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3
1159 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3
1160
1161 #endif /* SEAICE_DO_OPEN_WATER_GROWTH */
1162
1163 #ifdef ALLOW_DIAGNOSTICS
1164 DIAGarrayA(I,J) = maskC(I,J,kSurface,bi,bj)
1165 & * d_HEFFbyATMonOCN_open(I,J)
1166 #endif
1167 ENDDO
1168 ENDDO
1169
1170 #ifdef ALLOW_DIAGNOSTICS
1171 IF ( useDiagnostics ) THEN
1172 IF ( DIAGNOSTICS_IS_ON('SIdHbATO',myThid) ) THEN
1173 CALL DIAGNOSTICS_FILL(DIAGarrayA,
1174 & 'SIdHbATO',0,1,3,bi,bj,myThid)
1175 ENDIF
1176 ENDIF
1177 #endif /* ALLOW DIAGNOSTICS */
1178
1179 #endif /* SEAICE_GROWTH_LEGACY */
1180
1181 #ifdef ALLOW_SITRACER
1182 DO J=1,sNy
1183 DO I=1,sNx
1184 c needs to be here to allow use also with LEGACY branch
1185 SItrHEFF(I,J,bi,bj,4)=HEFF(I,J,bi,bj)
1186 ENDDO
1187 ENDDO
1188 #endif
1189
1190 C convert snow to ice if submerged.
1191 C =================================
1192
1193 #ifndef SEAICE_GROWTH_LEGACY
1194 C note: in legacy, this process is done at the end
1195 #ifdef ALLOW_SEAICE_FLOODING
1196 #ifdef ALLOW_AUTODIFF_TAMC
1197 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1198 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1199 #endif /* ALLOW_AUTODIFF_TAMC */
1200 IF ( SEAICEuseFlooding ) THEN
1201 DO J=1,sNy
1202 DO I=1,sNx
1203 hDraft = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
1204 & +HEFF(I,J,bi,bj)*SEAICE_rhoIce)/rhoConst
1205 tmpscal1 = MAX( 0. _d 0, hDraft - HEFF(I,J,bi,bj))
1206 d_HEFFbyFLOODING(I,J)=tmpscal1
1207 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)+d_HEFFbyFLOODING(I,J)
1208 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-
1209 & d_HEFFbyFLOODING(I,J)*ICE2SNOW
1210 ENDDO
1211 ENDDO
1212 ENDIF
1213 #endif /* ALLOW_SEAICE_FLOODING */
1214 #endif /* SEAICE_GROWTH_LEGACY */
1215
1216 C ===================================================================
1217 C ==========PART 4: determine ice cover fraction increments=========-
1218 C ===================================================================
1219
1220 #ifdef ALLOW_AUTODIFF_TAMC
1221 CADJ STORE d_HEFFbyATMonOCN = comlev1_bibj,key=iicekey,byte=isbyte
1222 CADJ STORE d_HEFFbyATMonOCN_cover = comlev1_bibj,key=iicekey,byte=isbyte
1223 CADJ STORE d_HEFFbyATMonOCN_open = comlev1_bibj,key=iicekey,byte=isbyte
1224 CADJ STORE d_HEFFbyOCNonICE = comlev1_bibj,key=iicekey,byte=isbyte
1225 CADJ STORE a_QbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte
1226 CADJ STORE heffActual = comlev1_bibj,key=iicekey,byte=isbyte
1227 CADJ STORE AREApreTH = comlev1_bibj,key=iicekey,byte=isbyte
1228 CADJ STORE HEFF(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1229 CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1230 CADJ STORE AREA(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1231 #endif /* ALLOW_AUTODIFF_TAMC */
1232
1233 DO J=1,sNy
1234 DO I=1,sNx
1235
1236 IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN
1237 HO_loc=HO_south
1238 ELSE
1239 HO_loc=HO
1240 ENDIF
1241
1242 C compute contraction/expansion from melt/growth
1243 #ifdef SEAICE_GROWTH_LEGACY
1244
1245 C compute heff after ice melt by ocn:
1246 tmpscal0=HEFF(I,J,bi,bj) - d_HEFFbyATMonOCN(I,J)
1247 C compute available heat left after snow melt by atm:
1248 tmpscal1= a_QbyATM_open(I,J)+a_QbyATM_cover(I,J)
1249 & - d_HSNWbyATMonSNW(I,J)/ICE2SNOW
1250 C (cannot melt more than all the ice)
1251 tmpscal2 = MAX(-tmpscal0,tmpscal1)
1252 tmpscal3 = MIN(ZERO,tmpscal2)
1253 #ifdef ALLOW_DIAGNOSTICS
1254 DIAGarray(I,J) = tmpscal2
1255 #endif
1256 C gain of new ice over open water
1257 tmpscal4=MAX(ZERO,a_QbyATM_open(I,J))
1258 C compute cover fraction tendency
1259 d_AREAbyATM(I,J)=tmpscal4/HO_loc+HALF*tmpscal3
1260 & *AREApreTH(I,J) /(tmpscal0+.00001 _d 0)
1261
1262 #else /* SEAICE_GROWTH_LEGACY */
1263
1264 # ifdef FENTY_AREA_EXPANSION_CONTRACTION
1265
1266 C the various volume tendency terms
1267 tmpscal1 = d_HEFFbyATMOnOCN_open(I,J)
1268 tmpscal2 = d_HEFFbyATMOnOCN_cover(I,J)
1269 tmpscal3 = d_HEFFbyOCNonICE(I,J)
1270
1271 C Part 1: expansion/contraction of ice cover in open-water areas.
1272 IF (tmpscal1 .GT. ZERO) then
1273 C Ice cover is all created from open ocean-water air-sea heat fluxes
1274 d_AREAbyATM(I,J)=tmpscal1/HO_loc
1275 ELSEIF (AREApreTH(I,J) .GT. ZERO) THEN
1276 C that can also act to remove ice cover.
1277 d_AREAbyATM(i,J) = HALF*tmpscal1*recip_heffActual(I,J)
1278 ENDIF
1279
1280 C Part 2: contraction from ice thinning from above
1281 if (tmpscal2 .LE. ZERO) d_AREAbyICE(I,J) =
1282 & HALF * tmpscal2 * recip_heffActual(I,J)
1283
1284 C Part 3: contraction from ice thinning from below
1285 if (tmpscal3 .LE.ZERO) d_AREAbyOCN(I,J) =
1286 & HALF * tmpscal3* recip_heffActual(I,J)
1287
1288 # else /* FENTY_AREA_EXPANSION_CONTRACTION */
1289
1290 # ifdef SEAICE_OCN_MELT_ACT_ON_AREA
1291 C ice cover reduction by joint OCN+ATM melt
1292 tmpscal3 = MIN( 0. _d 0 ,
1293 & d_HEFFbyATMonOCN(I,J)+d_HEFFbyOCNonICE(I,J) )
1294 # else
1295 C ice cover reduction by ATM melt only -- as in legacy code
1296 tmpscal3 = MIN( 0. _d 0 , d_HEFFbyATMonOCN(I,J) )
1297 # endif
1298 C gain of new ice over open water
1299
1300 # ifdef SEAICE_DO_OPEN_WATER_GROWTH
1301 C the one effectively used to increment HEFF
1302 tmpscal4 = MAX(ZERO,d_HEFFbyATMonOCN_open(I,J))
1303 # else
1304 C the virtual one -- as in legcy code
1305 tmpscal4 = MAX(ZERO,a_QbyATM_open(I,J))
1306 # endif
1307
1308 C compute cover fraction tendency
1309 d_AREAbyATM(I,J)=tmpscal4/HO_loc+
1310 & HALF*tmpscal3*recip_heffActual(I,J)
1311
1312 # endif /* FENTY_AREA_EXPANSION_CONTRACTION */
1313
1314 #endif /* SEAICE_GROWTH_LEGACY */
1315
1316 C apply tendency
1317 IF ( (HEFF(i,j,bi,bj).GT.0. _d 0).OR.
1318 & (HSNOW(i,j,bi,bj).GT.0. _d 0) ) THEN
1319 AREA(I,J,bi,bj)=max(0. _d 0, min( 1. _d 0, AREA(I,J,bi,bj)
1320 & +d_AREAbyATM(I,J) + d_AREAbyOCN(I,J) + d_AREAbyICE(I,J) ))
1321 ELSE
1322 AREA(I,J,bi,bj)=0. _d 0
1323 ENDIF
1324
1325 #ifdef ALLOW_DIAGNOSTICS
1326 DIAGarray(I,J) = d_AREAbyICE(I,J) + d_AREAbyATM(I,J)
1327 & + d_AREAbyOCN(I,J)
1328 #endif
1329
1330 ENDDO
1331 ENDDO
1332
1333 #ifdef ALLOW_DIAGNOSTICS
1334 IF ( useDiagnostics ) THEN
1335 IF ( DIAGNOSTICS_IS_ON('SIdA ',myThid) ) THEN
1336 CALL DIAGNOSTICS_FILL(DIAGarray,
1337 & 'SIdA ',0,1,3,bi,bj,myThid)
1338 ENDIF
1339 IF ( DIAGNOSTICS_IS_ON('SIdAbATO',myThid) ) THEN
1340 CALL DIAGNOSTICS_FILL(d_AREAbyATM,
1341 & 'SIdAbATO',0,1,3,bi,bj,myThid)
1342 ENDIF
1343 IF ( DIAGNOSTICS_IS_ON('SIdAbATC',myThid) ) THEN
1344 CALL DIAGNOSTICS_FILL(d_AREAbyICE,
1345 & 'SIdAbATC',0,1,3,bi,bj,myThid)
1346 ENDIF
1347 IF ( DIAGNOSTICS_IS_ON('SIdAbOCN',myThid) ) THEN
1348 CALL DIAGNOSTICS_FILL(d_AREAbyOCN,
1349 & 'SIdAbOCN',0,1,3,bi,bj,myThid)
1350 ENDIF
1351 ENDIF
1352 #endif
1353
1354 #ifdef ALLOW_AUTODIFF_TAMC
1355 #ifdef SEAICE_MODIFY_GROWTH_ADJ
1356 Cgf 'bulk' linearization of area=f(HEFF)
1357 if ( SEAICEadjMODE.GE.1 ) then
1358 DO J=1,sNy
1359 DO I=1,sNx
1360 C AREA(I,J,bi,bj) = 0.1 _d 0 * HEFF(I,J,bi,bj)
1361 AREA(I,J,bi,bj) = AREApreTH(I,J) + 0.1 _d 0 *
1362 & ( HEFF(I,J,bi,bj) - HEFFpreTH(I,J) )
1363 ENDDO
1364 ENDDO
1365 endif
1366 #endif
1367 #endif
1368
1369 #ifdef SEAICE_GROWTH_LEGACY
1370 #ifdef ALLOW_DIAGNOSTICS
1371 IF ( useDiagnostics ) THEN
1372 IF ( DIAGNOSTICS_IS_ON('SIfice ',myThid) ) THEN
1373 CALL DIAGNOSTICS_FILL(DIAGarray,'SIfice ',0,1,3,bi,bj,myThid)
1374 ENDIF
1375 ENDIF
1376 #endif
1377 #endif /* SEAICE_GROWTH_LEGACY */
1378
1379
1380 C ===================================================================
1381 C =============PART 5: determine ice salinity increments=============
1382 C ===================================================================
1383
1384 #ifndef SEAICE_VARIABLE_SALINITY
1385 # ifdef ALLOW_AUTODIFF_TAMC
1386 # ifdef ALLOW_SALT_PLUME
1387 CADJ STORE d_HEFFbyNEG = comlev1_bibj,key=iicekey,byte=isbyte
1388 CADJ STORE d_HEFFbyOCNonICE = comlev1_bibj,key=iicekey,byte=isbyte
1389 CADJ STORE d_HEFFbyATMonOCN = comlev1_bibj,key=iicekey,byte=isbyte
1390 CADJ STORE d_HEFFbyATMonOCN_open = comlev1_bibj,key=iicekey,byte=isbyte
1391 CADJ STORE d_HEFFbyATMonOCN_cover = comlev1_bibj,key=iicekey,byte=isbyte
1392 CADJ STORE d_HEFFbyFLOODING = comlev1_bibj,key=iicekey,byte=isbyte
1393 CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
1394 CADJ & key = iicekey, byte = isbyte
1395 # endif /* ALLOW_SALT_PLUME */
1396 # endif /* ALLOW_AUTODIFF_TAMC */
1397 DO J=1,sNy
1398 DO I=1,sNx
1399 Cgf note: flooding does count negatively
1400 tmpscal1 = d_HEFFbyNEG(I,J) + d_HEFFbyOCNonICE(I,J) +
1401 & d_HEFFbyATMonOCN(I,J) - d_HEFFbyFLOODING(I,J)
1402 tmpscal2 = tmpscal1 * SIsal0 * HEFFM(I,J,bi,bj)
1403 & /SEAICE_deltaTtherm * SEAICE_rhoIce
1404 saltFlux(I,J,bi,bj) = tmpscal2
1405 #ifdef ALLOW_SALT_PLUME
1406 tmpscal3 = tmpscal1*salt(I,j,kSurface,bi,bj)*HEFFM(I,J,bi,bj)
1407 & /SEAICE_deltaTtherm * SEAICE_rhoIce
1408 saltPlumeFlux(I,J,bi,bj) = MAX( tmpscal3-tmpscal2 , 0. _d 0)
1409 #endif /* ALLOW_SALT_PLUME */
1410 ENDDO
1411 ENDDO
1412 #endif
1413
1414 #ifdef ALLOW_ATM_TEMP
1415 #ifdef SEAICE_VARIABLE_SALINITY
1416
1417 #ifdef SEAICE_GROWTH_LEGACY
1418 # ifdef ALLOW_AUTODIFF_TAMC
1419 CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1420 # endif /* ALLOW_AUTODIFF_TAMC */
1421 DO J=1,sNy
1422 DO I=1,sNx
1423 C set HSALT = 0 if HSALT < 0 and compute salt to remove from ocean
1424 IF ( HSALT(I,J,bi,bj) .LT. 0.0 ) THEN
1425 saltFluxAdjust(I,J) = - HEFFM(I,J,bi,bj) *
1426 & HSALT(I,J,bi,bj) / SEAICE_deltaTtherm
1427 HSALT(I,J,bi,bj) = 0.0 _d 0
1428 ENDIF
1429 ENDDO
1430 ENDDO
1431 #endif /* SEAICE_GROWTH_LEGACY */
1432
1433 #ifdef ALLOW_AUTODIFF_TAMC
1434 CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1435 #endif /* ALLOW_AUTODIFF_TAMC */
1436
1437 DO J=1,sNy
1438 DO I=1,sNx
1439 C sum up the terms that affect the salt content of the ice pack
1440 tmpscal1=d_HEFFbyOCNonICE(I,J)+d_HEFFbyATMonOCN(I,J)
1441
1442 C recompute HEFF before thermodyncamic updates (which is not AREApreTH in legacy code)
1443 tmpscal2=HEFF(I,J,bi,bj)-tmpscal1-d_HEFFbyFLOODING(I,J)
1444 C tmpscal1 > 0 : m of sea ice that is created
1445 IF ( tmpscal1 .GE. 0.0 ) THEN
1446 saltFlux(I,J,bi,bj) =
1447 & HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm
1448 & *SIsalFRAC*salt(I,j,kSurface,bi,bj)
1449 & *tmpscal1*SEAICE_rhoIce
1450 #ifdef ALLOW_SALT_PLUME
1451 C saltPlumeFlux is defined only during freezing:
1452 saltPlumeFlux(I,J,bi,bj)=
1453 & HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm
1454 & *(1-SIsalFRAC)*salt(I,j,kSurface,bi,bj)
1455 & *tmpscal1*SEAICE_rhoIce
1456 C if SaltPlumeSouthernOcean=.FALSE. turn off salt plume in Southern Ocean
1457 IF ( .NOT. SaltPlumeSouthernOcean ) THEN
1458 IF ( YC(I,J,bi,bj) .LT. 0.0 _d 0 )
1459 & saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
1460 ENDIF
1461 #endif /* ALLOW_SALT_PLUME */
1462
1463 C tmpscal1 < 0 : m of sea ice that is melted
1464 ELSE
1465 saltFlux(I,J,bi,bj) =
1466 & HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm
1467 & *HSALT(I,J,bi,bj)
1468 & *tmpscal1/tmpscal2
1469 #ifdef ALLOW_SALT_PLUME
1470 saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
1471 #endif /* ALLOW_SALT_PLUME */
1472 ENDIF
1473 C update HSALT based on surface saltFlux
1474 HSALT(I,J,bi,bj) = HSALT(I,J,bi,bj) +
1475 & saltFlux(I,J,bi,bj) * SEAICE_deltaTtherm
1476 saltFlux(I,J,bi,bj) =
1477 & saltFlux(I,J,bi,bj) + saltFluxAdjust(I,J)
1478 #ifdef SEAICE_GROWTH_LEGACY
1479 C set HSALT = 0 if HEFF = 0 and compute salt to dump into ocean
1480 IF ( HEFF(I,J,bi,bj) .EQ. 0.0 ) THEN
1481 saltFlux(I,J,bi,bj) = saltFlux(I,J,bi,bj) -
1482 & HEFFM(I,J,bi,bj) * HSALT(I,J,bi,bj) /
1483 & SEAICE_deltaTtherm
1484 HSALT(I,J,bi,bj) = 0.0 _d 0
1485 #ifdef ALLOW_SALT_PLUME
1486 saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
1487 #endif /* ALLOW_SALT_PLUME */
1488 ENDIF
1489 #endif /* SEAICE_GROWTH_LEGACY */
1490 ENDDO
1491 ENDDO
1492 #endif /* SEAICE_VARIABLE_SALINITY */
1493 #endif /* ALLOW_ATM_TEMP */
1494
1495
1496 C =======================================================================
1497 C =====LEGACY PART 5.5: treat pathological cases, then do flooding ======
1498 C =======================================================================
1499
1500 #ifdef SEAICE_GROWTH_LEGACY
1501
1502 C treat values of ice cover fraction oustide
1503 C the [0 1] range, and other such issues.
1504 C ===========================================
1505
1506 Cgf note: this part cannot be heat and water conserving
1507
1508 #ifdef ALLOW_AUTODIFF_TAMC
1509 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
1510 CADJ & key = iicekey, byte = isbyte
1511 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
1512 CADJ & key = iicekey, byte = isbyte
1513 #endif /* ALLOW_AUTODIFF_TAMC */
1514 DO J=1,sNy
1515 DO I=1,sNx
1516 C NOW SET AREA(I,J,bi,bj)=0 WHERE NO ICE IS
1517 AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj)
1518 & ,HEFF(I,J,bi,bj)/.0001 _d 0)
1519 ENDDO
1520 ENDDO
1521
1522 #ifdef ALLOW_AUTODIFF_TAMC
1523 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
1524 CADJ & key = iicekey, byte = isbyte
1525 #endif /* ALLOW_AUTODIFF_TAMC */
1526 DO J=1,sNy
1527 DO I=1,sNx
1528 C NOW TRUNCATE AREA
1529 AREA(I,J,bi,bj)=MIN(ONE,AREA(I,J,bi,bj))
1530 ENDDO
1531 ENDDO
1532
1533 #ifdef ALLOW_AUTODIFF_TAMC
1534 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
1535 CADJ & key = iicekey, byte = isbyte
1536 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
1537 CADJ & key = iicekey, byte = isbyte
1538 #endif /* ALLOW_AUTODIFF_TAMC */
1539 DO J=1,sNy
1540 DO I=1,sNx
1541 AREA(I,J,bi,bj) = MAX(ZERO,AREA(I,J,bi,bj))
1542 HSNOW(I,J,bi,bj) = MAX(ZERO,HSNOW(I,J,bi,bj))
1543 AREA(I,J,bi,bj) = AREA(I,J,bi,bj)*HEFFM(I,J,bi,bj)
1544 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)*HEFFM(I,J,bi,bj)
1545 #ifdef SEAICE_CAP_HEFF
1546 C This is not energy conserving, but at least it conserves fresh water
1547 tmpscal0 = -MAX(HEFF(I,J,bi,bj)-MAX_HEFF,0. _d 0)
1548 d_HEFFbyNeg(I,J) = d_HEFFbyNeg(I,J) + tmpscal0
1549 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal0
1550 #endif /* SEAICE_CAP_HEFF */
1551 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj)
1552 ENDDO
1553 ENDDO
1554
1555 #ifdef ALLOW_DIAGNOSTICS
1556 IF ( useDiagnostics ) THEN
1557 IF ( DIAGNOSTICS_IS_ON('SIthdgrh',myThid) ) THEN
1558 DO J=1,sNy
1559 DO I=1,sNx
1560 tmparr1(I,J) = (HEFF(I,J,bi,bj)-HEFFpreTH(I,J))
1561 & /SEAICE_deltaTtherm
1562 ENDDO
1563 ENDDO
1564 CALL DIAGNOSTICS_FILL(tmparr1,'SIthdgrh',0,1,3,bi,bj,myThid)
1565 ENDIF
1566 ENDIF
1567 #endif /* ALLOW_DIAGNOSTICS */
1568
1569 C convert snow to ice if submerged.
1570 C =================================
1571
1572 #ifdef ALLOW_SEAICE_FLOODING
1573 IF ( SEAICEuseFlooding ) THEN
1574 DO J=1,sNy
1575 DO I=1,sNx
1576 hDraft = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
1577 & +HEFF(I,J,bi,bj)*SEAICE_rhoIce)/rhoConst
1578 tmpscal1 = MAX( 0. _d 0, hDraft - HEFF(I,J,bi,bj))
1579 d_HEFFbyFLOODING(I,J)=tmpscal1
1580 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)+d_HEFFbyFLOODING(I,J)
1581 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-
1582 & d_HEFFbyFLOODING(I,J)*ICE2SNOW
1583 ENDDO
1584 ENDDO
1585 #ifdef ALLOW_DIAGNOSTICS
1586 IF ( useDiagnostics ) THEN
1587 IF ( DIAGNOSTICS_IS_ON('SIsnwice',myThid) ) THEN
1588 DO J=1,sNy
1589 DO I=1,sNx
1590 tmparr1(I,J) = d_HEFFbyFLOODING(I,J)/SEAICE_deltaTtherm
1591 ENDDO
1592 ENDDO
1593 CALL DIAGNOSTICS_FILL(tmparr1,'SIsnwice',0,1,3,bi,bj,myThid)
1594 ENDIF
1595 ENDIF
1596 #endif /* ALLOW_DIAGNOSTICS */
1597 ENDIF
1598 #endif /* ALLOW_SEAICE_FLOODING */
1599
1600 #endif /* SEAICE_GROWTH_LEGACY */
1601
1602 #ifdef ALLOW_SITRACER
1603 DO J=1,sNy
1604 DO I=1,sNx
1605 c needs to be here to allow use also with LEGACY branch
1606 SItrHEFF(I,J,bi,bj,5)=HEFF(I,J,bi,bj)
1607 ENDDO
1608 ENDDO
1609 #endif
1610
1611 C ===================================================================
1612 C ===============PART 6: determine ice age increments================
1613 C ===================================================================
1614
1615 #ifdef SEAICE_AGE
1616 C Sources and sinks for sea ice age:
1617 C assume that a) freezing: new ice fraction forms with zero age
1618 C b) melting: ice fraction vanishes with current age
1619 DO J=1,sNy
1620 DO I=1,sNx
1621 C-- increase age as time passes
1622 IceAgeTr(i,j,bi,bj,1)=IceAgeTr(i,j,bi,bj,1)
1623 & +SEAICE_deltaTtherm*AREApreTH(i,j)
1624 IceAgeTr(i,j,bi,bj,2)=IceAgeTr(i,j,bi,bj,2)
1625 & +SEAICE_deltaTtherm*HEFFpreTH(i,j)
1626 C-- account for ice melt
1627 IF (AREApreTH(i,j).GT.AREA(i,j,bi,bj)) THEN
1628 IceAgeTr(i,j,bi,bj,1)=IceAgeTr(i,j,bi,bj,1)
1629 & *AREA(i,j,bi,bj)/AREApreTH(i,j)
1630 ENDIF
1631 IF (HEFFpreTH(i,j).GT.HEFF(i,j,bi,bj)) THEN
1632 IceAgeTr(i,j,bi,bj,2)=IceAgeTr(i,j,bi,bj,2)
1633 & *HEFF(i,j,bi,bj)/HEFFpreTH(i,j)
1634 ENDIF
1635 ENDDO
1636 ENDDO
1637
1638 C #ifdef ALLOW_DIAGNOSTICS
1639 C IF ( useDiagnostics ) THEN
1640 C IF ( DIAGNOSTICS_IS_ON('SIage03 ',myThid) ) THEN
1641 C DO J=1,sNy
1642 C DO I=1,sNx
1643 C DIAGarray(I,J) = IceAgeTr(i,j,bi,bj,1)/AREA(i,j,bi,bj)
1644 C ENDDO
1645 C ENDDO
1646 C CALL DIAGNOSTICS_FILL(DIAGarray,'SIage03 ',0,1,3,bi,bj,myThid)
1647 C ENDIF
1648 C IF ( DIAGNOSTICS_IS_ON('SIage04 ',myThid) ) THEN
1649 C DO J=1,sNy
1650 C DO I=1,sNx
1651 C DIAGarray(I,J) = IceAgeTr(i,j,bi,bj,2)/HEFF(i,j,bi,bj)
1652 C ENDDO
1653 C ENDDO
1654 C CALL DIAGNOSTICS_FILL(DIAGarray,'SIage04 ',0,1,3,bi,bj,myThid)
1655 C ENDIF
1656 C ENDIF
1657 C #endif /* ALLOW_DIAGNOSTICS */
1658
1659 #endif /* SEAICE_AGE */
1660
1661
1662 C ===================================================================
1663 C ==============PART 7: determine ocean model forcing================
1664 C ===================================================================
1665
1666 C compute net heat flux leaving/entering the ocean,
1667 C accounting for the part used in melt/freeze processes
1668 C =====================================================
1669
1670 DO J=1,sNy
1671 DO I=1,sNx
1672 QNET(I,J,bi,bj) = r_QbyATM_cover(I,J) + r_QbyATM_open(I,J)
1673 #ifndef SEAICE_GROWTH_LEGACY
1674 C in principle a_QSWbyATM_cover should always be included here, however
1675 C for backward compatibility it is left out of the LEGACY branch
1676 & + a_QSWbyATM_cover(I,J)
1677 #endif /* SEAICE_GROWTH_LEGACY */
1678 & - ( d_HEFFbyOCNonICE(I,J) +
1679 & d_HSNWbyOCNonSNW(I,J)/ICE2SNOW +
1680 & d_HEFFbyNEG(I,J) +
1681 & d_HSNWbyNEG(I,J)/ICE2SNOW )
1682 & * maskC(I,J,kSurface,bi,bj)
1683 QSW(I,J,bi,bj) = a_QSWbyATM_cover(I,J) + a_QSWbyATM_open(I,J)
1684 ENDDO
1685 ENDDO
1686
1687 #ifdef ALLOW_DIAGNOSTICS
1688 IF ( useDiagnostics ) THEN
1689 IF ( DIAGNOSTICS_IS_ON('SIqneto ',myThid) ) THEN
1690 DO J=1,sNy
1691 DO I=1,sNx
1692 DIAGarray(I,J) = r_QbyATM_open(I,J) * convertHI2Q
1693 ENDDO
1694 ENDDO
1695 CALL DIAGNOSTICS_FILL(DIAGarray,'SIqneto ',0,1,3,bi,bj,myThid)
1696 ENDIF
1697 IF ( DIAGNOSTICS_IS_ON('SIqneti ',myThid) ) THEN
1698 DO J=1,sNy
1699 DO I=1,sNx
1700 DIAGarray(I,J) = r_QbyATM_cover(I,J) * convertHI2Q
1701 ENDDO
1702 ENDDO
1703 CALL DIAGNOSTICS_FILL(DIAGarray,'SIqneti ',0,1,3,bi,bj,myThid)
1704 ENDIF
1705 ENDIF
1706 #endif /* ALLOW_DIAGNOSTICS */
1707
1708 C switch heat fluxes from 'effective' ice meters to W/m2
1709 C ======================================================
1710
1711 DO J=1,sNy
1712 DO I=1,sNx
1713 QNET(I,J,bi,bj) = QNET(I,J,bi,bj)*convertHI2Q
1714 QSW(I,J,bi,bj) = QSW(I,J,bi,bj)*convertHI2Q
1715 ENDDO
1716 ENDDO
1717
1718 C compute net fresh water flux leaving/entering
1719 C the ocean, accounting for fresh/salt water stocks.
1720 C ==================================================
1721
1722 #ifdef ALLOW_ATM_TEMP
1723 DO J=1,sNy
1724 DO I=1,sNx
1725 tmpscal1= d_HSNWbyATMonSNW(I,J)/ICE2SNOW
1726 & +d_HFRWbyRAIN(I,J)
1727 & +d_HSNWbyOCNonSNW(I,J)/ICE2SNOW
1728 & +d_HEFFbyOCNonICE(I,J)
1729 & +d_HEFFbyATMonOCN(I,J)
1730 & +d_HEFFbyNEG(I,J)
1731 & +d_HSNWbyNEG(I,J)/ICE2SNOW
1732 #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
1733 & +a_FWbySublim(I,J)
1734 #endif /* SEAICE_ADD_SUBLIMATION_TO_FWBUDGET */
1735 EmPmR(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
1736 & ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )
1737 & * ( ONE - AREApreTH(I,J) )
1738 #ifdef ALLOW_RUNOFF
1739 & - RUNOFF(I,J,bi,bj)
1740 #endif /* ALLOW_RUNOFF */
1741 & + tmpscal1*convertHI2PRECIP
1742 & )*rhoConstFresh
1743 ENDDO
1744 ENDDO
1745
1746 #ifdef ALLOW_MEAN_SFLUX_COST_CONTRIBUTION
1747 DO J=1,sNy
1748 DO I=1,sNx
1749 frWtrAtm(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
1750 & PRECIP(I,J,bi,bj)
1751 & - EVAP(I,J,bi,bj)
1752 & *( ONE - AREApreTH(I,J) )
1753 & + RUNOFF(I,J,bi,bj)
1754 & )*rhoConstFresh
1755 ENDDO
1756 ENDDO
1757 #endif
1758 #endif /* ALLOW_ATM_TEMP */
1759
1760 #ifdef SEAICE_DEBUG
1761 CALL PLOT_FIELD_XYRL( QSW,'Current QSW ', myIter, myThid )
1762 CALL PLOT_FIELD_XYRL( QNET,'Current QNET ', myIter, myThid )
1763 CALL PLOT_FIELD_XYRL( EmPmR,'Current EmPmR ', myIter, myThid )
1764 #endif /* SEAICE_DEBUG */
1765
1766 C Sea Ice Load on the sea surface.
1767 C =================================
1768
1769 #ifdef ALLOW_AUTODIFF_TAMC
1770 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1771 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1772 #endif /* ALLOW_AUTODIFF_TAMC */
1773
1774 IF ( useRealFreshWaterFlux ) THEN
1775 DO J=1,sNy
1776 DO I=1,sNx
1777 #ifdef SEAICE_CAP_ICELOAD
1778 tmpscal1 = HEFF(I,J,bi,bj)*SEAICE_rhoIce
1779 & + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
1780 tmpscal2 = min(tmpscal1,heffTooHeavy*rhoConst)
1781 #else
1782 tmpscal2 = HEFF(I,J,bi,bj)*SEAICE_rhoIce
1783 & + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
1784 #endif
1785 sIceLoad(i,j,bi,bj) = tmpscal2
1786 ENDDO
1787 ENDDO
1788 ENDIF
1789
1790 C close bi,bj loops
1791 ENDDO
1792 ENDDO
1793
1794 RETURN
1795 END
  ViewVC Help
Powered by ViewVC 1.1.22