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

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

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

revision 1.3 by torge, Wed Sep 26 17:50:17 2012 UTC revision 1.14 by torge, Mon Dec 10 22:19:49 2012 UTC
# Line 58  C     !FUNCTIONS: Line 58  C     !FUNCTIONS:
58    
59  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
60  C     === Local variables ===  C     === Local variables ===
 c ToM<<< debug seaice_growth  
 C     msgBuf      :: Informational/error message buffer  
       CHARACTER*(MAX_LEN_MBUF) msgBuf  
 c ToM>>>  
61  C  C
62  C unit/sign convention:  C unit/sign convention:
63  C    Within the thermodynamic computation all stocks, except HSNOW,  C    Within the thermodynamic computation all stocks, except HSNOW,
# Line 95  C     i,j,bi,bj :: Loop counters Line 91  C     i,j,bi,bj :: Loop counters
91        INTEGER i, j, bi, bj        INTEGER i, j, bi, bj
92  C     number of surface interface layer  C     number of surface interface layer
93        INTEGER kSurface        INTEGER kSurface
94    C     IT :: ice thickness category index (MULTICATEGORIES and ITD code)
95          INTEGER IT
96    C     msgBuf      :: Informational/error message buffer
97    #ifdef ALLOW_BALANCE_FLUXES
98          CHARACTER*(MAX_LEN_MBUF) msgBuf
99    #elif (defined (SEAICE_DEBUG))
100          CHARACTER*(MAX_LEN_MBUF) msgBuf
101          CHARACTER*12 msgBufForm
102    #endif
103  C     constants  C     constants
104          _RL pFac
105        _RL tempFrz, ICE2SNOW, SNOW2ICE        _RL tempFrz, ICE2SNOW, SNOW2ICE
106        _RL QI, QS, recip_QI        _RL QI, QS, recip_QI
107          _RL lhSublim
108    
109    C conversion factors to go from Q (W/m2) to HEFF (ice meters)
110          _RL convertQ2HI, convertHI2Q
111    C conversion factors to go from precip (m/s) unit to HEFF (ice meters)
112          _RL convertPRECIP2HI, convertHI2PRECIP
113    C     Factor by which we increase the upper ocean friction velocity (u*) when
114    C     ice is absent in a grid cell  (dimensionless)
115          _RL MixedLayerTurbulenceFactor
116    
117    C     wind speed square
118          _RL SPEED_SQ
119    
120    C     Regularization values squared
121          _RL area_reg_sq, hice_reg_sq
122    C     pathological cases thresholds
123          _RL heffTooHeavy
124    
125    C     Helper variables: reciprocal of some constants
126          _RL recip_multDim
127          _RL recip_deltaTtherm
128          _RL recip_rhoIce
129    C     local value (=1/HO or 1/HO_south)
130          _RL recip_HO
131    C     local value (=1/ice thickness)
132          _RL recip_HH
133    C     facilitate multi-category snow implementation
134          _RL pFacSnow
135    
136    C     temporary variables available for the various computations
137          _RL tmpscal0, tmpscal1, tmpscal2, tmpscal3, tmpscal4
138    #ifdef SEAICE_ITD
139          _RL tmpscal1itd(1:sNx,1:sNy), tmpscal2itd(1:sNx,1:sNy)
140          _RL tmpscal3itd(1:sNx,1:sNy)
141    #endif
142    
143    #ifdef ALLOW_SITRACER
144          INTEGER iTr
145    #ifdef ALLOW_DIAGNOSTICS
146          CHARACTER*8   diagName
147    #endif
148    #endif /* ALLOW_SITRACER */
149    #ifdef ALLOW_AUTODIFF_TAMC
150          INTEGER ilockey
151    #endif
152    
153    C==   local arrays ==
154  C--   TmixLoc        :: ocean surface/mixed-layer temperature (in K)  C--   TmixLoc        :: ocean surface/mixed-layer temperature (in K)
155        _RL TmixLoc       (1:sNx,1:sNy)        _RL TmixLoc       (1:sNx,1:sNy)
156    
157    C     actual ice thickness (with upper and lower limit)
158          _RL heffActual          (1:sNx,1:sNy)
159    C     actual snow thickness
160          _RL hsnowActual         (1:sNx,1:sNy)
161    C     actual ice thickness (with lower limit only) Reciprocal
162          _RL recip_heffActual    (1:sNx,1:sNy)
163    
164    C     AREA_PRE :: hold sea-ice fraction field before any seaice-thermo update
165          _RL AREApreTH           (1:sNx,1:sNy)
166          _RL HEFFpreTH           (1:sNx,1:sNy)
167          _RL HSNWpreTH           (1:sNx,1:sNy)
168    #ifdef SEAICE_ITD
169          _RL AREAITDpreTH        (1:sNx,1:sNy,1:nITD)
170          _RL HEFFITDpreTH        (1:sNx,1:sNy,1:nITD)
171          _RL HSNWITDpreTH        (1:sNx,1:sNy,1:nITD)
172          _RL areaFracFactor      (1:sNx,1:sNy,1:nITD)
173          _RL leadIceThickMin
174    #endif
175    
176    C     wind speed
177          _RL UG                  (1:sNx,1:sNy)
178    
179    C     temporary variables available for the various computations
180          _RL tmparr1             (1:sNx,1:sNy)
181    #ifdef SEAICE_VARIABLE_SALINITY
182          _RL saltFluxAdjust      (1:sNx,1:sNy)
183    #endif
184    
185          _RL ticeInMult          (1:sNx,1:sNy,MULTDIM)
186          _RL ticeOutMult         (1:sNx,1:sNy,MULTDIM)
187          _RL heffActualMult      (1:sNx,1:sNy,MULTDIM)
188          _RL hsnowActualMult     (1:sNx,1:sNy,MULTDIM)
189    #ifdef SEAICE_ITD
190          _RL recip_heffActualMult(1:sNx,1:sNy,MULTDIM)
191    #endif
192          _RL a_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
193          _RL a_QSWbyATMmult_cover(1:sNx,1:sNy,MULTDIM)
194          _RL a_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
195    #ifdef SEAICE_ITD
196          _RL r_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
197          _RL r_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
198    #endif
199    
200  C     a_QbyATM_cover :: available heat (in W/m^2) due to the interaction of  C     a_QbyATM_cover :: available heat (in W/m^2) due to the interaction of
201  C             the atmosphere and the ocean surface - for ice covered water  C             the atmosphere and the ocean surface - for ice covered water
202  C     a_QbyATM_open  :: same but for open water  C     a_QbyATM_open  :: same but for open water
# Line 115  C     a_QSWbyATM_open   - short wave hea Line 210  C     a_QSWbyATM_open   - short wave hea
210  C     a_QSWbyATM_cover  - short wave heat flux under ice in W/m^2  C     a_QSWbyATM_cover  - short wave heat flux under ice in W/m^2
211        _RL a_QSWbyATM_open     (1:sNx,1:sNy)        _RL a_QSWbyATM_open     (1:sNx,1:sNy)
212        _RL a_QSWbyATM_cover    (1:sNx,1:sNy)        _RL a_QSWbyATM_cover    (1:sNx,1:sNy)
213  C     a_QbyOCN :: available heat (in in W/m^2) due to the  C     a_QbyOCN :: available heat (in W/m^2) due to the
214  C             interaction of the ice pack and the ocean surface  C             interaction of the ice pack and the ocean surface
215  C     r_QbyOCN :: residual of a_QbyOCN after freezing/melting  C     r_QbyOCN :: residual of a_QbyOCN after freezing/melting
216  C             processes have been accounted for  C             processes have been accounted for
217        _RL a_QbyOCN            (1:sNx,1:sNy)        _RL a_QbyOCN            (1:sNx,1:sNy)
218        _RL r_QbyOCN            (1:sNx,1:sNy)        _RL r_QbyOCN            (1:sNx,1:sNy)
219    
220  C conversion factors to go from Q (W/m2) to HEFF (ice meters)  C     The change of mean ice thickness due to out-of-bounds values following
221        _RL convertQ2HI, convertHI2Q  C     sea ice dyhnamics
 C conversion factors to go from precip (m/s) unit to HEFF (ice meters)  
       _RL convertPRECIP2HI, convertHI2PRECIP  
   
 #ifdef ALLOW_DIAGNOSTICS  
 C ICE/SNOW stocks tendencies associated with the various melt/freeze processes  
       _RL d_AREAbyATM         (1:sNx,1:sNy)  
       _RL d_AREAbyOCN         (1:sNx,1:sNy)  
       _RL d_AREAbyICE         (1:sNx,1:sNy)  
 #endif  
   
 #ifdef SEAICE_ALLOW_AREA_RELAXATION  
 C     ICE/SNOW stocks tendency associated with relaxation towards observation  
       _RL d_AREAbyRLX         (1:sNx,1:sNy)  
 c     The change of mean ice thickness due to relaxation  
       _RL d_HEFFbyRLX         (1:sNx,1:sNy)  
 #endif  
   
 #ifdef SEAICE_ITD  
 c     The change of mean ice area due to out-of-bounds values following  
 c     sea ice dynamics  
       _RL d_AREAbyNEG         (1:sNx,1:sNy)  
 #endif  
 c     The change of mean ice thickness due to out-of-bounds values following  
 c     sea ice dynamics  
222        _RL d_HEFFbyNEG         (1:sNx,1:sNy)        _RL d_HEFFbyNEG         (1:sNx,1:sNy)
223    
224  c     The change of mean ice thickness due to turbulent ocean-sea ice heat fluxes  C     The change of mean ice thickness due to turbulent ocean-sea ice heat fluxes
225        _RL d_HEFFbyOCNonICE    (1:sNx,1:sNy)        _RL d_HEFFbyOCNonICE    (1:sNx,1:sNy)
226    
227  c     The sum of mean ice thickness increments due to atmospheric fluxes over the open water  C     The sum of mean ice thickness increments due to atmospheric fluxes over
228  c     fraction and ice-covered fractions of the grid cell  C     the open water fraction and ice-covered fractions of the grid cell
229        _RL d_HEFFbyATMonOCN    (1:sNx,1:sNy)        _RL d_HEFFbyATMonOCN    (1:sNx,1:sNy)
230  c     The change of mean ice thickness due to flooding by snow  C     The change of mean ice thickness due to flooding by snow
231        _RL d_HEFFbyFLOODING    (1:sNx,1:sNy)        _RL d_HEFFbyFLOODING    (1:sNx,1:sNy)
232    
233  c     The mean ice thickness increments due to atmospheric fluxes over the open water  C     The mean ice thickness increments due to atmospheric fluxes over the open
234  c     fraction and ice-covered fractions of the grid cell, respectively  C     water fraction and ice-covered fractions of the grid cell, respectively
235        _RL d_HEFFbyATMonOCN_open(1:sNx,1:sNy)        _RL d_HEFFbyATMonOCN_open(1:sNx,1:sNy)
236        _RL d_HEFFbyATMonOCN_cover(1:sNx,1:sNy)        _RL d_HEFFbyATMonOCN_cover(1:sNx,1:sNy)
237    
# Line 170  c     fraction and ice-covered fractions Line 241  c     fraction and ice-covered fractions
241        _RL d_HSNWbyRAIN        (1:sNx,1:sNy)        _RL d_HSNWbyRAIN        (1:sNx,1:sNy)
242    
243        _RL d_HFRWbyRAIN        (1:sNx,1:sNy)        _RL d_HFRWbyRAIN        (1:sNx,1:sNy)
244  C  
245  C     a_FWbySublim :: fresh water flux implied by latent heat of  C     a_FWbySublim :: fresh water flux implied by latent heat of
246  C                     sublimation to atmosphere, same sign convention  C                     sublimation to atmosphere, same sign convention
247  C                     as EVAP (positive upward)  C                     as EVAP (positive upward)
# Line 179  C                     as EVAP (positive Line 250  C                     as EVAP (positive
250        _RL d_HEFFbySublim      (1:sNx,1:sNy)        _RL d_HEFFbySublim      (1:sNx,1:sNy)
251        _RL d_HSNWbySublim      (1:sNx,1:sNy)        _RL d_HSNWbySublim      (1:sNx,1:sNy)
252    
253  #if (defined(SEAICE_CAP_SUBLIM) || defined(SEAICE_ITD))  #ifdef SEAICE_CAP_SUBLIM
254  C     The latent heat flux which will sublimate all snow and ice  C     The latent heat flux which will sublimate all snow and ice
255  C     over one time step  C     over one time step
256        _RL latentHeatFluxMax   (1:sNx,1:sNy)        _RL latentHeatFluxMax   (1:sNx,1:sNy)
257        _RL latentHeatFluxMaxMult  (1:sNx,1:sNy,MULTDIM)        _RL latentHeatFluxMaxMult(1:sNx,1:sNy,MULTDIM)
 #endif  
   
 C     actual ice thickness (with upper and lower limit)  
       _RL heffActual          (1:sNx,1:sNy)  
 C     actual snow thickness  
       _RL hsnowActual         (1:sNx,1:sNy)  
 C     actual ice thickness (with lower limit only) Reciprocal  
       _RL recip_heffActual    (1:sNx,1:sNy)  
 C     local value (=1/HO or 1/HO_south)  
       _RL recip_HO  
 C     local value (=1/ice thickness)  
       _RL recip_HH  
   
 C     AREA_PRE :: hold sea-ice fraction field before any seaice-thermo update  
       _RL AREApreTH           (1:sNx,1:sNy)  
       _RL HEFFpreTH           (1:sNx,1:sNy)  
       _RL HSNWpreTH           (1:sNx,1:sNy)  
 #ifdef SEAICE_ITD  
       _RL AREAITDpreTH        (1:sNx,1:sNy,1:nITD)  
       _RL HEFFITDpreTH        (1:sNx,1:sNy,1:nITD)  
       _RL HSNWITDpreTH        (1:sNx,1:sNy,1:nITD)  
       _RL areaFracFactor      (1:sNx,1:sNy,1:nITD)  
       _RL heffFracFactor      (1:sNx,1:sNy,1:nITD)  
 #endif  
   
 C     wind speed  
       _RL UG                  (1:sNx,1:sNy)  
 #ifdef ALLOW_ATM_WIND  
       _RL SPEED_SQ  
258  #endif  #endif
259    
260  C     Regularization values squared  #ifdef EXF_ALLOW_SEAICE_RELAX
261        _RL area_reg_sq, hice_reg_sq  C     ICE/SNOW stocks tendency associated with relaxation towards observation
262          _RL d_AREAbyRLX         (1:sNx,1:sNy)
263  C     pathological cases thresholds  C     The change of mean ice thickness due to relaxation
264        _RL heffTooHeavy        _RL d_HEFFbyRLX         (1:sNx,1:sNy)
   
       _RL lhSublim  
   
 C temporary variables available for the various computations  
       _RL tmpscal0, tmpscal1, tmpscal2, tmpscal3, tmpscal4  
       _RL tmparr1             (1:sNx,1:sNy)  
   
 #ifdef SEAICE_VARIABLE_SALINITY  
       _RL saltFluxAdjust      (1:sNx,1:sNy)  
265  #endif  #endif
266    
       INTEGER ilockey  
       INTEGER it  
267  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
268        INTEGER K        _RL d_HEFFbySublim_ITD         (1:sNx,1:sNy,1:nITD)
269  #endif        _RL d_HSNWbySublim_ITD         (1:sNx,1:sNy,1:nITD)
270        _RL pFac        _RL d_HEFFbyOCNonICE_ITD       (1:sNx,1:sNy,1:nITD)
271        _RL ticeInMult          (1:sNx,1:sNy,MULTDIM)        _RL d_HSNWbyATMonSNW_ITD       (1:sNx,1:sNy,1:nITD)
272        _RL ticeOutMult         (1:sNx,1:sNy,MULTDIM)        _RL d_HEFFbyATMonOCN_ITD       (1:sNx,1:sNy,1:nITD)
273        _RL heffActualMult      (1:sNx,1:sNy,MULTDIM)        _RL d_HEFFbyATMonOCN_cover_ITD (1:sNx,1:sNy,1:nITD)
274  #ifdef SEAICE_ITD        _RL d_HEFFbyATMonOCN_open_ITD  (1:sNx,1:sNy,1:nITD)
275        _RL hsnowActualMult     (1:sNx,1:sNy,MULTDIM)        _RL d_HSNWbyRAIN_ITD           (1:sNx,1:sNy,1:nITD)
276        _RL recip_heffActualMult(1:sNx,1:sNy,MULTDIM)        _RL d_HSNWbyOCNonSNW_ITD       (1:sNx,1:sNy,1:nITD)
277          _RL d_HEFFbyFLOODING_ITD       (1:sNx,1:sNy,1:nITD)
278  #endif  #endif
       _RL a_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)  
       _RL a_QSWbyATMmult_cover(1:sNx,1:sNy,MULTDIM)  
       _RL a_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)  
 #ifdef SEAICE_ITD  
       _RL r_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)  
       _RL r_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)  
 #endif  
 C     Helper variables: reciprocal of some constants  
       _RL recip_multDim  
       _RL recip_deltaTtherm  
       _RL recip_rhoIce  
279    
 C     Factor by which we increase the upper ocean friction velocity (u*) when  
 C     ice is absent in a grid cell  (dimensionless)  
       _RL MixedLayerTurbulenceFactor  
   
 #ifdef ALLOW_SITRACER  
       INTEGER iTr  
       CHARACTER*8   diagName  
 #endif  
280  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
281  c     Helper variables for diagnostics  C ICE/SNOW stocks tendencies associated with the various melt/freeze processes
282          _RL d_AREAbyATM         (1:sNx,1:sNy)
283          _RL d_AREAbyOCN         (1:sNx,1:sNy)
284          _RL d_AREAbyICE         (1:sNx,1:sNy)
285    C     Helper variables for diagnostics
286        _RL DIAGarrayA    (1:sNx,1:sNy)        _RL DIAGarrayA    (1:sNx,1:sNy)
287        _RL DIAGarrayB    (1:sNx,1:sNy)        _RL DIAGarrayB    (1:sNx,1:sNy)
288        _RL DIAGarrayC    (1:sNx,1:sNy)        _RL DIAGarrayC    (1:sNx,1:sNy)
289        _RL DIAGarrayD    (1:sNx,1:sNy)        _RL DIAGarrayD    (1:sNx,1:sNy)
290  #endif  #endif /* ALLOW_DIAGNOSTICS */
291    
292          _RL SItflux     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
293          _RL SIatmQnt    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
294          _RL SIatmFW     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
295    #ifdef ALLOW_BALANCE_FLUXES
296          _RL FWFsiTile(nSx,nSy)
297          _RL FWFsiGlob
298          _RL HFsiTile(nSx,nSy)
299          _RL HFsiGlob
300          _RL FWF2HFsiTile(nSx,nSy)
301          _RL FWF2HFsiGlob
302    #endif
303    
304  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
305    
# Line 286  C ====================================== Line 314  C ======================================
314        ENDIF        ENDIF
315    
316  C     avoid unnecessary divisions in loops  C     avoid unnecessary divisions in loops
317  #ifdef SEAICE_ITD  c#ifdef SEAICE_ITD
318  CToM  SEAICE_multDim = nITD (see SEAICE_SIZE.h and seaice_readparms.F)  CToM this is now set by MULTDIM = nITD in SEAICE_SIZE.h
319  #endif  C    (see SEAICE_SIZE.h and seaice_readparms.F)
320    c     SEAICE_multDim = nITD
321    c#endif
322        recip_multDim        = SEAICE_multDim        recip_multDim        = SEAICE_multDim
323        recip_multDim        = ONE / recip_multDim        recip_multDim        = ONE / recip_multDim
324  C     above/below: double/single precision calculation of recip_multDim  C     above/below: double/single precision calculation of recip_multDim
# Line 339  C conversion factors to go from precip ( Line 369  C conversion factors to go from precip (
369       &                       + act4*max1*max2*max3       &                       + act4*max1*max2*max3
370  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
371    
   
