/[MITgcm]/MITgcm_contrib/torge/itd/code/seaice_growth.F
ViewVC logotype

Annotation of /MITgcm_contrib/torge/itd/code/seaice_growth.F

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

Revision 1.9 - (hide annotations) (download)
Mon Oct 22 16:36:45 2012 UTC (12 years, 9 months ago) by torge
Branch: MAIN
Changes since 1.8: +76 -153 lines 
wrap all msgBuf write statements in SEAICE_DEBUG preprocessor option;
in seaice_model: add a constant divergence rate to mimic dynamics for the 1D verification experiment (not a permanent change)
1 torge 1.8 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_growth.F,v 1.162 2012/03/15 03:07:31 jmc Exp $
2 dimitri 1.1 C $Name: $
3    
4     #include "SEAICE_OPTIONS.h"
5     #ifdef ALLOW_EXF
6     # include "EXF_OPTIONS.h"
7     #endif
8    
9     CBOP
10     C !ROUTINE: SEAICE_GROWTH
11     C !INTERFACE:
12     SUBROUTINE SEAICE_GROWTH( myTime, myIter, myThid )
13     C !DESCRIPTION: \bv
14     C *==========================================================*
15     C | SUBROUTINE seaice_growth
16     C | o Updata ice thickness and snow depth
17     C *==========================================================*
18     C \ev
19    
20     C !USES:
21     IMPLICIT NONE
22     C === Global variables ===
23     #include "SIZE.h"
24     #include "EEPARAMS.h"
25     #include "PARAMS.h"
26     #include "DYNVARS.h"
27     #include "GRID.h"
28     #include "FFIELDS.h"
29     #include "SEAICE_SIZE.h"
30     #include "SEAICE_PARAMS.h"
31     #include "SEAICE.h"
32     #include "SEAICE_TRACER.h"
33     #ifdef ALLOW_EXF
34     # include "EXF_PARAM.h"
35     # include "EXF_FIELDS.h"
36     #endif
37     #ifdef ALLOW_SALT_PLUME
38     # include "SALT_PLUME.h"
39     #endif
40     #ifdef ALLOW_AUTODIFF_TAMC
41     # include "tamc.h"
42     #endif
43    
44     C !INPUT/OUTPUT PARAMETERS:
45     C === Routine arguments ===
46     C myTime :: Simulation time
47     C myIter :: Simulation timestep number
48     C myThid :: Thread no. that called this routine.
49     _RL myTime
50     INTEGER myIter, myThid
51     CEOP
52    
53     C !FUNCTIONS:
54     #ifdef ALLOW_DIAGNOSTICS
55     LOGICAL DIAGNOSTICS_IS_ON
56     EXTERNAL DIAGNOSTICS_IS_ON
57     #endif
58    
59     C !LOCAL VARIABLES:
60     C === Local variables ===
61 torge 1.9 #ifdef SEAICE_DEBUG
62 torge 1.3 c ToM<<< debug seaice_growth
63     C msgBuf :: Informational/error message buffer
64     CHARACTER*(MAX_LEN_MBUF) msgBuf
65     c ToM>>>
66 torge 1.9 #endif
67 dimitri 1.1 C
68     C unit/sign convention:
69     C Within the thermodynamic computation all stocks, except HSNOW,
70     C are in 'effective ice meters' units, and >0 implies more ice.
71     C This holds for stocks due to ocean and atmosphere heat,
72     C at the outset of 'PART 2: determine heat fluxes/stocks'
73     C and until 'PART 7: determine ocean model forcing'
74     C This strategy minimizes the need for multiplications/divisions
75     C by ice fraction, heat capacity, etc. The only conversions that
76     C occurs are for the HSNOW (in effective snow meters) and
77     C PRECIP (fresh water m/s).
78     C
79     C HEFF is effective Hice thickness (m3/m2)
80     C HSNOW is Heffective snow thickness (m3/m2)
81     C HSALT is Heffective salt content (g/m2)
82     C AREA is the seaice cover fraction (0<=AREA<=1)
83     C Q denotes heat stocks -- converted to ice stocks (m3/m2) early on
84     C
85     C For all other stocks/increments, such as d_HEFFbyATMonOCN
86     C or a_QbyATM_cover, the naming convention is as follows:
87     C The prefix 'a_' means available, the prefix 'd_' means delta
88     C (i.e. increment), and the prefix 'r_' means residual.
89     C The suffix '_cover' denotes a value for the ice covered fraction
90     C of the grid cell, whereas '_open' is for the open water fraction.
91     C The main part of the name states what ice/snow stock is concerned
92     C (e.g. QbyATM or HEFF), and how it is affected (e.g. d_HEFFbyATMonOCN
93     C is the increment of HEFF due to the ATMosphere extracting heat from the
94     C OCeaN surface, or providing heat to the OCeaN surface).
95    
96     C i,j,bi,bj :: Loop counters
97     INTEGER i, j, bi, bj
98     C number of surface interface layer
99     INTEGER kSurface
100     C constants
101     _RL tempFrz, ICE2SNOW, SNOW2ICE
102     _RL QI, QS, recip_QI
103    
104     C-- TmixLoc :: ocean surface/mixed-layer temperature (in K)
105     _RL TmixLoc (1:sNx,1:sNy)
106    
107     C a_QbyATM_cover :: available heat (in W/m^2) due to the interaction of
108     C the atmosphere and the ocean surface - for ice covered water
109     C a_QbyATM_open :: same but for open water
110     C r_QbyATM_cover :: residual of a_QbyATM_cover after freezing/melting processes
111     C r_QbyATM_open :: same but for open water
112     _RL a_QbyATM_cover (1:sNx,1:sNy)
113     _RL a_QbyATM_open (1:sNx,1:sNy)
114     _RL r_QbyATM_cover (1:sNx,1:sNy)
115     _RL r_QbyATM_open (1:sNx,1:sNy)
116     C a_QSWbyATM_open - short wave heat flux over ocean in W/m^2
117     C a_QSWbyATM_cover - short wave heat flux under ice in W/m^2
118     _RL a_QSWbyATM_open (1:sNx,1:sNy)
119     _RL a_QSWbyATM_cover (1:sNx,1:sNy)
120 torge 1.7 C a_QbyOCN :: available heat (in W/m^2) due to the
121 dimitri 1.1 C interaction of the ice pack and the ocean surface
122     C r_QbyOCN :: residual of a_QbyOCN after freezing/melting
123     C processes have been accounted for
124     _RL a_QbyOCN (1:sNx,1:sNy)
125     _RL r_QbyOCN (1:sNx,1:sNy)
126    
127     C conversion factors to go from Q (W/m2) to HEFF (ice meters)
128     _RL convertQ2HI, convertHI2Q
129     C conversion factors to go from precip (m/s) unit to HEFF (ice meters)
130     _RL convertPRECIP2HI, convertHI2PRECIP
131    
132     #ifdef ALLOW_DIAGNOSTICS
133     C ICE/SNOW stocks tendencies associated with the various melt/freeze processes
134     _RL d_AREAbyATM (1:sNx,1:sNy)
135     _RL d_AREAbyOCN (1:sNx,1:sNy)
136     _RL d_AREAbyICE (1:sNx,1:sNy)
137     #endif
138    
139     #ifdef SEAICE_ALLOW_AREA_RELAXATION
140     C ICE/SNOW stocks tendency associated with relaxation towards observation
141     _RL d_AREAbyRLX (1:sNx,1:sNy)
142     c The change of mean ice thickness due to relaxation
143     _RL d_HEFFbyRLX (1:sNx,1:sNy)
144     #endif
145    
146     c The change of mean ice thickness due to out-of-bounds values following
147 dimitri 1.2 c sea ice dynamics
148 dimitri 1.1 _RL d_HEFFbyNEG (1:sNx,1:sNy)
149    
150     c The change of mean ice thickness due to turbulent ocean-sea ice heat fluxes
151     _RL d_HEFFbyOCNonICE (1:sNx,1:sNy)
152    
153     c The sum of mean ice thickness increments due to atmospheric fluxes over the open water
154     c fraction and ice-covered fractions of the grid cell
155     _RL d_HEFFbyATMonOCN (1:sNx,1:sNy)
156     c The change of mean ice thickness due to flooding by snow
157     _RL d_HEFFbyFLOODING (1:sNx,1:sNy)
158    
159     c The mean ice thickness increments due to atmospheric fluxes over the open water
160     c fraction and ice-covered fractions of the grid cell, respectively
161     _RL d_HEFFbyATMonOCN_open(1:sNx,1:sNy)
162     _RL d_HEFFbyATMonOCN_cover(1:sNx,1:sNy)
163    
164     _RL d_HSNWbyNEG (1:sNx,1:sNy)
165     _RL d_HSNWbyATMonSNW (1:sNx,1:sNy)
166     _RL d_HSNWbyOCNonSNW (1:sNx,1:sNy)
167     _RL d_HSNWbyRAIN (1:sNx,1:sNy)
168    
169     _RL d_HFRWbyRAIN (1:sNx,1:sNy)
170     C
171     C a_FWbySublim :: fresh water flux implied by latent heat of
172     C sublimation to atmosphere, same sign convention
173     C as EVAP (positive upward)
174     _RL a_FWbySublim (1:sNx,1:sNy)
175     _RL r_FWbySublim (1:sNx,1:sNy)
176     _RL d_HEFFbySublim (1:sNx,1:sNy)
177     _RL d_HSNWbySublim (1:sNx,1:sNy)
178    
179 torge 1.5 #ifdef SEAICE_CAP_SUBLIM
180 dimitri 1.1 C The latent heat flux which will sublimate all snow and ice
181     C over one time step
182 torge 1.6 _RL latentHeatFluxMax (1:sNx,1:sNy)
183 torge 1.5 _RL latentHeatFluxMaxMult (1:sNx,1:sNy,MULTDIM)
184 dimitri 1.1 #endif
185    
186     C actual ice thickness (with upper and lower limit)
187     _RL heffActual (1:sNx,1:sNy)
188     C actual snow thickness
189     _RL hsnowActual (1:sNx,1:sNy)
190     C actual ice thickness (with lower limit only) Reciprocal
191     _RL recip_heffActual (1:sNx,1:sNy)
192     C local value (=1/HO or 1/HO_south)
193     _RL recip_HO
194     C local value (=1/ice thickness)
195     _RL recip_HH
196    
197     C AREA_PRE :: hold sea-ice fraction field before any seaice-thermo update
198     _RL AREApreTH (1:sNx,1:sNy)
199     _RL HEFFpreTH (1:sNx,1:sNy)
200     _RL HSNWpreTH (1:sNx,1:sNy)
201 dimitri 1.2 #ifdef SEAICE_ITD
202     _RL AREAITDpreTH (1:sNx,1:sNy,1:nITD)
203     _RL HEFFITDpreTH (1:sNx,1:sNy,1:nITD)
204     _RL HSNWITDpreTH (1:sNx,1:sNy,1:nITD)
205     _RL areaFracFactor (1:sNx,1:sNy,1:nITD)
206 torge 1.8 _RL leadIceThickMin
207 dimitri 1.2 #endif
208 dimitri 1.1
209     C wind speed
210     _RL UG (1:sNx,1:sNy)
211     #ifdef ALLOW_ATM_WIND
212     _RL SPEED_SQ
213     #endif
214    
215     C Regularization values squared
216     _RL area_reg_sq, hice_reg_sq
217    
218     C pathological cases thresholds
219     _RL heffTooHeavy
220    
221     _RL lhSublim
222    
223     C temporary variables available for the various computations
224     _RL tmpscal0, tmpscal1, tmpscal2, tmpscal3, tmpscal4
225     _RL tmparr1 (1:sNx,1:sNy)
226    
227     #ifdef SEAICE_VARIABLE_SALINITY
228     _RL saltFluxAdjust (1:sNx,1:sNy)
229     #endif
230    
231     INTEGER ilockey
232 torge 1.3 INTEGER it
233 dimitri 1.1 _RL pFac
234     _RL ticeInMult (1:sNx,1:sNy,MULTDIM)
235     _RL ticeOutMult (1:sNx,1:sNy,MULTDIM)
236     _RL heffActualMult (1:sNx,1:sNy,MULTDIM)
237 dimitri 1.2 #ifdef SEAICE_ITD
238     _RL hsnowActualMult (1:sNx,1:sNy,MULTDIM)
239     _RL recip_heffActualMult(1:sNx,1:sNy,MULTDIM)
240     #endif
241 dimitri 1.1 _RL a_QbyATMmult_cover (1:sNx,1:sNy,MULTDIM)
242     _RL a_QSWbyATMmult_cover(1:sNx,1:sNy,MULTDIM)
243     _RL a_FWbySublimMult (1:sNx,1:sNy,MULTDIM)
244 dimitri 1.2 #ifdef SEAICE_ITD
245     _RL r_QbyATMmult_cover (1:sNx,1:sNy,MULTDIM)
246     _RL r_FWbySublimMult (1:sNx,1:sNy,MULTDIM)
247     #endif
248 dimitri 1.1 C Helper variables: reciprocal of some constants
249     _RL recip_multDim
250     _RL recip_deltaTtherm
251     _RL recip_rhoIce
252    
253     C Factor by which we increase the upper ocean friction velocity (u*) when
254     C ice is absent in a grid cell (dimensionless)
255     _RL MixedLayerTurbulenceFactor
256    
257     #ifdef ALLOW_SITRACER
258     INTEGER iTr
259     CHARACTER*8 diagName
260     #endif
261     #ifdef ALLOW_DIAGNOSTICS
262     c Helper variables for diagnostics
263     _RL DIAGarrayA (1:sNx,1:sNy)
264     _RL DIAGarrayB (1:sNx,1:sNy)
265     _RL DIAGarrayC (1:sNx,1:sNy)
266     _RL DIAGarrayD (1:sNx,1:sNy)
267     #endif
268    
269    
270     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
271    
272     C ===================================================================
273     C =================PART 0: constants and initializations=============
274     C ===================================================================
275    
276     IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
277     kSurface = Nr
278     ELSE
279     kSurface = 1
280     ENDIF
281    
282     C avoid unnecessary divisions in loops
283 torge 1.8 c#ifdef SEAICE_ITD
284     CToM this is now set by MULTDIM = nITD in SEAICE_SIZE.h
285     C (see SEAICE_SIZE.h and seaice_readparms.F)
286     c SEAICE_multDim = nITD
287     c#endif
288 dimitri 1.1 recip_multDim = SEAICE_multDim
289     recip_multDim = ONE / recip_multDim
290     C above/below: double/single precision calculation of recip_multDim
291     c recip_multDim = 1./float(SEAICE_multDim)
292     recip_deltaTtherm = ONE / SEAICE_deltaTtherm
293     recip_rhoIce = ONE / SEAICE_rhoIce
294    
295     C Cutoff for iceload
296     heffTooHeavy=drF(kSurface) / 5. _d 0
297    
298     C RATIO OF SEA ICE DENSITY to SNOW DENSITY
299     ICE2SNOW = SEAICE_rhoIce/SEAICE_rhoSnow
300     SNOW2ICE = ONE / ICE2SNOW
301    
302     C HEAT OF FUSION OF ICE (J/m^3)
303     QI = SEAICE_rhoIce*SEAICE_lhFusion
304     recip_QI = ONE / QI
305     C HEAT OF FUSION OF SNOW (J/m^3)
306     QS = SEAICE_rhoSnow*SEAICE_lhFusion
307    
308     C ICE LATENT HEAT CONSTANT
309     lhSublim = SEAICE_lhEvap + SEAICE_lhFusion
310    
311     C regularization constants
312     area_reg_sq = SEAICE_area_reg * SEAICE_area_reg
313     hice_reg_sq = SEAICE_hice_reg * SEAICE_hice_reg
314    
315     C conversion factors to go from Q (W/m2) to HEFF (ice meters)
316     convertQ2HI=SEAICE_deltaTtherm/QI
317     convertHI2Q = ONE/convertQ2HI
318     C conversion factors to go from precip (m/s) unit to HEFF (ice meters)
319     convertPRECIP2HI=SEAICE_deltaTtherm*rhoConstFresh/SEAICE_rhoIce
320     convertHI2PRECIP = ONE/convertPRECIP2HI
321    
322     DO bj=myByLo(myThid),myByHi(myThid)
323     DO bi=myBxLo(myThid),myBxHi(myThid)
324    
325     #ifdef ALLOW_AUTODIFF_TAMC
326     act1 = bi - myBxLo(myThid)
327     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
328     act2 = bj - myByLo(myThid)
329     max2 = myByHi(myThid) - myByLo(myThid) + 1
330     act3 = myThid - 1
331     max3 = nTx*nTy
332     act4 = ikey_dynamics - 1
333     iicekey = (act1 + 1) + act2*max1
334     & + act3*max1*max2
335     & + act4*max1*max2*max3
336     #endif /* ALLOW_AUTODIFF_TAMC */
337    
338    
339     C array initializations
340     C =====================
341    
342     DO J=1,sNy
343     DO I=1,sNx
344     a_QbyATM_cover (I,J) = 0.0 _d 0
345     a_QbyATM_open(I,J) = 0.0 _d 0
346     r_QbyATM_cover (I,J) = 0.0 _d 0
347     r_QbyATM_open (I,J) = 0.0 _d 0
348    
349     a_QSWbyATM_open (I,J) = 0.0 _d 0
350     a_QSWbyATM_cover (I,J) = 0.0 _d 0
351    
352     a_QbyOCN (I,J) = 0.0 _d 0
353     r_QbyOCN (I,J) = 0.0 _d 0
354    
355     #ifdef ALLOW_DIAGNOSTICS
356     d_AREAbyATM(I,J) = 0.0 _d 0
357     d_AREAbyICE(I,J) = 0.0 _d 0
358     d_AREAbyOCN(I,J) = 0.0 _d 0
359     #endif
360    
361     #ifdef SEAICE_ALLOW_AREA_RELAXATION
362     d_AREAbyRLX(I,J) = 0.0 _d 0
363     d_HEFFbyRLX(I,J) = 0.0 _d 0
364     #endif
365    
366     d_HEFFbyNEG(I,J) = 0.0 _d 0
367     d_HEFFbyOCNonICE(I,J) = 0.0 _d 0
368     d_HEFFbyATMonOCN(I,J) = 0.0 _d 0
369     d_HEFFbyFLOODING(I,J) = 0.0 _d 0
370    
371     d_HEFFbyATMonOCN_open(I,J) = 0.0 _d 0
372     d_HEFFbyATMonOCN_cover(I,J) = 0.0 _d 0
373    
374     d_HSNWbyNEG(I,J) = 0.0 _d 0
375     d_HSNWbyATMonSNW(I,J) = 0.0 _d 0
376     d_HSNWbyOCNonSNW(I,J) = 0.0 _d 0
377     d_HSNWbyRAIN(I,J) = 0.0 _d 0
378     a_FWbySublim(I,J) = 0.0 _d 0
379     r_FWbySublim(I,J) = 0.0 _d 0
380     d_HEFFbySublim(I,J) = 0.0 _d 0
381     d_HSNWbySublim(I,J) = 0.0 _d 0
382     #ifdef SEAICE_CAP_SUBLIM
383     latentHeatFluxMax(I,J) = 0.0 _d 0
384     #endif
385     c
386     d_HFRWbyRAIN(I,J) = 0.0 _d 0
387    
388     tmparr1(I,J) = 0.0 _d 0
389    
390     #ifdef SEAICE_VARIABLE_SALINITY
391     saltFluxAdjust(I,J) = 0.0 _d 0
392     #endif
393     DO IT=1,SEAICE_multDim
394     ticeInMult(I,J,IT) = 0.0 _d 0
395     ticeOutMult(I,J,IT) = 0.0 _d 0
396     a_QbyATMmult_cover(I,J,IT) = 0.0 _d 0
397     a_QSWbyATMmult_cover(I,J,IT) = 0.0 _d 0
398     a_FWbySublimMult(I,J,IT) = 0.0 _d 0
399 torge 1.6 #ifdef SEAICE_CAP_SUBLIM
400     latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0
401     #endif
402 dimitri 1.2 #ifdef SEAICE_ITD
403     r_QbyATMmult_cover (I,J,IT) = 0.0 _d 0
404     r_FWbySublimMult(I,J,IT) = 0.0 _d 0
405     #endif
406 dimitri 1.1 ENDDO
407     ENDDO
408     ENDDO
409     #if (defined (ALLOW_MEAN_SFLUX_COST_CONTRIBUTION) || defined (ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION))
410     DO J=1-oLy,sNy+oLy
411     DO I=1-oLx,sNx+oLx
412     frWtrAtm(I,J,bi,bj) = 0.0 _d 0
413     ENDDO
414     ENDDO
415     #endif
416    
417    
418     C =====================================================================
419     C ===========PART 1: treat pathological cases (post advdiff)===========
420     C =====================================================================
421    
422     #ifdef SEAICE_GROWTH_LEGACY
423    
424     DO J=1,sNy
425     DO I=1,sNx
426     HEFFpreTH(I,J)=HEFFNM1(I,J,bi,bj)
427     HSNWpreTH(I,J)=HSNOW(I,J,bi,bj)
428     AREApreTH(I,J)=AREANM1(I,J,bi,bj)
429     d_HEFFbyNEG(I,J)=0. _d 0
430     d_HSNWbyNEG(I,J)=0. _d 0
431     #ifdef ALLOW_DIAGNOSTICS
432     DIAGarrayA(I,J) = AREANM1(I,J,bi,bj)
433     DIAGarrayB(I,J) = AREANM1(I,J,bi,bj)
434     DIAGarrayC(I,J) = HEFFNM1(I,J,bi,bj)
435     DIAGarrayD(I,J) = HSNOW(I,J,bi,bj)
436     #endif
437     ENDDO
438     ENDDO
439 dimitri 1.2 #ifdef SEAICE_ITD
440 torge 1.5 DO IT=1,nITD
441 dimitri 1.2 DO J=1,sNy
442     DO I=1,sNx
443 torge 1.5 HEFFITDpreTH(I,J,IT)=HEFFITD(I,J,IT,bi,bj)
444     HSNWITDpreTH(I,J,IT)=HSNOWITD(I,J,IT,bi,bj)
445     AREAITDpreTH(I,J,IT)=AREAITD(I,J,IT,bi,bj)
446 dimitri 1.2 ENDDO
447     ENDDO
448     ENDDO
449     #endif
450 dimitri 1.1
451     #else /* SEAICE_GROWTH_LEGACY */
452    
453     #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
454     Cgf no dependency through pathological cases treatment
455     IF ( SEAICEadjMODE.EQ.0 ) THEN
456     #endif
457    
458     #ifdef SEAICE_ALLOW_AREA_RELAXATION
459     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
460     CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
461     C 0) relax sea ice concentration towards observation
462     IF (SEAICE_tauAreaObsRelax .GT. 0.) THEN
463     DO J=1,sNy
464     DO I=1,sNx
465     C d_AREAbyRLX(i,j) = 0. _d 0
466     C d_HEFFbyRLX(i,j) = 0. _d 0
467     IF ( obsSIce(I,J,bi,bj).GT.AREA(I,J,bi,bj)) THEN
468     d_AREAbyRLX(i,j) =
469     & SEAICE_deltaTtherm/SEAICE_tauAreaObsRelax
470     & * (obsSIce(I,J,bi,bj) - AREA(I,J,bi,bj))
471     ENDIF
472     IF ( obsSIce(I,J,bi,bj).GT.0. _d 0 .AND.
473     & AREA(I,J,bi,bj).EQ.0. _d 0) THEN
474     C d_HEFFbyRLX(i,j) = 1. _d 1 * siEps * d_AREAbyRLX(i,j)
475     d_HEFFbyRLX(i,j) = 1. _d 1 * siEps
476     ENDIF
477 dimitri 1.2 #ifdef SEAICE_ITD
478     AREAITD(I,J,1,bi,bj) = AREAITD(I,J,1,bi,bj)
479     & + d_AREAbyRLX(i,j)
480     HEFFITD(I,J,1,bi,bj) = HEFFITD(I,J,1,bi,bj)
481     & + d_HEFFbyRLX(i,j)
482     #endif
483 dimitri 1.1 AREA(I,J,bi,bj) = AREA(I,J,bi,bj) + d_AREAbyRLX(i,j)
484     HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + d_HEFFbyRLX(i,j)
485     ENDDO
486     ENDDO
487     ENDIF
488     #endif /* SEAICE_ALLOW_AREA_RELAXATION */
489    
490     C 1) treat the case of negative values:
491    
492     #ifdef ALLOW_AUTODIFF_TAMC
493     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
494     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
495     CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
496     #endif /* ALLOW_AUTODIFF_TAMC */
497     DO J=1,sNy
498     DO I=1,sNx
499 dimitri 1.2 #ifdef SEAICE_ITD
500 torge 1.5 DO IT=1,nITD
501 dimitri 1.2 tmpscal2=0. _d 0
502     tmpscal3=0. _d 0
503 torge 1.5 tmpscal2=MAX(-HEFFITD(I,J,IT,bi,bj),0. _d 0)
504     HEFFITD(I,J,IT,bi,bj)=HEFFITD(I,J,IT,bi,bj)+tmpscal2
505 dimitri 1.2 d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
506 torge 1.5 tmpscal3=MAX(-HSNOWITD(I,J,IT,bi,bj),0. _d 0)
507     HSNOWITD(I,J,IT,bi,bj)=HSNOWITD(I,J,IT,bi,bj)+tmpscal3
508 dimitri 1.2 d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
509 torge 1.7 AREAITD(I,J,IT,bi,bj)=MAX(AREAITD(I,J,IT,bi,bj),0. _d 0)
510 dimitri 1.2 ENDDO
511 torge 1.3 CToM AREA, HEFF, and HSNOW will be updated at end of PART 1
512     C by calling SEAICE_ITD_SUM
513 dimitri 1.2 #else
514 dimitri 1.1 d_HEFFbyNEG(I,J)=MAX(-HEFF(I,J,bi,bj),0. _d 0)
515 torge 1.5 HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J)
516 dimitri 1.2 d_HSNWbyNEG(I,J)=MAX(-HSNOW(I,J,bi,bj),0. _d 0)
517 torge 1.5 HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J)
518 dimitri 1.2 AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),0. _d 0)
519 torge 1.3 #endif
520 dimitri 1.1 ENDDO
521     ENDDO
522    
523     C 1.25) treat the case of very thin ice:
524    
525     #ifdef ALLOW_AUTODIFF_TAMC
526     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
527     #endif /* ALLOW_AUTODIFF_TAMC */
528     DO J=1,sNy
529     DO I=1,sNx
530 dimitri 1.2 #ifdef SEAICE_ITD
531 torge 1.5 DO IT=1,nITD
532 torge 1.3 #endif
533     tmpscal2=0. _d 0
534     tmpscal3=0. _d 0
535     #ifdef SEAICE_ITD
536 torge 1.5 IF (HEFFITD(I,J,IT,bi,bj).LE.siEps) THEN
537     tmpscal2=-HEFFITD(I,J,IT,bi,bj)
538     tmpscal3=-HSNOWITD(I,J,IT,bi,bj)
539     TICES(I,J,IT,bi,bj)=celsius2K
540 torge 1.3 CToM TICE will be updated at end of Part 1 together with AREA and HEFF
541 dimitri 1.2 ENDIF
542 torge 1.5 HEFFITD(I,J,IT,bi,bj) =HEFFITD(I,J,IT,bi,bj) +tmpscal2
543     HSNOWITD(I,J,IT,bi,bj)=HSNOWITD(I,J,IT,bi,bj)+tmpscal3
544 dimitri 1.2 #else
545 dimitri 1.1 IF (HEFF(I,J,bi,bj).LE.siEps) THEN
546 torge 1.3 tmpscal2=-HEFF(I,J,bi,bj)
547     tmpscal3=-HSNOW(I,J,bi,bj)
548     TICE(I,J,bi,bj)=celsius2K
549     DO IT=1,SEAICE_multDim
550     TICES(I,J,IT,bi,bj)=celsius2K
551     ENDDO
552 dimitri 1.1 ENDIF
553     HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+tmpscal2
554 torge 1.3 HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3
555     #endif
556 dimitri 1.1 d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
557     d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
558 torge 1.3 #ifdef SEAICE_ITD
559     ENDDO
560 dimitri 1.2 #endif
561 dimitri 1.1 ENDDO
562     ENDDO
563    
564     C 1.5) treat the case of area but no ice/snow:
565    
566     #ifdef ALLOW_AUTODIFF_TAMC
567     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
568     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
569     #endif /* ALLOW_AUTODIFF_TAMC */
570     DO J=1,sNy
571     DO I=1,sNx
572 torge 1.3 #ifdef SEAICE_ITD
573 torge 1.5 DO IT=1,nITD
574     IF ((HEFFITD(I,J,IT,bi,bj).EQ.0. _d 0).AND.
575     & (HSNOWITD(I,J,IT,bi,bj).EQ.0. _d 0))
576     & AREAITD(I,J,IT,bi,bj)=0. _d 0
577 torge 1.3 ENDDO
578     #else
579 dimitri 1.1 IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.
580 torge 1.3 & (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0
581 dimitri 1.2 #endif
582 dimitri 1.1 ENDDO
583     ENDDO
584    
585     C 2) treat the case of very small area:
586    
587     #ifndef DISABLE_AREA_FLOOR
588     #ifdef ALLOW_AUTODIFF_TAMC
589     CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
590     #endif /* ALLOW_AUTODIFF_TAMC */
591     DO J=1,sNy
592     DO I=1,sNx
593 torge 1.3 #ifdef SEAICE_ITD
594 torge 1.5 DO IT=1,nITD
595     IF ((HEFFITD(I,J,IT,bi,bj).GT.0).OR.
596     & (HSNOWITD(I,J,IT,bi,bj).GT.0)) THEN
597 torge 1.3 CToM SEAICE_area_floor*nITD cannot be allowed to exceed 1
598     C hence use SEAICE_area_floor devided by nITD
599     C (or install a warning in e.g. seaice_readparms.F)
600 torge 1.5 AREAITD(I,J,IT,bi,bj)=
601     & MAX(AREAITD(I,J,IT,bi,bj),SEAICE_area_floor/float(nITD))
602 torge 1.3 ENDIF
603     ENDDO
604     #else
605 dimitri 1.1 IF ((HEFF(i,j,bi,bj).GT.0).OR.(HSNOW(i,j,bi,bj).GT.0)) THEN
606     AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),SEAICE_area_floor)
607 torge 1.3 ENDIF
608 dimitri 1.2 #endif
609 dimitri 1.1 ENDDO
610     ENDDO
611     #endif /* DISABLE_AREA_FLOOR */
612    
613     C 2.5) treat case of excessive ice cover, e.g., due to ridging:
614    
615 torge 1.3 CToM for SEAICE_ITD this case is treated in SEAICE_ITD_REDIST,
616     C which is called at end of PART 1 below
617     #ifndef SEAICE_ITD
618 dimitri 1.1 #ifdef ALLOW_AUTODIFF_TAMC
619     CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
620     #endif /* ALLOW_AUTODIFF_TAMC */
621     DO J=1,sNy
622     DO I=1,sNx
623     #ifdef ALLOW_DIAGNOSTICS
624     DIAGarrayA(I,J) = AREA(I,J,bi,bj)
625     #endif
626     #ifdef ALLOW_SITRACER
627     SItrAREA(I,J,bi,bj,1)=AREA(I,J,bi,bj)
628     #endif
629     AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),SEAICE_area_max)
630 torge 1.3 ENDDO
631     ENDDO
632 torge 1.5 #endif /* notSEAICE_ITD */
633 torge 1.3
634 dimitri 1.2 #ifdef SEAICE_ITD
635 torge 1.3 CToM catch up with items 1.25 and 2.5 involving category sums AREA and HEFF
636     DO J=1,sNy
637     DO I=1,sNx
638     C TICES was changed above (item 1.25), now update TICE as ice volume
639     C weighted average of TICES
640 torge 1.5 C also compute total of AREAITD (needed for finishing item 2.5, see below)
641 torge 1.3 tmpscal1 = 0. _d 0
642     tmpscal2 = 0. _d 0
643 torge 1.5 tmpscal3 = 0. _d 0
644     DO IT=1,nITD
645     tmpscal1=tmpscal1 + TICES(I,J,IT,bi,bj)*HEFFITD(I,J,IT,bi,bj)
646     tmpscal2=tmpscal2 + HEFFITD(I,J,IT,bi,bj)
647     tmpscal3=tmpscal3 + AREAITD(I,J,IT,bi,bj)
648 dimitri 1.2 ENDDO
649 torge 1.3 TICE(I,J,bi,bj)=tmpscal1/tmpscal2
650     C lines of item 2.5 that were omitted:
651     C in 2.5 these lines are executed before "ridging" is applied to AREA
652     C hence we execute them here before SEAICE_ITD_REDIST is called
653     C although this means that AREA has not been completely regularized
654     #ifdef ALLOW_DIAGNOSTICS
655 torge 1.5 DIAGarrayA(I,J) = tmpscal3
656 torge 1.3 #endif
657     #ifdef ALLOW_SITRACER
658 torge 1.5 SItrAREA(I,J,bi,bj,1)=tmpscal3
659 dimitri 1.2 #endif
660 dimitri 1.1 ENDDO
661     ENDDO
662    
663 torge 1.3 CToM finally make sure that all categories meet their thickness limits
664     C which includes ridging as in item 2.5
665     C and update AREA, HEFF, and HSNOW
666     CALL SEAICE_ITD_REDIST(bi, bj, myTime, myIter, myThid)
667     CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)
668    
669 torge 1.9 #ifdef SEAICE_DEBUG
670 torge 1.7 c ToM<<< debug seaice_growth
671 torge 1.8 WRITE(msgBuf,'(A,7F8.4)')
672 torge 1.7 & ' SEAICE_GROWTH: Heff increments 0, HEFFITD = ',
673 torge 1.8 & HEFFITD(1,1,:,bi,bj)
674 torge 1.7 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
675     & SQUEEZE_RIGHT , myThid)
676 torge 1.8 WRITE(msgBuf,'(A,7F8.4)')
677 torge 1.7 & ' SEAICE_GROWTH: Area increments 0, AREAITD = ',
678 torge 1.8 & AREAITD(1,1,:,bi,bj)
679 torge 1.7 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
680     & SQUEEZE_RIGHT , myThid)
681 torge 1.9 #endif
682 torge 1.7 #else
683 torge 1.9 #ifdef SEAICE_DEBUG
684 torge 1.8 WRITE(msgBuf,'(A,7F8.4)')
685 torge 1.7 & ' SEAICE_GROWTH: Heff increments 0, HEFF = ',
686 torge 1.8 & HEFF(1,1,bi,bj)
687 torge 1.7 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
688     & SQUEEZE_RIGHT , myThid)
689 torge 1.8 WRITE(msgBuf,'(A,7F8.4)')
690 torge 1.7 & ' SEAICE_GROWTH: Area increments 0, AREA = ',
691 torge 1.8 & AREA(1,1,bi,bj)
692 torge 1.7 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
693     & SQUEEZE_RIGHT , myThid)
694     c ToM>>>
695 torge 1.9 #endif
696 torge 1.7 #endif /* SEAICE_ITD */
697 dimitri 1.1 #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
698 torge 1.5 C end SEAICEadjMODE.EQ.0 statement:
699 dimitri 1.1 ENDIF
700     #endif
701    
702     C 3) store regularized values of heff, hsnow, area at the onset of thermo.
703     DO J=1,sNy
704     DO I=1,sNx
705     HEFFpreTH(I,J)=HEFF(I,J,bi,bj)
706     HSNWpreTH(I,J)=HSNOW(I,J,bi,bj)
707     AREApreTH(I,J)=AREA(I,J,bi,bj)
708     #ifdef ALLOW_DIAGNOSTICS
709     DIAGarrayB(I,J) = AREA(I,J,bi,bj)
710     DIAGarrayC(I,J) = HEFF(I,J,bi,bj)
711     DIAGarrayD(I,J) = HSNOW(I,J,bi,bj)
712     #endif
713     #ifdef ALLOW_SITRACER
714     SItrHEFF(I,J,bi,bj,1)=HEFF(I,J,bi,bj)
715     SItrAREA(I,J,bi,bj,2)=AREA(I,J,bi,bj)
716     #endif
717     ENDDO
718     ENDDO
719 dimitri 1.2 #ifdef SEAICE_ITD
720 torge 1.5 DO IT=1,nITD
721 dimitri 1.2 DO J=1,sNy
722     DO I=1,sNx
723 torge 1.5 HEFFITDpreTH(I,J,IT)=HEFFITD(I,J,IT,bi,bj)
724     HSNWITDpreTH(I,J,IT)=HSNOWITD(I,J,IT,bi,bj)
725     AREAITDpreTH(I,J,IT)=AREAITD(I,J,IT,bi,bj)
726 torge 1.3
727     C memorize areal and volume fraction of each ITD category
728 torge 1.8 IF (AREA(I,J,bi,bj) .GT. ZERO) THEN
729 torge 1.5 areaFracFactor(I,J,IT)=AREAITD(I,J,IT,bi,bj)/AREA(I,J,bi,bj)
730 torge 1.3 ELSE
731 torge 1.7 C if there's no ice, potential growth starts in 1st category
732     IF (IT .EQ. 1) THEN
733     areaFracFactor(I,J,IT)=ONE
734     ELSE
735     areaFracFactor(I,J,IT)=ZERO
736     ENDIF
737 torge 1.3 ENDIF
738     ENDDO
739     ENDDO
740     ENDDO
741     C prepare SItrHEFF to be computed as cumulative sum
742 torge 1.5 DO iTr=2,5
743 torge 1.3 DO J=1,sNy
744     DO I=1,sNx
745 torge 1.5 SItrHEFF(I,J,bi,bj,iTr)=ZERO
746 dimitri 1.2 ENDDO
747     ENDDO
748     ENDDO
749 torge 1.3 C prepare SItrAREA to be computed as cumulative sum
750     DO J=1,sNy
751     DO I=1,sNx
752     SItrAREA(I,J,bi,bj,3)=ZERO
753     ENDDO
754     ENDDO
755 dimitri 1.2 #endif
756 dimitri 1.1
757     C 4) treat sea ice salinity pathological cases
758     #ifdef SEAICE_VARIABLE_SALINITY
759     #ifdef ALLOW_AUTODIFF_TAMC
760     CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
761     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
762     #endif /* ALLOW_AUTODIFF_TAMC */
763     DO J=1,sNy
764     DO I=1,sNx
765     IF ( (HSALT(I,J,bi,bj) .LT. 0.0).OR.
766     & (HEFF(I,J,bi,bj) .EQ. 0.0) ) THEN
767     saltFluxAdjust(I,J) = - HEFFM(I,J,bi,bj) *
768     & HSALT(I,J,bi,bj) * recip_deltaTtherm
769     HSALT(I,J,bi,bj) = 0.0 _d 0
770     ENDIF
771     ENDDO
772     ENDDO
773     #endif /* SEAICE_VARIABLE_SALINITY */
774    
775     #endif /* SEAICE_GROWTH_LEGACY */
776    
777     #ifdef ALLOW_DIAGNOSTICS
778     IF ( useDiagnostics ) THEN
779     CALL DIAGNOSTICS_FILL(DIAGarrayA,'SIareaPR',0,1,3,bi,bj,myThid)
780     CALL DIAGNOSTICS_FILL(DIAGarrayB,'SIareaPT',0,1,3,bi,bj,myThid)
781     CALL DIAGNOSTICS_FILL(DIAGarrayC,'SIheffPT',0,1,3,bi,bj,myThid)
782     CALL DIAGNOSTICS_FILL(DIAGarrayD,'SIhsnoPT',0,1,3,bi,bj,myThid)
783     #ifdef ALLOW_SITRACER
784     DO iTr = 1, SItrNumInUse
785     WRITE(diagName,'(A4,I2.2,A2)') 'SItr',iTr,'PT'
786     IF (SItrMate(iTr).EQ.'HEFF') THEN
787     CALL DIAGNOSTICS_FRACT_FILL(
788     I SItracer(1-OLx,1-OLy,bi,bj,iTr),HEFF(1-OLx,1-OLy,bi,bj),
789     I ONE, 1, diagName,0,1,2,bi,bj,myThid )
790     ELSE
791     CALL DIAGNOSTICS_FRACT_FILL(
792     I SItracer(1-OLx,1-OLy,bi,bj,iTr),AREA(1-OLx,1-OLy,bi,bj),
793     I ONE, 1, diagName,0,1,2,bi,bj,myThid )
794     ENDIF
795     ENDDO
796     #endif /* ALLOW_SITRACER */
797     ENDIF
798     #endif /* ALLOW_DIAGNOSTICS */
799    
800     #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
801     Cgf no additional dependency of air-sea fluxes to ice
802     IF ( SEAICEadjMODE.GE.1 ) THEN
803     DO J=1,sNy
804     DO I=1,sNx
805     HEFFpreTH(I,J) = 0. _d 0
806     HSNWpreTH(I,J) = 0. _d 0
807     AREApreTH(I,J) = 0. _d 0
808     ENDDO
809     ENDDO
810 dimitri 1.2 #ifdef SEAICE_ITD
811 torge 1.5 DO IT=1,nITD
812 dimitri 1.2 DO J=1,sNy
813     DO I=1,sNx
814 torge 1.5 HEFFITDpreTH(I,J,IT) = 0. _d 0
815     HSNWITDpreTH(I,J,IT) = 0. _d 0
816     AREAITDpreTH(I,J,IT) = 0. _d 0
817 dimitri 1.2 ENDDO
818     ENDDO
819     ENDDO
820     #endif
821 dimitri 1.1 ENDIF
822     #endif
823    
824     #if (defined (ALLOW_MEAN_SFLUX_COST_CONTRIBUTION) || defined (ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION))
825     DO J=1,sNy
826     DO I=1,sNx
827     AREAforAtmFW(I,J,bi,bj) = AREApreTH(I,J)
828     ENDDO
829     ENDDO
830     #endif
831    
832     C 4) COMPUTE ACTUAL ICE/SNOW THICKNESS; USE MIN/MAX VALUES
833     C TO REGULARIZE SEAICE_SOLVE4TEMP/d_AREA COMPUTATIONS
834    
835     #ifdef ALLOW_AUTODIFF_TAMC
836     CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
837     CADJ STORE HEFFpreTH = comlev1_bibj, key = iicekey, byte = isbyte
838     CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte
839     #endif /* ALLOW_AUTODIFF_TAMC */
840 dimitri 1.2 #ifdef SEAICE_ITD
841 torge 1.5 DO IT=1,nITD
842 dimitri 1.2 #endif
843 dimitri 1.1 DO J=1,sNy
844     DO I=1,sNx
845    
846 dimitri 1.2 #ifdef SEAICE_ITD
847 torge 1.5 IF (HEFFITDpreTH(I,J,IT) .GT. ZERO) THEN
848 dimitri 1.2 #ifdef SEAICE_GROWTH_LEGACY
849 torge 1.3 tmpscal1 = MAX(SEAICE_area_reg/float(nITD),
850 torge 1.5 & AREAITDpreTH(I,J,IT))
851     hsnowActualMult(I,J,IT) = HSNWITDpreTH(I,J,IT)/tmpscal1
852     tmpscal2 = HEFFITDpreTH(I,J,IT)/tmpscal1
853     heffActualMult(I,J,IT) = MAX(tmpscal2,SEAICE_hice_reg)
854 dimitri 1.2 #else /* SEAICE_GROWTH_LEGACY */
855     cif regularize AREA with SEAICE_area_reg
856 torge 1.5 tmpscal1 = SQRT(AREAITDpreTH(I,J,IT) * AREAITDpreTH(I,J,IT)
857 dimitri 1.2 & + area_reg_sq)
858     cif heffActual calculated with the regularized AREA
859 torge 1.5 tmpscal2 = HEFFITDpreTH(I,J,IT) / tmpscal1
860 dimitri 1.2 cif regularize heffActual with SEAICE_hice_reg (add lower bound)
861 torge 1.5 heffActualMult(I,J,IT) = SQRT(tmpscal2 * tmpscal2
862 dimitri 1.2 & + hice_reg_sq)
863     cif hsnowActual calculated with the regularized AREA
864 torge 1.5 hsnowActualMult(I,J,IT) = HSNWITDpreTH(I,J,IT) / tmpscal1
865 dimitri 1.2 #endif /* SEAICE_GROWTH_LEGACY */
866     cif regularize the inverse of heffActual by hice_reg
867 torge 1.5 recip_heffActualMult(I,J,IT) = AREAITDpreTH(I,J,IT) /
868     & sqrt(HEFFITDpreTH(I,J,IT) * HEFFITDpreTH(I,J,IT)
869 dimitri 1.2 & + hice_reg_sq)
870     cif Do not regularize when HEFFpreTH = 0
871     ELSE
872 torge 1.5 heffActualMult(I,J,IT) = ZERO
873     hsnowActualMult(I,J,IT) = ZERO
874     recip_heffActualMult(I,J,IT) = ZERO
875 dimitri 1.2 ENDIF
876 torge 1.3 #else /* SEAICE_ITD */
877 dimitri 1.1 IF (HEFFpreTH(I,J) .GT. ZERO) THEN
878     #ifdef SEAICE_GROWTH_LEGACY
879     tmpscal1 = MAX(SEAICE_area_reg,AREApreTH(I,J))
880     hsnowActual(I,J) = HSNWpreTH(I,J)/tmpscal1
881     tmpscal2 = HEFFpreTH(I,J)/tmpscal1
882     heffActual(I,J) = MAX(tmpscal2,SEAICE_hice_reg)
883     #else /* SEAICE_GROWTH_LEGACY */
884     cif regularize AREA with SEAICE_area_reg
885     tmpscal1 = SQRT(AREApreTH(I,J)* AREApreTH(I,J) + area_reg_sq)
886     cif heffActual calculated with the regularized AREA
887     tmpscal2 = HEFFpreTH(I,J) / tmpscal1
888     cif regularize heffActual with SEAICE_hice_reg (add lower bound)
889     heffActual(I,J) = SQRT(tmpscal2 * tmpscal2 + hice_reg_sq)
890     cif hsnowActual calculated with the regularized AREA
891     hsnowActual(I,J) = HSNWpreTH(I,J) / tmpscal1
892     #endif /* SEAICE_GROWTH_LEGACY */
893     cif regularize the inverse of heffActual by hice_reg
894     recip_heffActual(I,J) = AREApreTH(I,J) /
895     & sqrt(HEFFpreTH(I,J)*HEFFpreTH(I,J) + hice_reg_sq)
896     cif Do not regularize when HEFFpreTH = 0
897     ELSE
898     heffActual(I,J) = ZERO
899     hsnowActual(I,J) = ZERO
900     recip_heffActual(I,J) = ZERO
901     ENDIF
902 torge 1.3 #endif /* SEAICE_ITD */
903 dimitri 1.1
904     ENDDO
905     ENDDO
906 dimitri 1.2 #ifdef SEAICE_ITD
907     ENDDO
908     #endif
909 dimitri 1.1
910     #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
911     CALL ZERO_ADJ_1D( sNx*sNy, heffActual, myThid)
912     CALL ZERO_ADJ_1D( sNx*sNy, hsnowActual, myThid)
913     CALL ZERO_ADJ_1D( sNx*sNy, recip_heffActual, myThid)
914     #endif
915    
916     #ifdef SEAICE_CAP_SUBLIM
917     C 5) COMPUTE MAXIMUM LATENT HEAT FLUXES FOR THE CURRENT ICE
918     C AND SNOW THICKNESS
919 dimitri 1.2 #ifdef SEAICE_ITD
920 torge 1.5 DO IT=1,nITD
921 dimitri 1.2 #endif
922 dimitri 1.1 DO J=1,sNy
923     DO I=1,sNx
924     c The latent heat flux over the sea ice which
925     c will sublimate all of the snow and ice over one time
926     c step (W/m^2)
927 dimitri 1.2 #ifdef SEAICE_ITD
928 torge 1.5 IF (HEFFITDpreTH(I,J,IT) .GT. ZERO) THEN
929     latentHeatFluxMaxMult(I,J,IT) = lhSublim*recip_deltaTtherm *
930     & (HEFFITDpreTH(I,J,IT)*SEAICE_rhoIce +
931     & HSNWITDpreTH(I,J,IT)*SEAICE_rhoSnow)
932     & /AREAITDpreTH(I,J,IT)
933 dimitri 1.2 ELSE
934 torge 1.5 latentHeatFluxMaxMult(I,J,IT) = ZERO
935 dimitri 1.2 ENDIF
936     #else
937 dimitri 1.1 IF (HEFFpreTH(I,J) .GT. ZERO) THEN
938     latentHeatFluxMax(I,J) = lhSublim * recip_deltaTtherm *
939     & (HEFFpreTH(I,J) * SEAICE_rhoIce +
940     & HSNWpreTH(I,J) * SEAICE_rhoSnow)/AREApreTH(I,J)
941     ELSE
942     latentHeatFluxMax(I,J) = ZERO
943     ENDIF
944 dimitri 1.2 #endif
945 dimitri 1.1 ENDDO
946     ENDDO
947 dimitri 1.2 #ifdef SEAICE_ITD
948     ENDDO
949     #endif
950 dimitri 1.1 #endif /* SEAICE_CAP_SUBLIM */
951    
952     C ===================================================================
953     C ================PART 2: determine heat fluxes/stocks===============
954     C ===================================================================
955    
956     C determine available heat due to the atmosphere -- for open water
957     C ================================================================
958    
959     DO j=1,sNy
960     DO i=1,sNx
961     C ocean surface/mixed layer temperature
962     TmixLoc(i,j) = theta(i,j,kSurface,bi,bj)+celsius2K
963     C wind speed from exf
964     UG(I,J) = MAX(SEAICE_EPS,wspeed(I,J,bi,bj))
965     ENDDO
966     ENDDO
967    
968     #ifdef ALLOW_AUTODIFF_TAMC
969     CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
970     CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
971     cCADJ STORE UG = comlev1_bibj, key = iicekey,byte=isbyte
972     cCADJ STORE TmixLoc = comlev1_bibj, key = iicekey,byte=isbyte
973     #endif /* ALLOW_AUTODIFF_TAMC */
974    
975     CALL SEAICE_BUDGET_OCEAN(
976     I UG,
977     I TmixLoc,
978     O a_QbyATM_open, a_QSWbyATM_open,
979     I bi, bj, myTime, myIter, myThid )
980    
981     C determine available heat due to the atmosphere -- for ice covered water
982     C =======================================================================
983    
984     #ifdef ALLOW_ATM_WIND
985     IF (useRelativeWind) THEN
986     C Compute relative wind speed over sea ice.
987     DO J=1,sNy
988     DO I=1,sNx
989     SPEED_SQ =
990     & (uWind(I,J,bi,bj)
991     & +0.5 _d 0*(uVel(i,j,kSurface,bi,bj)
992     & +uVel(i+1,j,kSurface,bi,bj))
993     & -0.5 _d 0*(uice(i,j,bi,bj)+uice(i+1,j,bi,bj)))**2
994     & +(vWind(I,J,bi,bj)
995     & +0.5 _d 0*(vVel(i,j,kSurface,bi,bj)
996     & +vVel(i,j+1,kSurface,bi,bj))
997     & -0.5 _d 0*(vice(i,j,bi,bj)+vice(i,j+1,bi,bj)))**2
998     IF ( SPEED_SQ .LE. SEAICE_EPS_SQ ) THEN
999     UG(I,J)=SEAICE_EPS
1000     ELSE
1001     UG(I,J)=SQRT(SPEED_SQ)
1002     ENDIF
1003     ENDDO
1004     ENDDO
1005     ENDIF
1006     #endif /* ALLOW_ATM_WIND */
1007    
1008     #ifdef ALLOW_AUTODIFF_TAMC
1009     CADJ STORE tice(:,:,bi,bj)
1010     CADJ & = comlev1_bibj, key = iicekey, byte = isbyte
1011     CADJ STORE hsnowActual = comlev1_bibj, key = iicekey, byte = isbyte
1012     CADJ STORE heffActual = comlev1_bibj, key = iicekey, byte = isbyte
1013     CADJ STORE UG = comlev1_bibj, key = iicekey, byte = isbyte
1014     CADJ STORE tices(:,:,:,bi,bj)
1015     CADJ & = comlev1_bibj, key = iicekey, byte = isbyte
1016     CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
1017     CADJ & key = iicekey, byte = isbyte
1018     #endif /* ALLOW_AUTODIFF_TAMC */
1019    
1020     C-- Start loop over multi-categories
1021 torge 1.3 #ifdef SEAICE_ITD
1022     CToM SEAICE_multDim = nITD (see SEAICE_SIZE.h and seaice_readparms.F)
1023     #endif
1024 dimitri 1.1 DO IT=1,SEAICE_multDim
1025     c homogeneous distribution between 0 and 2 x heffActual
1026 dimitri 1.2 #ifndef SEAICE_ITD
1027 dimitri 1.1 pFac = (2.0 _d 0*real(IT)-1.0 _d 0)*recip_multDim
1028 dimitri 1.2 #endif
1029 dimitri 1.1 DO J=1,sNy
1030     DO I=1,sNx
1031 dimitri 1.2 #ifndef SEAICE_ITD
1032     CToM for SEAICE_ITD heffActualMult and latentHeatFluxMaxMult are calculated above
1033     C (instead of heffActual and latentHeatFluxMax)
1034 dimitri 1.1 heffActualMult(I,J,IT)= heffActual(I,J)*pFac
1035     #ifdef SEAICE_CAP_SUBLIM
1036     latentHeatFluxMaxMult(I,J,IT) = latentHeatFluxMax(I,J)*pFac
1037     #endif
1038 dimitri 1.2 #endif
1039 dimitri 1.1 ticeInMult(I,J,IT)=TICES(I,J,IT,bi,bj)
1040     ticeOutMult(I,J,IT)=TICES(I,J,IT,bi,bj)
1041     TICE(I,J,bi,bj) = ZERO
1042     TICES(I,J,IT,bi,bj) = ZERO
1043     ENDDO
1044     ENDDO
1045     ENDDO
1046    
1047     #ifdef ALLOW_AUTODIFF_TAMC
1048     CADJ STORE heffActualMult = comlev1_bibj, key = iicekey, byte = isbyte
1049     CADJ STORE ticeInMult = comlev1_bibj, key = iicekey, byte = isbyte
1050     # ifdef SEAICE_CAP_SUBLIM
1051     CADJ STORE latentHeatFluxMaxMult
1052     CADJ & = comlev1_bibj, key = iicekey, byte = isbyte
1053     # endif
1054     CADJ STORE a_QbyATMmult_cover =
1055     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
1056     CADJ STORE a_QSWbyATMmult_cover =
1057     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
1058     CADJ STORE a_FWbySublimMult =
1059     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
1060     #endif /* ALLOW_AUTODIFF_TAMC */
1061    
1062     DO IT=1,SEAICE_multDim
1063     CALL SEAICE_SOLVE4TEMP(
1064 dimitri 1.2 #ifdef SEAICE_ITD
1065     I UG, heffActualMult(1,1,IT), hsnowActualMult(1,1,IT),
1066     #else
1067 dimitri 1.1 I UG, heffActualMult(1,1,IT), hsnowActual,
1068 dimitri 1.2 #endif
1069 dimitri 1.1 #ifdef SEAICE_CAP_SUBLIM
1070     I latentHeatFluxMaxMult(1,1,IT),
1071     #endif
1072     U ticeInMult(1,1,IT), ticeOutMult(1,1,IT),
1073     O a_QbyATMmult_cover(1,1,IT), a_QSWbyATMmult_cover(1,1,IT),
1074     O a_FWbySublimMult(1,1,IT),
1075     I bi, bj, myTime, myIter, myThid )
1076     ENDDO
1077    
1078     #ifdef ALLOW_AUTODIFF_TAMC
1079     CADJ STORE heffActualMult = comlev1_bibj, key = iicekey, byte = isbyte
1080     CADJ STORE ticeOutMult = comlev1_bibj, key = iicekey, byte = isbyte
1081     # ifdef SEAICE_CAP_SUBLIM
1082     CADJ STORE latentHeatFluxMaxMult
1083     CADJ & = comlev1_bibj, key = iicekey, byte = isbyte
1084     # endif
1085     CADJ STORE a_QbyATMmult_cover =
1086     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
1087     CADJ STORE a_QSWbyATMmult_cover =
1088     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
1089     CADJ STORE a_FWbySublimMult =
1090     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
1091     #endif /* ALLOW_AUTODIFF_TAMC */
1092    
1093     DO IT=1,SEAICE_multDim
1094     DO J=1,sNy
1095     DO I=1,sNx
1096     C update TICE & TICES
1097 dimitri 1.2 #ifdef SEAICE_ITD
1098     C calculate area weighted mean
1099 torge 1.3 C (although the ice's temperature relates to its energy content
1100 torge 1.7 C and hence should be averaged weighted by ice volume,
1101 torge 1.3 C the temperature here is a result of the fluxes through the ice surface
1102     C computed individually for each single category in SEAICE_SOLVE4TEMP
1103     C and hence is averaged area weighted [areaFracFactor])
1104 dimitri 1.2 TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
1105 torge 1.5 & + ticeOutMult(I,J,IT)*areaFracFactor(I,J,IT)
1106 dimitri 1.2 #else
1107 dimitri 1.1 TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
1108     & + ticeOutMult(I,J,IT)*recip_multDim
1109 dimitri 1.2 #endif
1110 dimitri 1.1 TICES(I,J,IT,bi,bj) = ticeOutMult(I,J,IT)
1111     C average over categories
1112 dimitri 1.2 #ifdef SEAICE_ITD
1113     C calculate area weighted mean
1114 torge 1.3 C (fluxes are per unit (ice surface) area and are thus area weighted)
1115 dimitri 1.2 a_QbyATM_cover (I,J) = a_QbyATM_cover(I,J)
1116 torge 1.5 & + a_QbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,IT)
1117 dimitri 1.2 a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
1118 torge 1.5 & + a_QSWbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,IT)
1119 dimitri 1.2 a_FWbySublim (I,J) = a_FWbySublim(I,J)
1120 torge 1.5 & + a_FWbySublimMult(I,J,IT)*areaFracFactor(I,J,IT)
1121 dimitri 1.2 #else
1122 dimitri 1.1 a_QbyATM_cover (I,J) = a_QbyATM_cover(I,J)
1123     & + a_QbyATMmult_cover(I,J,IT)*recip_multDim
1124     a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
1125     & + a_QSWbyATMmult_cover(I,J,IT)*recip_multDim
1126     a_FWbySublim (I,J) = a_FWbySublim(I,J)
1127     & + a_FWbySublimMult(I,J,IT)*recip_multDim
1128 dimitri 1.2 #endif
1129 dimitri 1.1 ENDDO
1130     ENDDO
1131     ENDDO
1132    
1133     #ifdef SEAICE_CAP_SUBLIM
1134     # ifdef ALLOW_DIAGNOSTICS
1135     DO J=1,sNy
1136     DO I=1,sNx
1137     c The actual latent heat flux realized by SOLVE4TEMP
1138     DIAGarrayA(I,J) = a_FWbySublim(I,J) * lhSublim
1139     ENDDO
1140     ENDDO
1141     cif The actual vs. maximum latent heat flux
1142     IF ( useDiagnostics ) THEN
1143     CALL DIAGNOSTICS_FILL(DIAGarrayA,
1144     & 'SIactLHF',0,1,3,bi,bj,myThid)
1145     CALL DIAGNOSTICS_FILL(latentHeatFluxMax,
1146     & 'SImaxLHF',0,1,3,bi,bj,myThid)
1147     ENDIF
1148     # endif /* ALLOW_DIAGNOSTICS */
1149     #endif /* SEAICE_CAP_SUBLIM */
1150    
1151     #ifdef ALLOW_AUTODIFF_TAMC
1152     CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
1153     CADJ STORE a_QbyATM_cover = comlev1_bibj, key = iicekey, byte = isbyte
1154     CADJ STORE a_QSWbyATM_cover= comlev1_bibj, key = iicekey, byte = isbyte
1155     CADJ STORE a_QbyATM_open = comlev1_bibj, key = iicekey, byte = isbyte
1156     CADJ STORE a_QSWbyATM_open = comlev1_bibj, key = iicekey, byte = isbyte
1157     CADJ STORE a_FWbySublim = comlev1_bibj, key = iicekey, byte = isbyte
1158     #endif /* ALLOW_AUTODIFF_TAMC */
1159    
1160     C switch heat fluxes from W/m2 to 'effective' ice meters
1161 torge 1.3 #ifdef SEAICE_ITD
1162 torge 1.5 DO IT=1,nITD
1163 torge 1.3 DO J=1,sNy
1164     DO I=1,sNx
1165 torge 1.5 a_QbyATMmult_cover(I,J,IT) = a_QbyATMmult_cover(I,J,IT)
1166     & * convertQ2HI * AREAITDpreTH(I,J,IT)
1167     a_QSWbyATMmult_cover(I,J,IT) = a_QSWbyATMmult_cover(I,J,IT)
1168     & * convertQ2HI * AREAITDpreTH(I,J,IT)
1169 torge 1.8 C and initialize r_QbyATMmult_cover
1170 torge 1.5 r_QbyATMmult_cover(I,J,IT)=a_QbyATMmult_cover(I,J,IT)
1171 torge 1.3 C Convert fresh water flux by sublimation to 'effective' ice meters.
1172     C Negative sublimation is resublimation and will be added as snow.
1173     #ifdef SEAICE_DISABLE_SUBLIM
1174 torge 1.5 a_FWbySublimMult(I,J,IT) = ZERO
1175 torge 1.3 #endif
1176 torge 1.5 a_FWbySublimMult(I,J,IT) = SEAICE_deltaTtherm*recip_rhoIce
1177     & * a_FWbySublimMult(I,J,IT)*AREAITDpreTH(I,J,IT)
1178     r_FWbySublimMult(I,J,IT)=a_FWbySublimMult(I,J,IT)
1179 torge 1.3 ENDDO
1180     ENDDO
1181     ENDDO
1182     DO J=1,sNy
1183     DO I=1,sNx
1184 torge 1.4 a_QbyATM_open(I,J) = a_QbyATM_open(I,J)
1185     & * convertQ2HI * ( ONE - AREApreTH(I,J) )
1186     a_QSWbyATM_open(I,J) = a_QSWbyATM_open(I,J)
1187     & * convertQ2HI * ( ONE - AREApreTH(I,J) )
1188 torge 1.3 C and initialize r_QbyATM_open
1189     r_QbyATM_open(I,J)=a_QbyATM_open(I,J)
1190     ENDDO
1191     ENDDO
1192     #else /* SEAICE_ITD */
1193 dimitri 1.1 DO J=1,sNy
1194     DO I=1,sNx
1195     a_QbyATM_cover(I,J) = a_QbyATM_cover(I,J)
1196     & * convertQ2HI * AREApreTH(I,J)
1197     a_QSWbyATM_cover(I,J) = a_QSWbyATM_cover(I,J)
1198     & * convertQ2HI * AREApreTH(I,J)
1199     a_QbyATM_open(I,J) = a_QbyATM_open(I,J)
1200     & * convertQ2HI * ( ONE - AREApreTH(I,J) )
1201     a_QSWbyATM_open(I,J) = a_QSWbyATM_open(I,J)
1202     & * convertQ2HI * ( ONE - AREApreTH(I,J) )
1203     C and initialize r_QbyATM_cover/r_QbyATM_open
1204     r_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)
1205     r_QbyATM_open(I,J)=a_QbyATM_open(I,J)
1206     C Convert fresh water flux by sublimation to 'effective' ice meters.
1207     C Negative sublimation is resublimation and will be added as snow.
1208     #ifdef SEAICE_DISABLE_SUBLIM
1209     cgf just for those who may need to omit this term to reproduce old results
1210     a_FWbySublim(I,J) = ZERO
1211 dimitri 1.2 #endif
1212 dimitri 1.1 a_FWbySublim(I,J) = SEAICE_deltaTtherm*recip_rhoIce
1213     & * a_FWbySublim(I,J)*AREApreTH(I,J)
1214     r_FWbySublim(I,J)=a_FWbySublim(I,J)
1215     ENDDO
1216     ENDDO
1217 torge 1.3 #endif /* SEAICE_ITD */
1218 dimitri 1.1
1219     #ifdef ALLOW_AUTODIFF_TAMC
1220     CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
1221     CADJ STORE a_QbyATM_cover = comlev1_bibj, key = iicekey, byte = isbyte
1222     CADJ STORE a_QSWbyATM_cover= comlev1_bibj, key = iicekey, byte = isbyte
1223     CADJ STORE a_QbyATM_open = comlev1_bibj, key = iicekey, byte = isbyte
1224     CADJ STORE a_QSWbyATM_open = comlev1_bibj, key = iicekey, byte = isbyte
1225     CADJ STORE a_FWbySublim = comlev1_bibj, key = iicekey, byte = isbyte
1226     CADJ STORE r_QbyATM_cover = comlev1_bibj, key = iicekey, byte = isbyte
1227     CADJ STORE r_QbyATM_open = comlev1_bibj, key = iicekey, byte = isbyte
1228     CADJ STORE r_FWbySublim = comlev1_bibj, key = iicekey, byte = isbyte
1229     #endif /* ALLOW_AUTODIFF_TAMC */
1230    
1231     #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
1232     Cgf no additional dependency through ice cover
1233     IF ( SEAICEadjMODE.GE.3 ) THEN
1234 dimitri 1.2 #ifdef SEAICE_ITD
1235 torge 1.5 DO IT=1,nITD
1236 dimitri 1.2 DO J=1,sNy
1237     DO I=1,sNx
1238 torge 1.5 a_QbyATMmult_cover(I,J,IT) = 0. _d 0
1239     r_QbyATMmult_cover(I,J,IT) = 0. _d 0
1240     a_QSWbyATMmult_cover(I,J,IT) = 0. _d 0
1241 dimitri 1.2 ENDDO
1242     ENDDO
1243     ENDDO
1244 torge 1.3 #else
1245     DO J=1,sNy
1246     DO I=1,sNx
1247     a_QbyATM_cover(I,J) = 0. _d 0
1248     r_QbyATM_cover(I,J) = 0. _d 0
1249     a_QSWbyATM_cover(I,J) = 0. _d 0
1250     ENDDO
1251     ENDDO
1252 dimitri 1.2 #endif
1253 dimitri 1.1 ENDIF
1254     #endif
1255    
1256     C determine available heat due to the ice pack tying the
1257     C underlying surface water temperature to freezing point
1258     C ======================================================
1259    
1260     #ifdef ALLOW_AUTODIFF_TAMC
1261     CADJ STORE theta(:,:,kSurface,bi,bj) = comlev1_bibj,
1262     CADJ & key = iicekey, byte = isbyte
1263     CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
1264     CADJ & key = iicekey, byte = isbyte
1265     #endif
1266    
1267     DO J=1,sNy
1268     DO I=1,sNx
1269     c FREEZING TEMP. OF SEA WATER (deg C)
1270     tempFrz = SEAICE_tempFrz0 +
1271     & SEAICE_dTempFrz_dS *salt(I,J,kSurface,bi,bj)
1272     c efficiency of turbulent fluxes : dependency to sign of THETA-TBC
1273     IF ( theta(I,J,kSurface,bi,bj) .GE. tempFrz ) THEN
1274     tmpscal1 = SEAICE_mcPheePiston
1275     ELSE
1276     tmpscal1 =SEAICE_frazilFrac*drF(kSurface)/SEAICE_deltaTtherm
1277     ENDIF
1278     c efficiency of turbulent fluxes : dependency to AREA (McPhee cases)
1279     IF ( (AREApreTH(I,J) .GT. 0. _d 0).AND.
1280     & (.NOT.SEAICE_mcPheeStepFunc) ) THEN
1281     MixedLayerTurbulenceFactor = ONE -
1282     & SEAICE_mcPheeTaper * AREApreTH(I,J)
1283     ELSEIF ( (AREApreTH(I,J) .GT. 0. _d 0).AND.
1284     & (SEAICE_mcPheeStepFunc) ) THEN
1285     MixedLayerTurbulenceFactor = ONE - SEAICE_mcPheeTaper
1286     ELSE
1287     MixedLayerTurbulenceFactor = ONE
1288     ENDIF
1289     c maximum turbulent flux, in ice meters
1290     tmpscal2= - (HeatCapacity_Cp*rhoConst * recip_QI)
1291     & * (theta(I,J,kSurface,bi,bj)-tempFrz)
1292     & * SEAICE_deltaTtherm * maskC(i,j,kSurface,bi,bj)
1293     c available turbulent flux
1294     a_QbyOCN(i,j) =
1295     & tmpscal1 * tmpscal2 * MixedLayerTurbulenceFactor
1296     r_QbyOCN(i,j) = a_QbyOCN(i,j)
1297     ENDDO
1298     ENDDO
1299    
1300     #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
1301     CALL ZERO_ADJ_1D( sNx*sNy, r_QbyOCN, myThid)
1302     #endif
1303    
1304    
1305     C ===================================================================
1306     C =========PART 3: determine effective thicknesses increments========
1307     C ===================================================================
1308    
1309     C compute snow/ice tendency due to sublimation
1310     C ============================================
1311    
1312     #ifdef ALLOW_AUTODIFF_TAMC
1313     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1314     CADJ STORE r_FWbySublim = comlev1_bibj,key=iicekey,byte=isbyte
1315     #endif /* ALLOW_AUTODIFF_TAMC */
1316 dimitri 1.2 #ifdef SEAICE_ITD
1317 torge 1.5 DO IT=1,nITD
1318 dimitri 1.2 #endif
1319 dimitri 1.1 DO J=1,sNy
1320     DO I=1,sNx
1321     C First sublimate/deposite snow
1322     tmpscal2 =
1323 dimitri 1.2 #ifdef SEAICE_ITD
1324 torge 1.5 & MAX(MIN(r_FWbySublimMult(I,J,IT),HSNOWITD(I,J,IT,bi,bj)
1325 dimitri 1.2 & *SNOW2ICE),ZERO)
1326     C accumulate change over ITD categories
1327 torge 1.5 d_HSNWbySublim(I,J) = d_HSNWbySublim(I,J) - tmpscal2
1328 dimitri 1.2 & *ICE2SNOW
1329 torge 1.5 HSNOWITD(I,J,IT,bi,bj) = HSNOWITD(I,J,IT,bi,bj) - tmpscal2
1330 dimitri 1.2 & *ICE2SNOW
1331 torge 1.5 r_FWbySublimMult(I,J,IT)= r_FWbySublimMult(I,J,IT) - tmpscal2
1332 dimitri 1.2 #else
1333 dimitri 1.1 & MAX(MIN(r_FWbySublim(I,J),HSNOW(I,J,bi,bj)*SNOW2ICE),ZERO)
1334     d_HSNWbySublim(I,J) = - tmpscal2 * ICE2SNOW
1335     HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) - tmpscal2*ICE2SNOW
1336     r_FWbySublim(I,J) = r_FWbySublim(I,J) - tmpscal2
1337 dimitri 1.2 #endif
1338 dimitri 1.1 ENDDO
1339     ENDDO
1340     #ifdef ALLOW_AUTODIFF_TAMC
1341     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1342     CADJ STORE r_FWbySublim = comlev1_bibj,key=iicekey,byte=isbyte
1343     #endif /* ALLOW_AUTODIFF_TAMC */
1344     DO J=1,sNy
1345     DO I=1,sNx
1346     C If anything is left, sublimate ice
1347     tmpscal2 =
1348 dimitri 1.2 #ifdef SEAICE_ITD
1349 torge 1.5 & MAX(MIN(r_FWbySublimMult(I,J,IT),HEFFITD(I,J,IT,bi,bj)),ZERO)
1350 torge 1.3 C accumulate change over ITD categories
1351 torge 1.8 d_HSNWbySublim(I,J) = d_HSNWbySublim(I,J) - tmpscal2
1352 torge 1.5 HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj) - tmpscal2
1353     r_FWbySublimMult(I,J,IT) = r_FWbySublimMult(I,J,IT) - tmpscal2
1354 dimitri 1.2 #else
1355 dimitri 1.1 & MAX(MIN(r_FWbySublim(I,J),HEFF(I,J,bi,bj)),ZERO)
1356     d_HEFFbySublim(I,J) = - tmpscal2
1357     HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) - tmpscal2
1358     r_FWbySublim(I,J) = r_FWbySublim(I,J) - tmpscal2
1359 dimitri 1.2 #endif
1360 dimitri 1.1 ENDDO
1361     ENDDO
1362     DO J=1,sNy
1363     DO I=1,sNx
1364     C If anything is left, it will be evaporated from the ocean rather than sublimated.
1365 dimitri 1.2 C Since a_QbyATM_cover was computed for sublimation, not simple evaporation, we need to
1366 dimitri 1.1 C remove the fusion part for the residual (that happens to be precisely r_FWbySublim).
1367 dimitri 1.2 #ifdef SEAICE_ITD
1368 torge 1.5 a_QbyATMmult_cover(I,J,IT) = a_QbyATMmult_cover(I,J,IT)
1369     & - r_FWbySublimMult(I,J,IT)
1370     r_QbyATMmult_cover(I,J,IT) = r_QbyATMmult_cover(I,J,IT)
1371     & - r_FWbySublimMult(I,J,IT)
1372 torge 1.8 #else
1373     a_QbyATM_cover(I,J) = a_QbyATM_cover(I,J)-r_FWbySublim(I,J)
1374     r_QbyATM_cover(I,J) = r_QbyATM_cover(I,J)-r_FWbySublim(I,J)
1375     #endif
1376 dimitri 1.2 ENDDO
1377     ENDDO
1378 torge 1.8 #ifdef SEAICE_ITD
1379 torge 1.5 C end IT loop
1380 dimitri 1.2 ENDDO
1381     #endif
1382 torge 1.9 #ifdef SEAICE_DEBUG
1383 torge 1.3 c ToM<<< debug seaice_growth
1384 torge 1.8 WRITE(msgBuf,'(A,7F8.4)')
1385 torge 1.4 #ifdef SEAICE_ITD
1386 torge 1.3 & ' SEAICE_GROWTH: Heff increments 1, HEFFITD = ',
1387 torge 1.8 & HEFFITD(1,1,:,bi,bj)
1388     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1389     & SQUEEZE_RIGHT , myThid)
1390     WRITE(msgBuf,'(A,7F8.4)')
1391     & ' SEAICE_GROWTH: Area increments 1, AREAITD = ',
1392     & AREAITD(1,1,:,bi,bj)
1393 torge 1.4 #else
1394     & ' SEAICE_GROWTH: Heff increments 1, HEFF = ',
1395 torge 1.8 & HEFF(1,1,bi,bj)
1396 torge 1.4 #endif
1397 torge 1.3 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1398     & SQUEEZE_RIGHT , myThid)
1399     c ToM>>>
1400 torge 1.9 #endif
1401 dimitri 1.1
1402     C compute ice thickness tendency due to ice-ocean interaction
1403     C ===========================================================
1404    
1405     #ifdef ALLOW_AUTODIFF_TAMC
1406     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1407     CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1408     #endif /* ALLOW_AUTODIFF_TAMC */
1409    
1410 dimitri 1.2 #ifdef SEAICE_ITD
1411 torge 1.5 DO IT=1,nITD
1412 dimitri 1.2 DO J=1,sNy
1413     DO I=1,sNx
1414 torge 1.7 C ice growth/melt due to ocean heat r_QbyOCN (W/m^2) is
1415     C equally distributed under the ice and hence weighted by
1416     C fractional area of each thickness category
1417 torge 1.5 tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,IT),
1418     & -HEFFITD(I,J,IT,bi,bj))
1419 torge 1.7 d_HEFFbyOCNonICE(I,J) = d_HEFFbyOCNonICE(I,J) + tmpscal1
1420     HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj) + tmpscal1
1421 torge 1.3 #ifdef ALLOW_SITRACER
1422     SItrHEFF(I,J,bi,bj,2) = SItrHEFF(I,J,bi,bj,2)
1423 torge 1.5 & + HEFFITD(I,J,IT,bi,bj)
1424 torge 1.3 #endif
1425 dimitri 1.2 ENDDO
1426     ENDDO
1427     ENDDO
1428 torge 1.7 DO J=1,sNy
1429     DO I=1,sNx
1430     r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)
1431     ENDDO
1432     ENDDO
1433 torge 1.3 #else /* SEAICE_ITD */
1434 dimitri 1.1 DO J=1,sNy
1435     DO I=1,sNx
1436     d_HEFFbyOCNonICE(I,J)=MAX(r_QbyOCN(i,j), -HEFF(I,J,bi,bj))
1437     r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)
1438     HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj) + d_HEFFbyOCNonICE(I,J)
1439     #ifdef ALLOW_SITRACER
1440     SItrHEFF(I,J,bi,bj,2)=HEFF(I,J,bi,bj)
1441     #endif
1442     ENDDO
1443     ENDDO
1444 torge 1.3 #endif /* SEAICE_ITD */
1445 torge 1.9 #ifdef SEAICE_DEBUG
1446 torge 1.3 c ToM<<< debug seaice_growth
1447 torge 1.8 WRITE(msgBuf,'(A,7F8.4)')
1448 torge 1.4 #ifdef SEAICE_ITD
1449 torge 1.3 & ' SEAICE_GROWTH: Heff increments 2, HEFFITD = ',
1450 torge 1.8 & HEFFITD(1,1,:,bi,bj)
1451     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1452     & SQUEEZE_RIGHT , myThid)
1453     WRITE(msgBuf,'(A,7F8.4)')
1454     & ' SEAICE_GROWTH: Area increments 2, AREAITD = ',
1455     & AREAITD(1,1,:,bi,bj)
1456 torge 1.4 #else
1457     & ' SEAICE_GROWTH: Heff increments 2, HEFF = ',
1458 torge 1.8 & HEFF(1,1,bi,bj)
1459 torge 1.4 #endif
1460 torge 1.3 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1461     & SQUEEZE_RIGHT , myThid)
1462     c ToM>>>
1463 torge 1.9 #endif
1464 dimitri 1.1
1465     C compute snow melt tendency due to snow-atmosphere interaction
1466     C ==================================================================
1467    
1468     #ifdef ALLOW_AUTODIFF_TAMC
1469     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1470     CADJ STORE r_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1471     #endif /* ALLOW_AUTODIFF_TAMC */
1472    
1473 dimitri 1.2 #ifdef SEAICE_ITD
1474 torge 1.5 DO IT=1,nITD
1475 dimitri 1.2 DO J=1,sNy
1476     DO I=1,sNx
1477     C Convert to standard units (meters of ice) rather than to meters
1478     C of snow. This appears to be more robust.
1479 torge 1.5 tmpscal1=MAX(r_QbyATMmult_cover(I,J,IT),
1480     & -HSNOWITD(I,J,IT,bi,bj)*SNOW2ICE)
1481     tmpscal2=MIN(tmpscal1,0. _d 0)
1482 dimitri 1.2 #ifdef SEAICE_MODIFY_GROWTH_ADJ
1483     Cgf no additional dependency through snow
1484 torge 1.5 IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1485 dimitri 1.2 #endif
1486 torge 1.5 d_HSNWbyATMonSNW(I,J) = d_HSNWbyATMonSNW(I,J)
1487     & + tmpscal2*ICE2SNOW
1488     HSNOWITD(I,J,IT,bi,bj)= HSNOWITD(I,J,IT,bi,bj)
1489     & + tmpscal2*ICE2SNOW
1490     r_QbyATMmult_cover(I,J,IT)=r_QbyATMmult_cover(I,J,IT)
1491     & - tmpscal2
1492 dimitri 1.2 ENDDO
1493     ENDDO
1494     ENDDO
1495 torge 1.3 #else /* SEAICE_ITD */
1496 dimitri 1.1 DO J=1,sNy
1497     DO I=1,sNx
1498     C Convert to standard units (meters of ice) rather than to meters
1499     C of snow. This appears to be more robust.
1500     tmpscal1=MAX(r_QbyATM_cover(I,J),-HSNOW(I,J,bi,bj)*SNOW2ICE)
1501     tmpscal2=MIN(tmpscal1,0. _d 0)
1502     #ifdef SEAICE_MODIFY_GROWTH_ADJ
1503     Cgf no additional dependency through snow
1504     IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1505     #endif
1506     d_HSNWbyATMonSNW(I,J)= tmpscal2*ICE2SNOW
1507 torge 1.3 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + tmpscal2*ICE2SNOW
1508 dimitri 1.1 r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2
1509     ENDDO
1510     ENDDO
1511 dimitri 1.2 #endif /* SEAICE_ITD */
1512 torge 1.9 #ifdef SEAICE_DEBUG
1513 torge 1.3 c ToM<<< debug seaice_growth
1514 torge 1.8 WRITE(msgBuf,'(A,7F8.4)')
1515 torge 1.4 #ifdef SEAICE_ITD
1516 torge 1.3 & ' SEAICE_GROWTH: Heff increments 3, HEFFITD = ',
1517 torge 1.8 & HEFFITD(1,1,:,bi,bj)
1518     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1519     & SQUEEZE_RIGHT , myThid)
1520     WRITE(msgBuf,'(A,7F8.4)')
1521     & ' SEAICE_GROWTH: Area increments 3, AREAITD = ',
1522     & AREAITD(1,1,:,bi,bj)
1523 torge 1.4 #else
1524     & ' SEAICE_GROWTH: Heff increments 3, HEFF = ',
1525 torge 1.8 & HEFF(1,1,bi,bj)
1526 torge 1.4 #endif
1527 torge 1.3 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1528     & SQUEEZE_RIGHT , myThid)
1529     c ToM>>>
1530 torge 1.9 #endif
1531 dimitri 1.1
1532     C compute ice thickness tendency due to the atmosphere
1533     C ====================================================
1534    
1535     #ifdef ALLOW_AUTODIFF_TAMC
1536     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1537     CADJ STORE r_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1538     #endif /* ALLOW_AUTODIFF_TAMC */
1539    
1540     Cgf note: this block is not actually tested by lab_sea
1541     Cgf where all experiments start in January. So even though
1542     Cgf the v1.81=>v1.82 revision would change results in
1543     Cgf warming conditions, the lab_sea results were not changed.
1544    
1545 dimitri 1.2 #ifdef SEAICE_ITD
1546 torge 1.5 DO IT=1,nITD
1547 dimitri 1.2 DO J=1,sNy
1548     DO I=1,sNx
1549     #ifdef SEAICE_GROWTH_LEGACY
1550 torge 1.5 tmpscal2 = MAX(-HEFFITD(I,J,IT,bi,bj),
1551     & r_QbyATMmult_cover(I,J,IT))
1552 dimitri 1.2 #else
1553 torge 1.5 tmpscal2 = MAX(-HEFFITD(I,J,IT,bi,bj),
1554     & r_QbyATMmult_cover(I,J,IT)
1555 dimitri 1.2 c Limit ice growth by potential melt by ocean
1556 torge 1.5 & + AREAITDpreTH(I,J,IT) * r_QbyOCN(I,J))
1557 dimitri 1.2 #endif /* SEAICE_GROWTH_LEGACY */
1558     d_HEFFbyATMonOCN_cover(I,J) = d_HEFFbyATMonOCN_cover(I,J)
1559     & + tmpscal2
1560     d_HEFFbyATMonOCN(I,J) = d_HEFFbyATMonOCN(I,J)
1561     & + tmpscal2
1562 torge 1.8 r_QbyATMmult_cover(I,J,IT) = r_QbyATMmult_cover(I,J,IT)
1563 dimitri 1.2 & - tmpscal2
1564 torge 1.5 HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj) + tmpscal2
1565 torge 1.3
1566     #ifdef ALLOW_SITRACER
1567     SItrHEFF(I,J,bi,bj,3) = SItrHEFF(I,J,bi,bj,3)
1568 torge 1.5 & + HEFFITD(I,J,IT,bi,bj)
1569 torge 1.3 #endif
1570 dimitri 1.2 ENDDO
1571     ENDDO
1572     ENDDO
1573 torge 1.3 #else /* SEAICE_ITD */
1574 dimitri 1.1 DO J=1,sNy
1575     DO I=1,sNx
1576    
1577     #ifdef SEAICE_GROWTH_LEGACY
1578     tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J))
1579     #else
1580     tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J)+
1581     c Limit ice growth by potential melt by ocean
1582     & AREApreTH(I,J) * r_QbyOCN(I,J))
1583     #endif /* SEAICE_GROWTH_LEGACY */
1584    
1585     d_HEFFbyATMonOCN_cover(I,J)=tmpscal2
1586     d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal2
1587     r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)-tmpscal2
1588 torge 1.3 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal2
1589 dimitri 1.1
1590     #ifdef ALLOW_SITRACER
1591     SItrHEFF(I,J,bi,bj,3)=HEFF(I,J,bi,bj)
1592     #endif
1593 torge 1.3 ENDDO
1594     ENDDO
1595     #endif /* SEAICE_ITD */
1596 torge 1.9 #ifdef SEAICE_DEBUG
1597 torge 1.3 c ToM<<< debug seaice_growth
1598 torge 1.8 WRITE(msgBuf,'(A,7F8.4)')
1599 torge 1.4 #ifdef SEAICE_ITD
1600 torge 1.3 & ' SEAICE_GROWTH: Heff increments 4, HEFFITD = ',
1601 torge 1.8 & HEFFITD(1,1,:,bi,bj)
1602     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1603     & SQUEEZE_RIGHT , myThid)
1604     WRITE(msgBuf,'(A,7F8.4)')
1605     & ' SEAICE_GROWTH: Area increments 4, AREAITD = ',
1606     & AREAITD(1,1,:,bi,bj)
1607 torge 1.4 #else
1608     & ' SEAICE_GROWTH: Heff increments 4, HEFF = ',
1609 torge 1.8 & HEFF(1,1,bi,bj)
1610 torge 1.4 #endif
1611 torge 1.3 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1612     & SQUEEZE_RIGHT , myThid)
1613     c ToM>>>
1614 torge 1.9 #endif
1615 dimitri 1.1
1616     C attribute precip to fresh water or snow stock,
1617     C depending on atmospheric conditions.
1618     C =================================================
1619     #ifdef ALLOW_ATM_TEMP
1620     #ifdef ALLOW_AUTODIFF_TAMC
1621     CADJ STORE a_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1622     CADJ STORE PRECIP(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1623     CADJ STORE AREApreTH = comlev1_bibj,key=iicekey,byte=isbyte
1624     #endif /* ALLOW_AUTODIFF_TAMC */
1625     DO J=1,sNy
1626     DO I=1,sNx
1627     C possible alternatives to the a_QbyATM_cover criterium
1628     c IF (TICE(I,J,bi,bj) .LT. TMIX) THEN
1629     c IF (atemp(I,J,bi,bj) .LT. celsius2K) THEN
1630     IF ( a_QbyATM_cover(I,J).GE. 0. _d 0 ) THEN
1631     C add precip as snow
1632     d_HFRWbyRAIN(I,J)=0. _d 0
1633     d_HSNWbyRAIN(I,J)=convertPRECIP2HI*ICE2SNOW*
1634     & PRECIP(I,J,bi,bj)*AREApreTH(I,J)
1635     ELSE
1636     C add precip to the fresh water bucket
1637     d_HFRWbyRAIN(I,J)=-convertPRECIP2HI*
1638     & PRECIP(I,J,bi,bj)*AREApreTH(I,J)
1639     d_HSNWbyRAIN(I,J)=0. _d 0
1640     ENDIF
1641     ENDDO
1642     ENDDO
1643 dimitri 1.2 #ifdef SEAICE_ITD
1644 torge 1.5 DO IT=1,nITD
1645 dimitri 1.2 DO J=1,sNy
1646     DO I=1,sNx
1647 torge 1.5 HSNOWITD(I,J,IT,bi,bj) = HSNOWITD(I,J,IT,bi,bj)
1648     & + d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,IT)
1649 dimitri 1.2 ENDDO
1650     ENDDO
1651     ENDDO
1652 torge 1.3 #else
1653     DO J=1,sNy
1654     DO I=1,sNx
1655     HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyRAIN(I,J)
1656     ENDDO
1657     ENDDO
1658 dimitri 1.2 #endif
1659 dimitri 1.1 Cgf note: this does not affect air-sea heat flux,
1660     Cgf since the implied air heat gain to turn
1661     Cgf rain to snow is not a surface process.
1662     #endif /* ALLOW_ATM_TEMP */
1663 torge 1.9 #ifdef SEAICE_DEBUG
1664 torge 1.3 c ToM<<< debug seaice_growth
1665 torge 1.8 WRITE(msgBuf,'(A,7F8.4)')
1666 torge 1.4 #ifdef SEAICE_ITD
1667 torge 1.3 & ' SEAICE_GROWTH: Heff increments 5, HEFFITD = ',
1668 torge 1.8 & HEFFITD(1,1,:,bi,bj)
1669     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1670     & SQUEEZE_RIGHT , myThid)
1671     WRITE(msgBuf,'(A,7F8.4)')
1672     & ' SEAICE_GROWTH: Area increments 5, AREAITD = ',
1673     & AREAITD(1,1,:,bi,bj)
1674 torge 1.4 #else
1675     & ' SEAICE_GROWTH: Heff increments 5, HEFF = ',
1676 torge 1.8 & HEFF(1,1,bi,bj)
1677 torge 1.4 #endif
1678 torge 1.3 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1679     & SQUEEZE_RIGHT , myThid)
1680     c ToM>>>
1681 torge 1.9 #endif
1682 dimitri 1.1
1683     C compute snow melt due to heat available from ocean.
1684     C =================================================================
1685    
1686     Cgf do we need to keep this comment and cpp bracket?
1687     Cph( very sensitive bit here by JZ
1688     #ifndef SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING
1689     #ifdef ALLOW_AUTODIFF_TAMC
1690     CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1691     CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1692     #endif /* ALLOW_AUTODIFF_TAMC */
1693 dimitri 1.2
1694     #ifdef SEAICE_ITD
1695 torge 1.5 DO IT=1,nITD
1696 dimitri 1.2 DO J=1,sNy
1697     DO I=1,sNx
1698 torge 1.5 tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,IT),
1699     & -HSNOWITD(I,J,IT,bi,bj))
1700 dimitri 1.2 tmpscal2=MIN(tmpscal1,0. _d 0)
1701     #ifdef SEAICE_MODIFY_GROWTH_ADJ
1702     Cgf no additional dependency through snow
1703     if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1704     #endif
1705 torge 1.3 d_HSNWbyOCNonSNW(I,J) = d_HSNWbyOCNonSNW(I,J) + tmpscal2
1706     r_QbyOCN(I,J)=r_QbyOCN(I,J) - tmpscal2*SNOW2ICE
1707 torge 1.5 HSNOWITD(I,J,IT,bi,bj) = HSNOWITD(I,J,IT,bi,bj) + tmpscal2
1708 dimitri 1.2 ENDDO
1709     ENDDO
1710     ENDDO
1711 torge 1.3 #else /* SEAICE_ITD */
1712 dimitri 1.1 DO J=1,sNy
1713     DO I=1,sNx
1714     tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW, -HSNOW(I,J,bi,bj))
1715     tmpscal2=MIN(tmpscal1,0. _d 0)
1716     #ifdef SEAICE_MODIFY_GROWTH_ADJ
1717     Cgf no additional dependency through snow
1718     if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1719     #endif
1720     d_HSNWbyOCNonSNW(I,J) = tmpscal2
1721     r_QbyOCN(I,J)=r_QbyOCN(I,J)
1722     & -d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
1723     HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+d_HSNWbyOCNonSNW(I,J)
1724     ENDDO
1725     ENDDO
1726 torge 1.3 #endif /* SEAICE_ITD */
1727 dimitri 1.1 #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */
1728     Cph)
1729 torge 1.9 #ifdef SEAICE_DEBUG
1730 torge 1.3 c ToM<<< debug seaice_growth
1731 torge 1.8 WRITE(msgBuf,'(A,7F8.4)')
1732 torge 1.4 #ifdef SEAICE_ITD
1733 torge 1.3 & ' SEAICE_GROWTH: Heff increments 6, HEFFITD = ',
1734 torge 1.8 & HEFFITD(1,1,:,bi,bj)
1735     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1736     & SQUEEZE_RIGHT , myThid)
1737     WRITE(msgBuf,'(A,7F8.4)')
1738     & ' SEAICE_GROWTH: Area increments 6, AREAITD = ',
1739     & AREAITD(1,1,:,bi,bj)
1740 torge 1.4 #else
1741     & ' SEAICE_GROWTH: Heff increments 6, HEFF = ',
1742 torge 1.8 & HEFF(1,1,bi,bj)
1743 torge 1.4 #endif
1744 torge 1.3 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1745     & SQUEEZE_RIGHT , myThid)
1746     c ToM>>>
1747 torge 1.9 #endif
1748 dimitri 1.1
1749     C gain of new ice over open water
1750     C ===============================
1751     #ifdef ALLOW_AUTODIFF_TAMC
1752     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1753     CADJ STORE r_QbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte
1754     CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1755     CADJ STORE a_QSWbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1756     CADJ STORE a_QSWbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte
1757     #endif /* ALLOW_AUTODIFF_TAMC */
1758    
1759     DO J=1,sNy
1760     DO I=1,sNx
1761     c Initial ice growth is triggered by open water
1762     c heat flux overcoming potential melt by ocean
1763     tmpscal1=r_QbyATM_open(I,J)+r_QbyOCN(i,j) *
1764     & (1.0 _d 0 - AREApreTH(I,J))
1765     c Penetrative shortwave flux beyond first layer
1766     c that is therefore not available to ice growth/melt
1767     tmpscal2=SWFracB * a_QSWbyATM_open(I,J)
1768     C impose -HEFF as the maxmum melting if SEAICE_doOpenWaterMelt
1769     C or 0. otherwise (no melting if not SEAICE_doOpenWaterMelt)
1770     tmpscal3=facOpenGrow*MAX(tmpscal1-tmpscal2,
1771     & -HEFF(I,J,bi,bj)*facOpenMelt)*HEFFM(I,J,bi,bj)
1772     d_HEFFbyATMonOCN_open(I,J)=tmpscal3
1773     d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3
1774     r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3
1775 dimitri 1.2 #ifdef SEAICE_ITD
1776 torge 1.9 C open water area fraction
1777     tmpscal0 = ONE-AREApreTH(I,J)
1778     C determine thickness of new ice
1779     ctomC considering the entire open water area to refreeze
1780     ctom tmpscal1 = tmpscal3/tmpscal0
1781     C considering a minimum lead ice thickness of 10 cm
1782     C WATCH that leadIceThickMin is smaller that Hlimit(1)!
1783     leadIceThickMin = 0.1
1784     tmpscal1 = MAX(leadIceThickMin,tmpscal3/tmpscal0)
1785     C adjust ice area fraction covered by new ice
1786     tmpscal0 = tmpscal3/tmpscal1
1787     C then add new ice volume to appropriate thickness category
1788     DO IT=1,nITD
1789     IF (tmpscal1.LT.Hlimit(IT)) THEN
1790     HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj) + tmpscal3
1791     tmpscal3=ZERO
1792     C not sure if AREAITD should be changed here since AREA is incremented
1793     C in PART 4 below in non-itd code
1794     C in this scenario all open water is covered by new ice instantaneously,
1795     C i.e. no delayed lead closing is concidered such as is associated with
1796     C Hibler's h_0 parameter
1797     AREAITD(I,J,IT,bi,bj) = AREAITD(I,J,IT,bi,bj)
1798     & + tmpscal0
1799     tmpscal0=ZERO
1800     ENDIF
1801     ENDDO
1802 torge 1.3 #else
1803     HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3
1804 dimitri 1.2 #endif
1805 dimitri 1.1 ENDDO
1806     ENDDO
1807    
1808     #ifdef ALLOW_SITRACER
1809 torge 1.3 #ifdef SEAICE_ITD
1810 torge 1.5 DO IT=1,nITD
1811 torge 1.3 DO J=1,sNy
1812     DO I=1,sNx
1813     c needs to be here to allow use also with LEGACY branch
1814     SItrHEFF(I,J,bi,bj,4) = SItrHEFF(I,J,bi,bj,4)
1815 torge 1.5 & + HEFFITD(I,J,IT,bi,bj)
1816 torge 1.3 ENDDO
1817     ENDDO
1818     ENDDO
1819     #else
1820 dimitri 1.1 DO J=1,sNy
1821     DO I=1,sNx
1822     c needs to be here to allow use also with LEGACY branch
1823     SItrHEFF(I,J,bi,bj,4)=HEFF(I,J,bi,bj)
1824     ENDDO
1825     ENDDO
1826 torge 1.3 #endif
1827 dimitri 1.1 #endif /* ALLOW_SITRACER */
1828 torge 1.9 #ifdef SEAICE_DEBUG
1829 torge 1.3 c ToM<<< debug seaice_growth
1830 torge 1.8 WRITE(msgBuf,'(A,7F8.4)')
1831 torge 1.4 #ifdef SEAICE_ITD
1832 torge 1.3 & ' SEAICE_GROWTH: Heff increments 7, HEFFITD = ',
1833 torge 1.8 & HEFFITD(1,1,:,bi,bj)
1834     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1835     & SQUEEZE_RIGHT , myThid)
1836     WRITE(msgBuf,'(A,7F8.4)')
1837     & ' SEAICE_GROWTH: Area increments 7, AREAITD = ',
1838     & AREAITD(1,1,:,bi,bj)
1839 torge 1.4 #else
1840     & ' SEAICE_GROWTH: Heff increments 7, HEFF = ',
1841 torge 1.8 & HEFF(1,1,bi,bj)
1842     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1843     & SQUEEZE_RIGHT , myThid)
1844     WRITE(msgBuf,'(A,7F8.4)')
1845     & ' SEAICE_GROWTH: Area increments 7, AREA = ',
1846     & AREA(1,1,bi,bj)
1847 torge 1.4 #endif
1848 torge 1.3 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1849     & SQUEEZE_RIGHT , myThid)
1850     c ToM>>>
1851 torge 1.9 #endif
1852 dimitri 1.1
1853     C convert snow to ice if submerged.
1854     C =================================
1855    
1856     #ifndef SEAICE_GROWTH_LEGACY
1857     C note: in legacy, this process is done at the end
1858     #ifdef ALLOW_AUTODIFF_TAMC
1859     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1860     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1861     #endif /* ALLOW_AUTODIFF_TAMC */
1862     IF ( SEAICEuseFlooding ) THEN
1863 dimitri 1.2 #ifdef SEAICE_ITD
1864 torge 1.5 DO IT=1,nITD
1865 dimitri 1.2 DO J=1,sNy
1866     DO I=1,sNx
1867 torge 1.5 tmpscal0 = (HSNOWITD(I,J,IT,bi,bj)*SEAICE_rhoSnow
1868     & + HEFFITD(I,J,IT,bi,bj) *SEAICE_rhoIce)
1869     & *recip_rhoConst
1870     tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFFITD(I,J,IT,bi,bj))
1871     d_HEFFbyFLOODING(I,J) = d_HEFFbyFLOODING(I,J) + tmpscal1
1872     HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj) + tmpscal1
1873     HSNOWITD(I,J,IT,bi,bj)= HSNOWITD(I,J,IT,bi,bj) - tmpscal1
1874 dimitri 1.2 & * ICE2SNOW
1875     ENDDO
1876     ENDDO
1877     ENDDO
1878     #else
1879 dimitri 1.1 DO J=1,sNy
1880     DO I=1,sNx
1881     tmpscal0 = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
1882     & +HEFF(I,J,bi,bj)*SEAICE_rhoIce)*recip_rhoConst
1883     tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFF(I,J,bi,bj))
1884     d_HEFFbyFLOODING(I,J)=tmpscal1
1885 torge 1.3 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)+d_HEFFbyFLOODING(I,J)
1886     HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-
1887     & d_HEFFbyFLOODING(I,J)*ICE2SNOW
1888 dimitri 1.2 ENDDO
1889     ENDDO
1890     #endif
1891 dimitri 1.1 ENDIF
1892     #endif /* SEAICE_GROWTH_LEGACY */
1893 torge 1.9 #ifdef SEAICE_DEBUG
1894 torge 1.3 c ToM<<< debug seaice_growth
1895 torge 1.8 WRITE(msgBuf,'(A,7F8.4)')
1896 torge 1.4 #ifdef SEAICE_ITD
1897 torge 1.3 & ' SEAICE_GROWTH: Heff increments 8, HEFFITD = ',
1898 torge 1.8 & HEFFITD(1,1,:,bi,bj)
1899 torge 1.7 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1900     & SQUEEZE_RIGHT , myThid)
1901 torge 1.8 WRITE(msgBuf,'(A,7F8.4)')
1902 torge 1.7 & ' SEAICE_GROWTH: Area increments 8, AREAITD = ',
1903 torge 1.8 & AREAITD(1,1,:,bi,bj)
1904 torge 1.4 #else
1905     & ' SEAICE_GROWTH: Heff increments 8, HEFF = ',
1906 torge 1.8 & HEFF(1,1,bi,bj)
1907 torge 1.7 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1908     & SQUEEZE_RIGHT , myThid)
1909 torge 1.8 WRITE(msgBuf,'(A,7F8.4)')
1910 torge 1.7 & ' SEAICE_GROWTH: Area increments 8, AREA = ',
1911 torge 1.8 & AREA(1,1,bi,bj)
1912 torge 1.4 #endif
1913 torge 1.3 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1914     & SQUEEZE_RIGHT , myThid)
1915     c ToM>>>
1916 torge 1.9 #endif
1917 dimitri 1.1
1918     C ===================================================================
1919     C ==========PART 4: determine ice cover fraction increments=========-
1920     C ===================================================================
1921    
1922     #ifdef ALLOW_AUTODIFF_TAMC
1923     CADJ STORE d_HEFFbyATMonOCN = comlev1_bibj,key=iicekey,byte=isbyte
1924     CADJ STORE d_HEFFbyATMonOCN_cover = comlev1_bibj,key=iicekey,byte=isbyte
1925     CADJ STORE d_HEFFbyATMonOCN_open = comlev1_bibj,key=iicekey,byte=isbyte
1926     CADJ STORE d_HEFFbyOCNonICE = comlev1_bibj,key=iicekey,byte=isbyte
1927     CADJ STORE recip_heffActual = comlev1_bibj,key=iicekey,byte=isbyte
1928     CADJ STORE d_hsnwbyatmonsnw = comlev1_bibj,key=iicekey,byte=isbyte
1929     cph(
1930     cphCADJ STORE d_AREAbyATM = comlev1_bibj,key=iicekey,byte=isbyte
1931     cphCADJ STORE d_AREAbyICE = comlev1_bibj,key=iicekey,byte=isbyte
1932     cphCADJ STORE d_AREAbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1933     cph)
1934     CADJ STORE a_QbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte
1935     CADJ STORE heffActual = comlev1_bibj,key=iicekey,byte=isbyte
1936     CADJ STORE AREApreTH = comlev1_bibj,key=iicekey,byte=isbyte
1937     CADJ STORE HEFF(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1938     CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1939     CADJ STORE AREA(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1940     #endif /* ALLOW_AUTODIFF_TAMC */
1941    
1942 torge 1.9 #ifdef SEAICE_ITD
1943     C replaces Hibler '79 scheme and lead closing parameter
1944     C because ITD accounts explicitly for lead openings and
1945     C different melt rates due to varying ice thickness
1946     C
1947     C only consider ice area loss due to total ice thickness loss;
1948     C ice area gain due to freezing of open water is handled above
1949     C under "gain of new ice over open water"
1950     C
1951     C does not account for lateral melt of ice floes
1952     C
1953     C AREAITD is incremented in section "gain of new ice over open water" above
1954     C
1955     DO IT=1,nITD
1956     DO J=1,sNy
1957     DO I=1,sNx
1958     IF (HEFFITD(I,J,IT,bi,bj).LE.ZERO) THEN
1959     AREAITD(I,J,IT,bi,bj)=ZERO
1960     ENDIF
1961     #ifdef ALLOW_SITRACER
1962     SItrAREA(I,J,bi,bj,3) = SItrAREA(I,J,bi,bj,3)
1963     & + AREAITD(I,J,IT,bi,bj)
1964     #endif /* ALLOW_SITRACER */
1965     ENDDO
1966     ENDDO
1967     ENDDO
1968     #else /* SEAICE_ITD */
1969 dimitri 1.1 DO J=1,sNy
1970     DO I=1,sNx
1971    
1972     IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN
1973     recip_HO=1. _d 0 / HO_south
1974     ELSE
1975     recip_HO=1. _d 0 / HO
1976     ENDIF
1977     #ifdef SEAICE_GROWTH_LEGACY
1978     tmpscal0=HEFF(I,J,bi,bj) - d_HEFFbyATMonOCN(I,J)
1979     recip_HH = AREApreTH(I,J) /(tmpscal0+.00001 _d 0)
1980     #else
1981     recip_HH = recip_heffActual(I,J)
1982     #endif
1983    
1984     C gain of ice over open water : computed from
1985     C (SEAICE_areaGainFormula.EQ.1) from growth by ATM
1986     C (SEAICE_areaGainFormula.EQ.2) from predicted growth by ATM
1987     IF (SEAICE_areaGainFormula.EQ.1) THEN
1988     tmpscal4 = MAX(ZERO,d_HEFFbyATMonOCN_open(I,J))
1989     ELSE
1990     tmpscal4=MAX(ZERO,a_QbyATM_open(I,J))
1991     ENDIF
1992    
1993     C loss of ice cover by melting : computed from
1994     C (SEAICE_areaLossFormula.EQ.1) from all but only melt conributions by ATM and OCN
1995     C (SEAICE_areaLossFormula.EQ.2) from net melt-growth>0 by ATM and OCN
1996     C (SEAICE_areaLossFormula.EQ.3) from predicted melt by ATM
1997     IF (SEAICE_areaLossFormula.EQ.1) THEN
1998     tmpscal3 = MIN( 0. _d 0 , d_HEFFbyATMonOCN_cover(I,J) )
1999     & + MIN( 0. _d 0 , d_HEFFbyATMonOCN_open(I,J) )
2000     & + MIN( 0. _d 0 , d_HEFFbyOCNonICE(I,J) )
2001     ELSEIF (SEAICE_areaLossFormula.EQ.2) THEN
2002     tmpscal3 = MIN( 0. _d 0 , d_HEFFbyATMonOCN_cover(I,J)
2003     & + d_HEFFbyATMonOCN_open(I,J) + d_HEFFbyOCNonICE(I,J) )
2004     ELSE
2005     C compute heff after ice melt by ocn:
2006     tmpscal0=HEFF(I,J,bi,bj) - d_HEFFbyATMonOCN(I,J)
2007     C compute available heat left after snow melt by atm:
2008     tmpscal1= a_QbyATM_open(I,J)+a_QbyATM_cover(I,J)
2009     & - d_HSNWbyATMonSNW(I,J)*SNOW2ICE
2010     C could not melt more than all the ice
2011     tmpscal2 = MAX(-tmpscal0,tmpscal1)
2012     tmpscal3 = MIN(ZERO,tmpscal2)
2013     ENDIF
2014    
2015     C apply tendency
2016     IF ( (HEFF(i,j,bi,bj).GT.0. _d 0).OR.
2017     & (HSNOW(i,j,bi,bj).GT.0. _d 0) ) THEN
2018     AREA(I,J,bi,bj)=MAX(0. _d 0,
2019     & MIN( SEAICE_area_max, AREA(I,J,bi,bj)
2020     & + recip_HO*tmpscal4+HALF*recip_HH*tmpscal3 ))
2021     ELSE
2022     AREA(I,J,bi,bj)=0. _d 0
2023     ENDIF
2024     #ifdef ALLOW_SITRACER
2025     SItrAREA(I,J,bi,bj,3)=AREA(I,J,bi,bj)
2026     #endif /* ALLOW_SITRACER */
2027     #ifdef ALLOW_DIAGNOSTICS
2028     d_AREAbyATM(I,J)=
2029     & recip_HO*MAX(ZERO,d_HEFFbyATMonOCN_open(I,J))
2030     & +HALF*recip_HH*MIN(0. _d 0,d_HEFFbyATMonOCN_open(I,J))
2031     d_AREAbyICE(I,J)=
2032     & HALF*recip_HH*MIN(0. _d 0,d_HEFFbyATMonOCN_cover(I,J))
2033     d_AREAbyOCN(I,J)=
2034     & HALF*recip_HH*MIN( 0. _d 0,d_HEFFbyOCNonICE(I,J) )
2035     #endif /* ALLOW_DIAGNOSTICS */
2036     ENDDO
2037     ENDDO
2038 torge 1.9 #endif /* SEAICE_ITD */
2039 dimitri 1.1
2040     #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
2041     Cgf 'bulk' linearization of area=f(HEFF)
2042     IF ( SEAICEadjMODE.GE.1 ) THEN
2043 dimitri 1.2 #ifdef SEAICE_ITD
2044 torge 1.5 DO IT=1,nITD
2045 dimitri 1.2 DO J=1,sNy
2046     DO I=1,sNx
2047 torge 1.5 AREAITD(I,J,IT,bi,bj) = AREAITDpreTH(I,J,IT) + 0.1 _d 0 *
2048     & ( HEFFITD(I,J,IT,bi,bj) - HEFFITDpreTH(I,J,IT) )
2049 dimitri 1.2 ENDDO
2050     ENDDO
2051     ENDDO
2052     #else
2053 dimitri 1.1 DO J=1,sNy
2054     DO I=1,sNx
2055     C AREA(I,J,bi,bj) = 0.1 _d 0 * HEFF(I,J,bi,bj)
2056     AREA(I,J,bi,bj) = AREApreTH(I,J) + 0.1 _d 0 *
2057     & ( HEFF(I,J,bi,bj) - HEFFpreTH(I,J) )
2058     ENDDO
2059     ENDDO
2060 dimitri 1.2 #endif
2061 dimitri 1.1 ENDIF
2062     #endif
2063 torge 1.3 #ifdef SEAICE_ITD
2064     C check categories for consistency with limits after growth/melt
2065     CALL SEAICE_ITD_REDIST(bi, bj, myTime,myIter,myThid)
2066     C finally update total AREA, HEFF, HSNOW
2067     CALL SEAICE_ITD_SUM(bi, bj, myTime,myIter,myThid)
2068     #endif
2069 dimitri 1.1
2070     C ===================================================================
2071     C =============PART 5: determine ice salinity increments=============
2072     C ===================================================================
2073    
2074     #ifndef SEAICE_VARIABLE_SALINITY
2075     # if (defined ALLOW_AUTODIFF_TAMC && defined ALLOW_SALT_PLUME)
2076     CADJ STORE d_HEFFbyNEG = comlev1_bibj,key=iicekey,byte=isbyte
2077     CADJ STORE d_HEFFbyOCNonICE = comlev1_bibj,key=iicekey,byte=isbyte
2078     CADJ STORE d_HEFFbyATMonOCN = comlev1_bibj,key=iicekey,byte=isbyte
2079     CADJ STORE d_HEFFbyATMonOCN_open = comlev1_bibj,key=iicekey,byte=isbyte
2080     CADJ STORE d_HEFFbyATMonOCN_cover = comlev1_bibj,key=iicekey,byte=isbyte
2081     CADJ STORE d_HEFFbyFLOODING = comlev1_bibj,key=iicekey,byte=isbyte
2082     CADJ STORE d_HEFFbySublim = comlev1_bibj,key=iicekey,byte=isbyte
2083     CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
2084     CADJ & key = iicekey, byte = isbyte
2085     # endif /* ALLOW_AUTODIFF_TAMC and ALLOW_SALT_PLUME */
2086     DO J=1,sNy
2087     DO I=1,sNx
2088     tmpscal1 = d_HEFFbyNEG(I,J) + d_HEFFbyOCNonICE(I,J) +
2089     & d_HEFFbyATMonOCN(I,J) + d_HEFFbyFLOODING(I,J)
2090     & + d_HEFFbySublim(I,J)
2091     #ifdef SEAICE_ALLOW_AREA_RELAXATION
2092     + d_HEFFbyRLX(I,J)
2093     #endif
2094     tmpscal2 = tmpscal1 * SEAICE_salt0 * HEFFM(I,J,bi,bj)
2095     & * recip_deltaTtherm * SEAICE_rhoIce
2096     saltFlux(I,J,bi,bj) = tmpscal2
2097     #ifdef ALLOW_SALT_PLUME
2098     tmpscal3 = tmpscal1*salt(I,J,kSurface,bi,bj)*HEFFM(I,J,bi,bj)
2099     & * recip_deltaTtherm * SEAICE_rhoIce
2100     saltPlumeFlux(I,J,bi,bj) = MAX( tmpscal3-tmpscal2 , 0. _d 0)
2101     & *SPsalFRAC
2102     #endif /* ALLOW_SALT_PLUME */
2103     ENDDO
2104     ENDDO
2105     #endif /* ndef SEAICE_VARIABLE_SALINITY */
2106    
2107     #ifdef SEAICE_VARIABLE_SALINITY
2108    
2109     #ifdef SEAICE_GROWTH_LEGACY
2110     # ifdef ALLOW_AUTODIFF_TAMC
2111     CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
2112     # endif /* ALLOW_AUTODIFF_TAMC */
2113     DO J=1,sNy
2114     DO I=1,sNx
2115     C set HSALT = 0 if HSALT < 0 and compute salt to remove from ocean
2116     IF ( HSALT(I,J,bi,bj) .LT. 0.0 ) THEN
2117     saltFluxAdjust(I,J) = - HEFFM(I,J,bi,bj) *
2118     & HSALT(I,J,bi,bj) * recip_deltaTtherm
2119     HSALT(I,J,bi,bj) = 0.0 _d 0
2120     ENDIF
2121     ENDDO
2122     ENDDO
2123     #endif /* SEAICE_GROWTH_LEGACY */
2124    
2125     #ifdef ALLOW_AUTODIFF_TAMC
2126     CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
2127     #endif /* ALLOW_AUTODIFF_TAMC */
2128    
2129     DO J=1,sNy
2130     DO I=1,sNx
2131     C sum up the terms that affect the salt content of the ice pack
2132     tmpscal1=d_HEFFbyOCNonICE(I,J)+d_HEFFbyATMonOCN(I,J)
2133    
2134     C recompute HEFF before thermodynamic updates (which is not AREApreTH in legacy code)
2135     tmpscal2=HEFF(I,J,bi,bj)-tmpscal1-d_HEFFbyFLOODING(I,J)
2136     C tmpscal1 > 0 : m of sea ice that is created
2137     IF ( tmpscal1 .GE. 0.0 ) THEN
2138     saltFlux(I,J,bi,bj) =
2139     & HEFFM(I,J,bi,bj)*recip_deltaTtherm
2140     & *SEAICE_saltFrac*salt(I,J,kSurface,bi,bj)
2141     & *tmpscal1*SEAICE_rhoIce
2142     #ifdef ALLOW_SALT_PLUME
2143     C saltPlumeFlux is defined only during freezing:
2144     saltPlumeFlux(I,J,bi,bj)=
2145     & HEFFM(I,J,bi,bj)*recip_deltaTtherm
2146     & *(ONE-SEAICE_saltFrac)*salt(I,J,kSurface,bi,bj)
2147     & *tmpscal1*SEAICE_rhoIce
2148     & *SPsalFRAC
2149     C if SaltPlumeSouthernOcean=.FALSE. turn off salt plume in Southern Ocean
2150     IF ( .NOT. SaltPlumeSouthernOcean ) THEN
2151     IF ( YC(I,J,bi,bj) .LT. 0.0 _d 0 )
2152     & saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
2153     ENDIF
2154     #endif /* ALLOW_SALT_PLUME */
2155    
2156     C tmpscal1 < 0 : m of sea ice that is melted
2157     ELSE
2158     saltFlux(I,J,bi,bj) =
2159     & HEFFM(I,J,bi,bj)*recip_deltaTtherm
2160     & *HSALT(I,J,bi,bj)
2161     & *tmpscal1/tmpscal2
2162     #ifdef ALLOW_SALT_PLUME
2163     saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
2164     #endif /* ALLOW_SALT_PLUME */
2165     ENDIF
2166     C update HSALT based on surface saltFlux
2167     HSALT(I,J,bi,bj) = HSALT(I,J,bi,bj) +
2168     & saltFlux(I,J,bi,bj) * SEAICE_deltaTtherm
2169     saltFlux(I,J,bi,bj) =
2170     & saltFlux(I,J,bi,bj) + saltFluxAdjust(I,J)
2171     #ifdef SEAICE_GROWTH_LEGACY
2172     C set HSALT = 0 if HEFF = 0 and compute salt to dump into ocean
2173     IF ( HEFF(I,J,bi,bj) .EQ. 0.0 ) THEN
2174     saltFlux(I,J,bi,bj) = saltFlux(I,J,bi,bj) -
2175     & HEFFM(I,J,bi,bj) * HSALT(I,J,bi,bj) * recip_deltaTtherm
2176     HSALT(I,J,bi,bj) = 0.0 _d 0
2177     #ifdef ALLOW_SALT_PLUME
2178     saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
2179     #endif /* ALLOW_SALT_PLUME */
2180     ENDIF
2181     #endif /* SEAICE_GROWTH_LEGACY */
2182     ENDDO
2183     ENDDO
2184     #endif /* SEAICE_VARIABLE_SALINITY */
2185    
2186    
2187     C =======================================================================
2188     C ==LEGACY PART 6 (LEGACY) treat pathological cases, then do flooding ===
2189     C =======================================================================
2190    
2191     #ifdef SEAICE_GROWTH_LEGACY
2192    
2193     C treat values of ice cover fraction oustide
2194     C the [0 1] range, and other such issues.
2195     C ===========================================
2196    
2197     Cgf note: this part cannot be heat and water conserving
2198    
2199     #ifdef ALLOW_AUTODIFF_TAMC
2200     CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
2201     CADJ & key = iicekey, byte = isbyte
2202     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
2203     CADJ & key = iicekey, byte = isbyte
2204     #endif /* ALLOW_AUTODIFF_TAMC */
2205     DO J=1,sNy
2206     DO I=1,sNx
2207     C NOW SET AREA(I,J,bi,bj)=0 WHERE THERE IS NO ICE
2208     CML replaced "/.0001 _d 0" by "*1. _d 4", 1e-4 is probably
2209     CML meant to be something like a minimum thickness
2210     AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),HEFF(I,J,bi,bj)*1. _d 4)
2211     ENDDO
2212     ENDDO
2213    
2214     #ifdef ALLOW_AUTODIFF_TAMC
2215     CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
2216     CADJ & key = iicekey, byte = isbyte
2217     #endif /* ALLOW_AUTODIFF_TAMC */
2218     DO J=1,sNy
2219     DO I=1,sNx
2220     C NOW TRUNCATE AREA
2221     AREA(I,J,bi,bj)=MIN(ONE,AREA(I,J,bi,bj))
2222     ENDDO
2223     ENDDO
2224    
2225     #ifdef ALLOW_AUTODIFF_TAMC
2226     CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
2227     CADJ & key = iicekey, byte = isbyte
2228     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
2229     CADJ & key = iicekey, byte = isbyte
2230     #endif /* ALLOW_AUTODIFF_TAMC */
2231     DO J=1,sNy
2232     DO I=1,sNx
2233     AREA(I,J,bi,bj) = MAX(ZERO,AREA(I,J,bi,bj))
2234     HSNOW(I,J,bi,bj) = MAX(ZERO,HSNOW(I,J,bi,bj))
2235     AREA(I,J,bi,bj) = AREA(I,J,bi,bj)*HEFFM(I,J,bi,bj)
2236     HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)*HEFFM(I,J,bi,bj)
2237     #ifdef SEAICE_CAP_HEFF
2238     C This is not energy conserving, but at least it conserves fresh water
2239     tmpscal0 = -MAX(HEFF(I,J,bi,bj)-MAX_HEFF,0. _d 0)
2240     d_HEFFbyNeg(I,J) = d_HEFFbyNeg(I,J) + tmpscal0
2241     HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal0
2242     #endif /* SEAICE_CAP_HEFF */
2243     HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj)
2244     ENDDO
2245     ENDDO
2246    
2247     C convert snow to ice if submerged.
2248     C =================================
2249    
2250     IF ( SEAICEuseFlooding ) THEN
2251     DO J=1,sNy
2252     DO I=1,sNx
2253     tmpscal0 = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
2254     & +HEFF(I,J,bi,bj)*SEAICE_rhoIce)*recip_rhoConst
2255     tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFF(I,J,bi,bj))
2256     d_HEFFbyFLOODING(I,J)=tmpscal1
2257     HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)+d_HEFFbyFLOODING(I,J)
2258     HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-
2259     & d_HEFFbyFLOODING(I,J)*ICE2SNOW
2260     ENDDO
2261     ENDDO
2262     ENDIF
2263    
2264     #endif /* SEAICE_GROWTH_LEGACY */
2265    
2266     #ifdef ALLOW_SITRACER
2267     DO J=1,sNy
2268     DO I=1,sNx
2269     c needs to be here to allow use also with LEGACY branch
2270     SItrHEFF(I,J,bi,bj,5)=HEFF(I,J,bi,bj)
2271     ENDDO
2272     ENDDO
2273     #endif /* ALLOW_SITRACER */
2274    
2275     C ===================================================================
2276     C ==============PART 7: determine ocean model forcing================
2277     C ===================================================================
2278    
2279     C compute net heat flux leaving/entering the ocean,
2280     C accounting for the part used in melt/freeze processes
2281     C =====================================================
2282    
2283 torge 1.8 #ifdef SEAICE_ITD
2284     C compute total of "mult" fluxes for ocean forcing
2285     DO J=1,sNy
2286     DO I=1,sNx
2287     a_QbyATM_cover(I,J) = 0.0 _d 0
2288     r_QbyATM_cover(I,J) = 0.0 _d 0
2289     a_QSWbyATM_cover(I,J) = 0.0 _d 0
2290     r_FWbySublim(I,J) = 0.0 _d 0
2291     ENDDO
2292     ENDDO
2293     DO IT=1,nITD
2294     DO J=1,sNy
2295     DO I=1,sNx
2296     cToM if fluxes in W/m^2 then
2297     c a_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)
2298     c & + a_QbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
2299     c r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)
2300     c & + r_QbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
2301     c a_QSWbyATM_cover(I,J)=a_QSWbyATM_cover(I,J)
2302     c & + a_QSWbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
2303     c r_FWbySublim(I,J)=r_FWbySublim(I,J)
2304     c & + r_FWbySublimMult(I,J,IT) * areaFracFactor(I,J,IT)
2305     cToM if fluxes in effective ice meters, i.e. ice volume per area, then
2306     a_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)
2307     & + a_QbyATMmult_cover(I,J,IT)
2308     r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)
2309     & + r_QbyATMmult_cover(I,J,IT)
2310     a_QSWbyATM_cover(I,J)=a_QSWbyATM_cover(I,J)
2311     & + a_QSWbyATMmult_cover(I,J,IT)
2312     r_FWbySublim(I,J)=r_FWbySublim(I,J)
2313     & + r_FWbySublimMult(I,J,IT)
2314     ENDDO
2315     ENDDO
2316     ENDDO
2317     #endif
2318    
2319 dimitri 1.1 #ifdef ALLOW_AUTODIFF_TAMC
2320     CADJ STORE d_hsnwbyneg = comlev1_bibj,key=iicekey,byte=isbyte
2321     CADJ STORE d_hsnwbyocnonsnw = comlev1_bibj,key=iicekey,byte=isbyte
2322     #endif /* ALLOW_AUTODIFF_TAMC */
2323    
2324     DO J=1,sNy
2325     DO I=1,sNx
2326     QNET(I,J,bi,bj) = r_QbyATM_cover(I,J) + r_QbyATM_open(I,J)
2327     #ifndef SEAICE_GROWTH_LEGACY
2328     C in principle a_QSWbyATM_cover should always be included here, however
2329     C for backward compatibility it is left out of the LEGACY branch
2330     & + a_QSWbyATM_cover(I,J)
2331     #endif /* SEAICE_GROWTH_LEGACY */
2332     & - ( d_HEFFbyOCNonICE(I,J) +
2333     & d_HSNWbyOCNonSNW(I,J)*SNOW2ICE +
2334     & d_HEFFbyNEG(I,J) +
2335     #ifdef SEAICE_ALLOW_AREA_RELAXATION
2336     & d_HEFFbyRLX(I,J) +
2337     #endif
2338     & d_HSNWbyNEG(I,J)*SNOW2ICE )
2339     & * maskC(I,J,kSurface,bi,bj)
2340     QSW(I,J,bi,bj) = a_QSWbyATM_cover(I,J) + a_QSWbyATM_open(I,J)
2341     ENDDO
2342     ENDDO
2343    
2344     C switch heat fluxes from 'effective' ice meters to W/m2
2345     C ======================================================
2346    
2347     DO J=1,sNy
2348     DO I=1,sNx
2349     QNET(I,J,bi,bj) = QNET(I,J,bi,bj)*convertHI2Q
2350     QSW(I,J,bi,bj) = QSW(I,J,bi,bj)*convertHI2Q
2351     ENDDO
2352     ENDDO
2353    
2354     #ifndef SEAICE_DISABLE_HEATCONSFIX
2355     C treat advective heat flux by ocean to ice water exchange (at 0decC)
2356     C ===================================================================
2357     # ifdef ALLOW_AUTODIFF_TAMC
2358     CADJ STORE d_HEFFbyNEG = comlev1_bibj,key=iicekey,byte=isbyte
2359     CADJ STORE d_HEFFbyOCNonICE = comlev1_bibj,key=iicekey,byte=isbyte
2360     CADJ STORE d_HEFFbyATMonOCN = comlev1_bibj,key=iicekey,byte=isbyte
2361     CADJ STORE d_HSNWbyNEG = comlev1_bibj,key=iicekey,byte=isbyte
2362     CADJ STORE d_HSNWbyOCNonSNW = comlev1_bibj,key=iicekey,byte=isbyte
2363     CADJ STORE d_HSNWbyATMonSNW = comlev1_bibj,key=iicekey,byte=isbyte
2364     CADJ STORE theta(:,:,kSurface,bi,bj) = comlev1_bibj,
2365     CADJ & key = iicekey, byte = isbyte
2366     # endif /* ALLOW_AUTODIFF_TAMC */
2367     IF ( SEAICEheatConsFix ) THEN
2368     c Unlike for evap and precip, the temperature of gained/lost
2369     c ocean liquid water due to melt/freeze of solid water cannot be chosen
2370     c to be e.g. the ocean SST. It must be done at 0degC. The fix below anticipates
2371     c on external_forcing_surf.F and applies the correction to QNET.
2372     IF ((convertFW2Salt.EQ.-1.).OR.(temp_EvPrRn.EQ.UNSET_RL)) THEN
2373     c I leave alone the exotic case when onvertFW2Salt.NE.-1 and temp_EvPrRn.NE.UNSET_RL and
2374     c the small error of the synchronous time stepping case (see external_forcing_surf.F).
2375     DO J=1,sNy
2376     DO I=1,sNx
2377     #ifdef ALLOW_DIAGNOSTICS
2378     c store unaltered QNET for diagnostic purposes
2379     DIAGarrayA(I,J)=QNET(I,J,bi,bj)
2380     #endif
2381     c compute the ocean water going to ice/snow, in precip units
2382     tmpscal3=rhoConstFresh*maskC(I,J,kSurface,bi,bj)*
2383     & ( d_HSNWbyATMonSNW(I,J)*SNOW2ICE
2384     & + d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
2385     & + d_HEFFbyOCNonICE(I,J) + d_HEFFbyATMonOCN(I,J)
2386     & + d_HEFFbyNEG(I,J) + d_HSNWbyNEG(I,J)*SNOW2ICE )
2387     & * convertHI2PRECIP
2388     c factor in the heat content that external_forcing_surf.F
2389     c will associate with EMPMR, and remove it from QNET, so that
2390     c melt/freez water is in effect consistently gained/lost at 0degC
2391     IF (temp_EvPrRn.NE.UNSET_RL) THEN
2392     QNET(I,J,bi,bj)=QNET(I,J,bi,bj) - tmpscal3*
2393     & HeatCapacity_Cp * temp_EvPrRn
2394     ELSE
2395     QNET(I,J,bi,bj)=QNET(I,J,bi,bj) - tmpscal3*
2396     & HeatCapacity_Cp * theta(I,J,kSurface,bi,bj)
2397     ENDIF
2398     #ifdef ALLOW_DIAGNOSTICS
2399     c back out the eventual TFLUX adjustement and fill diag
2400     DIAGarrayA(I,J)=QNET(I,J,bi,bj)-DIAGarrayA(I,J)
2401     #endif
2402     ENDDO
2403     ENDDO
2404     ENDIF
2405     #ifdef ALLOW_DIAGNOSTICS
2406     CALL DIAGNOSTICS_FILL(DIAGarrayA,
2407     & 'SIaaflux',0,1,3,bi,bj,myThid)
2408     #endif
2409     ENDIF
2410     #endif /* ndef SEAICE_DISABLE_HEATCONSFIX */
2411    
2412     C compute net fresh water flux leaving/entering
2413     C the ocean, accounting for fresh/salt water stocks.
2414     C ==================================================
2415    
2416     #ifdef ALLOW_ATM_TEMP
2417     DO J=1,sNy
2418     DO I=1,sNx
2419     tmpscal1= d_HSNWbyATMonSNW(I,J)*SNOW2ICE
2420     & +d_HFRWbyRAIN(I,J)
2421     & +d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
2422     & +d_HEFFbyOCNonICE(I,J)
2423     & +d_HEFFbyATMonOCN(I,J)
2424     & +d_HEFFbyNEG(I,J)
2425     #ifdef SEAICE_ALLOW_AREA_RELAXATION
2426     & +d_HEFFbyRLX(I,J)
2427     #endif
2428     & +d_HSNWbyNEG(I,J)*SNOW2ICE
2429     C If r_FWbySublim>0, then it is evaporated from ocean.
2430     & +r_FWbySublim(I,J)
2431     EmPmR(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
2432     & ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )
2433     & * ( ONE - AREApreTH(I,J) )
2434     #ifdef ALLOW_RUNOFF
2435     & - RUNOFF(I,J,bi,bj)
2436     #endif /* ALLOW_RUNOFF */
2437     & + tmpscal1*convertHI2PRECIP
2438     & )*rhoConstFresh
2439     ENDDO
2440     ENDDO
2441    
2442     #ifdef ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION
2443     C--
2444     DO J=1,sNy
2445     DO I=1,sNx
2446     frWtrAtm(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
2447     & PRECIP(I,J,bi,bj)
2448     & - EVAP(I,J,bi,bj)*( ONE - AREApreTH(I,J) )
2449     # ifdef ALLOW_RUNOFF
2450     & + RUNOFF(I,J,bi,bj)
2451     # endif /* ALLOW_RUNOFF */
2452     & )*rhoConstFresh
2453     # ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
2454     & - a_FWbySublim(I,J)*AREApreTH(I,J)
2455     # endif /* SEAICE_ADD_SUBLIMATION_TO_FWBUDGET */
2456     ENDDO
2457     ENDDO
2458     C--
2459     #else /* ndef ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION */
2460     C--
2461     # ifdef ALLOW_MEAN_SFLUX_COST_CONTRIBUTION
2462     DO J=1,sNy
2463     DO I=1,sNx
2464     frWtrAtm(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
2465     & PRECIP(I,J,bi,bj)
2466     & - EVAP(I,J,bi,bj)
2467     & *( ONE - AREApreTH(I,J) )
2468     # ifdef ALLOW_RUNOFF
2469     & + RUNOFF(I,J,bi,bj)
2470     # endif /* ALLOW_RUNOFF */
2471     & )*rhoConstFresh
2472     & - a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm
2473     ENDDO
2474     ENDDO
2475     # endif
2476     C--
2477     #endif /* ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION */
2478    
2479     #endif /* ALLOW_ATM_TEMP */
2480    
2481     #ifdef SEAICE_DEBUG
2482     CALL PLOT_FIELD_XYRL( QSW,'Current QSW ', myIter, myThid )
2483     CALL PLOT_FIELD_XYRL( QNET,'Current QNET ', myIter, myThid )
2484     CALL PLOT_FIELD_XYRL( EmPmR,'Current EmPmR ', myIter, myThid )
2485     #endif /* SEAICE_DEBUG */
2486    
2487     C Sea Ice Load on the sea surface.
2488     C =================================
2489    
2490     #ifdef ALLOW_AUTODIFF_TAMC
2491     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
2492     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
2493     #endif /* ALLOW_AUTODIFF_TAMC */
2494    
2495     IF ( useRealFreshWaterFlux ) THEN
2496     DO J=1,sNy
2497     DO I=1,sNx
2498     #ifdef SEAICE_CAP_ICELOAD
2499     tmpscal1 = HEFF(I,J,bi,bj)*SEAICE_rhoIce
2500     & + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
2501     tmpscal2 = MIN(tmpscal1,heffTooHeavy*rhoConst)
2502     #else
2503     tmpscal2 = HEFF(I,J,bi,bj)*SEAICE_rhoIce
2504     & + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
2505     #endif
2506     sIceLoad(i,j,bi,bj) = tmpscal2
2507     ENDDO
2508     ENDDO
2509     ENDIF
2510    
2511     C ===================================================================
2512     C ======================PART 8: diagnostics==========================
2513     C ===================================================================
2514    
2515     #ifdef ALLOW_DIAGNOSTICS
2516     IF ( useDiagnostics ) THEN
2517     tmpscal1=1. _d 0 * recip_deltaTtherm
2518     CALL DIAGNOSTICS_SCALE_FILL(a_QbyATM_cover,
2519     & tmpscal1,1,'SIaQbATC',0,1,3,bi,bj,myThid)
2520     CALL DIAGNOSTICS_SCALE_FILL(a_QbyATM_open,
2521     & tmpscal1,1,'SIaQbATO',0,1,3,bi,bj,myThid)
2522     CALL DIAGNOSTICS_SCALE_FILL(a_QbyOCN,
2523     & tmpscal1,1,'SIaQbOCN',0,1,3,bi,bj,myThid)
2524     CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyOCNonICE,
2525     & tmpscal1,1,'SIdHbOCN',0,1,3,bi,bj,myThid)
2526     CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyATMonOCN_cover,
2527     & tmpscal1,1,'SIdHbATC',0,1,3,bi,bj,myThid)
2528     CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyATMonOCN_open,
2529     & tmpscal1,1,'SIdHbATO',0,1,3,bi,bj,myThid)
2530     CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyFLOODING,
2531     & tmpscal1,1,'SIdHbFLO',0,1,3,bi,bj,myThid)
2532     CALL DIAGNOSTICS_SCALE_FILL(d_HSNWbyOCNonSNW,
2533     & tmpscal1,1,'SIdSbOCN',0,1,3,bi,bj,myThid)
2534     CALL DIAGNOSTICS_SCALE_FILL(d_HSNWbyATMonSNW,
2535     & tmpscal1,1,'SIdSbATC',0,1,3,bi,bj,myThid)
2536     CALL DIAGNOSTICS_SCALE_FILL(d_AREAbyATM,
2537     & tmpscal1,1,'SIdAbATO',0,1,3,bi,bj,myThid)
2538     CALL DIAGNOSTICS_SCALE_FILL(d_AREAbyICE,
2539     & tmpscal1,1,'SIdAbATC',0,1,3,bi,bj,myThid)
2540     CALL DIAGNOSTICS_SCALE_FILL(d_AREAbyOCN,
2541     & tmpscal1,1,'SIdAbOCN',0,1,3,bi,bj,myThid)
2542     CALL DIAGNOSTICS_SCALE_FILL(r_QbyATM_open,
2543     & convertHI2Q,1, 'SIqneto ',0,1,3,bi,bj,myThid)
2544     CALL DIAGNOSTICS_SCALE_FILL(r_QbyATM_cover,
2545     & convertHI2Q,1, 'SIqneti ',0,1,3,bi,bj,myThid)
2546     C three that actually need intermediate storage
2547     DO J=1,sNy
2548     DO I=1,sNx
2549     DIAGarrayA(I,J) = maskC(I,J,kSurface,bi,bj)
2550     & * d_HSNWbyRAIN(I,J)*SEAICE_rhoSnow*recip_deltaTtherm
2551     DIAGarrayB(I,J) = AREA(I,J,bi,bj)-AREApreTH(I,J)
2552     ENDDO
2553     ENDDO
2554     CALL DIAGNOSTICS_FILL(DIAGarrayA,
2555     & 'SIsnPrcp',0,1,3,bi,bj,myThid)
2556     CALL DIAGNOSTICS_SCALE_FILL(DIAGarrayB,
2557     & tmpscal1,1,'SIdA ',0,1,3,bi,bj,myThid)
2558     #ifdef ALLOW_ATM_TEMP
2559     DO J=1,sNy
2560     DO I=1,sNx
2561     CML If I consider the atmosphere above the ice, the surface flux
2562     CML which is relevant for the air temperature dT/dt Eq
2563     CML accounts for sensible and radiation (with different treatment
2564     CML according to wave-length) fluxes but not for "latent heat flux",
2565     CML since it does not contribute to heating the air.
2566     CML So this diagnostic is only good for heat budget calculations within
2567     CML the ice-ocean system.
2568     DIAGarrayA(I,J) = maskC(I,J,kSurface,bi,bj)*convertHI2Q*(
2569     #ifndef SEAICE_GROWTH_LEGACY
2570     & a_QSWbyATM_cover(I,J) +
2571     #endif /* SEAICE_GROWTH_LEGACY */
2572     & a_QbyATM_cover(I,J) + a_QbyATM_open(I,J) )
2573     C
2574     DIAGarrayB(I,J) = maskC(I,J,kSurface,bi,bj) *
2575     & a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm
2576     C
2577     DIAGarrayC(I,J) = maskC(I,J,kSurface,bi,bj)*(
2578     & PRECIP(I,J,bi,bj)
2579     & - EVAP(I,J,bi,bj)*( ONE - AREApreTH(I,J) )
2580     #ifdef ALLOW_RUNOFF
2581     & + RUNOFF(I,J,bi,bj)
2582     #endif /* ALLOW_RUNOFF */
2583     & )*rhoConstFresh
2584     & - a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm
2585     ENDDO
2586     ENDDO
2587     CALL DIAGNOSTICS_FILL(DIAGarrayA,
2588     & 'SIatmQnt',0,1,3,bi,bj,myThid)
2589     CALL DIAGNOSTICS_FILL(DIAGarrayB,
2590     & 'SIfwSubl',0,1,3,bi,bj,myThid)
2591     CALL DIAGNOSTICS_FILL(DIAGarrayC,
2592     & 'SIatmFW ',0,1,3,bi,bj,myThid)
2593     C
2594     DO J=1,sNy
2595     DO I=1,sNx
2596     C the actual Freshwater flux of sublimated ice, >0 decreases ice
2597     DIAGarrayA(I,J) = maskC(I,J,kSurface,bi,bj)
2598     & * (a_FWbySublim(I,J)-r_FWbySublim(I,J))
2599     & * SEAICE_rhoIce * recip_deltaTtherm
2600     c the residual Freshwater flux of sublimated ice
2601     DIAGarrayC(I,J) = maskC(I,J,kSurface,bi,bj)
2602     & * r_FWbySublim(I,J)
2603     & * SEAICE_rhoIce * recip_deltaTtherm
2604     C the latent heat flux
2605     tmpscal1= EVAP(I,J,bi,bj)*( ONE - AREApreTH(I,J) )
2606     & + r_FWbySublim(I,J)*convertHI2PRECIP
2607     tmpscal2= ( a_FWbySublim(I,J)-r_FWbySublim(I,J) )
2608     & * convertHI2PRECIP
2609     tmpscal3= SEAICE_lhEvap+SEAICE_lhFusion
2610     DIAGarrayB(I,J) = -maskC(I,J,kSurface,bi,bj)*rhoConstFresh
2611     & * ( tmpscal1*SEAICE_lhEvap + tmpscal2*tmpscal3 )
2612     ENDDO
2613     ENDDO
2614     CALL DIAGNOSTICS_FILL(DIAGarrayA,'SIacSubl',0,1,3,bi,bj,myThid)
2615     CALL DIAGNOSTICS_FILL(DIAGarrayC,'SIrsSubl',0,1,3,bi,bj,myThid)
2616     CALL DIAGNOSTICS_FILL(DIAGarrayB,'SIhl ',0,1,3,bi,bj,myThid)
2617    
2618     DO J=1,sNy
2619     DO I=1,sNx
2620     c compute ice/snow water going to atm, in precip units
2621     tmpscal1 = rhoConstFresh*maskC(I,J,kSurface,bi,bj)
2622     & * convertHI2PRECIP * ( - d_HSNWbyRAIN(I,J)*SNOW2ICE
2623     & + a_FWbySublim(I,J) - r_FWbySublim(I,J) )
2624     c compute ocean water going to atm, in precip units
2625     tmpscal2=rhoConstFresh*maskC(I,J,kSurface,bi,bj)*
2626     & ( ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )
2627     & * ( ONE - AREApreTH(I,J) )
2628     #ifdef ALLOW_RUNOFF
2629     & - RUNOFF(I,J,bi,bj)
2630     #endif /* ALLOW_RUNOFF */
2631     & + ( d_HFRWbyRAIN(I,J) + r_FWbySublim(I,J) )
2632     & *convertHI2PRECIP )
2633     c factor in the advected specific energy (referenced to 0 for 0deC liquid water)
2634     tmpscal1= - tmpscal1*
2635     & ( -SEAICE_lhFusion + HeatCapacity_Cp * ZERO )
2636     IF (temp_EvPrRn.NE.UNSET_RL) THEN
2637     tmpscal2= - tmpscal2*
2638     & ( ZERO + HeatCapacity_Cp * temp_EvPrRn )
2639     ELSE
2640     tmpscal2= - tmpscal2*
2641     & ( ZERO + HeatCapacity_Cp * theta(I,J,kSurface,bi,bj) )
2642     ENDIF
2643     c add to SIatmQnt, leading to SItflux, which is analogous to TFLUX
2644     DIAGarrayA(I,J)=maskC(I,J,kSurface,bi,bj)*convertHI2Q*(
2645     #ifndef SEAICE_GROWTH_LEGACY
2646     & a_QSWbyATM_cover(I,J) +
2647     #endif
2648     & a_QbyATM_cover(I,J) + a_QbyATM_open(I,J) )
2649     & -tmpscal1-tmpscal2
2650     ENDDO
2651     ENDDO
2652     CALL DIAGNOSTICS_FILL(DIAGarrayA,
2653     & 'SItflux ',0,1,3,bi,bj,myThid)
2654     #endif /* ALLOW_ATM_TEMP */
2655    
2656     ENDIF
2657     #endif /* ALLOW_DIAGNOSTICS */
2658    
2659     C close bi,bj loops
2660     ENDDO
2661     ENDDO
2662    
2663     RETURN
2664     END
  ViewVC Help
Powered by ViewVC 1.1.22