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