372  C array initializations  C array initializations
373  C =====================  C =====================
374    
# Line 362  C ===================== Line 391  C =====================
391            d_AREAbyOCN(I,J)           = 0.0 _d 0            d_AREAbyOCN(I,J)           = 0.0 _d 0
392  #endif  #endif
393    
394  #ifdef SEAICE_ALLOW_AREA_RELAXATION  #ifdef EXF_ALLOW_SEAICE_RELAX
395            d_AREAbyRLX(I,J)           = 0.0 _d 0            d_AREAbyRLX(I,J)           = 0.0 _d 0
396            d_HEFFbyRLX(I,J)           = 0.0 _d 0            d_HEFFbyRLX(I,J)           = 0.0 _d 0
397  #endif  #endif
398    
 #ifdef SEAICE_ITD  
           d_AREAbyNEG(I,J)           = 0.0 _d 0  
 #endif  
399            d_HEFFbyNEG(I,J)           = 0.0 _d 0            d_HEFFbyNEG(I,J)           = 0.0 _d 0
400            d_HEFFbyOCNonICE(I,J)      = 0.0 _d 0            d_HEFFbyOCNonICE(I,J)      = 0.0 _d 0
401            d_HEFFbyATMonOCN(I,J)      = 0.0 _d 0            d_HEFFbyATMonOCN(I,J)      = 0.0 _d 0
# Line 389  C ===================== Line 415  C =====================
415  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
416            latentHeatFluxMax(I,J)     = 0.0 _d 0            latentHeatFluxMax(I,J)     = 0.0 _d 0
417  #endif  #endif
 c  
418            d_HFRWbyRAIN(I,J)          = 0.0 _d 0            d_HFRWbyRAIN(I,J)          = 0.0 _d 0
   
419            tmparr1(I,J)               = 0.0 _d 0            tmparr1(I,J)               = 0.0 _d 0
   
420  #ifdef SEAICE_VARIABLE_SALINITY  #ifdef SEAICE_VARIABLE_SALINITY
421            saltFluxAdjust(I,J)        = 0.0 _d 0            saltFluxAdjust(I,J)        = 0.0 _d 0
422  #endif  #endif
# Line 403  c Line 426  c
426              a_QbyATMmult_cover(I,J,IT)    = 0.0 _d 0              a_QbyATMmult_cover(I,J,IT)    = 0.0 _d 0
427              a_QSWbyATMmult_cover(I,J,IT)  = 0.0 _d 0              a_QSWbyATMmult_cover(I,J,IT)  = 0.0 _d 0
428              a_FWbySublimMult(I,J,IT)      = 0.0 _d 0              a_FWbySublimMult(I,J,IT)      = 0.0 _d 0
429  #ifdef SEAICE_ITD  #ifdef SEAICE_CAP_SUBLIM
             r_QbyATMmult_cover (I,J,IT)   = 0.0 _d 0  
             r_FWbySublimMult(I,J,IT)      = 0.0 _d 0  
 #endif  
 #if (defined(SEAICE_CAP_SUBLIM) || defined(SEAICE_ITD))  
430              latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0              latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0
431  #endif  #endif
432    #ifdef SEAICE_ITD
433                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    #endif
446            ENDDO            ENDDO
447           ENDDO           ENDDO
448          ENDDO          ENDDO
# Line 421  c Line 454  c
454          ENDDO          ENDDO
455  #endif  #endif
456    
   
457  C =====================================================================  C =====================================================================
458  C ===========PART 1: treat pathological cases (post advdiff)===========  C ===========PART 1: treat pathological cases (post advdiff)===========
459  C =====================================================================  C =====================================================================
# Line 444  C ====================================== Line 476  C ======================================
476           ENDDO           ENDDO
477          ENDDO          ENDDO
478  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
479          DO K=1,nITD          DO IT=1,nITD
480           DO J=1,sNy           DO J=1,sNy
481            DO I=1,sNx            DO I=1,sNx
482             HEFFITDpreTH(I,J,K)=HEFFITD(I,J,K,bi,bj)             HEFFITDpreTH(I,J,IT)=HEFFITD(I,J,IT,bi,bj)
483             HSNWITDpreTH(I,J,K)=HSNOWITD(I,J,K,bi,bj)             HSNWITDpreTH(I,J,IT)=HSNOWITD(I,J,IT,bi,bj)
484             AREAITDpreTH(I,J,K)=AREAITD(I,J,K,bi,bj)             AREAITDpreTH(I,J,IT)=AREAITD(I,J,IT,bi,bj)
485            ENDDO            ENDDO
486           ENDDO           ENDDO
487          ENDDO          ENDDO
# Line 462  Cgf no dependency through pathological c Line 494  Cgf no dependency through pathological c
494          IF ( SEAICEadjMODE.EQ.0 ) THEN          IF ( SEAICEadjMODE.EQ.0 ) THEN
495  #endif  #endif
496    
497  #ifdef SEAICE_ALLOW_AREA_RELAXATION  #ifdef EXF_ALLOW_SEAICE_RELAX
498  CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte  CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
499  CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte  CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
500  C 0) relax sea ice concentration towards observation  C 0) relax sea ice concentration towards observation
# Line 482  C           d_HEFFbyRLX(i,j) = 1. _d 1 * Line 514  C           d_HEFFbyRLX(i,j) = 1. _d 1 *
514              d_HEFFbyRLX(i,j) = 1. _d 1 * siEps              d_HEFFbyRLX(i,j) = 1. _d 1 * siEps
515             ENDIF             ENDIF
516  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
517              AREAITD(I,J,1,bi,bj) = AREAITD(I,J,1,bi,bj)              AREAITD(I,J,1,bi,bj) = AREAITD(I,J,1,bi,bj)
518       &                           +  d_AREAbyRLX(i,j)       &                           +  d_AREAbyRLX(i,j)
519              HEFFITD(I,J,1,bi,bj) = HEFFITD(I,J,1,bi,bj)              HEFFITD(I,J,1,bi,bj) = HEFFITD(I,J,1,bi,bj)
520       &                           +  d_HEFFbyRLX(i,j)       &                           +  d_HEFFbyRLX(i,j)
# Line 492  C           d_HEFFbyRLX(i,j) = 1. _d 1 * Line 524  C           d_HEFFbyRLX(i,j) = 1. _d 1 *
524           ENDDO           ENDDO
525          ENDDO          ENDDO
526          ENDIF          ENDIF
527  #endif /* SEAICE_ALLOW_AREA_RELAXATION */  #endif /* EXF_ALLOW_SEAICE_RELAX */
528    
529  C 1) treat the case of negative values:  C 1) treat the case of negative values:
530    
# Line 501  CADJ STORE heff(:,:,bi,bj)  = comlev1_bi Line 533  CADJ STORE heff(:,:,bi,bj)  = comlev1_bi
533  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
534  CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte  CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
535  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
536    #ifdef SEAICE_ITD
537            DO IT=1,nITD
538    #endif
539          DO J=1,sNy          DO J=1,sNy
540           DO I=1,sNx           DO I=1,sNx
541  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
542            DO K=1,nITD            tmpscal2=0. _d 0
543             tmpscal1=0. _d 0            tmpscal3=0. _d 0
544             tmpscal2=0. _d 0            tmpscal2=MAX(-HEFFITD(I,J,IT,bi,bj),0. _d 0)
545             tmpscal3=0. _d 0            HEFFITD(I,J,IT,bi,bj)=HEFFITD(I,J,IT,bi,bj)+tmpscal2
546             tmpscal2=MAX(-HEFFITD(I,J,K,bi,bj),0. _d 0)            d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
547             HEFFITD(I,J,K,bi,bj)=HEFFITD(I,J,K,bi,bj)+tmpscal2            tmpscal3=MAX(-HSNOWITD(I,J,IT,bi,bj),0. _d 0)
548             d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2            HSNOWITD(I,J,IT,bi,bj)=HSNOWITD(I,J,IT,bi,bj)+tmpscal3
549             tmpscal3=MAX(-HSNOWITD(I,J,K,bi,bj),0. _d 0)            d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
550             HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal3            AREAITD(I,J,IT,bi,bj)=MAX(AREAITD(I,J,IT,bi,bj),0. _d 0)
            d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3  
            tmpscal1=MAX(-AREAITD(I,J,K,bi,bj),0. _d 0)  
            AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj)+tmpscal1  
            d_AREAbyNEG(I,J)=d_AREAbyNEG(I,J)+tmpscal1  
           ENDDO  
551  CToM      AREA, HEFF, and HSNOW will be updated at end of PART 1  CToM      AREA, HEFF, and HSNOW will be updated at end of PART 1
552  C         by calling SEAICE_ITD_SUM  C         by calling SEAICE_ITD_SUM
553  #else  #else
554            d_HEFFbyNEG(I,J)=MAX(-HEFF(I,J,bi,bj),0. _d 0)            d_HEFFbyNEG(I,J)=MAX(-HEFF(I,J,bi,bj),0. _d 0)
           d_HSNWbyNEG(I,J)=MAX(-HSNOW(I,J,bi,bj),0. _d 0)  
           AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),0. _d 0)  
555            HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J)            HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J)
556              d_HSNWbyNEG(I,J)=MAX(-HSNOW(I,J,bi,bj),0. _d 0)
557            HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J)            HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J)
558              AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),0. _d 0)
559  #endif  #endif
560           ENDDO           ENDDO
561          ENDDO          ENDDO
562    #ifdef SEAICE_ITD
563            ENDDO
564    #endif
565    
566  C 1.25) treat the case of very thin ice:  C 1.25) treat the case of very thin ice:
567    
568  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
569  CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte  CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
570  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
         DO J=1,sNy  
          DO I=1,sNx  
571  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
572            DO K=1,nITD          DO IT=1,nITD
573  #endif  #endif
574            tmpscal2=0. _d 0          DO J=1,sNy
575            tmpscal3=0. _d 0           DO I=1,sNx
576              tmpscal2=0. _d 0
577              tmpscal3=0. _d 0
578  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
579             IF (HEFFITD(I,J,K,bi,bj).LE.siEps) THEN            IF (HEFFITD(I,J,IT,bi,bj).LE.siEps) THEN
580              tmpscal2=-HEFFITD(I,J,K,bi,bj)             tmpscal2=-HEFFITD(I,J,IT,bi,bj)
581              tmpscal3=-HSNOWITD(I,J,K,bi,bj)             tmpscal3=-HSNOWITD(I,J,IT,bi,bj)
582              TICES(I,J,K,bi,bj)=celsius2K             TICES(I,J,IT,bi,bj)=celsius2K
583  CToM        TICE will be updated at end of Part 1 together with AREA and HEFF  CToM       TICE will be updated at end of Part 1 together with AREA and HEFF
584             ENDIF            ENDIF
585             HEFFITD(I,J,K,bi,bj) =HEFFITD(I,J,K,bi,bj) +tmpscal2            HEFFITD(I,J,IT,bi,bj) =HEFFITD(I,J,IT,bi,bj) +tmpscal2
586             HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal3            HSNOWITD(I,J,IT,bi,bj)=HSNOWITD(I,J,IT,bi,bj)+tmpscal3
587  #else  #else
588            IF (HEFF(I,J,bi,bj).LE.siEps) THEN            IF (HEFF(I,J,bi,bj).LE.siEps) THEN
589             tmpscal2=-HEFF(I,J,bi,bj)             tmpscal2=-HEFF(I,J,bi,bj)
# Line 565  CToM        TICE will be updated at end Line 598  CToM        TICE will be updated at end
598  #endif  #endif
599            d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2            d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
600            d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3            d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
 #ifdef SEAICE_ITD  
           ENDDO  
 #endif  
601           ENDDO           ENDDO
602          ENDDO          ENDDO
603    #ifdef SEAICE_ITD
604            ENDDO
605    #endif
606    
607  C 1.5) treat the case of area but no ice/snow:  C 1.5) treat the case of area but no ice/snow:
608    
# Line 577  C 1.5) treat the case of area but no ice Line 610  C 1.5) treat the case of area but no ice
610  CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte  CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
611  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
612  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
613    #ifdef SEAICE_ITD
614            DO IT=1,nITD
615    #endif
616          DO J=1,sNy          DO J=1,sNy
617           DO I=1,sNx           DO I=1,sNx
618  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
619            DO K=1,nITD            IF ((HEFFITD(I,J,IT,bi,bj).EQ.0. _d 0).AND.
620             IF ((HEFFITD(i,j,k,bi,bj).EQ.0. _d 0).AND.       &        (HSNOWITD(I,J,IT,bi,bj).EQ.0. _d 0))
621       &         (HSNOWITD(i,j,k,bi,bj).EQ.0. _d 0))       &     AREAITD(I,J,IT,bi,bj)=0. _d 0
      &      AREAITD(I,J,K,bi,bj)=0. _d 0  
           ENDDO  
622  #else  #else
623            IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.            IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.
624       &        (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0       &        (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0
625  #endif  #endif
626           ENDDO           ENDDO
627          ENDDO          ENDDO
628    #ifdef SEAICE_ITD
629            ENDDO
630    #endif
631    
632  C 2) treat the case of very small area:  C 2) treat the case of very small area:
633    
# Line 598  C 2) treat the case of very small area: Line 635  C 2) treat the case of very small area:
635  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
636  CADJ STORE area(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte  CADJ STORE area(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
637  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
638    #ifdef SEAICE_ITD
639            DO IT=1,nITD
640    #endif
641          DO J=1,sNy          DO J=1,sNy
642           DO I=1,sNx           DO I=1,sNx
643  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
644            DO K=1,nITD            IF ((HEFFITD(I,J,IT,bi,bj).GT.0).OR.
645             IF ((HEFFITD(i,j,k,bi,bj).GT.0).OR.       &        (HSNOWITD(I,J,IT,bi,bj).GT.0)) THEN
646       &         (HSNOWITD(i,j,k,bi,bj).GT.0)) THEN  CToM       SEAICE_area_floor*nITD cannot be allowed to exceed 1
647  CToM        SEAICE_area_floor*nITD cannot be allowed to exceed 1  C          hence use SEAICE_area_floor devided by nITD
648  C           hence use SEAICE_area_floor devided by nITD  C          (or install a warning in e.g. seaice_readparms.F)
649  C           (or install a warning in e.g. seaice_readparms.F)             AREAITD(I,J,IT,bi,bj)=
650              AREAITD(I,J,K,bi,bj)=       &        MAX(AREAITD(I,J,IT,bi,bj),SEAICE_area_floor/float(nITD))
651       &         MAX(AREAITD(I,J,K,bi,bj),SEAICE_area_floor/float(nITD))            ENDIF
            ENDIF  
           ENDDO  
652  #else  #else
653            IF ((HEFF(i,j,bi,bj).GT.0).OR.(HSNOW(i,j,bi,bj).GT.0)) THEN            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)             AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),SEAICE_area_floor)
# Line 618  C           (or install a warning in e.g Line 656  C           (or install a warning in e.g
656  #endif  #endif
657           ENDDO           ENDDO
658          ENDDO          ENDDO
659    #ifdef SEAICE_ITD
660            ENDDO
661    #endif
662  #endif /* DISABLE_AREA_FLOOR */  #endif /* DISABLE_AREA_FLOOR */
663    
664  C 2.5) treat case of excessive ice cover, e.g., due to ridging:  C 2.5) treat case of excessive ice cover, e.g., due to ridging:
# Line 639  CADJ STORE area(:,:,bi,bj)  = comlev1_bi Line 680  CADJ STORE area(:,:,bi,bj)  = comlev1_bi
680            AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),SEAICE_area_max)            AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),SEAICE_area_max)
681           ENDDO           ENDDO
682          ENDDO          ENDDO
683  #endif /* SEAICE_ITD */  #endif /* notSEAICE_ITD */
684    
685  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
686  CToM catch up with items 1.25 and 2.5 involving category sums AREA and HEFF  CToM catch up with items 1.25 and 2.5 involving category sums AREA and HEFF
687  C    first, update AREA and HEFF:          DO IT=1,nITD
688          CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)           DO J=1,sNy
689  C            DO I=1,sNx
         DO J=1,sNy  
          DO I=1,sNx  
