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