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