690  C    TICES was changed above (item 1.25), now update TICE as ice volume  C    TICES was changed above (item 1.25), now update TICE as ice volume
691  C     weighted average of TICES  C     weighted average of TICES
692            tmpscal1 = 0. _d 0  C    also compute total of AREAITD (needed for finishing item 2.5, see below)
693            tmpscal2 = 0. _d 0             IF (IT .eq. 1) THEN
694            DO K=1,nITD              tmpscal1itd(i,j) = 0. _d 0
695             tmpscal1=tmpscal1 + TICES(I,J,K,bi,bj)*HEFFITD(I,J,K,bi,bj)              tmpscal2itd(i,j) = 0. _d 0
696             tmpscal2=tmpscal2 + HEFFITD(I,J,K,bi,bj)              tmpscal3itd(i,j) = 0. _d 0
697            ENDDO             ENDIF
698            TICE(I,J,bi,bj)=tmpscal1/tmpscal2             tmpscal1itd(i,j)=tmpscal1itd(i,j) + TICES(I,J,IT,bi,bj)
699  C    lines of item 2.5 that were omitted:       &                                       * 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    C    lines of item 2.5 that were omitted:
705  C    in 2.5 these lines are executed before "ridging" is applied to AREA  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  C    hence we execute them here before SEAICE_ITD_REDIST is called
707  C    although this means that AREA has not been completely regularized  C    although this means that AREA has not been completely regularized
708  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
709            DIAGarrayA(I,J) = AREA(I,J,bi,bj)              DIAGarrayA(I,J) = tmpscal3itd(i,j)
710  #endif  #endif
711  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
712            SItrAREA(I,J,bi,bj,1)=AREA(I,J,bi,bj)              SItrAREA(I,J,bi,bj,1)=tmpscal3itd(i,j)
713  #endif  #endif
714               ENDIF
715              ENDDO
716           ENDDO           ENDDO
717          ENDDO          ENDDO
718    
# Line 675  C    which includes ridging as in item 2 Line 721  C    which includes ridging as in item 2
721  C    and update AREA, HEFF, and HSNOW  C    and update AREA, HEFF, and HSNOW
722          CALL SEAICE_ITD_REDIST(bi, bj, myTime, myIter, myThid)          CALL SEAICE_ITD_REDIST(bi, bj, myTime, myIter, myThid)
723          CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)          CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)
724    #endif /* SEAICE_ITD */
725    
726    #ifdef SEAICE_DEBUG
727    #ifdef SEAICE_ITD
728            WRITE(msgBufForm,'(A,I2,A)') '(A,',nITD,'F14.10)'
729    #else
730            WRITE(msgBufForm,'(A,A)') '(A,  F14.10)'
731    #endif
732            WRITE(msgBuf,msgBufForm)
733         &    ' SEAICE_GROWTH: Heff increments 0, HEFF = ',
734    #ifdef SEAICE_ITD
735         &     HEFFITD(1,1,:,bi,bj)
736    #else
737         &     HEFF(1,1,bi,bj)
738  #endif  #endif
739            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
740         &    SQUEEZE_RIGHT , myThid)
741            WRITE(msgBuf,msgBufForm)
742         &    ' SEAICE_GROWTH: Area increments 0, AREA = ',
743    #ifdef SEAICE_ITD
744         &     AREAITD(1,1,:,bi,bj)
745    #else
746         &     AREA(1,1,bi,bj)
747    #endif
748            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
749         &    SQUEEZE_RIGHT , myThid)
750    #endif /* SEAICE_DEBUG */
751    
752  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
753  C        ENDIF SEAICEadjMODE.EQ.0  C        end SEAICEadjMODE.EQ.0 statement:
754          ENDIF          ENDIF
755  #endif  #endif
756    
# Line 700  C 3) store regularized values of heff, h Line 772  C 3) store regularized values of heff, h
772           ENDDO           ENDDO
773          ENDDO          ENDDO
774  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
775          DO K=1,nITD          DO IT=1,nITD
776           DO J=1,sNy           DO J=1,sNy
777            DO I=1,sNx            DO I=1,sNx
778             HEFFITDpreTH(I,J,K)=HEFFITD(I,J,K,bi,bj)             HEFFITDpreTH(I,J,IT)=HEFFITD(I,J,IT,bi,bj)
779             HSNWITDpreTH(I,J,K)=HSNOWITD(I,J,K,bi,bj)             HSNWITDpreTH(I,J,IT)=HSNOWITD(I,J,IT,bi,bj)
780             AREAITDpreTH(I,J,K)=AREAITD(I,J,K,bi,bj)             AREAITDpreTH(I,J,IT)=AREAITD(I,J,IT,bi,bj)
781    
782  C memorize areal and volume fraction of each ITD category  C memorize areal and volume fraction of each ITD category
783             IF (AREA(I,J,bi,bj).GT.0) THEN             IF (AREA(I,J,bi,bj) .GT. ZERO) THEN
784              areaFracFactor(I,J,K)=AREAITD(I,J,K,bi,bj)/AREA(I,J,bi,bj)              areaFracFactor(I,J,IT)=AREAITD(I,J,IT,bi,bj)/AREA(I,J,bi,bj)
785             ELSE             ELSE
786              areaFracFactor(I,J,K)=ZERO  C           if there is no ice, potential growth starts in 1st category
787             ENDIF              IF (IT .EQ. 1) THEN
788             IF (HEFF(I,J,bi,bj).GT.0) THEN               areaFracFactor(I,J,IT)=ONE
789              heffFracFactor(I,J,K)=HEFFITD(I,J,K,bi,bj)/HEFF(I,J,bi,bj)              ELSE
790             ELSE               areaFracFactor(I,J,IT)=ZERO
791              heffFracFactor(I,J,K)=ZERO              ENDIF
792             ENDIF             ENDIF
793            ENDDO            ENDDO
794           ENDDO           ENDDO
795          ENDDO          ENDDO
796    #ifdef ALLOW_SITRACER
797  C prepare SItrHEFF to be computed as cumulative sum  C prepare SItrHEFF to be computed as cumulative sum
798          DO K=2,5          DO iTr=2,5
799           DO J=1,sNy           DO J=1,sNy
800            DO I=1,sNx            DO I=1,sNx
801             SItrHEFF(I,J,bi,bj,K)=ZERO             SItrHEFF(I,J,bi,bj,iTr)=ZERO
802            ENDDO            ENDDO
803           ENDDO           ENDDO
804          ENDDO          ENDDO
# Line 736  C prepare SItrAREA to be computed as cum Line 809  C prepare SItrAREA to be computed as cum
809           ENDDO           ENDDO
810          ENDDO          ENDDO
811  #endif  #endif
812    #endif /* SEAICE_ITD */
813    
814  C 4) treat sea ice salinity pathological cases  C 4) treat sea ice salinity pathological cases
815  #ifdef SEAICE_VARIABLE_SALINITY  #ifdef SEAICE_VARIABLE_SALINITY
# Line 791  Cgf no additional dependency of air-sea Line 865  Cgf no additional dependency of air-sea
865            ENDDO            ENDDO
866           ENDDO           ENDDO
867  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
868           DO K=1,nITD           DO IT=1,nITD
869            DO J=1,sNy            DO J=1,sNy
870             DO I=1,sNx             DO I=1,sNx
871              HEFFITDpreTH(I,J,K) = 0. _d 0              HEFFITDpreTH(I,J,IT) = 0. _d 0
872              HSNWITDpreTH(I,J,K) = 0. _d 0              HSNWITDpreTH(I,J,IT) = 0. _d 0
873              AREAITDpreTH(I,J,K) = 0. _d 0              AREAITDpreTH(I,J,IT) = 0. _d 0
874             ENDDO             ENDDO
875            ENDDO            ENDDO
876           ENDDO           ENDDO
# Line 821  CADJ STORE HEFFpreTH = comlev1_bibj, key Line 895  CADJ STORE HEFFpreTH = comlev1_bibj, key
895  CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte  CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte
896  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
897  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
898          DO K=1,nITD          DO IT=1,nITD
899  #endif  #endif
900          DO J=1,sNy          DO J=1,sNy
901           DO I=1,sNx           DO I=1,sNx
902    
903  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
904            IF (HEFFITDpreTH(I,J,K) .GT. ZERO) THEN            IF (HEFFITDpreTH(I,J,IT) .GT. ZERO) THEN
905  #ifdef SEAICE_GROWTH_LEGACY  #ifdef SEAICE_GROWTH_LEGACY
906              tmpscal1          = MAX(SEAICE_area_reg/float(nITD),              tmpscal1          = MAX(SEAICE_area_reg/float(nITD),
907       &                              AREAITDpreTH(I,J,K))       &                              AREAITDpreTH(I,J,IT))
908              hsnowActualMult(I,J,K) = HSNWITDpreTH(I,J,K)/tmpscal1              hsnowActualMult(I,J,IT) = HSNWITDpreTH(I,J,IT)/tmpscal1
909              tmpscal2               = HEFFITDpreTH(I,J,K)/tmpscal1              tmpscal2               = HEFFITDpreTH(I,J,IT)/tmpscal1
910              heffActualMult(I,J,K)  = MAX(tmpscal2,SEAICE_hice_reg)              heffActualMult(I,J,IT)  = MAX(tmpscal2,SEAICE_hice_reg)
911  #else /* SEAICE_GROWTH_LEGACY */  #else /* SEAICE_GROWTH_LEGACY */
912  cif        regularize AREA with SEAICE_area_reg  cif        regularize AREA with SEAICE_area_reg
913             tmpscal1 = SQRT(AREAITDpreTH(I,J,K) * AREAITDpreTH(I,J,K)             tmpscal1 = SQRT(AREAITDpreTH(I,J,IT) * AREAITDpreTH(I,J,IT)
914       &                     + area_reg_sq)       &                     + area_reg_sq)
915  cif        heffActual calculated with the regularized AREA  cif        heffActual calculated with the regularized AREA
916             tmpscal2 = HEFFITDpreTH(I,J,K) / tmpscal1             tmpscal2 = HEFFITDpreTH(I,J,IT) / tmpscal1
917  cif        regularize heffActual with SEAICE_hice_reg (add lower bound)  cif        regularize heffActual with SEAICE_hice_reg (add lower bound)
918             heffActualMult(I,J,K) = SQRT(tmpscal2 * tmpscal2             heffActualMult(I,J,IT) = SQRT(tmpscal2 * tmpscal2
919       &                                  + hice_reg_sq)       &                                  + hice_reg_sq)
920  cif        hsnowActual calculated with the regularized AREA  cif        hsnowActual calculated with the regularized AREA
921             hsnowActualMult(I,J,K) = HSNWITDpreTH(I,J,K) / tmpscal1             hsnowActualMult(I,J,IT) = HSNWITDpreTH(I,J,IT) / tmpscal1
922  #endif /* SEAICE_GROWTH_LEGACY */  #endif /* SEAICE_GROWTH_LEGACY */
923  cif        regularize the inverse of heffActual by hice_reg  cif        regularize the inverse of heffActual by hice_reg
924             recip_heffActualMult(I,J,K)  = AREAITDpreTH(I,J,K) /             recip_heffActualMult(I,J,IT)  = AREAITDpreTH(I,J,IT) /
925       &                 sqrt(HEFFITDpreTH(I,J,K) * HEFFITDpreTH(I,J,K)       &                 sqrt(HEFFITDpreTH(I,J,IT) * HEFFITDpreTH(I,J,IT)
926       &                      + hice_reg_sq)       &                      + hice_reg_sq)
927  cif       Do not regularize when HEFFpreTH = 0  cif       Do not regularize when HEFFpreTH = 0
928            ELSE            ELSE
929              heffActualMult(I,J,K) = ZERO              heffActualMult(I,J,IT) = ZERO
930              hsnowActualMult(I,J,K) = ZERO              hsnowActualMult(I,J,IT) = ZERO
931              recip_heffActualMult(I,J,K)  = ZERO              recip_heffActualMult(I,J,IT)  = ZERO
932            ENDIF            ENDIF
933  #else /* SEAICE_ITD */  #else /* SEAICE_ITD */
934            IF (HEFFpreTH(I,J) .GT. ZERO) THEN            IF (HEFFpreTH(I,J) .GT. ZERO) THEN
# Line 864  cif       Do not regularize when HEFFpre Line 938  cif       Do not regularize when HEFFpre
938              tmpscal2         = HEFFpreTH(I,J)/tmpscal1              tmpscal2         = HEFFpreTH(I,J)/tmpscal1
939              heffActual(I,J)  = MAX(tmpscal2,SEAICE_hice_reg)              heffActual(I,J)  = MAX(tmpscal2,SEAICE_hice_reg)
940  #else /* SEAICE_GROWTH_LEGACY */  #else /* SEAICE_GROWTH_LEGACY */
941  cif        regularize AREA with SEAICE_area_reg  Cif        regularize AREA with SEAICE_area_reg
942             tmpscal1 = SQRT(AREApreTH(I,J)* AREApreTH(I,J) + area_reg_sq)             tmpscal1 = SQRT(AREApreTH(I,J)* AREApreTH(I,J) + area_reg_sq)
943  cif        heffActual calculated with the regularized AREA  Cif        heffActual calculated with the regularized AREA
944             tmpscal2 = HEFFpreTH(I,J) / tmpscal1             tmpscal2 = HEFFpreTH(I,J) / tmpscal1
945  cif        regularize heffActual with SEAICE_hice_reg (add lower bound)  Cif        regularize heffActual with SEAICE_hice_reg (add lower bound)
946             heffActual(I,J) = SQRT(tmpscal2 * tmpscal2 + hice_reg_sq)             heffActual(I,J) = SQRT(tmpscal2 * tmpscal2 + hice_reg_sq)
947  cif        hsnowActual calculated with the regularized AREA  Cif        hsnowActual calculated with the regularized AREA
948             hsnowActual(I,J) = HSNWpreTH(I,J) / tmpscal1             hsnowActual(I,J) = HSNWpreTH(I,J) / tmpscal1
949  #endif /* SEAICE_GROWTH_LEGACY */  #endif /* SEAICE_GROWTH_LEGACY */
950  cif        regularize the inverse of heffActual by hice_reg  Cif        regularize the inverse of heffActual by hice_reg
951             recip_heffActual(I,J)  = AREApreTH(I,J) /             recip_heffActual(I,J)  = AREApreTH(I,J) /
952       &                 sqrt(HEFFpreTH(I,J)*HEFFpreTH(I,J) + hice_reg_sq)       &                 sqrt(HEFFpreTH(I,J)*HEFFpreTH(I,J) + hice_reg_sq)
953  cif       Do not regularize when HEFFpreTH = 0  Cif       Do not regularize when HEFFpreTH = 0
954            ELSE            ELSE
955              heffActual(I,J) = ZERO              heffActual(I,J) = ZERO
956              hsnowActual(I,J) = ZERO              hsnowActual(I,J) = ZERO
# Line 900  cif       Do not regularize when HEFFpre Line 974  cif       Do not regularize when HEFFpre
974  C 5) COMPUTE MAXIMUM LATENT HEAT FLUXES FOR THE CURRENT ICE  C 5) COMPUTE MAXIMUM LATENT HEAT FLUXES FOR THE CURRENT ICE
975  C    AND SNOW THICKNESS  C    AND SNOW THICKNESS
976  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
977          DO K=1,nITD          DO IT=1,nITD
978  #endif  #endif
979          DO J=1,sNy          DO J=1,sNy
980           DO I=1,sNx           DO I=1,sNx
981  c        The latent heat flux over the sea ice which  C        The latent heat flux over the sea ice which
982  c        will sublimate all of the snow and ice over one time  C        will sublimate all of the snow and ice over one time
983  c        step (W/m^2)  C        step (W/m^2)
984  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
985            IF (HEFFITDpreTH(I,J,K) .GT. ZERO) THEN            IF (HEFFITDpreTH(I,J,IT) .GT. ZERO) THEN
986              latentHeatFluxMaxMult(I,J,K) = lhSublim*recip_deltaTtherm *              latentHeatFluxMaxMult(I,J,IT) = lhSublim*recip_deltaTtherm *
987       &         (HEFFITDpreTH(I,J,K)*SEAICE_rhoIce +       &         (HEFFITDpreTH(I,J,IT)*SEAICE_rhoIce +
988       &          HSNWITDpreTH(I,J,K)*SEAICE_rhoSnow)/AREAITDpreTH(I,J,K)       &          HSNWITDpreTH(I,J,IT)*SEAICE_rhoSnow)
989         &         /AREAITDpreTH(I,J,IT)
990            ELSE            ELSE
991              latentHeatFluxMaxMult(I,J,K) = ZERO              latentHeatFluxMaxMult(I,J,IT) = ZERO
992            ENDIF            ENDIF
993  #else  #else
994            IF (HEFFpreTH(I,J) .GT. ZERO) THEN            IF (HEFFpreTH(I,J) .GT. ZERO) THEN
# Line 963  cCADJ STORE TmixLoc = comlev1_bibj, key Line 1038  cCADJ STORE TmixLoc = comlev1_bibj, key
1038  C determine available heat due to the atmosphere -- for ice covered water  C determine available heat due to the atmosphere -- for ice covered water
1039  C =======================================================================  C =======================================================================
1040    
1041  #ifdef ALLOW_ATM_WIND          IF (useRelativeWind.AND.useAtmWind) THEN
         IF (useRelativeWind) THEN  
1042  C     Compute relative wind speed over sea ice.  C     Compute relative wind speed over sea ice.
1043           DO J=1,sNy           DO J=1,sNy
1044            DO I=1,sNx            DO I=1,sNx
# Line 985  C     Compute relative wind speed over s Line 1059  C     Compute relative wind speed over s
1059            ENDDO            ENDDO
1060           ENDDO           ENDDO
1061          ENDIF          ENDIF
 #endif /* ALLOW_ATM_WIND */  
1062    
1063  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
1064  CADJ STORE tice(:,:,bi,bj)  CADJ STORE tice(:,:,bi,bj)
# Line 1001  CADJ &                       key = iicek Line 1074  CADJ &                       key = iicek
1074    
1075  C--   Start loop over multi-categories  C--   Start loop over multi-categories
1076  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1077  CToM  SEAICE_multDim = nITD (see SEAICE_SIZE.h and seaice_readparms.F)          DO IT=1,nITD
 #endif  
         DO IT=1,SEAICE_multDim  
 c        homogeneous distribution between 0 and 2 x heffActual  
 #ifndef SEAICE_ITD  
          pFac = (2.0 _d 0*real(IT)-1.0 _d 0)*recip_multDim  
 #endif  
1078           DO J=1,sNy           DO J=1,sNy
1079            DO I=1,sNx            DO I=1,sNx
1080  #ifndef SEAICE_ITD  CToM for SEAICE_ITD heffActualMult and latentHeatFluxMaxMult are calculated above
 CToM for SEAICE_ITD heffActualMult and latentHeatFluxMaxMult are calculated above  
1081  C    (instead of heffActual and latentHeatFluxMax)  C    (instead of heffActual and latentHeatFluxMax)
1082               ticeInMult(I,J,IT)=TICES(I,J,IT,bi,bj)
1083               ticeOutMult(I,J,IT)=TICES(I,J,IT,bi,bj)
1084               TICE(I,J,bi,bj) = ZERO
1085               TICES(I,J,IT,bi,bj) = ZERO
1086              ENDDO
1087             ENDDO
1088            ENDDO
1089    #else
1090            DO IT=1,SEAICE_multDim
1091    C        homogeneous distribution between 0 and 2 x heffActual
1092             pFac = (2.0 _d 0*IT - 1.0 _d 0)*recip_multDim
1093             pFacSnow = 1. _d 0
1094             IF ( SEAICE_useMultDimSnow ) pFacSnow=pFac
1095             DO J=1,sNy
1096              DO I=1,sNx
1097             heffActualMult(I,J,IT)= heffActual(I,J)*pFac             heffActualMult(I,J,IT)= heffActual(I,J)*pFac
1098               hsnowActualMult(I,J,IT)=hsnowActual(I,J)*pFacSnow
1099  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
1100             latentHeatFluxMaxMult(I,J,IT) = latentHeatFluxMax(I,J)*pFac             latentHeatFluxMaxMult(I,J,IT) = latentHeatFluxMax(I,J)*pFac
1101  #endif  #endif
 #endif  
1102             ticeInMult(I,J,IT)=TICES(I,J,IT,bi,bj)             ticeInMult(I,J,IT)=TICES(I,J,IT,bi,bj)
1103             ticeOutMult(I,J,IT)=TICES(I,J,IT,bi,bj)             ticeOutMult(I,J,IT)=TICES(I,J,IT,bi,bj)
1104             TICE(I,J,bi,bj) = ZERO             TICE(I,J,bi,bj) = ZERO
# Line 1025  C    (instead of heffActual and latentHe Line 1106  C    (instead of heffActual and latentHe
1106            ENDDO            ENDDO
1107           ENDDO           ENDDO
1108          ENDDO          ENDDO
1109    #endif
1110    
1111  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
1112  CADJ STORE heffActualMult = comlev1_bibj, key = iicekey, byte = isbyte  CADJ STORE heffActualMult = comlev1_bibj, key = iicekey, byte = isbyte
1113    CADJ STORE hsnowActualMult= comlev1_bibj, key = iicekey, byte = isbyte
1114  CADJ STORE ticeInMult     = comlev1_bibj, key = iicekey, byte = isbyte  CADJ STORE ticeInMult     = comlev1_bibj, key = iicekey, byte = isbyte
1115  # ifdef SEAICE_CAP_SUBLIM  # ifdef SEAICE_CAP_SUBLIM
1116  CADJ STORE latentHeatFluxMaxMult  CADJ STORE latentHeatFluxMaxMult
# Line 1043  CADJ &     comlev1_bibj, key = iicekey, Line 1126  CADJ &     comlev1_bibj, key = iicekey,
1126    
1127          DO IT=1,SEAICE_multDim          DO IT=1,SEAICE_multDim
1128           CALL SEAICE_SOLVE4TEMP(           CALL SEAICE_SOLVE4TEMP(
 #ifdef SEAICE_ITD  
1129       I        UG, heffActualMult(1,1,IT), hsnowActualMult(1,1,IT),       I        UG, heffActualMult(1,1,IT), hsnowActualMult(1,1,IT),
 #else  
      I        UG, heffActualMult(1,1,IT), hsnowActual,  
 #endif  
1130  #ifdef SEAICE_CAP_SUBLIM  #ifdef SEAICE_CAP_SUBLIM
1131       I        latentHeatFluxMaxMult(1,1,IT),       I        latentHeatFluxMaxMult(1,1,IT),
1132  #endif  #endif
1133       U        ticeInMult(1,1,IT), ticeOutMult(1,1,IT),       U        ticeInMult(1,1,IT), ticeOutMult(1,1,IT),
1134       O        a_QbyATMmult_cover(1,1,IT), a_QSWbyATMmult_cover(1,1,IT),       O        a_QbyATMmult_cover(1,1,IT),
1135         O        a_QSWbyATMmult_cover(1,1,IT),
1136       O        a_FWbySublimMult(1,1,IT),       O        a_FWbySublimMult(1,1,IT),
1137       I        bi, bj, myTime, myIter, myThid )       I        bi, bj, myTime, myIter, myThid )
1138          ENDDO          ENDDO
1139    
1140  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
1141  CADJ STORE heffActualMult = comlev1_bibj, key = iicekey, byte = isbyte  CADJ STORE heffActualMult = comlev1_bibj, key = iicekey, byte = isbyte
1142    CADJ STORE hsnowActualMult= comlev1_bibj, key = iicekey, byte = isbyte
1143  CADJ STORE ticeOutMult    = comlev1_bibj, key = iicekey, byte = isbyte  CADJ STORE ticeOutMult    = comlev1_bibj, key = iicekey, byte = isbyte
1144  # ifdef SEAICE_CAP_SUBLIM  # ifdef SEAICE_CAP_SUBLIM
1145  CADJ STORE latentHeatFluxMaxMult  CADJ STORE latentHeatFluxMaxMult
# Line 1078  CADJ &     comlev1_bibj, key = iicekey, Line 1159  CADJ &     comlev1_bibj, key = iicekey,
1159  C     update TICE & TICES  C     update TICE & TICES
1160  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1161  C     calculate area weighted mean  C     calculate area weighted mean
1162  C     (although the ice's temperature relates to its energy content  C     (although the ice temperature relates to its energy content
1163  C      and hence should be averaged weighted by ice volume [heffFracFactor],  C      and hence should be averaged weighted by ice volume,
1164  C      the temperature here is a result of the fluxes through the ice surface  C      the temperature here is a result of the fluxes through the ice surface
1165  C      computed individually for each single category in SEAICE_SOLVE4TEMP  C      computed individually for each single category in SEAICE_SOLVE4TEMP
1166  C      and hence is averaged area weighted [areaFracFactor])  C      and hence is averaged area weighted [areaFracFactor])
1167             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
1168       &          +  ticeOutMult(I,J,IT)*areaFracFactor(I,J,K)       &          +  ticeOutMult(I,J,IT)*areaFracFactor(I,J,IT)
1169  #else  #else
1170             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
1171       &          +  ticeOutMult(I,J,IT)*recip_multDim       &          +  ticeOutMult(I,J,IT)*recip_multDim
# Line 1095  C     average over categories Line 1176  C     average over categories
1176  C     calculate area weighted mean  C     calculate area weighted mean
1177  C     (fluxes are per unit (ice surface) area and are thus area weighted)  C     (fluxes are per unit (ice surface) area and are thus area weighted)
1178             a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)             a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)
1179       &          + a_QbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,K)       &          + a_QbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,IT)
1180             a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)             a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
1181       &          + a_QSWbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,K)       &          + a_QSWbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,IT)
1182             a_FWbySublim     (I,J) = a_FWbySublim(I,J)             a_FWbySublim     (I,J) = a_FWbySublim(I,J)
1183       &          + a_FWbySublimMult(I,J,IT)*areaFracFactor(I,J,K)       &          + a_FWbySublimMult(I,J,IT)*areaFracFactor(I,J,IT)
1184  #else  #else
1185             a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)             a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)
1186       &          + a_QbyATMmult_cover(I,J,IT)*recip_multDim       &          + a_QbyATMmult_cover(I,J,IT)*recip_multDim
# Line 1116  C     (fluxes are per unit (ice surface) Line 1197  C     (fluxes are per unit (ice surface)
1197  # ifdef ALLOW_DIAGNOSTICS  # ifdef ALLOW_DIAGNOSTICS
1198          DO J=1,sNy          DO J=1,sNy
1199           DO I=1,sNx           DO I=1,sNx
1200  c          The actual latent heat flux realized by SOLVE4TEMP  C          The actual latent heat flux realized by SOLVE4TEMP
1201             DIAGarrayA(I,J) = a_FWbySublim(I,J) * lhSublim             DIAGarrayA(I,J) = a_FWbySublim(I,J) * lhSublim
1202           ENDDO           ENDDO
1203          ENDDO          ENDDO
1204  cif     The actual vs. maximum latent heat flux  Cif     The actual vs. maximum latent heat flux
1205          IF ( useDiagnostics ) THEN          IF ( useDiagnostics ) THEN
1206            CALL DIAGNOSTICS_FILL(DIAGarrayA,            CALL DIAGNOSTICS_FILL(DIAGarrayA,
1207       &     'SIactLHF',0,1,3,bi,bj,myThid)       &     'SIactLHF',0,1,3,bi,bj,myThid)
# Line 1141  CADJ STORE a_FWbySublim    = comlev1_bib Line 1222  CADJ STORE a_FWbySublim    = comlev1_bib
1222    
1223  C switch heat fluxes from W/m2 to 'effective' ice meters  C switch heat fluxes from W/m2 to 'effective' ice meters
1224  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1225          DO K=1,nITD          DO IT=1,nITD
1226           DO J=1,sNy           DO J=1,sNy
1227            DO I=1,sNx            DO I=1,sNx
1228             a_QbyATMmult_cover(I,J,K)   = a_QbyATMmult_cover(I,J,K)             a_QbyATMmult_cover(I,J,IT)   = a_QbyATMmult_cover(I,J,IT)
1229       &          * convertQ2HI * AREAITDpreTH(I,J,K)       &          * convertQ2HI * AREAITDpreTH(I,J,IT)
1230             a_QSWbyATMmult_cover(I,J,K) = a_QSWbyATMmult_cover(I,J,K)             a_QSWbyATMmult_cover(I,J,IT) = a_QSWbyATMmult_cover(I,J,IT)
1231       &          * convertQ2HI * AREAITDpreTH(I,J,K)       &          * convertQ2HI * AREAITDpreTH(I,J,IT)
1232  C and initialize r_QbyATM_cover  C and initialize r_QbyATMmult_cover
1233             r_QbyATMmult_cover(I,J,K)=a_QbyATMmult_cover(I,J,K)             r_QbyATMmult_cover(I,J,IT)=a_QbyATMmult_cover(I,J,IT)
1234  C     Convert fresh water flux by sublimation to 'effective' ice meters.  C     Convert fresh water flux by sublimation to 'effective' ice meters.
1235  C     Negative sublimation is resublimation and will be added as snow.  C     Negative sublimation is resublimation and will be added as snow.
1236  #ifdef SEAICE_DISABLE_SUBLIM  #ifdef SEAICE_DISABLE_SUBLIM
1237             a_FWbySublimMult(I,J,K) = ZERO             a_FWbySublimMult(I,J,IT) = ZERO
1238  #endif  #endif
1239             a_FWbySublimMult(I,J,K) = SEAICE_deltaTtherm*recip_rhoIce             a_FWbySublimMult(I,J,IT) = SEAICE_deltaTtherm*recip_rhoIce
1240       &            * a_FWbySublimMult(I,J,K)*AREAITDpreTH(I,J,K)       &            * a_FWbySublimMult(I,J,IT)*AREAITDpreTH(I,J,IT)
1241             r_FWbySublimMult(I,J,K)=a_FWbySublimMult(I,J,K)             r_FWbySublimMult(I,J,IT)=a_FWbySublimMult(I,J,IT)
1242            ENDDO            ENDDO
1243           ENDDO           ENDDO
1244          ENDDO          ENDDO
1245          DO J=1,sNy          DO J=1,sNy
1246           DO I=1,sNx           DO I=1,sNx
1247              a_QbyATM_open(I,J)    = a_QbyATM_open(I,J)
1248         &         * convertQ2HI * ( ONE - AREApreTH(I,J) )
1249              a_QSWbyATM_open(I,J)  = a_QSWbyATM_open(I,J)
1250         &         * convertQ2HI * ( ONE - AREApreTH(I,J) )
1251  C and initialize r_QbyATM_open  C and initialize r_QbyATM_open
1252            r_QbyATM_open(I,J)=a_QbyATM_open(I,J)            r_QbyATM_open(I,J)=a_QbyATM_open(I,J)
1253           ENDDO           ENDDO
# Line 1184  C and initialize r_QbyATM_cover/r_QbyATM Line 1269  C and initialize r_QbyATM_cover/r_QbyATM
1269  C     Convert fresh water flux by sublimation to 'effective' ice meters.  C     Convert fresh water flux by sublimation to 'effective' ice meters.
1270  C     Negative sublimation is resublimation and will be added as snow.  C     Negative sublimation is resublimation and will be added as snow.
1271  #ifdef SEAICE_DISABLE_SUBLIM  #ifdef SEAICE_DISABLE_SUBLIM
1272  cgf just for those who may need to omit this term to reproduce old results  Cgf just for those who may need to omit this term to reproduce old results
1273            a_FWbySublim(I,J) = ZERO            a_FWbySublim(I,J) = ZERO
1274  #endif  #endif /* SEAICE_DISABLE_SUBLIM */
1275            a_FWbySublim(I,J) = SEAICE_deltaTtherm*recip_rhoIce            a_FWbySublim(I,J) = SEAICE_deltaTtherm*recip_rhoIce
1276       &           * a_FWbySublim(I,J)*AREApreTH(I,J)       &           * a_FWbySublim(I,J)*AREApreTH(I,J)
1277            r_FWbySublim(I,J)=a_FWbySublim(I,J)            r_FWbySublim(I,J)=a_FWbySublim(I,J)
# Line 1210  CADJ STORE r_FWbySublim    = comlev1_bib Line 1295  CADJ STORE r_FWbySublim    = comlev1_bib
1295  Cgf no additional dependency through ice cover  Cgf no additional dependency through ice cover
1296          IF ( SEAICEadjMODE.GE.3 ) THEN          IF ( SEAICEadjMODE.GE.3 ) THEN
1297  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1298           DO K=1,nITD           DO IT=1,nITD
1299            DO J=1,sNy            DO J=1,sNy
1300             DO I=1,sNx             DO I=1,sNx
1301              a_QbyATMmult_cover(I,J,K)   = 0. _d 0              a_QbyATMmult_cover(I,J,IT)   = 0. _d 0
1302              r_QbyATMmult_cover(I,J,K)   = 0. _d 0              r_QbyATMmult_cover(I,J,IT)   = 0. _d 0
1303              a_QSWbyATMmult_cover(I,J,K) = 0. _d 0              a_QSWbyATMmult_cover(I,J,IT) = 0. _d 0
1304             ENDDO             ENDDO
1305            ENDDO            ENDDO
1306           ENDDO           ENDDO
1307  #else  #else
1308           DO J=1,sNy           DO J=1,sNy
1309            DO I=1,sNx            DO I=1,sNx
# Line 1244  CADJ &                       key = iicek Line 1329  CADJ &                       key = iicek
1329    
1330          DO J=1,sNy          DO J=1,sNy
1331           DO I=1,sNx           DO I=1,sNx
1332  c         FREEZING TEMP. OF SEA WATER (deg C)  C         FREEZING TEMP. OF SEA WATER (deg C)
1333            tempFrz = SEAICE_tempFrz0 +            tempFrz = SEAICE_tempFrz0 +
1334       &              SEAICE_dTempFrz_dS *salt(I,J,kSurface,bi,bj)       &              SEAICE_dTempFrz_dS *salt(I,J,kSurface,bi,bj)
1335  c efficiency of turbulent fluxes : dependency to sign of THETA-TBC  C efficiency of turbulent fluxes : dependency to sign of THETA-TBC
1336            IF ( theta(I,J,kSurface,bi,bj) .GE. tempFrz ) THEN            IF ( theta(I,J,kSurface,bi,bj) .GE. tempFrz ) THEN
1337             tmpscal1 = SEAICE_mcPheePiston             tmpscal1 = SEAICE_mcPheePiston
1338            ELSE            ELSE
1339             tmpscal1 =SEAICE_frazilFrac*drF(kSurface)/SEAICE_deltaTtherm             tmpscal1 =SEAICE_frazilFrac*drF(kSurface)/SEAICE_deltaTtherm
1340            ENDIF            ENDIF
1341  c efficiency of turbulent fluxes : dependency to AREA (McPhee cases)  C efficiency of turbulent fluxes : dependency to AREA (McPhee cases)
1342            IF ( (AREApreTH(I,J) .GT. 0. _d 0).AND.            IF ( (AREApreTH(I,J) .GT. 0. _d 0).AND.
1343       &         (.NOT.SEAICE_mcPheeStepFunc) ) THEN       &         (.NOT.SEAICE_mcPheeStepFunc) ) THEN
1344             MixedLayerTurbulenceFactor = ONE -             MixedLayerTurbulenceFactor = ONE -
# Line 1264  c efficiency of turbulent fluxes : depen Line 1349  c efficiency of turbulent fluxes : depen
1349            ELSE            ELSE
1350             MixedLayerTurbulenceFactor = ONE             MixedLayerTurbulenceFactor = ONE
1351            ENDIF            ENDIF
1352  c maximum turbulent flux, in ice meters  C maximum turbulent flux, in ice meters
1353            tmpscal2= - (HeatCapacity_Cp*rhoConst * recip_QI)            tmpscal2= - (HeatCapacity_Cp*rhoConst * recip_QI)
1354       &         * (theta(I,J,kSurface,bi,bj)-tempFrz)       &         * (theta(I,J,kSurface,bi,bj)-tempFrz)
1355       &         * SEAICE_deltaTtherm * maskC(i,j,kSurface,bi,bj)       &         * SEAICE_deltaTtherm * maskC(i,j,kSurface,bi,bj)
1356  c available turbulent flux  C available turbulent flux
1357            a_QbyOCN(i,j) =            a_QbyOCN(i,j) =
1358       &         tmpscal1 * tmpscal2 * MixedLayerTurbulenceFactor       &         tmpscal1 * tmpscal2 * MixedLayerTurbulenceFactor
1359            r_QbyOCN(i,j) = a_QbyOCN(i,j)            r_QbyOCN(i,j) = a_QbyOCN(i,j)
 ctm  
           if (i.eq.20 .and. j.eq.20) then  
             print *, HeatCapacity_Cp  
             print *, rhoConst  
             print *, recip_QI  
             print *, theta(20,20,kSurface,bi,bj)  
             print *, tempFrz  
             print *, SEAICE_deltaTtherm  
             print *, maskC(20,20,kSurface,bi,bj)  
             print *, tmpscal2  
             print *, a_QbyOCN(20,20)  
           endif  
 ctm  
1360           ENDDO           ENDDO
1361          ENDDO          ENDDO
1362    
# Line 1292  ctm Line 1364  ctm
1364          CALL ZERO_ADJ_1D( sNx*sNy, r_QbyOCN, myThid)          CALL ZERO_ADJ_1D( sNx*sNy, r_QbyOCN, myThid)
1365  #endif  #endif
1366    
   
1367  C ===================================================================  C ===================================================================
1368  C =========PART 3: determine effective thicknesses increments========  C =========PART 3: determine effective thicknesses increments========
1369  C ===================================================================  C ===================================================================
# Line 1305  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 1376  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
1376  CADJ STORE r_FWbySublim     = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE r_FWbySublim     = comlev1_bibj,key=iicekey,byte=isbyte
1377  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1378  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1379          DO K=1,nITD          DO IT=1,nITD
1380  #endif  #endif
1381          DO J=1,sNy          DO J=1,sNy
1382           DO I=1,sNx           DO I=1,sNx
1383  C     First sublimate/deposite snow  C     First sublimate/deposite snow
1384            tmpscal2 =            tmpscal2 =
1385  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1386       &     MAX(MIN(r_FWbySublimMult(I,J,K),HSNOWITD(I,J,K,bi,bj)       &     MAX(MIN(r_FWbySublimMult(I,J,IT),HSNOWITD(I,J,IT,bi,bj)
1387       &             *SNOW2ICE),ZERO)       &             *SNOW2ICE),ZERO)
1388              d_HSNWbySublim_ITD(I,J,IT) = - tmpscal2 * ICE2SNOW
1389  C         accumulate change over ITD categories  C         accumulate change over ITD categories
1390            d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)     - tmpscal2            d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)      - tmpscal2
      &                                                       *ICE2SNOW  
           HSNOWITD(I,J,K,bi,bj)   = HSNOWITD(I,J,K,bi,bj)   - tmpscal2  
1391       &                                                       *ICE2SNOW       &                                                       *ICE2SNOW
1392            r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - tmpscal2            r_FWbySublimMult(I,J,IT)= r_FWbySublimMult(I,J,IT) - tmpscal2
 C         keep total up to date, too  
           r_FWbySublim(I,J)       = r_FWbySublim(I,J)       - tmpscal2  
1393  #else  #else
1394       &     MAX(MIN(r_FWbySublim(I,J),HSNOW(I,J,bi,bj)*SNOW2ICE),ZERO)       &     MAX(MIN(r_FWbySublim(I,J),HSNOW(I,J,bi,bj)*SNOW2ICE),ZERO)
1395            d_HSNWbySublim(I,J) = - tmpscal2 * ICE2SNOW            d_HSNWbySublim(I,J) = - tmpscal2 * ICE2SNOW
# Line 1339  CADJ STORE r_FWbySublim    = comlev1_bib Line 1407  CADJ STORE r_FWbySublim    = comlev1_bib
1407  C     If anything is left, sublimate ice  C     If anything is left, sublimate ice
1408            tmpscal2 =            tmpscal2 =
1409  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1410       &     MAX(MIN(r_FWbySublimMult(I,J,K),HEFFITD(I,J,K,bi,bj)),ZERO)       &     MAX(MIN(r_FWbySublimMult(I,J,IT),HEFFITD(I,J,IT,bi,bj)),ZERO)
1411              d_HEFFbySublim_ITD(I,J,IT) = - tmpscal2
1412  C         accumulate change over ITD categories  C         accumulate change over ITD categories
1413            d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)     - tmpscal2            d_HEFFbySublim(I,J)      = d_HEFFbySublim(I,J)      - tmpscal2
1414            HEFFITD(I,J,K,bi,bj)    = HEFFITD(I,J,K,bi,bj)    - tmpscal2            r_FWbySublimMult(I,J,IT) = r_FWbySublimMult(I,J,IT) - tmpscal2
           r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - tmpscal2  
 C         keep total up to date, too  
           r_FWbySublim(I,J)       = r_FWbySublim(I,J)       - tmpscal2  
1415  #else  #else
1416       &     MAX(MIN(r_FWbySublim(I,J),HEFF(I,J,bi,bj)),ZERO)       &     MAX(MIN(r_FWbySublim(I,J),HEFF(I,J,bi,bj)),ZERO)
1417            d_HEFFbySublim(I,J) = - tmpscal2            d_HEFFbySublim(I,J) = - tmpscal2
# Line 1360  C     If anything is left, it will be ev Line 1426  C     If anything is left, it will be ev
1426  C     Since a_QbyATM_cover was computed for sublimation, not simple evaporation, we need to  C     Since a_QbyATM_cover was computed for sublimation, not simple evaporation, we need to
1427  C     remove the fusion part for the residual (that happens to be precisely r_FWbySublim).  C     remove the fusion part for the residual (that happens to be precisely r_FWbySublim).
1428  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1429            a_QbyATMmult_cover(I,J,K) = a_QbyATMmult_cover(I,J,K)            a_QbyATMmult_cover(I,J,IT) = a_QbyATMmult_cover(I,J,IT)
1430       &                              - r_FWbySublimMult(I,J,K)       &                               - r_FWbySublimMult(I,J,IT)
1431            r_QbyATMmult_cover(I,J,K) = r_QbyATMmult_cover(I,J,K)            r_QbyATMmult_cover(I,J,IT) = r_QbyATMmult_cover(I,J,IT)
1432       &                              - r_FWbySublimMult(I,J,K)       &                               - r_FWbySublimMult(I,J,IT)
1433           ENDDO  #else
         ENDDO  
 C       end K loop  
         ENDDO  
 C       then update totals        
         DO J=1,sNy  
          DO I=1,sNx  
 #endif  
1434            a_QbyATM_cover(I,J) = a_QbyATM_cover(I,J)-r_FWbySublim(I,J)            a_QbyATM_cover(I,J) = a_QbyATM_cover(I,J)-r_FWbySublim(I,J)
1435            r_QbyATM_cover(I,J) = r_QbyATM_cover(I,J)-r_FWbySublim(I,J)            r_QbyATM_cover(I,J) = r_QbyATM_cover(I,J)-r_FWbySublim(I,J)
1436    #endif
1437           ENDDO           ENDDO
1438          ENDDO          ENDDO
1439    #ifdef SEAICE_ITD
1440    C       end IT loop
1441            ENDDO
1442    #endif
1443    #ifdef SEAICE_DEBUG
1444  c ToM<<< debug seaice_growth  c ToM<<< debug seaice_growth
1445          WRITE(msgBuf,'(A,7F6.2)')          WRITE(msgBuf,msgBufForm)
1446       &    ' SEAICE_GROWTH: Heff increments 1, HEFFITD = ',       &    ' SEAICE_GROWTH: Hsnow increments 1, d_HSNWySublim = ',
1447       &     HEFFITD(20,20,:,bi,bj)  #ifdef SEAICE_ITD
1448         &     d_HSNWbySublim_ITD(1,1,:)
1449    #else
1450         &     d_HSNWbySublim(1,1)
1451    #endif
1452            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1453         &    SQUEEZE_RIGHT , myThid)
1454            WRITE(msgBuf,msgBufForm)
1455         &    ' SEAICE_GROWTH: Heff increments 1, d_HEFFbySublim = ',
1456    #ifdef SEAICE_ITD
1457         &     d_HEFFbySublim_ITD(1,1,:)
1458    #else
1459         &     d_HEFFbySublim(1,1)
1460    #endif
1461          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1462       &    SQUEEZE_RIGHT , myThid)       &    SQUEEZE_RIGHT , myThid)
1463  c ToM>>>  c ToM>>>
1464    #endif /* SEAICE_DEBUG */
1465    
1466  C compute ice thickness tendency due to ice-ocean interaction  C compute ice thickness tendency due to ice-ocean interaction
1467  C ===========================================================  C ===========================================================
# Line 1393  CADJ STORE r_QbyOCN = comlev1_bibj,key=i Line 1472  CADJ STORE r_QbyOCN = comlev1_bibj,key=i
1472  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1473    
1474  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1475          DO K=1,nITD          DO IT=1,nITD
1476           DO J=1,sNy           DO J=1,sNy
1477            DO I=1,sNx            DO I=1,sNx
1478  C          ice growth/melt due to ocean heat is equally distributed under the ice  C          ice growth/melt due to ocean heat r_QbyOCN (W/m^2) is
1479  C          and hence weighted by fractional area of each thickness category  C          equally distributed under the ice and hence weighted by
1480             tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,K),  C          fractional area of each thickness category
1481       &                               -HEFFITD(I,J,K,bi,bj))             tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,IT),
1482             d_HEFFbyOCNonICE(I,J)= d_HEFFbyOCNonICE(I,J) + tmpscal1       &                               -HEFFITD(I,J,IT,bi,bj))
1483             r_QbyOCN(I,J)        = r_QbyOCN(I,J)         - tmpscal1             d_HEFFbyOCNonICE_ITD(I,J,IT)=tmpscal1
1484             HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj)  + tmpscal1             d_HEFFbyOCNonICE(I,J) = d_HEFFbyOCNonICE(I,J) + tmpscal1
1485              ENDDO
1486             ENDDO
1487            ENDDO
1488  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
1489             SItrHEFF(I,J,bi,bj,2) = SItrHEFF(I,J,bi,bj,2)          DO J=1,sNy
1490       &                           + HEFFITD(I,J,K,bi,bj)           DO I=1,sNx
1491              SItrHEFF(I,J,bi,bj,2) = HEFFpreTH(I,J)
1492         &                          + d_HEFFbySublim(I,J)
1493         &                          + d_HEFFbyOCNonICE(I,J)
1494             ENDDO
1495            ENDDO
1496  #endif  #endif
1497            ENDDO          DO J=1,sNy
1498             DO I=1,sNx
1499              r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)
1500           ENDDO           ENDDO
1501          ENDDO          ENDDO
 c ToM<<< debug seaice_growth  
         WRITE(msgBuf,'(A,7F9.6)')  
      &    ' SEAICE_GROWTH: d_HEFFbyOCNonICE w/ITD: ',  
      &     d_HEFFbyOCNonICE(20,20)  
         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,  
      &    SQUEEZE_RIGHT , myThid)  
 c ToM>>>  
1502  #else /* SEAICE_ITD */  #else /* SEAICE_ITD */
1503          DO J=1,sNy          DO J=1,sNy
1504           DO I=1,sNx           DO I=1,sNx
# Line 1428  c ToM>>> Line 1510  c ToM>>>
1510  #endif  #endif
1511           ENDDO           ENDDO
1512          ENDDO          ENDDO
 c ToM<<< debug seaice_growth  
         WRITE(msgBuf,'(A,7F9.6)')  
      &    ' SEAICE_GROWTH: d_HEFFbyOCNonICE w/o ITD: ',  
      &     d_HEFFbyOCNonICE(20,20)  
         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,  
      &    SQUEEZE_RIGHT , myThid)  
 c ToM>>>  
1513  #endif /* SEAICE_ITD */  #endif /* SEAICE_ITD */
1514    #ifdef SEAICE_DEBUG
1515  c ToM<<< debug seaice_growth  c ToM<<< debug seaice_growth
1516          WRITE(msgBuf,'(A,7F6.2)')          WRITE(msgBuf,msgBufForm)
1517       &    ' SEAICE_GROWTH: Heff increments 2, HEFFITD = ',       &    ' SEAICE_GROWTH: Heff increments 2, d_HEFFbyOCNonICE = ',
1518       &     HEFFITD(20,20,:,bi,bj)  #ifdef SEAICE_ITD
1519         &     d_HEFFbyOCNonICE_ITD(1,1,:)
1520    #else
1521         &     d_HEFFbyOCNonICE(1,1)
1522    #endif
1523          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1524       &    SQUEEZE_RIGHT , myThid)       &    SQUEEZE_RIGHT , myThid)
1525  c ToM>>>  c ToM>>>
1526    #endif /* SEAICE_DEBUG */
1527    
1528  C compute snow melt tendency due to snow-atmosphere interaction  C compute snow melt tendency due to snow-atmosphere interaction
1529  C ==================================================================  C ==================================================================
# Line 1453  CADJ STORE r_QbyATM_cover = comlev1_bibj Line 1534  CADJ STORE r_QbyATM_cover = comlev1_bibj
1534  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1535    
1536  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1537          DO K=1,nITD          DO IT=1,nITD
1538           DO J=1,sNy           DO J=1,sNy
1539            DO I=1,sNx            DO I=1,sNx
1540  C     Convert to standard units (meters of ice) rather than to meters  C     Convert to standard units (meters of ice) rather than to meters
1541  C     of snow. This appears to be more robust.  C     of snow. This appears to be more robust.
1542            tmpscal1=MAX(r_QbyATMmult_cover(I,J,K),-HSNOWITD(I,J,K,bi,bj)             tmpscal1=MAX(r_QbyATMmult_cover(I,J,IT),
1543       &                                           *SNOW2ICE)       &                  -HSNOWITD(I,J,IT,bi,bj)*SNOW2ICE)
1544            tmpscal2=MIN(tmpscal1,0. _d 0)             tmpscal2=MIN(tmpscal1,0. _d 0)
1545  #ifdef SEAICE_MODIFY_GROWTH_ADJ  #ifdef SEAICE_MODIFY_GROWTH_ADJ
1546  Cgf no additional dependency through snow  Cgf no additional dependency through snow
1547            IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0             IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1548  #endif  #endif
1549            d_HSNWbyATMonSNW(I,J)=d_HSNWbyATMonSNW(I,J)+tmpscal2*ICE2SNOW             d_HSNWbyATMonSNW_ITD(I,J,IT) = tmpscal2*ICE2SNOW
1550            HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal2*ICE2SNOW             d_HSNWbyATMonSNW(I,J) = d_HSNWbyATMonSNW(I,J)
1551            r_QbyATMmult_cover(I,J,K)=r_QbyATMmult_cover(I,J,K) - tmpscal2       &                           + tmpscal2*ICE2SNOW
1552  C         keep the total up to date, too             r_QbyATMmult_cover(I,J,IT)=r_QbyATMmult_cover(I,J,IT)
1553            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2       &                           - tmpscal2
1554            ENDDO            ENDDO
1555           ENDDO           ENDDO
1556          ENDDO          ENDDO
1557  #else /* SEAICE_ITD */  #else /* SEAICE_ITD */
1558          DO J=1,sNy          DO J=1,sNy
1559           DO I=1,sNx           DO I=1,sNx
# Line 1490  Cgf no additional dependency through sno Line 1571  Cgf no additional dependency through sno
1571           ENDDO           ENDDO
1572          ENDDO          ENDDO
1573  #endif /* SEAICE_ITD */  #endif /* SEAICE_ITD */
1574    #ifdef SEAICE_DEBUG
1575  c ToM<<< debug seaice_growth  c ToM<<< debug seaice_growth
1576          WRITE(msgBuf,'(A,7F6.2)')          WRITE(msgBuf,msgBufForm)
1577       &    ' SEAICE_GROWTH: Heff increments 3, HEFFITD = ',       &    ' SEAICE_GROWTH: Hsnow increments 3, d_HSNWbyATMonSNW = ',
1578       &     HEFFITD(20,20,:,bi,bj)  #ifdef SEAICE_ITD
1579         &    d_HSNWbyATMonSNW_ITD(1,1,:)
1580    #else
1581         &    d_HSNWbyATMonSNW(1,1)
1582    #endif
1583          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1584       &    SQUEEZE_RIGHT , myThid)       &    SQUEEZE_RIGHT , myThid)
1585  c ToM>>>  c ToM>>>
1586    #endif /* SEAICE_DEBUG */
1587    
1588  C compute ice thickness tendency due to the atmosphere  C compute ice thickness tendency due to the atmosphere
1589  C ====================================================  C ====================================================
# Line 1512  Cgf the v1.81=>v1.82 revision would chan Line 1599  Cgf the v1.81=>v1.82 revision would chan
1599  Cgf warming conditions, the lab_sea results were not changed.  Cgf warming conditions, the lab_sea results were not changed.
1600    
1601  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1602          DO K=1,nITD          DO IT=1,nITD
1603           DO J=1,sNy           DO J=1,sNy
1604            DO I=1,sNx            DO I=1,sNx
1605               tmpscal1 = HEFFITDpreTH(I,J,IT)
1606         &              + d_HEFFbySublim_ITD(I,J,IT)
1607         &              + d_HEFFbyOCNonICE_ITD(I,J,IT)
1608  #ifdef SEAICE_GROWTH_LEGACY  #ifdef SEAICE_GROWTH_LEGACY
1609             tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K))             tmpscal2 = MAX(-tmpscal1,
1610         &                     r_QbyATMmult_cover(I,J,IT))
1611  #else  #else
1612             tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K)             tmpscal2 = MAX(-tmpscal1,
1613         &                     r_QbyATMmult_cover(I,J,IT)
1614  c         Limit ice growth by potential melt by ocean  c         Limit ice growth by potential melt by ocean
1615       &     + AREAITDpreTH(I,J,K) * r_QbyOCN(I,J))       &              + AREAITDpreTH(I,J,IT) * r_QbyOCN(I,J))
1616  #endif /* SEAICE_GROWTH_LEGACY */  #endif /* SEAICE_GROWTH_LEGACY */
1617               d_HEFFbyATMonOCN_cover_ITD(I,J,IT) = tmpscal2
1618             d_HEFFbyATMonOCN_cover(I,J) = d_HEFFbyATMonOCN_cover(I,J)             d_HEFFbyATMonOCN_cover(I,J) = d_HEFFbyATMonOCN_cover(I,J)
1619       &                                 + tmpscal2       &                                 + tmpscal2
1620               d_HEFFbyATMonOCN_ITD(I,J,IT) = d_HEFFbyATMonOCN_ITD(I,J,IT)
1621         &                                 + tmpscal2
1622             d_HEFFbyATMonOCN(I,J)       = d_HEFFbyATMonOCN(I,J)             d_HEFFbyATMonOCN(I,J)       = d_HEFFbyATMonOCN(I,J)
1623       &                                 + tmpscal2       &                                 + tmpscal2
1624             r_QbyATM_cover(I,J)         = r_QbyATM_cover(I,J)             r_QbyATMmult_cover(I,J,IT)  = r_QbyATMmult_cover(I,J,IT)
1625       &                                 - tmpscal2       &                                 - tmpscal2
1626             HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal2            ENDDO
1627             ENDDO
1628            ENDDO
1629  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
1630             SItrHEFF(I,J,bi,bj,3) = SItrHEFF(I,J,bi,bj,3)          DO J=1,sNy
1631       &                           + HEFFITD(I,J,K,bi,bj)           DO I=1,sNx
1632              SItrHEFF(I,J,bi,bj,3) = SItrHEFF(I,J,bi,bj,2)
1633         &                          + d_HEFFbyATMonOCN_cover(I,J)
1634             ENDDO
1635            ENDDO
1636  #endif  #endif
           ENDDO  
          ENDDO  
         ENDDO  
1637  #else /* SEAICE_ITD */  #else /* SEAICE_ITD */
1638          DO J=1,sNy          DO J=1,sNy
1639           DO I=1,sNx           DO I=1,sNx
# Line 1545  c         Limit ice growth by potential Line 1642  c         Limit ice growth by potential
1642            tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J))            tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J))
1643  #else  #else
1644            tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J)+            tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J)+
1645  c         Limit ice growth by potential melt by ocean  C         Limit ice growth by potential melt by ocean
1646       &        AREApreTH(I,J) * r_QbyOCN(I,J))       &        AREApreTH(I,J) * r_QbyOCN(I,J))
1647  #endif /* SEAICE_GROWTH_LEGACY */  #endif /* SEAICE_GROWTH_LEGACY */
1648    
# Line 1557  c         Limit ice growth by potential Line 1654  c         Limit ice growth by potential
1654  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
1655            SItrHEFF(I,J,bi,bj,3)=HEFF(I,J,bi,bj)            SItrHEFF(I,J,bi,bj,3)=HEFF(I,J,bi,bj)
1656  #endif  #endif
1657           ENDDO           ENDDO
1658          ENDDO          ENDDO
1659  #endif /* SEAICE_ITD */  #endif /* SEAICE_ITD */
1660    #ifdef SEAICE_DEBUG
1661  c ToM<<< debug seaice_growth  c ToM<<< debug seaice_growth
1662          WRITE(msgBuf,'(A,7F6.2)')          WRITE(msgBuf,msgBufForm)
1663       &    ' SEAICE_GROWTH: Heff increments 4, HEFFITD = ',       &   ' SEAICE_GROWTH: Heff increments 4, d_HEFFbyATMonOCN_cover = ',
1664       &     HEFFITD(20,20,:,bi,bj)  #ifdef SEAICE_ITD
1665         &     d_HEFFbyATMonOCN_cover_ITD(1,1,:)
1666    #else
1667         &     d_HEFFbyATMonOCN_cover(1,1)
1668    #endif
1669            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1670         &    SQUEEZE_RIGHT , myThid)
1671            WRITE(msgBuf,msgBufForm)
1672         &    ' SEAICE_GROWTH: Heff increments 4, d_HEFFbyATMonOCN = ',
1673    #ifdef SEAICE_ITD
1674         &     d_HEFFbyATMonOCN_ITD(1,1,:)
1675    #else
1676         &     d_HEFFbyATMonOCN(1,1)
1677    #endif
1678          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1679       &    SQUEEZE_RIGHT , myThid)       &    SQUEEZE_RIGHT , myThid)
1680  c ToM>>>  c ToM>>>
1681    #endif /* SEAICE_DEBUG */
1682    
1683  C attribute precip to fresh water or snow stock,  C add snow precipitation to HSNOW.
 C depending on atmospheric conditions.  
1684  C =================================================  C =================================================
1685  #ifdef ALLOW_ATM_TEMP  #ifdef ALLOW_ATM_TEMP
1686  #ifdef ALLOW_AUTODIFF_TAMC  # ifdef ALLOW_AUTODIFF_TAMC
1687  CADJ STORE a_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE a_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1688  CADJ STORE PRECIP(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE PRECIP(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1689  CADJ STORE AREApreTH = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE AREApreTH = comlev1_bibj,key=iicekey,byte=isbyte
1690  #endif /* ALLOW_AUTODIFF_TAMC */  # endif /* ALLOW_AUTODIFF_TAMC */
1691          DO J=1,sNy          IF ( snowPrecipFile .NE. ' ' ) THEN
1692           DO I=1,sNx  C add snowPrecip to HSNOW
1693             DO J=1,sNy
1694              DO I=1,sNx
1695               d_HSNWbyRAIN(I,J) = convertPRECIP2HI * ICE2SNOW *
1696         &          snowPrecip(i,j,bi,bj) * AREApreTH(I,J)
1697               d_HFRWbyRAIN(I,J) = -convertPRECIP2HI *
1698         &          ( PRECIP(I,J,bi,bj) - snowPrecip(I,J,bi,bj) ) *
1699         &          AREApreTH(I,J)
1700               HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyRAIN(I,J)
1701              ENDDO
1702             ENDDO
1703            ELSE
1704    C attribute precip to fresh water or snow stock,
1705    C depending on atmospheric conditions.
1706             DO J=1,sNy
1707              DO I=1,sNx
1708  C possible alternatives to the a_QbyATM_cover criterium  C possible alternatives to the a_QbyATM_cover criterium
1709  c          IF (TICE(I,J,bi,bj) .LT. TMIX) THEN  c          IF (TICE(I,J,bi,bj) .LT. TMIX) THEN
1710  c          IF (atemp(I,J,bi,bj) .LT. celsius2K) THEN  c          IF (atemp(I,J,bi,bj) .LT. celsius2K) THEN
1711            IF ( a_QbyATM_cover(I,J).GE. 0. _d 0 ) THEN             IF ( a_QbyATM_cover(I,J).GE. 0. _d 0 ) THEN
1712  C           add precip as snow  C           add precip as snow
1713              d_HFRWbyRAIN(I,J)=0. _d 0              d_HFRWbyRAIN(I,J)=0. _d 0
1714              d_HSNWbyRAIN(I,J)=convertPRECIP2HI*ICE2SNOW*              d_HSNWbyRAIN(I,J)=convertPRECIP2HI*ICE2SNOW*
1715       &            PRECIP(I,J,bi,bj)*AREApreTH(I,J)       &            PRECIP(I,J,bi,bj)*AREApreTH(I,J)
1716            ELSE             ELSE
1717  C           add precip to the fresh water bucket  C           add precip to the fresh water bucket
1718              d_HFRWbyRAIN(I,J)=-convertPRECIP2HI*              d_HFRWbyRAIN(I,J)=-convertPRECIP2HI*
1719       &            PRECIP(I,J,bi,bj)*AREApreTH(I,J)       &            PRECIP(I,J,bi,bj)*AREApreTH(I,J)
1720              d_HSNWbyRAIN(I,J)=0. _d 0              d_HSNWbyRAIN(I,J)=0. _d 0
1721            ENDIF             ENDIF
1722           ENDDO           ENDDO
1723          ENDDO          ENDDO
1724  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1725          DO K=1,nITD          DO IT=1,nITD
1726           DO J=1,sNy           DO J=1,sNy
1727            DO I=1,sNx            DO I=1,sNx
1728             HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj)             d_HSNWbyRAIN_ITD(I,J,IT)
1729       &     + d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,K)       &     = d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,IT)
1730            ENDDO            ENDDO
1731           ENDDO           ENDDO
1732          ENDDO          ENDDO
1733  #else  #else
1734          DO J=1,sNy          DO J=1,sNy
1735           DO I=1,sNx           DO I=1,sNx
# Line 1614  C           add precip to the fresh wate Line 1740  C           add precip to the fresh wate
1740  Cgf note: this does not affect air-sea heat flux,  Cgf note: this does not affect air-sea heat flux,
1741  Cgf since the implied air heat gain to turn  Cgf since the implied air heat gain to turn
1742  Cgf rain to snow is not a surface process.  Cgf rain to snow is not a surface process.
1743    C end of IF statement snowPrecipFile:
1744            ENDIF
1745  #endif /* ALLOW_ATM_TEMP */  #endif /* ALLOW_ATM_TEMP */
1746    #ifdef SEAICE_DEBUG
1747  c ToM<<< debug seaice_growth  c ToM<<< debug seaice_growth
1748          WRITE(msgBuf,'(A,7F6.2)')          WRITE(msgBuf,msgBufForm)
1749       &    ' SEAICE_GROWTH: Heff increments 5, HEFFITD = ',       &    ' SEAICE_GROWTH: Hsnow increments 5, d_HSNWbyRAIN = ',
1750       &     HEFFITD(20,20,:,bi,bj)  #ifdef SEAICE_ITD
1751         &     d_HSNWbyRAIN_ITD(1,1,:)
1752    #else
1753         &     d_HSNWbyRAIN(1,1)
1754    #endif
1755          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1756       &    SQUEEZE_RIGHT , myThid)       &    SQUEEZE_RIGHT , myThid)
1757  c ToM>>>  c ToM>>>
1758    #endif /* SEAICE_DEBUG */
1759    
1760  C compute snow melt due to heat available from ocean.  C compute snow melt due to heat available from ocean.
1761  C =================================================================  C =================================================================
# Line 1635  CADJ STORE r_QbyOCN = comlev1_bibj,key=i Line 1769  CADJ STORE r_QbyOCN = comlev1_bibj,key=i
1769  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1770    
1771  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1772          DO K=1,nITD          DO IT=1,nITD
1773           DO J=1,sNy           DO J=1,sNy
1774            DO I=1,sNx            DO I=1,sNx
1775             tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,K),             tmpscal4 = HSNWITDpreTH(I,J,IT)
1776       &                  -HSNOWITD(I,J,K,bi,bj))       &              + d_HSNWbySublim_ITD(I,J,IT)
1777         &              + d_HSNWbyATMonSNW_ITD(I,J,IT)
1778         &              + d_HSNWbyRAIN_ITD(I,J,IT)
1779               tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,IT),
1780         &                  -tmpscal4)
1781             tmpscal2=MIN(tmpscal1,0. _d 0)             tmpscal2=MIN(tmpscal1,0. _d 0)
1782  #ifdef SEAICE_MODIFY_GROWTH_ADJ  #ifdef SEAICE_MODIFY_GROWTH_ADJ
1783  Cgf no additional dependency through snow  Cgf no additional dependency through snow
1784             if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0             if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1785  #endif  #endif
1786               d_HSNWbyOCNonSNW_ITD(I,J,IT) = tmpscal2
1787             d_HSNWbyOCNonSNW(I,J) = d_HSNWbyOCNonSNW(I,J) + tmpscal2             d_HSNWbyOCNonSNW(I,J) = d_HSNWbyOCNonSNW(I,J) + tmpscal2
1788             r_QbyOCN(I,J)=r_QbyOCN(I,J) - tmpscal2*SNOW2ICE             r_QbyOCN(I,J)=r_QbyOCN(I,J) - tmpscal2*SNOW2ICE
            HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) + tmpscal2  
1789            ENDDO            ENDDO
1790           ENDDO           ENDDO
1791          ENDDO          ENDDO
# Line 1669  Cgf no additional dependency through sno Line 1807  Cgf no additional dependency through sno
1807  #endif /* SEAICE_ITD */  #endif /* SEAICE_ITD */
1808  #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */  #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */
1809  Cph)  Cph)
1810    #ifdef SEAICE_DEBUG
1811  c ToM<<< debug seaice_growth  c ToM<<< debug seaice_growth
1812          WRITE(msgBuf,'(A,7F6.2)')          WRITE(msgBuf,msgBufForm)
1813       &    ' SEAICE_GROWTH: Heff increments 6, HEFFITD = ',       &    ' SEAICE_GROWTH: Hsnow increments 6, d_HSNWbyOCNonSNW = ',
1814       &     HEFFITD(20,20,:,bi,bj)  #ifdef SEAICE_ITD
1815         &     d_HSNWbyOCNonSNW_ITD(1,1,:)
1816    #else
1817         &     d_HSNWbyOCNonSNW(1,1)
1818    #endif
1819          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1820       &    SQUEEZE_RIGHT , myThid)       &    SQUEEZE_RIGHT , myThid)
1821  c ToM>>>  c ToM>>>
1822    #endif /* SEAICE_DEBUG */
1823    
1824  C gain of new ice over open water  C gain of new ice over open water
1825  C ===============================  C ===============================
# Line 1689  CADJ STORE a_QSWbyATM_open = comlev1_bib Line 1833  CADJ STORE a_QSWbyATM_open = comlev1_bib
1833    
1834          DO J=1,sNy          DO J=1,sNy
1835           DO I=1,sNx           DO I=1,sNx
1836  c           Initial ice growth is triggered by open water  #ifdef SEAICE_ITD
1837  c           heat flux overcoming potential melt by ocean  C         HEFF will be updated at the end of PART 3,
1838    C         hence sum of tendencies so far is needed
1839              tmpscal4 = HEFFpreTH(I,J)
1840         &             + d_HEFFbySublim(I,J)
1841         &             + d_HEFFbyOCNonICE(I,J)
1842         &             + d_HEFFbyATMonOCN(I,J)
1843    #else
1844    C         HEFF is updated step by step throughout seaice_growth
1845              tmpscal4 = HEFF(I,J,bi,bj)
1846    #endif
1847    C           Initial ice growth is triggered by open water
1848    C           heat flux overcoming potential melt by ocean
1849              tmpscal1=r_QbyATM_open(I,J)+r_QbyOCN(i,j) *              tmpscal1=r_QbyATM_open(I,J)+r_QbyOCN(i,j) *
1850       &            (1.0 _d 0 - AREApreTH(I,J))       &            (1.0 _d 0 - AREApreTH(I,J))
1851  c           Penetrative shortwave flux beyond first layer  C           Penetrative shortwave flux beyond first layer
1852  c           that is therefore not available to ice growth/melt  C           that is therefore not available to ice growth/melt
1853              tmpscal2=SWFracB * a_QSWbyATM_open(I,J)              tmpscal2=SWFracB * a_QSWbyATM_open(I,J)
1854  C           impose -HEFF as the maxmum melting if SEAICE_doOpenWaterMelt  C           impose -HEFF as the maxmum melting if SEAICE_doOpenWaterMelt
1855  C           or 0. otherwise (no melting if not SEAICE_doOpenWaterMelt)  C           or 0. otherwise (no melting if not SEAICE_doOpenWaterMelt)
1856              tmpscal3=facOpenGrow*MAX(tmpscal1-tmpscal2,              tmpscal3=facOpenGrow*MAX(tmpscal1-tmpscal2,
1857       &           -HEFF(I,J,bi,bj)*facOpenMelt)*HEFFM(I,J,bi,bj)       &           -tmpscal4*facOpenMelt)*HEFFM(I,J,bi,bj)
1858    #ifdef SEAICE_ITD
1859    C         ice growth in open water adds to first category
1860              d_HEFFbyATMonOCN_open_ITD(I,J,1)=tmpscal3
1861              d_HEFFbyATMonOCN_ITD(I,J,1)     =d_HEFFbyATMonOCN_ITD(I,J,1)
1862         &                                    +tmpscal3
1863    #endif
1864            d_HEFFbyATMonOCN_open(I,J)=tmpscal3            d_HEFFbyATMonOCN_open(I,J)=tmpscal3
1865            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3
1866            r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3            r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3
 #ifdef SEAICE_ITD  
 C         open water area fraction  
           tmpscal0 = ONE-AREApreTH(I,J)  
 C         determine thickness of new ice  
 C         considering the entire open water area to refreeze  
           tmpscal1 = tmpscal3/tmpscal0  
 C         then add new ice volume to appropriate thickness category  
           DO K=1,nITD  
            IF (tmpscal1.LT.Hlimit(K)) THEN  
             HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal3  
             tmpscal3=ZERO  
 C not sure if AREAITD should be changed here since AREA is incremented  
 C   in PART 4 below in non-itd code  
 C in this scenario all open water is covered by new ice instantaneously,  
 C   i.e. no delayed lead closing is concidered such as is associated with  
 C   Hibler's h_0 parameter  
             AREAITD(I,J,K,bi,bj) = AREAITD(I,J,K,bi,bj)  
      &                           + tmpscal0  
             tmpscal0=ZERO  
            ENDIF  
           ENDDO  
 #else  
1867            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3
 #endif  
1868           ENDDO           ENDDO
1869          ENDDO          ENDDO
1870    
1871  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
 #ifdef SEAICE_ITD  
         DO K=1,nITD  
          DO J=1,sNy  
           DO I=1,sNx  
 c needs to be here to allow use also with LEGACY branch  
            SItrHEFF(I,J,bi,bj,4) = SItrHEFF(I,J,bi,bj,4)  
      &                           + HEFFITD(I,J,K,bi,bj)  
           ENDDO  
          ENDDO  
         ENDDO  
 #else  
1872          DO J=1,sNy          DO J=1,sNy
1873           DO I=1,sNx           DO I=1,sNx
1874  c needs to be here to allow use also with LEGACY branch  C needs to be here to allow use also with LEGACY branch
1875    #ifdef SEAICE_ITD
1876              SItrHEFF(I,J,bi,bj,4)=SItrHEFF(I,J,bi,bj,3)
1877         &                         +d_HEFFbyATMonOCN_open(I,J)
1878    #else
1879            SItrHEFF(I,J,bi,bj,4)=HEFF(I,J,bi,bj)            SItrHEFF(I,J,bi,bj,4)=HEFF(I,J,bi,bj)
1880    #endif
1881           ENDDO           ENDDO
1882          ENDDO          ENDDO
 #endif  
1883  #endif /* ALLOW_SITRACER */  #endif /* ALLOW_SITRACER */
1884    #ifdef SEAICE_DEBUG
1885  c ToM<<< debug seaice_growth  c ToM<<< debug seaice_growth
1886          WRITE(msgBuf,'(A,7F6.2)')          WRITE(msgBuf,msgBufForm)
1887       &    ' SEAICE_GROWTH: Heff increments 7, HEFFITD = ',       &    ' SEAICE_GROWTH: Heff increments 7, d_HEFFbyATMonOCN_open = ',
1888       &     HEFFITD(20,20,:,bi,bj)  #ifdef SEAICE_ITD
1889         &     d_HEFFbyATMonOCN_open_ITD(1,1,:)
1890    #else
1891         &     d_HEFFbyATMonOCN_open(1,1)
1892    #endif
1893            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1894         &    SQUEEZE_RIGHT , myThid)
1895            WRITE(msgBuf,msgBufForm)
1896         &    ' SEAICE_GROWTH: Heff increments 7, d_HEFFbyATMonOCN = ',
1897    #ifdef SEAICE_ITD
1898         &     d_HEFFbyATMonOCN_ITD(1,1,:)
1899    #else
1900         &     d_HEFFbyATMonOCN(1,1)
1901    #endif
1902          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1903       &    SQUEEZE_RIGHT , myThid)       &    SQUEEZE_RIGHT , myThid)
1904  c ToM>>>  c ToM>>>
1905    #endif /* SEAICE_DEBUG */
1906    
1907  C convert snow to ice if submerged.  C convert snow to ice if submerged.
1908  C =================================  C =================================
# Line 1769  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 1915  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
1915  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
1916          IF ( SEAICEuseFlooding ) THEN          IF ( SEAICEuseFlooding ) THEN
1917  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
1918           DO K=1,nITD           DO IT=1,nITD
1919            DO J=1,sNy            DO J=1,sNy
1920             DO I=1,sNx             DO I=1,sNx
1921              tmpscal0 = (HSNOWITD(I,J,K,bi,bj)*SEAICE_rhoSnow              tmpscal3 = HEFFITDpreTH(I,J,IT)
1922       &               +HEFFITD(I,J,K,bi,bj)*SEAICE_rhoIce)*recip_rhoConst       &               + d_HEFFbySublim_ITD(I,J,IT)
1923              tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFFITD(I,J,K,bi,bj))       &               + d_HEFFbyOCNonICE_ITD(I,J,IT)
1924              d_HEFFbyFLOODING(I,J) = d_HEFFbyFLOODING(I,J) + tmpscal1       &               + d_HEFFbyATMonOCN_ITD(I,J,IT)
1925              HEFFITD(I,J,K,bi,bj)  = HEFFITD(I,J,K,bi,bj)  + tmpscal1              tmpscal4 = HSNWITDpreTH(I,J,IT)
1926              HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) - tmpscal1       &               + d_HSNWbySublim_ITD(I,J,IT)
1927       &                            * ICE2SNOW       &               + d_HSNWbyATMonSNW_ITD(I,J,IT)
1928             ENDDO       &               + d_HSNWbyRAIN_ITD(I,J,IT)
1929            ENDDO              tmpscal0 = (tmpscal4*SEAICE_rhoSnow
1930           ENDDO       &               +  tmpscal3*SEAICE_rhoIce)
1931         &               * recip_rhoConst
1932                tmpscal1 = MAX( 0. _d 0, tmpscal0 - tmpscal3)
1933                d_HEFFbyFLOODING_ITD(I,J,IT) = tmpscal1
1934                d_HEFFbyFLOODING(I,J) = d_HEFFbyFLOODING(I,J)  + tmpscal1
1935               ENDDO
1936              ENDDO
1937             ENDDO
1938  #else  #else
1939           DO J=1,sNy           DO J=1,sNy
1940            DO I=1,sNx            DO I=1,sNx
# Line 1792  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 1945  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
1945             HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)+d_HEFFbyFLOODING(I,J)             HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)+d_HEFFbyFLOODING(I,J)
1946             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-
1947       &                           d_HEFFbyFLOODING(I,J)*ICE2SNOW       &                           d_HEFFbyFLOODING(I,J)*ICE2SNOW
1948            ENDDO            ENDDO
1949           ENDDO           ENDDO
1950  #endif  #endif
1951          ENDIF          ENDIF
1952  #endif /* SEAICE_GROWTH_LEGACY */  #endif /* SEAICE_GROWTH_LEGACY */
1953    #ifdef SEAICE_DEBUG
1954    c ToM<<< debug seaice_growth
1955            WRITE(msgBuf,msgBufForm)
1956         &    ' SEAICE_GROWTH: Heff increments 8, d_HEFFbyFLOODING = ',
1957    #ifdef SEAICE_ITD
1958         &     d_HEFFbyFLOODING_ITD(1,1,:)
1959    #else
1960         &     d_HEFFbyFLOODING(1,1)
1961    #endif
1962            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1963         &    SQUEEZE_RIGHT , myThid)
1964    c ToM>>>
1965    #endif /* SEAICE_DEBUG */
1966    #ifdef SEAICE_ITD
1967    C apply ice and snow thickness changes
1968    C =================================================================
1969             DO IT=1,nITD
1970              DO J=1,sNy
1971               DO I=1,sNx
1972                HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj)
1973         &                            + d_HEFFbySublim_ITD(I,J,IT)
1974         &                            + d_HEFFbyOCNonICE_ITD(I,J,IT)
1975         &                            + d_HEFFbyATMonOCN_ITD(I,J,IT)
1976         &                            + d_HEFFbyFLOODING_ITD(I,J,IT)
1977                HSNOWITD(I,J,IT,bi,bj) = HSNOWITD(I,J,IT,bi,bj)
1978         &                            + d_HSNWbySublim_ITD(I,J,IT)
1979         &                            + d_HSNWbyATMonSNW_ITD(I,J,IT)
1980         &                            + d_HSNWbyRAIN_ITD(I,J,IT)
1981         &                            + d_HSNWbyOCNonSNW_ITD(I,J,IT)
1982         &                            - d_HEFFbyFLOODING_ITD(I,J,IT)
1983         &                            * ICE2SNOW
1984               ENDDO
1985              ENDDO
1986             ENDDO
1987    #endif
1988  c ToM<<< debug seaice_growth  c ToM<<< debug seaice_growth
1989          WRITE(msgBuf,'(A,7F6.2)')          WRITE(msgBuf,msgBufForm)
1990       &    ' SEAICE_GROWTH: Heff increments 8, HEFFITD = ',       &    ' SEAICE_GROWTH: Heff increments 9, HEFF = ',
1991       &     HEFFITD(20,20,:,bi,bj)  #ifdef SEAICE_ITD
1992         &     HEFFITD(1,1,:,bi,bj)
1993    #else
1994         &     HEFF(1,1,bi,bj)
1995    #endif
1996            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1997         &    SQUEEZE_RIGHT , myThid)
1998            WRITE(msgBuf,msgBufForm)
1999         &    ' SEAICE_GROWTH: Area increments 9, AREA = ',
2000    #ifdef SEAICE_ITD
2001         &     AREAITD(1,1,:,bi,bj)
2002    #else
2003         &     AREA(1,1,bi,bj)
2004    #endif
2005          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
2006       &    SQUEEZE_RIGHT , myThid)       &    SQUEEZE_RIGHT , myThid)
2007  c ToM>>>  c ToM>>>
# Line 1830  CADJ STORE AREA(:,:,bi,bj) = comlev1_bib Line 2031  CADJ STORE AREA(:,:,bi,bj) = comlev1_bib
2031  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
2032    
2033  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
2034  C       replaces Hibler '79 scheme and lead closing parameter  C--     account for lateral ice growth and melt only in thinnest category
2035  C       because ITD accounts explicitly for lead openings and  C--     use HEFF, ARE, HSNOW, etc. temporarily for 1st category
2036  C       different melt rates due to varying ice thickness  C       (this way we can use same code for ITD and non-ITD case)
2037  C          DO J=1,sNy
2038  C       only consider ice area loss due to total ice thickness loss           DO I=1,sNx
2039  C       ice area gain due to freezing of open water as handled above            HEFF(I,J,bi,bj)=HEFFITD(I,J,1,bi,bj)
2040  C       under "gain of new ice over open water"            AREA(I,J,bi,bj)=AREAITD(I,J,1,bi,bj)
2041  C            HSNOW(I,J,bi,bj)=HSNOWITD(I,J,1,bi,bj)
2042  C       does not account for lateral melt of ice floes            HEFFpreTH(I,J)=HEFFITDpreTH(I,J,1)
2043  C            AREApreTH(I,J)=AREAITDpreTH(I,J,1)
2044  C        AREAITD is incremented in section "gain of new ice over open water" above            recip_heffActual(I,J)=recip_heffActualMult(I,J,1)
2045  C           ENDDO
2046          DO K=1,nITD          ENDDO
2047           DO J=1,sNy  C       all other categories only experience basal growth or melt,
2048            DO I=1,sNx  C       i.e. change sin AREA only occur when all ice in a category is melted
2049             IF (HEFFITD(I,J,K,bi,bj).LE.ZERO) THEN          IF (nITD .gt. 1) THEN
2050              AREAITD(I,J,K,bi,bj)=ZERO           DO IT=2,nITD
2051             ENDIF            DO J=1,sNy
2052  #ifdef ALLOW_SITRACER             DO I=1,sNx
2053             SItrAREA(I,J,bi,bj,3) = SItrAREA(I,J,bi,bj,3)              IF (HEFFITD(I,J,IT,bi,bj).LE.ZERO) THEN
2054       &                           + AREAITD(I,J,K,bi,bj)               AREAITD(I,J,IT,bi,bj)=ZERO
2055  #endif /* ALLOW_SITRACER */              ENDIF
2056            ENDDO             ENDDO
2057           ENDDO            ENDDO
2058          ENDDO           ENDDO
2059  #else /* SEAICE_ITD */          ENDIF
2060    #endif
2061          DO J=1,sNy          DO J=1,sNy
2062           DO I=1,sNx           DO I=1,sNx
2063    
# Line 1925  C apply tendency Line 2127  C apply tendency
2127  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
2128           ENDDO           ENDDO
2129          ENDDO          ENDDO
2130  #endif /* SEAICE_ITD */  #ifdef SEAICE_ITD
2131    C       transfer 1st category values back into ITD variables
2132            DO J=1,sNy
2133             DO I=1,sNx
2134              HEFFITD(I,J,1,bi,bj)=HEFF(I,J,bi,bj)
2135              AREAITD(I,J,1,bi,bj)=AREA(I,J,bi,bj)
2136              HSNOWITD(I,J,1,bi,bj)=HSNOW(I,J,bi,bj)
2137             ENDDO
2138            ENDDO
2139    #endif
2140    
2141  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
2142  Cgf 'bulk' linearization of area=f(HEFF)  Cgf 'bulk' linearization of area=f(HEFF)
2143          IF ( SEAICEadjMODE.GE.1 ) THEN          IF ( SEAICEadjMODE.GE.1 ) THEN
2144  #ifdef SEAICE_ITD  #ifdef SEAICE_ITD
2145           DO K=1,nITD           DO IT=1,nITD
2146            DO J=1,sNy            DO J=1,sNy
2147             DO I=1,sNx             DO I=1,sNx
2148              AREAITD(I,J,K,bi,bj) = AREAITDpreTH(I,J,K) + 0.1 _d 0 *              AREAITD(I,J,IT,bi,bj) = AREAITDpreTH(I,J,IT) + 0.1 _d 0 *
2149       &               ( HEFFITD(I,J,K,bi,bj) - HEFFITDpreTH(I,J,K) )       &               ( HEFFITD(I,J,IT,bi,bj) - HEFFITDpreTH(I,J,IT) )
2150             ENDDO             ENDDO
2151            ENDDO            ENDDO
2152           ENDDO           ENDDO
# Line 1978  CADJ &                       key = iicek Line 2189  CADJ &                       key = iicek
2189            tmpscal1 = d_HEFFbyNEG(I,J) + d_HEFFbyOCNonICE(I,J) +            tmpscal1 = d_HEFFbyNEG(I,J) + d_HEFFbyOCNonICE(I,J) +
2190       &               d_HEFFbyATMonOCN(I,J) + d_HEFFbyFLOODING(I,J)       &               d_HEFFbyATMonOCN(I,J) + d_HEFFbyFLOODING(I,J)
2191       &             + d_HEFFbySublim(I,J)       &             + d_HEFFbySublim(I,J)
2192  #ifdef SEAICE_ALLOW_AREA_RELAXATION  #ifdef EXF_ALLOW_SEAICE_RELAX
2193                     + d_HEFFbyRLX(I,J)       &             + d_HEFFbyRLX(I,J)
2194  #endif  #endif
2195            tmpscal2 = tmpscal1 * SEAICE_salt0 * HEFFM(I,J,bi,bj)            tmpscal2 = tmpscal1 * SEAICE_salt0 * HEFFM(I,J,bi,bj)
2196       &            * recip_deltaTtherm * SEAICE_rhoIce       &            * recip_deltaTtherm * SEAICE_rhoIce
# Line 2073  C set HSALT = 0 if HEFF = 0 and compute Line 2284  C set HSALT = 0 if HEFF = 0 and compute
2284          ENDDO          ENDDO
2285  #endif /* SEAICE_VARIABLE_SALINITY */  #endif /* SEAICE_VARIABLE_SALINITY */
2286    
   
2287  C =======================================================================  C =======================================================================
2288  C ==LEGACY PART 6 (LEGACY) treat pathological cases, then do flooding ===  C ==LEGACY PART 6 (LEGACY) treat pathological cases, then do flooding ===
2289  C =======================================================================  C =======================================================================
# Line 2156  C ================================= Line 2366  C =================================
2366  #ifdef ALLOW_SITRACER  #ifdef ALLOW_SITRACER
2367          DO J=1,sNy          DO J=1,sNy
2368           DO I=1,sNx           DO I=1,sNx
2369  c needs to be here to allow use also with LEGACY branch  C needs to be here to allow use also with LEGACY branch
2370            SItrHEFF(I,J,bi,bj,5)=HEFF(I,J,bi,bj)            SItrHEFF(I,J,bi,bj,5)=HEFF(I,J,bi,bj)
2371           ENDDO           ENDDO
2372          ENDDO          ENDDO
# Line 2170  C compute net heat flux leaving/entering Line 2380  C compute net heat flux leaving/entering
2380  C accounting for the part used in melt/freeze processes  C accounting for the part used in melt/freeze processes
2381  C =====================================================  C =====================================================
2382    
2383    #ifdef SEAICE_ITD
2384    C compute total of "mult" fluxes for ocean forcing
2385            DO J=1,sNy
2386             DO I=1,sNx
2387              a_QbyATM_cover(I,J)   = 0.0 _d 0
2388              r_QbyATM_cover(I,J)   = 0.0 _d 0
2389              a_QSWbyATM_cover(I,J) = 0.0 _d 0
2390              r_FWbySublim(I,J)     = 0.0 _d 0
2391             ENDDO
2392            ENDDO
2393            DO IT=1,nITD
2394             DO J=1,sNy
2395              DO I=1,sNx
2396    cToM if fluxes in W/m^2 then
2397    c           a_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)
2398    c     &      + a_QbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
2399    c           r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)
2400    c     &      + r_QbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
2401    c           a_QSWbyATM_cover(I,J)=a_QSWbyATM_cover(I,J)
2402    c     &      + a_QSWbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
2403    c           r_FWbySublim(I,J)=r_FWbySublim(I,J)
2404    c     &      + r_FWbySublimMult(I,J,IT) * areaFracFactor(I,J,IT)
2405    cToM if fluxes in effective ice meters, i.e. ice volume per area, then
2406               a_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)
2407         &      + a_QbyATMmult_cover(I,J,IT)
2408               r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)
2409         &      + r_QbyATMmult_cover(I,J,IT)
2410               a_QSWbyATM_cover(I,J)=a_QSWbyATM_cover(I,J)
2411         &      + a_QSWbyATMmult_cover(I,J,IT)
2412               r_FWbySublim(I,J)=r_FWbySublim(I,J)
2413         &      + r_FWbySublimMult(I,J,IT)
2414              ENDDO
2415             ENDDO
2416            ENDDO
2417    #endif
2418    
2419  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
2420  CADJ STORE d_hsnwbyneg = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE d_hsnwbyneg = comlev1_bibj,key=iicekey,byte=isbyte
2421  CADJ STORE d_hsnwbyocnonsnw = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE d_hsnwbyocnonsnw = comlev1_bibj,key=iicekey,byte=isbyte
# Line 2183  C in principle a_QSWbyATM_cover should a Line 2429  C in principle a_QSWbyATM_cover should a
2429  C for backward compatibility it is left out of the LEGACY branch  C for backward compatibility it is left out of the LEGACY branch
2430       &         +   a_QSWbyATM_cover(I,J)       &         +   a_QSWbyATM_cover(I,J)
2431  #endif /* SEAICE_GROWTH_LEGACY */  #endif /* SEAICE_GROWTH_LEGACY */
2432       &         - ( d_HEFFbyOCNonICE(I,J) +       &         - ( d_HEFFbyOCNonICE(I,J)
2433       &             d_HSNWbyOCNonSNW(I,J)*SNOW2ICE +       &           + d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
2434       &             d_HEFFbyNEG(I,J) +       &           + d_HEFFbyNEG(I,J)
2435  #ifdef SEAICE_ALLOW_AREA_RELAXATION  #ifdef EXF_ALLOW_SEAICE_RELAX
2436       &             d_HEFFbyRLX(I,J) +       &           + d_HEFFbyRLX(I,J)
2437  #endif  #endif
2438       &             d_HSNWbyNEG(I,J)*SNOW2ICE )       &           + d_HSNWbyNEG(I,J)*SNOW2ICE
2439       &         * maskC(I,J,kSurface,bi,bj)       &           - convertPRECIP2HI *
2440         &             snowPrecip(i,j,bi,bj) * (ONE-AREApreTH(I,J))
2441         &           ) * maskC(I,J,kSurface,bi,bj)
2442             ENDDO
2443            ENDDO
2444            DO J=1,sNy
2445             DO I=1,sNx
2446            QSW(I,J,bi,bj)  = a_QSWbyATM_cover(I,J) + a_QSWbyATM_open(I,J)            QSW(I,J,bi,bj)  = a_QSWbyATM_cover(I,J) + a_QSWbyATM_open(I,J)
2447           ENDDO           ENDDO
2448          ENDDO          ENDDO
# Line 2218  CADJ STORE d_HSNWbyATMonSNW = comlev1_bi Line 2470  CADJ STORE d_HSNWbyATMonSNW = comlev1_bi
2470  CADJ STORE theta(:,:,kSurface,bi,bj) = comlev1_bibj,  CADJ STORE theta(:,:,kSurface,bi,bj) = comlev1_bibj,
2471  CADJ &                       key = iicekey, byte = isbyte  CADJ &                       key = iicekey, byte = isbyte
2472  # endif /* ALLOW_AUTODIFF_TAMC */  # endif /* ALLOW_AUTODIFF_TAMC */
2473        IF ( SEAICEheatConsFix ) THEN  cgf Unlike for evap and precip, the temperature of gained/lost
2474  c Unlike for evap and precip, the temperature of gained/lost  C ocean liquid water due to melt/freeze of solid water cannot be chosen
2475  c ocean liquid water due to melt/freeze of solid water cannot be chosen  C arbitrarily to be e.g. the ocean SST. Indeed the present seaice model
2476  c to be e.g. the ocean SST. It must be done at 0degC. The fix below anticipates  C implies a constant ice temperature of 0degC. If melt/freeze water is exchanged
2477  c on external_forcing_surf.F and applies the correction to QNET.  C at a different temperature, it leads to a loss of conservation in the
2478        IF ((convertFW2Salt.EQ.-1.).OR.(temp_EvPrRn.EQ.UNSET_RL)) THEN  C ocean+ice system. While this is mostly a serious issue in the
2479  c I leave alone the exotic case when onvertFW2Salt.NE.-1 and temp_EvPrRn.NE.UNSET_RL and  C real fresh water + non linear free surface framework, a mismatch
2480  c the small error of the synchronous time stepping case (see external_forcing_surf.F).  C between ice and ocean boundary condition can result in all cases.
2481    C Below we therefore anticipate on external_forcing_surf.F
2482    C to diagnoze and/or apply the correction to QNET.
2483          DO J=1,sNy          DO J=1,sNy
2484           DO I=1,sNx           DO I=1,sNx
2485  #ifdef ALLOW_DIAGNOSTICS  C ocean water going to ice/snow, in precip units
2486  c store unaltered QNET for diagnostic purposes             tmpscal3=rhoConstFresh*maskC(I,J,kSurface,bi,bj)*(
            DIAGarrayA(I,J)=QNET(I,J,bi,bj)  
 #endif  
 c compute the ocean water going to ice/snow, in precip units  
            tmpscal3=rhoConstFresh*maskC(I,J,kSurface,bi,bj)*  
2487       &       ( d_HSNWbyATMonSNW(I,J)*SNOW2ICE       &       ( d_HSNWbyATMonSNW(I,J)*SNOW2ICE
2488       &       + d_HSNWbyOCNonSNW(I,J)*SNOW2ICE       &       + d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
2489       &       + d_HEFFbyOCNonICE(I,J) + d_HEFFbyATMonOCN(I,J)       &       + d_HEFFbyOCNonICE(I,J) + d_HEFFbyATMonOCN(I,J)
2490       &       + d_HEFFbyNEG(I,J) + d_HSNWbyNEG(I,J)*SNOW2ICE )       &       + d_HEFFbyNEG(I,J) + d_HSNWbyNEG(I,J)*SNOW2ICE )
2491       &       * convertHI2PRECIP       &       * convertHI2PRECIP
2492  c factor in the heat content that external_forcing_surf.F       &       - snowPrecip(i,j,bi,bj) * (ONE-AREApreTH(I,J)) )
2493  c will associate with EMPMR, and remove it from QNET, so that  C factor in the heat content as done in external_forcing_surf.F
2494  c melt/freez water is in effect consistently gained/lost at 0degC        IF ( (temp_EvPrRn.NE.UNSET_RL).AND.useRealFreshWaterFlux
2495             IF (temp_EvPrRn.NE.UNSET_RL) THEN       &     .AND.(nonlinFreeSurf.NE.0) ) THEN
2496               QNET(I,J,bi,bj)=QNET(I,J,bi,bj) - tmpscal3*               tmpscal1 = - tmpscal3*
2497       &         HeatCapacity_Cp * temp_EvPrRn       &         HeatCapacity_Cp * temp_EvPrRn
2498             ELSE        ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL).AND.useRealFreshWaterFlux
2499               QNET(I,J,bi,bj)=QNET(I,J,bi,bj) - tmpscal3*       &         .AND.(nonlinFreeSurf.NE.0) ) THEN
2500                 tmpscal1 = - tmpscal3*
2501       &         HeatCapacity_Cp * theta(I,J,kSurface,bi,bj)       &         HeatCapacity_Cp * theta(I,J,kSurface,bi,bj)
2502             ENDIF        ELSEIF ( (temp_EvPrRn.NE.UNSET_RL) ) THEN
2503               tmpscal1 = - tmpscal3*HeatCapacity_Cp*
2504         &       ( temp_EvPrRn - theta(I,J,kSurface,bi,bj) )
2505          ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL) ) THEN
2506               tmpscal1 = ZERO
2507          ENDIF
2508  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
2509  c back out the eventual TFLUX adjustement and fill diag  C in all cases, diagnoze the boundary condition mismatch to SIaaflux
2510             DIAGarrayA(I,J)=QNET(I,J,bi,bj)-DIAGarrayA(I,J)             DIAGarrayA(I,J)=tmpscal1
2511  #endif  #endif
2512    C remove the mismatch when real fresh water is exchanged (at 0degC here)
2513          IF ( useRealFreshWaterFlux.AND.(nonlinFreeSurf.GT.0).AND.
2514         &   SEAICEheatConsFix ) QNET(I,J,bi,bj)=QNET(I,J,bi,bj)+tmpscal1
2515           ENDDO           ENDDO
2516          ENDDO          ENDDO
       ENDIF  
2517  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
2518          CALL DIAGNOSTICS_FILL(DIAGarrayA,          CALL DIAGNOSTICS_FILL(DIAGarrayA,
2519       &      'SIaaflux',0,1,3,bi,bj,myThid)       &      'SIaaflux',0,1,3,bi,bj,myThid)
2520  #endif  #endif
       ENDIF  
2521  #endif /* ndef SEAICE_DISABLE_HEATCONSFIX */  #endif /* ndef SEAICE_DISABLE_HEATCONSFIX */
2522    
2523    C compute the net heat flux, incl. adv. by water, entering ocean+ice
2524    C ===================================================================
2525             DO J=1,sNy
2526              DO I=1,sNx
2527    cgf 1) SIatmQnt (analogous to qnet; excl. adv. by water exch.)
2528    CML If I consider the atmosphere above the ice, the surface flux
2529    CML which is relevant for the air temperature dT/dt Eq
2530    CML accounts for sensible and radiation (with different treatment
2531    CML according to wave-length) fluxes but not for "latent heat flux",
2532    CML since it does not contribute to heating the air.
2533    CML So this diagnostic is only good for heat budget calculations within
2534    CML the ice-ocean system.
2535               SIatmQnt(I,J,bi,bj) =
2536         &            maskC(I,J,kSurface,bi,bj)*convertHI2Q*(
2537    #ifndef SEAICE_GROWTH_LEGACY
2538         &            a_QSWbyATM_cover(I,J) +
2539    #endif /* SEAICE_GROWTH_LEGACY */
2540         &            a_QbyATM_cover(I,J) + a_QbyATM_open(I,J) )
2541    cgf 2) SItflux (analogous to tflux; includes advection by water
2542    C             exchanged between atmosphere and ocean+ice)
2543    C solid water going to atm, in precip units
2544               tmpscal1 = rhoConstFresh*maskC(I,J,kSurface,bi,bj)
2545         &       * convertHI2PRECIP * ( - d_HSNWbyRAIN(I,J)*SNOW2ICE
2546         &       + a_FWbySublim(I,J) - r_FWbySublim(I,J) )
2547    C liquid water going to atm, in precip units
2548               tmpscal2=rhoConstFresh*maskC(I,J,kSurface,bi,bj)*
2549         &       ( ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )
2550         &         * ( ONE - AREApreTH(I,J) )
2551    #ifdef ALLOW_RUNOFF
2552         &         - RUNOFF(I,J,bi,bj)
2553    #endif /* ALLOW_RUNOFF */
2554         &         + ( d_HFRWbyRAIN(I,J) + r_FWbySublim(I,J) )
2555         &         *convertHI2PRECIP )
2556    C In real fresh water flux + nonlinFS, we factor in the advected specific
2557    C energy (referenced to 0 for 0deC liquid water). In virtual salt flux or
2558    C linFS, rain/evap get a special treatment (see external_forcing_surf.F).
2559               tmpscal1= - tmpscal1*
2560         &       ( -SEAICE_lhFusion + HeatCapacity_Cp * ZERO )
2561          IF ( (temp_EvPrRn.NE.UNSET_RL).AND.useRealFreshWaterFlux
2562         &     .AND.(nonlinFreeSurf.NE.0) ) THEN
2563               tmpscal2= - tmpscal2*
2564         &       ( ZERO + HeatCapacity_Cp * temp_EvPrRn )
2565          ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL).AND.useRealFreshWaterFlux
2566         &        .AND.(nonlinFreeSurf.NE.0) ) THEN
2567               tmpscal2= - tmpscal2*
2568         &       ( ZERO + HeatCapacity_Cp * theta(I,J,kSurface,bi,bj) )
2569          ELSEIF ( (temp_EvPrRn.NE.UNSET_RL) ) THEN
2570               tmpscal2= - tmpscal2*HeatCapacity_Cp*
2571         &       ( temp_EvPrRn - theta(I,J,kSurface,bi,bj) )
2572          ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL) ) THEN
2573               tmpscal2= ZERO
2574          ENDIF
2575               SItflux(I,J,bi,bj)=
2576         &            SIatmQnt(I,J,bi,bj)-tmpscal1-tmpscal2
2577              ENDDO
2578             ENDDO
2579    
2580  C compute net fresh water flux leaving/entering  C compute net fresh water flux leaving/entering
2581  C the ocean, accounting for fresh/salt water stocks.  C the ocean, accounting for fresh/salt water stocks.
2582  C ==================================================  C ==================================================
# Line 2276  C ====================================== Line 2590  C ======================================
2590       &             +d_HEFFbyOCNonICE(I,J)       &             +d_HEFFbyOCNonICE(I,J)
2591       &             +d_HEFFbyATMonOCN(I,J)       &             +d_HEFFbyATMonOCN(I,J)
2592       &             +d_HEFFbyNEG(I,J)       &             +d_HEFFbyNEG(I,J)
2593  #ifdef SEAICE_ALLOW_AREA_RELAXATION  #ifdef EXF_ALLOW_SEAICE_RELAX
2594       &             +d_HEFFbyRLX(I,J)       &             +d_HEFFbyRLX(I,J)
2595  #endif  #endif
2596       &             +d_HSNWbyNEG(I,J)*SNOW2ICE       &             +d_HSNWbyNEG(I,J)*SNOW2ICE
# Line 2290  C     If r_FWbySublim>0, then it is evap Line 2604  C     If r_FWbySublim>0, then it is evap
2604  #endif /* ALLOW_RUNOFF */  #endif /* ALLOW_RUNOFF */
2605       &         + tmpscal1*convertHI2PRECIP       &         + tmpscal1*convertHI2PRECIP
2606       &         )*rhoConstFresh       &         )*rhoConstFresh
2607    c and the flux leaving/entering the ocean+ice
2608               SIatmFW(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
2609         &          EVAP(I,J,bi,bj)*( ONE - AREApreTH(I,J) )
2610         &          - PRECIP(I,J,bi,bj)
2611    #ifdef ALLOW_RUNOFF
2612         &          - RUNOFF(I,J,bi,bj)
2613    #endif /* ALLOW_RUNOFF */
2614         &           )*rhoConstFresh
2615         &     + a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm
2616    
2617           ENDDO           ENDDO
2618          ENDDO          ENDDO
2619    
# Line 2362  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi Line 2686  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
2686           ENDDO           ENDDO
2687          ENDIF          ENDIF
2688    
2689    #ifdef ALLOW_BALANCE_FLUXES
2690    C Compute tile integrals of heat/fresh water fluxes to/from atm.
2691    C ==============================================================
2692           FWFsiTile(bi,bj) = 0. _d 0
2693           IF ( balanceEmPmR ) THEN
2694            DO j=1,sNy
2695             DO i=1,sNx
2696              FWFsiTile(bi,bj) =
2697         &      FWFsiTile(bi,bj) + SIatmFW(i,j,bi,bj)
2698         &      * rA(i,j,bi,bj) * maskInC(i,j,bi,bj)
2699             ENDDO
2700            ENDDO
2701           ENDIF
2702    c to translate global mean FWF adjustements (see below) we may need :
2703           FWF2HFsiTile(bi,bj) = 0. _d 0
2704           IF ( balanceEmPmR.AND.(temp_EvPrRn.EQ.UNSET_RL) ) THEN
2705            DO j=1,sNy
2706             DO i=1,sNx
2707              FWF2HFsiTile(bi,bj) = FWF2HFsiTile(bi,bj) +
2708         &      HeatCapacity_Cp * theta(I,J,kSurface,bi,bj)
2709         &      * rA(i,j,bi,bj) * maskInC(i,j,bi,bj)
2710             ENDDO
2711            ENDDO
2712           ENDIF
2713           HFsiTile(bi,bj) = 0. _d 0
2714           IF ( balanceQnet ) THEN
2715            DO j=1,sNy
2716             DO i=1,sNx
2717              HFsiTile(bi,bj) =
2718         &      HFsiTile(bi,bj) + SItflux(i,j,bi,bj)
2719         &      * rA(i,j,bi,bj) * maskInC(i,j,bi,bj)
2720             ENDDO
2721            ENDDO
2722           ENDIF
2723    #endif
2724    
2725  C ===================================================================  C ===================================================================
2726  C ======================PART 8: diagnostics==========================  C ======================PART 8: diagnostics==========================
2727  C ===================================================================  C ===================================================================
# Line 2412  C three that actually need intermediate Line 2772  C three that actually need intermediate
2772  #ifdef ALLOW_ATM_TEMP  #ifdef ALLOW_ATM_TEMP
2773           DO J=1,sNy           DO J=1,sNy
2774            DO I=1,sNx            DO I=1,sNx
 CML If I consider the atmosphere above the ice, the surface flux  
 CML which is relevant for the air temperature dT/dt Eq  
 CML accounts for sensible and radiation (with different treatment  
 CML according to wave-length) fluxes but not for "latent heat flux",  
 CML since it does not contribute to heating the air.  
 CML So this diagnostic is only good for heat budget calculations within  
 CML the ice-ocean system.  
            DIAGarrayA(I,J) = maskC(I,J,kSurface,bi,bj)*convertHI2Q*(  
 #ifndef SEAICE_GROWTH_LEGACY  
      &            a_QSWbyATM_cover(I,J) +  
 #endif /* SEAICE_GROWTH_LEGACY */  
      &            a_QbyATM_cover(I,J) + a_QbyATM_open(I,J) )  
 C  
2775             DIAGarrayB(I,J) = maskC(I,J,kSurface,bi,bj) *             DIAGarrayB(I,J) = maskC(I,J,kSurface,bi,bj) *
2776       &       a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm       &       a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm
 C  
            DIAGarrayC(I,J) = maskC(I,J,kSurface,bi,bj)*(  
      &          PRECIP(I,J,bi,bj)  
      &          - EVAP(I,J,bi,bj)*( ONE - AREApreTH(I,J) )  
 #ifdef ALLOW_RUNOFF  
      &          + RUNOFF(I,J,bi,bj)  
 #endif /* ALLOW_RUNOFF */  
      &           )*rhoConstFresh  
      &     - a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm  
2777            ENDDO            ENDDO
2778           ENDDO           ENDDO
          CALL DIAGNOSTICS_FILL(DIAGarrayA,  
      &      'SIatmQnt',0,1,3,bi,bj,myThid)  
2779           CALL DIAGNOSTICS_FILL(DIAGarrayB,           CALL DIAGNOSTICS_FILL(DIAGarrayB,
2780       &      'SIfwSubl',0,1,3,bi,bj,myThid)       &      'SIfwSubl',0,1,3,bi,bj,myThid)
          CALL DIAGNOSTICS_FILL(DIAGarrayC,  
      &      'SIatmFW ',0,1,3,bi,bj,myThid)  
2781  C  C
2782           DO J=1,sNy           DO J=1,sNy
2783            DO I=1,sNx            DO I=1,sNx
# Line 2451  C the actual Freshwater flux of sublimat Line 2785  C the actual Freshwater flux of sublimat
2785             DIAGarrayA(I,J) = maskC(I,J,kSurface,bi,bj)             DIAGarrayA(I,J) = maskC(I,J,kSurface,bi,bj)
2786       &       * (a_FWbySublim(I,J)-r_FWbySublim(I,J))       &       * (a_FWbySublim(I,J)-r_FWbySublim(I,J))
2787       &       * SEAICE_rhoIce * recip_deltaTtherm       &       * SEAICE_rhoIce * recip_deltaTtherm
2788  c the residual Freshwater flux of sublimated ice  C the residual Freshwater flux of sublimated ice
2789             DIAGarrayC(I,J) = maskC(I,J,kSurface,bi,bj)             DIAGarrayC(I,J) = maskC(I,J,kSurface,bi,bj)
2790       &       * r_FWbySublim(I,J)       &       * r_FWbySublim(I,J)
2791       &       * SEAICE_rhoIce * recip_deltaTtherm       &       * SEAICE_rhoIce * recip_deltaTtherm
# Line 2468  C the latent heat flux Line 2802  C the latent heat flux
2802           CALL DIAGNOSTICS_FILL(DIAGarrayA,'SIacSubl',0,1,3,bi,bj,myThid)           CALL DIAGNOSTICS_FILL(DIAGarrayA,'SIacSubl',0,1,3,bi,bj,myThid)
2803           CALL DIAGNOSTICS_FILL(DIAGarrayC,'SIrsSubl',0,1,3,bi,bj,myThid)           CALL DIAGNOSTICS_FILL(DIAGarrayC,'SIrsSubl',0,1,3,bi,bj,myThid)
2804           CALL DIAGNOSTICS_FILL(DIAGarrayB,'SIhl    ',0,1,3,bi,bj,myThid)           CALL DIAGNOSTICS_FILL(DIAGarrayB,'SIhl    ',0,1,3,bi,bj,myThid)
   
          DO J=1,sNy  
           DO I=1,sNx  
 c compute ice/snow water going to atm, in precip units  
            tmpscal1 = rhoConstFresh*maskC(I,J,kSurface,bi,bj)  
      &       * convertHI2PRECIP * ( - d_HSNWbyRAIN(I,J)*SNOW2ICE  
      &       + a_FWbySublim(I,J) - r_FWbySublim(I,J) )  
 c compute ocean water going to atm, in precip units  
            tmpscal2=rhoConstFresh*maskC(I,J,kSurface,bi,bj)*  
      &       ( ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )  
      &         * ( ONE - AREApreTH(I,J) )  
 #ifdef ALLOW_RUNOFF  
      &         - RUNOFF(I,J,bi,bj)  
 #endif /* ALLOW_RUNOFF */  
      &         + ( d_HFRWbyRAIN(I,J) + r_FWbySublim(I,J) )  
      &         *convertHI2PRECIP )  
 c factor in the advected specific energy (referenced to 0 for 0deC liquid water)  
            tmpscal1= - tmpscal1*  
      &       ( -SEAICE_lhFusion + HeatCapacity_Cp * ZERO )  
            IF (temp_EvPrRn.NE.UNSET_RL) THEN  
            tmpscal2= - tmpscal2*  
      &       ( ZERO + HeatCapacity_Cp * temp_EvPrRn )  
            ELSE  
            tmpscal2= - tmpscal2*  
      &       ( ZERO + HeatCapacity_Cp * theta(I,J,kSurface,bi,bj) )  
            ENDIF  
 c add to SIatmQnt, leading to SItflux, which is analogous to TFLUX  
            DIAGarrayA(I,J)=maskC(I,J,kSurface,bi,bj)*convertHI2Q*(  
 #ifndef SEAICE_GROWTH_LEGACY  
      &            a_QSWbyATM_cover(I,J) +  
 #endif  
      &            a_QbyATM_cover(I,J) + a_QbyATM_open(I,J) )  
      &            -tmpscal1-tmpscal2  
           ENDDO  
          ENDDO  
          CALL DIAGNOSTICS_FILL(DIAGarrayA,  
      &      'SItflux ',0,1,3,bi,bj,myThid)  
2805  #endif /* ALLOW_ATM_TEMP */  #endif /* ALLOW_ATM_TEMP */
2806    
2807          ENDIF          ENDIF
# Line 2514  C close bi,bj loops Line 2811  C close bi,bj loops
2811         ENDDO         ENDDO
2812        ENDDO        ENDDO
2813    
2814    
2815    C ===================================================================
2816    C =========PART 9: HF/FWF global integrals and balancing=============
2817    C ===================================================================
2818    
2819    #ifdef ALLOW_BALANCE_FLUXES
2820    
2821    c 1)  global sums
2822    # ifdef ALLOW_AUTODIFF_TAMC
2823    CADJ STORE FWFsiTile = comlev1, key=ikey_dynamics, kind=isbyte
2824    CADJ STORE HFsiTile = comlev1, key=ikey_dynamics, kind=isbyte
2825    CADJ STORE FWF2HFsiTile = comlev1, key=ikey_dynamics, kind=isbyte
2826    # endif /* ALLOW_AUTODIFF_TAMC */
2827          FWFsiGlob=0. _d 0
2828          IF ( balanceEmPmR )
2829         &   CALL GLOBAL_SUM_TILE_RL( FWFsiTile, FWFsiGlob, myThid )
2830          FWF2HFsiGlob=0. _d 0
2831          IF ( balanceEmPmR.AND.(temp_EvPrRn.EQ.UNSET_RL) ) THEN
2832             CALL GLOBAL_SUM_TILE_RL(FWF2HFsiTile, FWF2HFsiGlob, myThid)
2833          ELSEIF ( balanceEmPmR ) THEN
2834             FWF2HFsiGlob=HeatCapacity_Cp * temp_EvPrRn * globalArea
2835          ENDIF
2836          HFsiGlob=0. _d 0
2837          IF ( balanceQnet )
2838         &   CALL GLOBAL_SUM_TILE_RL( HFsiTile, HFsiGlob, myThid )
2839    
2840    c 2) global means
2841    c mean SIatmFW
2842          tmpscal0=FWFsiGlob / globalArea
2843    c corresponding mean advection by atm to ocean+ice water exchange
2844    c        (if mean SIatmFW was removed uniformely from ocean)
2845          tmpscal1=FWFsiGlob / globalArea * FWF2HFsiGlob / globalArea
2846    c mean SItflux (before potential adjustement due to SIatmFW)
2847          tmpscal2=HFsiGlob / globalArea
2848    c mean SItflux (after potential adjustement due to SIatmFW)
2849          IF ( balanceEmPmR ) tmpscal2=tmpscal2-tmpscal1
2850    
2851    c 3) balancing adjustments
2852          IF ( balanceEmPmR ) THEN
2853          DO bj=myByLo(myThid),myByHi(myThid)
2854           DO bi=myBxLo(myThid),myBxHi(myThid)
2855            DO j=1-OLy,sNy+OLy
2856             DO i=1-OLx,sNx+OLx
2857                empmr(i,j,bi,bj) = empmr(i,j,bi,bj) - tmpscal0
2858                SIatmFW(i,j,bi,bj) = SIatmFW(i,j,bi,bj) - tmpscal0
2859    c           adjust SItflux consistently
2860                IF ( (temp_EvPrRn.NE.UNSET_RL).AND.
2861         &        useRealFreshWaterFlux.AND.(nonlinFreeSurf.NE.0) ) THEN
2862                tmpscal1=
2863         &       ( ZERO + HeatCapacity_Cp * temp_EvPrRn )
2864                ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL).AND.
2865         &        useRealFreshWaterFlux.AND.(nonlinFreeSurf.NE.0) ) THEN
2866                tmpscal1=
2867         &       ( ZERO + HeatCapacity_Cp * theta(I,J,kSurface,bi,bj) )
2868                ELSEIF ( (temp_EvPrRn.NE.UNSET_RL) ) THEN
2869                tmpscal1=
2870         &       HeatCapacity_Cp*(temp_EvPrRn - theta(I,J,kSurface,bi,bj))
2871                ELSE
2872                tmpscal1=ZERO
2873                ENDIF
2874                SItflux(i,j,bi,bj) = SItflux(i,j,bi,bj) - tmpscal0*tmpscal1
2875    c           no qnet or tflux adjustement is needed
2876             ENDDO
2877            ENDDO
2878           ENDDO
2879          ENDDO
2880          IF ( balancePrintMean ) THEN
2881           _BEGIN_MASTER( myThid )
2882           WRITE(msgBuf,'(a,a,e24.17)') 'rm Global mean of ',
2883         &      'SIatmFW = ', tmpscal0
2884           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
2885         &      SQUEEZE_RIGHT , myThid)
2886           _END_MASTER( myThid )
2887          ENDIF
2888          ENDIF
2889          IF ( balanceQnet ) THEN
2890          DO bj=myByLo(myThid),myByHi(myThid)
2891           DO bi=myBxLo(myThid),myBxHi(myThid)
2892            DO j=1-OLy,sNy+OLy
2893             DO i=1-OLx,sNx+OLx
2894                SItflux(i,j,bi,bj) = SItflux(i,j,bi,bj) - tmpscal2
2895                qnet(i,j,bi,bj) = qnet(i,j,bi,bj) - tmpscal2
2896                SIatmQnt(i,j,bi,bj) = SIatmQnt(i,j,bi,bj) - tmpscal2
2897             ENDDO
2898            ENDDO
2899           ENDDO
2900          ENDDO
2901          IF ( balancePrintMean ) THEN
2902           _BEGIN_MASTER( myThid )
2903           WRITE(msgBuf,'(a,a,e24.17)') 'rm Global mean of ',
2904         &      'SItflux = ', tmpscal2
2905           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
2906         &      SQUEEZE_RIGHT , myThid)
2907           _END_MASTER( myThid )
2908          ENDIF
2909          ENDIF
2910    #endif /* ALLOW_BALANCE_FLUXES */
2911    
2912    #ifdef ALLOW_DIAGNOSTICS
2913    c these diags need to be done outside of the bi,bj loop so that
2914    c we may do potential global mean adjustement to them consistently.
2915             CALL DIAGNOSTICS_FILL(SItflux,
2916         &      'SItflux ',0,1,0,1,1,myThid)
2917             CALL DIAGNOSTICS_FILL(SIatmQnt,
2918         &      'SIatmQnt',0,1,0,1,1,myThid)
2919    c SIatmFW follows the same convention as empmr -- SIatmFW diag does not
2920             tmpscal1= - 1. _d 0
2921             CALL DIAGNOSTICS_SCALE_FILL(SIatmFW,
2922         &      tmpscal1,1,'SIatmFW ',0,1,0,1,1,myThid)
2923    #endif /* ALLOW_DIAGNOSTICS */
2924    
2925        RETURN        RETURN
2926        END        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.22