Back to home page

MITgcm

 
 

    


File indexing completed on 2025-07-08 05:11:23 UTC

view on githubraw file Latest commit 00c7090d on 2025-07-07 16:10:22 UTC
dc54d31829 Ian *0001 #include "SEAICE_OPTIONS.h"
                0002 #ifdef ALLOW_EXF
                0003 # include "EXF_OPTIONS.h"
                0004 #endif
                0005 #ifdef ALLOW_SALT_PLUME
                0006 # include "SALT_PLUME_OPTIONS.h"
                0007 #endif
                0008 #ifdef ALLOW_AUTODIFF
                0009 # include "AUTODIFF_OPTIONS.h"
                0010 #endif
                0011 
2383f7d4e9 Mart*0012 #undef SEAICE_GROWTH_ADX_STORE_MORE
                0013 
dc54d31829 Ian *0014 CBOP
                0015 C     !ROUTINE: SEAICE_GROWTH_ADX
                0016 C     !INTERFACE:
4dd39c50d9 Mart*0017       SUBROUTINE SEAICE_GROWTH_ADX( myTime, myIter, myThid )
dc54d31829 Ian *0018 C     !DESCRIPTION: \bv
                0019 C     *==========================================================*
                0020 C     | SUBROUTINE seaice_growth_adx
4dd39c50d9 Mart*0021 C     | o Adjointable seaice thermodynamic code based on
                0022 C     |   Ian Fentys code with modifications for nonlinear free
dc54d31829 Ian *0023 C     |   surface budget closure by A.Bigdeli.
                0024 C     *==========================================================*
                0025 C     \ev
                0026 
                0027 C     !USES:
                0028       IMPLICIT NONE
                0029 C     === Global variables ===
                0030 #include "SIZE.h"
                0031 #include "EEPARAMS.h"
                0032 #include "PARAMS.h"
                0033 #include "DYNVARS.h"
                0034 #include "GRID.h"
                0035 #include "FFIELDS.h"
                0036 #include "SEAICE_SIZE.h"
                0037 #include "SEAICE_PARAMS.h"
                0038 #include "SEAICE.h"
                0039 #include "SEAICE_TRACER.h"
                0040 #ifdef ALLOW_EXF
                0041 # include "EXF_PARAM.h"
                0042 # include "EXF_FIELDS.h"
                0043 #endif
                0044 #ifdef ALLOW_SALT_PLUME
                0045 # include "SALT_PLUME.h"
                0046 #endif
                0047 #ifdef ALLOW_AUTODIFF_TAMC
                0048 # include "tamc.h"
                0049 #endif
                0050 
                0051 C     !INPUT/OUTPUT PARAMETERS:
                0052 C     === Routine arguments ===
                0053 C     myTime :: Simulation time
                0054 C     myIter :: Simulation timestep number
                0055 C     myThid :: Thread no. that called this routine.
                0056       _RL myTime
                0057       INTEGER myIter, myThid
                0058 CEOP
                0059 
5337b47ce7 Jean*0060 #ifdef SEAICE_USE_GROWTH_ADX
dc54d31829 Ian *0061 #if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP)
                0062 C     !FUNCTIONS:
                0063 #ifdef ALLOW_DIAGNOSTICS
                0064       LOGICAL  DIAGNOSTICS_IS_ON
                0065       EXTERNAL DIAGNOSTICS_IS_ON
                0066 #endif
                0067 
                0068 C     !LOCAL VARIABLES:
                0069 C     === Local variables ===
                0070 C
                0071 C unit/sign convention:
                0072 
                0073 C HEFF is effective Hice thickness (m3/m2)
                0074 C HSNOW is Heffective snow thickness (m3/m2)
                0075 C HSALT is Heffective salt content (g/m2)
                0076 C AREA is the seaice cover fraction (0<=AREA<=1)
                0077 
                0078 C     i,j,bi,bj :: Loop counters
                0079       INTEGER i, j, bi, bj
00c7090dc0 Mart*0080 C     number of surface interface layer, and second interface layer
                0081       INTEGER kSurface, kSrfS
dc54d31829 Ian *0082       INTEGER kSurface
                0083 C     IT :: ice thickness category index (MULTICATEGORIES and ITD code)
                0084       INTEGER IT
                0085 C     msgBuf      :: Informational/error message buffer
914da317ef jm-c 0086 c     CHARACTER*(MAX_LEN_MBUF) msgBuf
dc54d31829 Ian *0087 C     constants
                0088       _RL tempFrz, ICE2SNOW, SNOW2ICE, surf_theta
914da317ef jm-c 0089       _RL QI, QS
dc54d31829 Ian *0090       _RL lhSublim
                0091 
                0092 C     Regularization values squared
                0093       _RL area_reg_sq, hice_reg_sq
                0094 
                0095 C     pathological cases thresholds
                0096       _RL heffTooHeavy
                0097 
                0098 C     Helper variables: reciprocal of some constants
                0099       _RL recip_multDim
                0100       _RL recip_deltaTtherm
                0101       _RL recip_rhoIce
4dd39c50d9 Mart*0102       _RL recip_QI
                0103       _RL recip_QS
                0104       _RL recip_HO
                0105       _RL recip_HO_south
                0106       _RL rhoFresh2rhoSnow
                0107       _RL rhoSnow2rhoFresh
                0108       _RL rhoIce2rhoFresh
dc54d31829 Ian *0109 
                0110 C     facilitate multi-category snow implementation
                0111       _RL pFac, pFacSnow
                0112 
                0113 C     temporary variables available for the various computations
914da317ef jm-c 0114       _RL tmpscal0, tmpscal1, tmpscal2, tmpscal3
dc54d31829 Ian *0115 
7c50f07931 Mart*0116 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0117 C     tkey :: tape key (depends on tiles)
                0118       INTEGER tkey
7c50f07931 Mart*0119 #endif
                0120 
dc54d31829 Ian *0121 C==   local arrays ==
                0122 C--   TmixLoc        :: ocean surface/mixed-layer temperature (in K)
                0123       _RL TmixLoc       (1:sNx,1:sNy)
                0124 
                0125 C     actual ice thickness (with upper and lower limit)
                0126       _RL hiceActual          (1:sNx,1:sNy)
                0127 C     actual snow thickness
                0128       _RL hsnowActual         (1:sNx,1:sNy)
                0129 C     actual ice thickness (with lower limit only) Reciprocal
                0130       _RL recip_hiceActual    (1:sNx,1:sNy)
                0131 
                0132 C     AREA_PRE :: hold sea-ice fraction field before any seaice-thermo update
                0133       _RL AREApreTH           (1:sNx,1:sNy)
                0134       _RL HEFFpreTH           (1:sNx,1:sNy)
                0135       _RL HSNWpreTH           (1:sNx,1:sNy)
                0136 
                0137 C     wind speed
                0138       _RL UG                  (1:sNx,1:sNy)
                0139 
                0140 C     temporary variables available for the various computations
                0141       _RL ticeInMult          (1:sNx,1:sNy,nITD)
                0142       _RL ticeOutMult         (1:sNx,1:sNy,nITD)
                0143       _RL hiceActualMult      (1:sNx,1:sNy,nITD)
                0144       _RL hsnowActualMult     (1:sNx,1:sNy,nITD)
                0145 
                0146       _RL F_io_net_mult     (1:sNx,1:sNy,nITD)
                0147       _RL F_ia_net_mult     (1:sNx,1:sNy,nITD)
                0148       _RL F_ia_mult     (1:sNx,1:sNy,nITD)
                0149       _RL QSWI_mult     (1:sNx,1:sNy,nITD)
a4e168e012 antn*0150 
                0151 C Variables related to FWsublim are dummies required for call
                0152 C to seaice_solve4temp subroutine and are not used afterwards
                0153 C in this routine (by choice).
dc54d31829 Ian *0154       _RL FWsublim_mult     (1:sNx,1:sNy,nITD)
a4e168e012 antn*0155       _RL FWsublim          (1:sNx,1:sNy)
dc54d31829 Ian *0156 
                0157 #ifdef ALLOW_SALT_PLUME
                0158 c     d(HEFF)/dt from heat fluxes in the open water fraction of the grid cell
                0159       _RL IceGrowthRateInLeads          (1:sNx,1:sNy)
                0160 
                0161 c     The fraction of salt released in leads by new ice production there
                0162 c     which is to be sent to the salt plume package
                0163       _RL leadPlumeFraction             (1:sNx,1:sNy)
                0164 #endif
                0165 
                0166       _RL QSWO_BELOW_FIRST_LAYER (1:sNx,1:sNy)
                0167       _RL QSWO_IN_FIRST_LAYER (1:sNx,1:sNy)
                0168       _RL QSWO (1:sNx,1:sNy)
                0169       _RL QSWI (1:sNx,1:sNy)
                0170 
                0171 C     Sea ice growth rates (m/s)
                0172       _RL ActualNewTotalVolumeChange     (1:sNx,1:sNy)
                0173       _RL ActualNewTotalSnowMelt     (1:sNx,1:sNy)
                0174 
                0175       _RL NetExistingIceGrowthRate      (1:sNx,1:sNy)
                0176       _RL IceGrowthRateUnderExistingIce (1:sNx,1:sNy)
                0177       _RL IceGrowthRateFromSurface      (1:sNx,1:sNy)
                0178       _RL IceGrowthRateOpenWater        (1:sNx,1:sNy)
                0179       _RL IceGrowthRateMixedLayer       (1:sNx,1:sNy)
                0180 
                0181       _RL EnergyInNewTotalIceVolume     (1:sNx,1:sNy)
                0182       _RL NetEnergyFluxOutOfOcean       (1:sNx,1:sNy)
                0183 
d73d15cdef antn*0184 #ifdef ALLOW_DIAGNOSTICS
a4e168e012 antn*0185 C Local arrays for diagnostics
4dd39c50d9 Mart*0186       _RL SItflux     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
a4e168e012 antn*0187       _RL SIeprflx    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
4dd39c50d9 Mart*0188       _RL SIatmFW     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
d73d15cdef antn*0189 #endif
dc54d31829 Ian *0190 CAB these two terms are placeholders for extra bits from advection,
                0191 C   used in budget
4dd39c50d9 Mart*0192       _RL SIheffNeg   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0193       _RL SIhsnwNeg   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
dc54d31829 Ian *0194 CAB
d73d15cdef antn*0195       _RL d_HEFFbyFLOODING              (1:sNx,1:sNy)
dc54d31829 Ian *0196 
                0197 C     The energy taken out of the ocean which is not converted
                0198 C     to sea ice (Joules)
                0199       _RL ResidualEnergyOutOfOcean      (1:sNx,1:sNy)
                0200 
                0201 C     Snow accumulation rate over ice (m/s)
                0202       _RL SnowAccRateOverIce            (1:sNx,1:sNy)
                0203 
                0204 C     Total snow accumulation over ice (m)
                0205       _RL SnowAccOverIce                (1:sNx,1:sNy)
                0206 
                0207 C     The precipitation rate over the ice which goes immediately into the ocean
                0208       _RL PrecipRateOverIceSurfaceToSea (1:sNx,1:sNy)
                0209 
                0210 C     The potential snow melt rate if all snow surface heat flux convergences
                0211 C     goes to melting snow (m/s)
                0212       _RL PotSnowMeltRateFromSurf       (1:sNx,1:sNy)
                0213 
                0214 C     The potential thickness of snow which could be melted by snow surface
                0215 C     heat flux convergence (m)
                0216       _RL PotSnowMeltFromSurf           (1:sNx,1:sNy)
                0217 
                0218 C     The actual snow melt rate due to snow surface  heat flux convergence
                0219       _RL SnowMeltRateFromSurface       (1:sNx,1:sNy)
                0220 
                0221 C     The actual surface heat flux convergence used to melt snow (W/m^2)
                0222       _RL SurfHeatFluxConvergToSnowMelt (1:sNx,1:sNy)
                0223 
                0224 C     The actual thickness of snow to be melted by snow surface
                0225 C     heat flux convergence (m)
                0226       _RL SnowMeltFromSurface           (1:sNx,1:sNy)
                0227 
                0228 C     The freshwater contribution to the ocean from melting snow (m)
                0229       _RL FreshwaterContribFromSnowMelt (1:sNx,1:sNy)
                0230 
                0231 C     The freshwater contribution to (from) the ocean from melting (growing) ice (m)
                0232       _RL FreshwaterContribFromIce      (1:sNx,1:sNy)
                0233 
                0234 C     S_a : d(AREA)/dt
                0235       _RL S_a                           (1:sNx,1:sNy)
                0236 
                0237 C     S_h : d(HEFF)/dt
                0238       _RL S_h                           (1:sNx,1:sNy)
                0239 
                0240 C     S_hsnow : d(HSNOW)/dt
                0241       _RL S_hsnow                       (1:sNx,1:sNy)
                0242 
                0243 C     F_ia  - sea ice/snow surface heat flux with atmosphere (W/m^2)
                0244 C       F_ia > 0, heat loss to atmosphere
                0245 C       F_ia < 0, atmospheric heat flux convergence (ice/snow surface melt)
                0246       _RL F_ia                          (1:sNx,1:sNy)
                0247 
                0248 C     F_ia_net - the net heat flux divergence at the sea ice/snow surface
                0249 C                including sea ice conductive fluxes and atmospheric fluxes (W/m^2)
                0250 C       F_ia_net = 0, sea ice/snow surface energy balance condition met
                0251 C                     upward conductive fluxes balance surface heat loss
                0252 C       F_ia_net < 0, net heat flux convergence at ice/snow surface
                0253 C                     zero conductive fluxes and net atmospheric convergence
                0254       _RL F_ia_net                      (1:sNx,1:sNy)
                0255 
                0256 C     F_ia_net - the net heat flux divergence at the sea ice/snow surface
                0257 C                before snow is melted with any convergence (W/m^2)
                0258 C        F_ia_net < 0, some snow, if present, will melt
                0259       _RL F_ia_net_before_snow          (1:sNx,1:sNy)
                0260 
                0261 C     F_io_net - the net upward conductive heat flux through the ice+snow system
                0262 C                realized at the sea ice/snow surface
                0263 C        F_io_net > 0, heat conducting upward from ice base --> basal thickening
                0264 C        F_io_net = 0, no upward heat conduction
                0265 C                      ice/snow surface temperature > SEAICE_freeze)
                0266       _RL F_io_net                      (1:sNx,1:sNy)
                0267 
                0268 C     F_ao  - heat flux from atmosphere to ocean (W/m^2)
                0269 C        F_ao > 0
                0270       _RL F_ao                          (1:sNx,1:sNy)
                0271 
                0272 C     F_mi - heat flux from ocean to the ice (W/m^2)
                0273       _RL F_mi                          (1:sNx,1:sNy)
                0274 
                0275 C     S_a_from_IGROW : d(AREA)/dt [from ice growth rate from open water fluxes]
                0276       _RL S_a_IGROW                     (1:sNx,1:sNy)
                0277 
                0278 C     S_a_from_IGROW : d(AREA)/dt [from ice growth rate from ocean-ice fluxes]
                0279       _RL S_a_IGRML                     (1:sNx,1:sNy)
                0280 
                0281 C     S_a_from_IGROW : d(AREA)/dt [from ice growth rate from ice-atm fluxes]
                0282       _RL S_a_IGRNE                     (1:sNx,1:sNy)
                0283 
                0284 C     Factor by which we increase the upper ocean friction velocity (u*) when
                0285 C     ice is absent in a grid cell  (dimensionless)
4dd39c50d9 Mart*0286       _RL MLTF
dc54d31829 Ian *0287 
                0288 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0289 
                0290 C ===================================================================
                0291 C =================PART 0: constants and initializations=============
                0292 C ===================================================================
                0293 
0320e25227 Mart*0294       IF ( usingPCoords ) THEN
00c7090dc0 Mart*0295        kSurface = Nr
                0296        kSrfS    = Nr
dc54d31829 Ian *0297       ELSE
00c7090dc0 Mart*0298        kSurface = 1
                0299        kSrfS    = 2
dc54d31829 Ian *0300       ENDIF
                0301 
                0302 C     Cutoff for iceload
                0303       heffTooHeavy=drF(kSurface) / 5. _d 0
0320e25227 Mart*0304       IF ( usingPCoords )
                0305      &     heffTooHeavy = heffTooHeavy * recip_rhoConst * recip_gravity
dc54d31829 Ian *0306 
                0307 C     RATIO OF SEA ICE DENSITY to SNOW DENSITY
                0308       ICE2SNOW     = SEAICE_rhoIce/SEAICE_rhoSnow
                0309       SNOW2ICE     = ONE / ICE2SNOW
                0310 
                0311 C  Heat Flux * QI = [W] [m^3/kg] [kg/J] = [m^3/s] or [m/s] per unit area
                0312       QI           = ONE/(SEAICE_rhoIce*SEAICE_lhFusion)
                0313 C  Heat Flux * QS = [W] [m^3/kg] [kg/J] = [m^3/s] or [m/s] per unit area
                0314       QS           = ONE/(SEAICE_rhoSnow*SEAICE_lhFusion)
                0315 
                0316 C     ICE LATENT HEAT CONSTANT
                0317       lhSublim = SEAICE_lhEvap + SEAICE_lhFusion
                0318 
4dd39c50d9 Mart*0319 C     avoid unnecessary divisions in loops
                0320       recip_multDim     = SEAICE_multDim
                0321       recip_multDim     = ONE / recip_multDim
                0322 C     above/below: double/single precision calculation of recip_multDim
                0323 c     recip_multDim     = 1./float(SEAICE_multDim)
                0324       recip_deltaTtherm = ONE / SEAICE_deltaTtherm
                0325       recip_rhoIce      = ONE / SEAICE_rhoIce
                0326       recip_QI          = ONE / QI
                0327       recip_QS          = ONE / QS
                0328       recip_HO          = ONE / HO
                0329       recip_HO_south    = ONE / HO_south
                0330       rhoFresh2rhoSnow  = rhoConstFresh / SEAICE_rhoSnow
                0331       rhoSnow2rhoFresh  = ONE / rhoFresh2rhoSnow
                0332       rhoIce2rhoFresh   = SEAICE_rhoIce / rhoConstFresh
                0333 
dc54d31829 Ian *0334 C     regularization constants
                0335       area_reg_sq = SEAICE_area_reg * SEAICE_area_reg
                0336       hice_reg_sq = SEAICE_hice_reg * SEAICE_hice_reg
                0337 
                0338 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0339       DO bj=myByLo(myThid),myByHi(myThid)
                0340        DO bi=myBxLo(myThid),myBxHi(myThid)
                0341 
                0342 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0343         tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
4dd39c50d9 Mart*0344 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
edb6656069 Mart*0345 CADJ STORE area (:,:,bi,bj) = comlev1_bibj, key = tkey, byte = isbyte
                0346 CADJ STORE heff (:,:,bi,bj) = comlev1_bibj, key = tkey, byte = isbyte
                0347 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = tkey, byte = isbyte
                0348 CADJ STORE qnet (:,:,bi,bj) = comlev1_bibj, key = tkey, byte = isbyte
                0349 CADJ STORE qsw  (:,:,bi,bj) = comlev1_bibj, key = tkey, byte = isbyte
                0350 CADJ STORE tices(:,:,:,bi,bj)=comlev1_bibj, key = tkey, byte = isbyte
2383f7d4e9 Mart*0351 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0352 #endif
dc54d31829 Ian *0353 
                0354 C array initializations
                0355 C =====================
                0356 
                0357         DO J=1,sNy
                0358          DO I=1,sNx
                0359 
                0360 C NEW VARIABLE NAMES
                0361           NetExistingIceGrowthRate(I,J) = 0.0 _d 0
                0362           IceGrowthRateOpenWater(I,J)   = 0.0 _d 0
                0363           IceGrowthRateFromSurface(I,J) = 0.0 _d 0
                0364           IceGrowthRateMixedLayer(I,J)  = 0.0 _d 0
                0365 
                0366           EnergyInNewTotalIceVolume(I,J) = 0.0 _d 0
                0367           NetEnergyFluxOutOfOcean(I,J)   = 0.0 _d 0
                0368           ResidualEnergyOutOfOcean(I,J)  = 0.0 _d 0
                0369 
                0370           PrecipRateOverIceSurfaceToSea(I,J) = 0.0 _d 0
                0371 
                0372           DO IT=1,SEAICE_multDim
                0373             ticeInMult(I,J,IT)            = 0.0 _d 0
                0374             ticeOutMult(I,J,IT)           = 0.0 _d 0
                0375 
                0376             F_io_net_mult(I,J,IT)  = 0.0 _d 0
                0377             F_ia_net_mult(I,J,IT)  = 0.0 _d 0
                0378             F_ia_mult(I,J,IT)      = 0.0 _d 0
                0379 
                0380             QSWI_mult(I,J,IT)      = 0.0 _d 0
                0381             FWsublim_mult(I,J,IT)  = 0.0 _d 0
                0382           ENDDO
                0383 
                0384           F_io_net(I,J) = 0.0 _d 0
                0385           F_ia_net(I,J) = 0.0 _d 0
                0386           F_ia(I,J)     = 0.0 _d 0
                0387 
                0388           QSWI(I,J) = 0.0 _d 0
                0389           FWsublim(I,J) = 0.0 _d 0
                0390 
                0391           QSWO_BELOW_FIRST_LAYER (I,J) = 0.0 _d 0
                0392           QSWO_IN_FIRST_LAYER (I,J) = 0.0 _d 0
                0393           QSWO(I,J) = 0.0 _d 0
                0394 
                0395           ActualNewTotalVolumeChange(I,J) = 0.0 _d 0
                0396           ActualNewTotalSnowMelt(I,J)     = 0.0 _d 0
                0397 
                0398           SnowAccOverIce(I,J) = 0.0 _d 0
                0399           SnowAccRateOverIce(I,J) = 0.0 _d 0
                0400 
                0401           PotSnowMeltRateFromSurf(I,J)    = 0.0 _d 0
                0402           PotSnowMeltFromSurf(I,J)        = 0.0 _d 0
                0403           SnowMeltRateFromSurface(I,J)    = 0.0 _d 0
                0404           SurfHeatFluxConvergToSnowMelt(I,J) = 0.0 _d 0
                0405           SnowMeltFromSurface(I,J)       = 0.0 _d 0
                0406 
                0407           FreshwaterContribFromSnowMelt(I,J) = 0.0 _d 0
                0408           FreshwaterContribFromIce(I,J) = 0.0 _d 0
                0409 
                0410           S_a(I,J)       = 0.0 _d 0
                0411           S_a_IGROW(I,J) = 0.0 _d 0
                0412           S_a_IGRML(I,J) = 0.0 _d 0
                0413           S_a_IGRNE(I,J) = 0.0 _d 0
                0414 
                0415           S_h(I,J)       = 0.0 _d 0
                0416           S_hsnow(I,J)   = 0.0 _d 0
                0417 
d73d15cdef antn*0418           d_HEFFbyFLOODING(I,J) = 0.0 _d 0
dc54d31829 Ian *0419 #ifdef ALLOW_SALT_PLUME
                0420           IceGrowthRateInLeads            (I,J) = 0. _d 0
                0421           leadPlumeFraction               (I,J) = 0. _d 0
                0422           saltPlumeFlux             (I,J,bi,bj) = 0. _d 0
                0423 #endif
af61e5eb16 Mart*0424 #ifdef ALLOW_DIAGNOSTICS
a4e168e012 antn*0425           SIeprflx(I,J,bi,bj) = 0. _d 0
af61e5eb16 Mart*0426 #endif
dc54d31829 Ian *0427          ENDDO
                0428         ENDDO
                0429 
                0430 C =====================================================================
                0431 C  PART 1: Store ice and snow state on onset + regularize actual
                0432 C               snow and ice thickness
                0433 C =====================================================================
                0434 CAB
                0435         DO J=1,sNy
                0436          DO I=1,sNx
                0437          SIheffNeg(I,J,bi,bj)=d_HEFFbyNEG(I,J,bi,bj)*SINegFac
                0438          SIhsnwNeg(I,J,bi,bj)=d_HSNWbyNEG(I,J,bi,bj)*SINegFac
                0439          ENDDO
                0440         ENDDO
                0441 
                0442         DO J=1,sNy
                0443          DO I=1,sNx
                0444 
4dd39c50d9 Mart*0445           tmpscal0 = HEFF(I,J,bi,bj)
                0446           HEFF(I,J,bi,bj) = MAX(ZERO, tmpscal0)
                0447           tmpscal1 = AREA(I,J,bi,bj)
                0448           AREA(I,J,bi,bj) = MAX(ZERO, tmpscal1)
dc54d31829 Ian *0449 
                0450           IF (HEFF(I,J,bi,bj) .LE. ZERO) then
                0451              AREA(I,J, bi,bj)  = ZERO
                0452              HSNOW(I,J, bi,bj) = ZERO
                0453           ELSEIF (AREA(I,J,bi,bj) .LE. ZERO) then
                0454              HEFF(I,J,bi,bj)  = ZERO
                0455              HSNOW(I,J,bi,bj) = ZERO
                0456           ENDIF
                0457 
                0458           HEFFpreTH(I,J) = HEFF(I,J,bi,bj)
                0459           HSNWpreTH(I,J) = HSNOW(I,J,bi,bj)
                0460           AREApreTH(I,J) = AREA(I,J,bi,bj)
                0461 
                0462          ENDDO
                0463         ENDDO
                0464 
                0465 #ifdef ALLOW_DIAGNOSTICS
                0466         IF ( useDiagnostics ) THEN
d73d15cdef antn*0467          CALL DIAGNOSTICS_FILL(AREA, 'SIareaPT',0,1,1,bi,bj,myThid)
                0468          CALL DIAGNOSTICS_FILL(HEFF, 'SIheffPT',0,1,1,bi,bj,myThid)
                0469          CALL DIAGNOSTICS_FILL(HSNOW,'SIhsnoPT',0,1,1,bi,bj,myThid)
dc54d31829 Ian *0470         ENDIF
                0471 #endif /* ALLOW_DIAGNOSTICS */
                0472 
2383f7d4e9 Mart*0473 #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_GROWTH_ADX_STORE_MORE)
dc54d31829 Ian *0474 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
edb6656069 Mart*0475 CADJ STORE AREApreTH = comlev1_bibj, key = tkey, byte = isbyte
                0476 CADJ STORE HEFFpreTH = comlev1_bibj, key = tkey, byte = isbyte
                0477 CADJ STORE HSNWpreTH = comlev1_bibj, key = tkey, byte = isbyte
dc54d31829 Ian *0478 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
2383f7d4e9 Mart*0479 #endif /* ALLOW_AUTODIFF_TAMC */
dc54d31829 Ian *0480 
                0481 C     COMPUTE ACTUAL ICE/SNOW THICKNESS; USE MIN/MAX VALUES
                0482 C     TO REGULARIZE SEAICE_SOLVE4TEMP/d_AREA COMPUTATIONS
                0483 
                0484         DO J=1,sNy
                0485          DO I=1,sNx
                0486           IF (HEFFpreTH(I,J) .GT. ZERO) THEN
                0487 c          regularize AREA with SEAICE_area_reg
                0488            tmpscal1 = SQRT(AREApreTH(I,J)* AREApreTH(I,J) + area_reg_sq)
                0489 
                0490 c          hiceActual calculated with the regularized AREA
                0491            tmpscal2 = HEFFpreTH(I,J) / tmpscal1
                0492 
                0493 c          regularize hiceActual with SEAICE_hice_reg (add lower bound)
                0494 c           hiceActual(I,J) = SQRT(tmpscal2 * tmpscal2 + hice_reg_sq)
                0495            hiceActual(I,J) = MAX(0.05 _d 0, tmpscal2)
                0496 
                0497 c          hsnowActual calculated with the regularized AREA (no lower bound)
                0498 c          hsnowActual(I,J) = HSNWpreTH(I,J) / tmpscal1
4dd39c50d9 Mart*0499 c          actually I do not think we need to regularize this.
dc54d31829 Ian *0500            hsnowActual(I,J) = HSNWpreTH(I,J) / AREApreTH(I,J)
                0501 
                0502 c          regularize the inverse of hiceActual by hice_reg
                0503            recip_hiceActual(I,J)  = AREApreTH(I,J) /
                0504      &         sqrt(HEFFpreTH(I,J)*HEFFpreTH(I,J) + hice_reg_sq)
                0505 
                0506 c         Do not regularize when HEFFpreTH = 0
                0507           ELSE
                0508            hiceActual (I,J)       = ZERO
                0509            hsnowActual(I,J)       = ZERO
                0510            recip_hiceActual(I,J)  = ZERO
                0511           ENDIF
                0512 
                0513          ENDDO
                0514         ENDDO
                0515 
                0516 C =============================================================================
                0517 C Part 2: Precipitation as snow or rain over ice
                0518 C =============================================================================
                0519 
                0520         DO J=1,sNy
                0521          DO I=1,sNx
                0522 c        if we have ice and the temperature of the ice is below the freezing point
                0523 c        then the precip falls and accumulates as snow
                0524           IF (( AREApreTH(I,J)     .GT. ZERO) .AND.
                0525      &        ( TICES(I,J,1,bi,bj) .LT. celsius2k) ) THEN
                0526 
                0527 c use either prescribed snowfall or PRECIP rate
                0528            IF ( snowPrecipFile .NE. ' ' ) THEN
                0529 c    rate of snow accumulation in m/s over ice
                0530 c                      y [m/s] \approx 1.0 [kg/m^3] / 0.9 [m^3/kg] * x [m/s]
4dd39c50d9 Mart*0531              SnowAccRateOverIce(I,J) = rhoFresh2rhoSnow *
dc54d31829 Ian *0532      &               snowPrecip(i,j,bi,bj)
                0533 
                0534            ELSE
4dd39c50d9 Mart*0535              SnowAccRateOverIce(I,J) = rhoFresh2rhoSnow *
dc54d31829 Ian *0536      &           PRECIP(i,j,bi,bj)
                0537 
                0538            ENDIF
                0539 
                0540            PrecipRateOverIceSurfaceToSea(I,J) = ZERO
                0541 
                0542           ELSE
                0543 c            The snow/ice surface is not frozen (wet) so the precipitation
                0544 c            remains wet and runs into the ocean
                0545              SnowAccRateOverIce(I,J) = ZERO
                0546              PrecipRateOverIceSurfaceToSea(I,J) = PRECIP(i,j,bi,bj)
                0547           ENDIF
                0548 
                0549 c  actual change in snow thickness due to precipitation in units of
                0550 c  mean snow thickness
                0551           SnowAccOverIce(I,J) =
                0552      &     SnowAccRateOverIce(I,J) * SEAICE_deltaTtherm * AREApreTH(I,J)
                0553 
                0554 c I,J
                0555          ENDDO
                0556         ENDDO
                0557 
                0558 C =============================================================================
                0559 C FIND WIND SPEED
                0560 C =============================================================================
                0561 
                0562         DO j=1,sNy
                0563          DO i=1,sNx
                0564 C ocean surface/mixed layer temperature
                0565           TmixLoc(i,j) = theta(i,j,kSurface,bi,bj) + celsius2K
                0566 C wind speed from exf
                0567           UG(I,J) = MAX(SEAICE_EPS, wspeed(I,J,bi,bj))
                0568          ENDDO
                0569         ENDDO
2383f7d4e9 Mart*0570 #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_GROWTH_ADX_STORE_MORE)
                0571 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
edb6656069 Mart*0572 CADJ STORE UG                 = comlev1_bibj,key=tkey,byte=isbyte
                0573 CADJ STORE SnowAccRateOverIce = comlev1_bibj,key=tkey,byte=isbyte
2383f7d4e9 Mart*0574 CADJ STORE PrecipRateOverIceSurfaceToSea
edb6656069 Mart*0575 CADJ &                        = comlev1_bibj,key=tkey,byte=isbyte
2383f7d4e9 Mart*0576 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0577 #endif /* ALLOW_AUTODIFF_TAMC */
dc54d31829 Ian *0578 
                0579 C =============================================================================
                0580 c           Retrieve the air-sea heat and shortwave radiative fluxes
                0581 C =============================================================================
                0582 
                0583         CALL SEAICE_BUDGET_OCEAN(
                0584      I       UG,
                0585      I       TmixLoc,
                0586      O       F_ao, QSWO,
                0587      I       bi, bj, myTime, myIter, myThid )
                0588 
2383f7d4e9 Mart*0589 #if (!defined SEAICE_EXTERNAL_FLUXES && defined ALLOW_AUTODIFF_TAMC)
                0590 C     In this case it might be better save the output of
                0591 C     seaice_budget_ocean rather than recompute the bulkformulae.
                0592 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
edb6656069 Mart*0593 CADJ STORE f_ao = comlev1_bibj, key = tkey, byte = isbyte
                0594 CADJ STORE qswo = comlev1_bibj, key = tkey, byte = isbyte
2383f7d4e9 Mart*0595 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0596 #endif
dc54d31829 Ian *0597 C =============================================================================
                0598 C      Calc air-sea fluxes in the uppermost grid cell
                0599 C =============================================================================
                0600 
8e32c48b8f Mart*0601 C--   Not all of the sw radiation is absorbed in the uppermost ocean
00c7090dc0 Mart*0602 C     grid cell layer. Only that fraction which converges in the
                0603 C     uppermost ocean grid cell is used to melt ice.
                0604 C     SWFrac3D - the fraction of incoming sw radiation absorbed in the
                0605 C     uppermost ocean grid cell (calculated in initialise_fixed.F)
8e32c48b8f Mart*0606 
00c7090dc0 Mart*0607 C The contribution of shortwave heating is
                0608 C not included without #define SHORTWAVE_HEATING
dc54d31829 Ian *0609 #ifdef SHORTWAVE_HEATING
00c7090dc0 Mart*0610         IF ( selectPenetratingSW .GE. 1 ) THEN
                0611          DO J=1,sNy
                0612           DO I=1,sNx
                0613            tmpscal0 = SWFrac3D(I,J,kSrfS,bi,bj)
                0614            QSWO_BELOW_FIRST_LAYER(I,J)= QSWO(I,J)*tmpscal0
                0615            QSWO_IN_FIRST_LAYER(I,J)   = QSWO(I,J)*(ONE - tmpscal0)
                0616           ENDDO
                0617          ENDDO
                0618 C     This is not  necessary, as both variables  are already initialised
                0619 C     to zero in the beginning of the routine.
                0620 c        ELSE
                0621 c#endif
                0622 c         DO J=1,sNy
                0623 c          DO I=1,sNx
                0624 c           QSWO_BELOW_FIRST_LAYER(I,J)= ZERO
                0625 c           QSWO_IN_FIRST_LAYER(I,J)   = ZERO
                0626 c          ENDDO
                0627 c         ENDDO
                0628 c#ifdef SHORTWAVE_HEATING
                0629         ENDIF
dc54d31829 Ian *0630 #endif
00c7090dc0 Mart*0631         DO J=1,sNy
                0632          DO I=1,sNx
                0633           IceGrowthRateOpenWater(I,J) = QI *
                0634      &         (F_ao(I,J) - QSWO(I,J) + QSWO_IN_FIRST_LAYER(I,J))
dc54d31829 Ian *0635          ENDDO
                0636         ENDDO
                0637 
                0638 #ifdef ALLOW_AUTODIFF_TAMC
2383f7d4e9 Mart*0639 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
4dd39c50d9 Mart*0640 CADJ STORE salt(:,:,kSurface,bi,bj)
edb6656069 Mart*0641 CADJ &                 = comlev1_bibj, key = tkey, byte = isbyte
dc54d31829 Ian *0642 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
2383f7d4e9 Mart*0643 #endif /* ALLOW_AUTODIFF_TAMC */
dc54d31829 Ian *0644 
                0645 C =============================================================================
                0646 c  calculate heat fluxes within ice (conduction), F_io, and across the
                0647 c  ice/atmosphere interface, F_ia
                0648 C =============================================================================
                0649 
                0650 C--   Start loop over multi-categories
                0651         DO IT=1,SEAICE_multDim
                0652 
                0653          DO J=1,sNy
                0654           DO I=1,sNx
                0655 c record prior ice surface temperatures
                0656            ticeInMult(I,J,IT)  = TICES(I,J,IT,bi,bj)
                0657            ticeOutMult(I,J,IT) = TICES(I,J,IT,bi,bj)
                0658            TICES(I,J,IT,bi,bj) = ZERO
                0659           ENDDO
                0660          ENDDO
                0661 
                0662 c set relative thickness of ice categories
                0663          pFac = (2.0 _d 0*IT - 1.0 _d 0)*recip_multDim
                0664          pFacSnow = 1. _d 0
                0665 
                0666 c find actual snow and ice thickness within categories categories
                0667          IF ( SEAICE_useMultDimSnow ) pFacSnow=pFac
                0668 
                0669          DO J=1,sNy
                0670           DO I=1,sNx
                0671             hiceActualMult(I,J,IT)   = hiceActual(I,J) *pFac
                0672             hsnowActualMult(I,J,IT)  = hsnowActual(I,J)*pFacSnow
                0673           ENDDO
                0674          ENDDO
4dd39c50d9 Mart*0675 C     multDim-loop
                0676         ENDDO
dc54d31829 Ian *0677 
2383f7d4e9 Mart*0678 #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_GROWTH_ADX_STORE_MORE)
dc54d31829 Ian *0679 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
edb6656069 Mart*0680 CADJ STORE hiceActualMult = comlev1_bibj, key = tkey, byte = isbyte
                0681 CADJ STORE hsnowActualMult= comlev1_bibj, key = tkey, byte = isbyte
dc54d31829 Ian *0682 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
2383f7d4e9 Mart*0683 #endif /* ALLOW_AUTODIFF_TAMC */
dc54d31829 Ian *0684 
                0685 c =======================================================================
                0686 c  find calculate heat fluxes within ice (conduction) and across the
                0687 c  ice/atmosphere interface for each thickness category
                0688 c =======================================================================
                0689 
                0690         DO IT=1,SEAICE_multDim
4dd39c50d9 Mart*0691          CALL SEAICE_SOLVE4TEMP(
dc54d31829 Ian *0692      I        UG, hiceActualMult(1,1,IT), hsnowActualMult(1,1,IT),
4dd39c50d9 Mart*0693 #ifdef SEAICE_CAP_SUBLIM
                0694          This error is put here intentionally, because SEAICE_CAP_SUBLIM
                0695          cannot be defined together with SEAICE_USE_GROWTH_ADX
                0696      I        latentHeatFluxMaxMult(1,1,IT),
                0697 #endif
                0698      I        ticeInMult(1,1,IT),
                0699      O        ticeOutMult(1,1,IT),
dc54d31829 Ian *0700      O        F_io_net_mult(1,1,IT),
                0701      O        F_ia_net_mult(1,1,IT),
                0702      O        F_ia_mult(1,1,IT),
                0703      O        QSWI_mult(1,1,IT),
                0704      O        FWsublim_mult(1,1,IT),
                0705      I        bi, bj, myTime, myIter, myThid )
                0706         ENDDO
                0707 
                0708 c =======================================================================
                0709 c  record the ice surface temperature in each category
                0710 c  and find the average of fluxes across each category
                0711 
                0712         DO IT=1,SEAICE_multDim
                0713          DO J=1,sNy
                0714           DO I=1,sNx
                0715 
                0716 C     update TICES
                0717             TICES(I,J,IT,bi,bj) = ticeOutMult(I,J,IT)
                0718 
                0719             F_io_net(I,J) = F_io_net(I,J) +
                0720      &        F_io_net_mult(I,J,IT)*recip_multDim
                0721 
                0722             F_ia_net(I,J) = F_ia_net(I,J) +
                0723      &        F_ia_net_mult(I,J,IT)*recip_multDim
                0724 
                0725             F_ia(I,J)  = F_ia(I,J)     +
                0726      &        F_ia_mult(I,J,IT)*recip_multDim
                0727 
                0728             QSWI(I,J)  = QSWI(I,J) + QSWI_mult(I,J,IT)*recip_multDim
                0729 
                0730             FWsublim(I,J) = FWsublim(I,J) +
                0731      &         FWsublim_mult(I,J,IT)*recip_multDim
                0732 
                0733           ENDDO
                0734          ENDDO
                0735         ENDDO
2383f7d4e9 Mart*0736 #ifdef ALLOW_AUTODIFF_TAMC
                0737 C     This store is needed to avoid recomputing seaice_solve4temp in AD-mode,
                0738 C     because F_ia_net is modified in the block below.
                0739 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
edb6656069 Mart*0740 CADJ STORE F_ia_net  = comlev1_bibj, key = tkey, byte = isbyte
2383f7d4e9 Mart*0741 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0742 #endif /* ALLOW_AUTODIFF_TAMC */
dc54d31829 Ian *0743 
                0744         DO J=1,sNy
                0745          DO I=1,sNx
                0746 c          If there is heat flux convergence at the snow surface,
4dd39c50d9 Mart*0747 c          use that energy to melt snow before melting ice.  It is
dc54d31829 Ian *0748 c          possible that some snow will remain after melting,
                0749 c          which will drive F_ia_net to zero, or that all of the
                0750 c          snow will be melted, leaving a nonzero F_ia_net to melt
                0751 c          some ice.
                0752           F_ia_net_before_snow(I,J) = F_ia_net(I,J)
                0753 
                0754 c         Only continue if there is snow and ice in the cell
                0755           IF (AREApreTH(I,J) .LE. ZERO) THEN
                0756             IceGrowthRateUnderExistingIce(I,J) = 0. _d 0
                0757             IceGrowthRateFromSurface(I,J)      = 0. _d 0
                0758             NetExistingIceGrowthRate(I,J)      = 0. _d 0
                0759           ELSE
                0760 c           The growth rate (m/s) beneath existing ice is given by the upward
                0761 c           ocean-ice conductive flux, F_io_net, and QI.
                0762             IceGrowthRateUnderExistingIce(I,J) = F_io_net(I,J)*QI
                0763 
                0764 c           The rate at which snow is melted (m/s) because of surface
                0765 c           heat flux convergence.  Note, during snow melt, F_ia_net must
                0766 c           be negative (implying convergence) to make PSMRFW is positive
                0767             PotSnowMeltRateFromSurf(I,J) = - F_ia_net(I,J)*QS
                0768 
                0769 c           This is the depth of snow (m) that would be melted in one dt
                0770             PotSnowMeltFromSurf(I,J)  =
                0771      &        PotSnowMeltRateFromSurf(I,J) * SEAICE_deltaTtherm
                0772 
                0773 c           If we can melt MORE than is actually there, then the melt
                0774 c           rate is reduced so that only that which is there
                0775 c           is melted during the time step.  In this case, not all of the
                0776 c           heat flux convergence at the surface is used to melt snow.
                0777 c           Any remaining energy will melt ice.
                0778 
                0779 c           SurfHeatFluxConvergToSnowMelt is the part of the total heat
                0780 c           flux convergence which melts snow.
                0781 
                0782 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
4dd39c50d9 Mart*0783 CCHECK: HSNOW ACTUAL SHOULD NOT BE REGULARIZED FOR THIS
dc54d31829 Ian *0784             IF (PotSnowMeltFromSurf(I,J) .GE. hsnowActual(I,J)) THEN
                0785 
                0786 c           Snow melt and melt rate [m] (actual snow thickness)
                0787               SnowMeltFromSurface(I,J)     = hsnowActual(I,J)
                0788 
                0789               SnowMeltRateFromSurface(I,J) =
4dd39c50d9 Mart*0790      &           SnowMeltFromSurface(I,J) * recip_deltaTtherm
dc54d31829 Ian *0791 
                0792               SurfHeatFluxConvergToSnowMelt(I,J) =
4dd39c50d9 Mart*0793      &          - hsnowActual(I,J)*recip_QS*recip_deltaTtherm
dc54d31829 Ian *0794            ELSE
                0795 c             In this case there will be snow remaining after melting.
                0796 c             All of the surface heat convergence will be redirected to
                0797 c             this effort.
                0798               SnowMeltFromSurface(I,J) = PotSnowMeltFromSurf(I,J)
                0799 
                0800               SnowMeltRateFromSurface(I,J) =PotSnowMeltRateFromSurf(I,J)
                0801 
                0802               SurfHeatFluxConvergToSnowMelt(I,J) = F_ia_net(I,J)
                0803 
                0804            ENDIF
                0805 
                0806 c          Reduce the heat flux convergence available to melt surface
                0807 c          ice by the amount used to melt snow
                0808            F_ia_net(I,J) = F_ia_net(I,J) -
                0809      &         SurfHeatFluxConvergToSnowMelt(I,J)
                0810 
                0811            IceGrowthRateFromSurface(I,J) = F_ia_net(I,J) * QI
                0812 
                0813 c          The total growth rate (m/s) of the existing ice - the rate of
                0814 c          new ice accretion at the base less the rate due to surface melt
                0815            NetExistingIceGrowthRate(I,J) =
                0816      &            IceGrowthRateUnderExistingIce(I,J) +
                0817      &            IceGrowthRateFromSurface(I,J)
                0818           ENDIF
                0819 
                0820          ENDDO
                0821         ENDDO
2383f7d4e9 Mart*0822 #ifdef ALLOW_AUTODIFF_TAMC
                0823 C     Store F_io/ia_net and QSWI to avoid calling seaice_solve4temp in AD-mode
                0824 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
edb6656069 Mart*0825 CADJ STORE F_io_net      = comlev1_bibj, key = tkey, byte = isbyte
                0826 CADJ STORE F_ia_net      = comlev1_bibj, key = tkey, byte = isbyte
                0827 CADJ STORE QSWI          = comlev1_bibj, key = tkey, byte = isbyte
2383f7d4e9 Mart*0828 CADJ STORE theta(:,:,kSurface,bi,bj)
edb6656069 Mart*0829 CADJ &                   = comlev1_bibj, key = tkey, byte = isbyte
2383f7d4e9 Mart*0830 CADJ STORE salt(:,:,kSurface,bi,bj)
edb6656069 Mart*0831 CADJ &                   = comlev1_bibj, key = tkey, byte = isbyte
2383f7d4e9 Mart*0832 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0833 #endif /* ALLOW_AUTODIFF_TAMC */
4dd39c50d9 Mart*0834 c     inflection point
                0835         tmpscal0 = 0.4 _d 0
                0836 c     steepness/inflection point
                0837         tmpscal1 = 7.0 _d 0/tmpscal0
                0838 c     avoid recomputing constant coefficients
                0839         tmpscal2 = STANTON_NUMBER*USTAR_BASE*rhoConst*HeatCapacity_Cp
2383f7d4e9 Mart*0840 C     Calculate the heat fluxes from the ocean to the sea ice
dc54d31829 Ian *0841         DO J=1,sNy
                0842          DO I=1,sNx
                0843 c
4dd39c50d9 Mart*0844 c     Bound the ocean temperature to be at or above the freezing point.
dc54d31829 Ian *0845           tempFrz = SEAICE_tempFrz0 +
                0846      &              SEAICE_dTempFrz_dS * salt(I,J,kSurface,bi,bj)
                0847 
                0848           surf_theta = max(theta(I,J,kSurface,bi,bj), tempFrz)
                0849 
4dd39c50d9 Mart*0850 C     MCPHEE_TAPER_FAC = 12.5 (CONSTANT PARAMETER)
                0851           MLTF = ONE + (MCPHEE_TAPER_FAC - ONE)
                0852      &         / ( ONE + EXP( (AREApreTH(I,J) - tmpscal0)*tmpscal1 ) )
dc54d31829 Ian *0853 
4dd39c50d9 Mart*0854 C          IF (AREApreTH(I,J) .GT. ZERO) THEN
                0855 CC     If ice is present, MixedLayerTurbulenceFactor = 1.0, else 12.50
                0856 C           MLTF = ONE
                0857 C          ELSE
                0858 C           MLTF = MCPHEE_TAPER_FAC
                0859 C          ENDIF
dc54d31829 Ian *0860 
4dd39c50d9 Mart*0861           F_mi(I,J) = - tmpscal2 * (surf_theta - tempFrz) * MLTF
dc54d31829 Ian *0862 
                0863           IceGrowthRateMixedLayer (I,J) = F_mi(I,J) * QI
                0864 
                0865          ENDDO
                0866         ENDDO
                0867 #ifdef ALLOW_AUTODIFF_TAMC
2383f7d4e9 Mart*0868 C     The first two store directives are also needed to avoid calling
                0869 C     seaice_solve4temp in AD-mode.
dc54d31829 Ian *0870 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
2383f7d4e9 Mart*0871 CADJ STORE NetExistingIceGrowthRate      = comlev1_bibj,
edb6656069 Mart*0872 CADJ &     key = tkey, byte = isbyte
2383f7d4e9 Mart*0873 CADJ STORE IceGrowthRateMixedLayer       = comlev1_bibj,
edb6656069 Mart*0874 CADJ &     key = tkey, byte = isbyte
2383f7d4e9 Mart*0875 # ifdef SEAICE_GROWTH_ADX_STORE_MORE
                0876 CADJ STORE SnowMeltRateFromSurface       = comlev1_bibj,
edb6656069 Mart*0877 CADJ &     key = tkey, byte = isbyte
2383f7d4e9 Mart*0878 CADJ STORE IceGrowthRateUnderExistingIce = comlev1_bibj,
edb6656069 Mart*0879 CADJ &     key = tkey, byte = isbyte
2383f7d4e9 Mart*0880 CADJ STORE IceGrowthRateFromSurface      = comlev1_bibj,
edb6656069 Mart*0881 CADJ &     key = tkey, byte = isbyte
2383f7d4e9 Mart*0882 # endif
                0883 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0884 #endif /* ALLOW_AUTODIFF_TAMC */
dc54d31829 Ian *0885 
                0886 C       CALCULATE THICKNESS DERIVATIVES of ice (dhdt) and snow (dhsdt)
                0887         DO J=1,sNy
                0888          DO I=1,sNx
                0889 
                0890           S_h(I,J) =
                0891      &       NetExistingIceGrowthRate(I,J) *  AREApreTH(I,J)
                0892      &     + IceGrowthRateOpenWater(I,J)   * (ONE - AREApreTH(I,J))
                0893      &     + IceGrowthRateMixedLayer(I,J)
                0894 
                0895 c         Both the accumulation and melt rates are in terms
                0896 c         of actual snow thickness.  As with ice, multiplying
                0897 c         with area converts to mean snow thickness.
                0898           S_hsnow(I,J) =  AREApreTH(I,J) *
                0899      &      ( SnowAccRateOverIce(I,J) - SnowMeltRateFromSurface(I,J))
                0900 
                0901          ENDDO
                0902         ENDDO
2383f7d4e9 Mart*0903 #ifdef ALLOW_AUTODIFF_TAMC
                0904 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0905 C     storing S_h would save one line of recomputation so do not do it
edb6656069 Mart*0906 cnostore CADJ STORE S_h = comlev1_bibj, key = tkey, byte = isbyte
2383f7d4e9 Mart*0907 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0908 #endif /* ALLOW_AUTODIFF_TAMC */
dc54d31829 Ian *0909 
                0910 #ifdef ALLOW_SALT_PLUME
4dd39c50d9 Mart*0911 # ifdef SALT_PLUME_IN_LEADS
                0912 C     avoid repetitive division by a constant
                0913         tmpscal0 = 5. _d 0 / (ONE - SPinflectionPoint)
                0914 C     Now that we know the thickness tendency terms, we can calculate
                0915 C     the saltPlumeFlux
dc54d31829 Ian *0916         DO J=1,sNy
                0917          DO I=1,sNx
                0918 
4dd39c50d9 Mart*0919 c     It is assumed that ice production in leads can generate
                0920 c     salt plumes. The fraction of the salt sent to the plume package
                0921 c     from ice produced in leads (defined as the open water
                0922 c     fraction) is a nonlinear function of ice concentration, AREA.
dc54d31829 Ian *0923 c
4dd39c50d9 Mart*0924 c     Specifically, function is a logistic curve (sigmoid) with a range and
                0925 c     domain {0,1}. The function, f(AREA), has a single free parameter,
                0926 c     SEAICE_plumeInflectionPoint, the inflection point of the curve.
                0927 c     By construction, the function has the following properties:
                0928 c     f(1) \approx 1.0
                0929 c     f(SEAICE_plumeInflectionPoint) = 0.5
                0930 c     f(0) \approx 0.0 (when SEAICE_plumeInflectionPoint \geq 0.5)
                0931 c     f(0) > 0.0 (when SEAICE_plumeInflectionPoint < 0.5)
dc54d31829 Ian *0932 c
4dd39c50d9 Mart*0933 c     As AREA --> 1, the open water fraction occurs in narrow leads,
                0934 c     new ice production become spatially non-uniform, and the assumptions
                0935 c     motivating KPP no longer hold. To treat overturning in a more
                0936 c     physically realistic way, the salt produced in the leads should
                0937 c     be sent to depth via the plume package. To assure only narrow leads
                0938 c     generate plumes, choose a SEAICE_plumeInflectionPoint of > 0.8.
                0939 
                0940 c     Ensure that there is already ice present or that the total ice
                0941 c     ice tendency term is positive.  We do not want to release
                0942 c     salt if sea ice is not established in the cell.
                0943 
                0944           IF (AREApreTH(I,J) .GT. ZERO .OR. S_h(I,J) .GT. ZERO ) THEN
                0945 
                0946            leadPlumeFraction(I,J) = ONE
                0947      &          / ( ONE + EXP( ( SPinflectionPoint - AREApreTH(I,J) )
                0948      &                         * tmpscal0 )
                0949      &            )
                0950 
                0951 c     Only consider positive ice growth rate in leads for salt production
                0952            IceGrowthRateInLeads(I,J) = max( ZERO,
                0953      &          (ONE - AREApreTH(I,J)) * IceGrowthRateOpenWater(I,J))
                0954 
                0955            saltPlumeFlux(I,J,bi,bj) = leadPlumeFraction(I,J) *
                0956      &          HEFFM(I,J,bi,bj)*IceGrowthRateInLeads(I,J)*
                0957      &          ICE2WATR*rhoConstFresh*
                0958      &          (salt(I,J,kSurface,bi,bj) - SEAICE_salt0)
                0959           ELSE
                0960 
                0961            saltPlumeFlux(I,J,bi,bj) = ZERO
                0962 
                0963           ENDIF
dc54d31829 Ian *0964 
                0965          ENDDO
                0966         ENDDO
                0967 #endif /* SALT_PLUME_IN_LEADS */
                0968 #endif /* ALLOW_SALT_PLUME */
                0969 
                0970 c       Caculate dA/dt (S_a)
                0971         DO J=1,sNy
                0972          DO I=1,sNx
                0973 
                0974           S_a(I,J) =  0. _d 0
                0975 
4dd39c50d9 Mart*0976 c     Caculate the ice area growth rate from the open water fluxes.
                0977 c     First, determine whether the open water growth rate is positive or
                0978 c     negative.  If positive, make sure that ice is present or that the
                0979 c     net ice thickness growth rate is positive before extending ice cover
dc54d31829 Ian *0980 
4dd39c50d9 Mart*0981 c     this is the geometric term: Area/(2*Heff),
                0982 c     with Area/Heff regularized as Area/(Heff^2 + epsilon^2)
                0983 c     Area/(Heff^2 + ep^2) = recip_hiceActual
                0984 c     epsilon = SEAICE_hice_reg
dc54d31829 Ian *0985 
4dd39c50d9 Mart*0986           tmpscal0 = 0.5 _d 0 * recip_hiceActual(I,J)
dc54d31829 Ian *0987 
                0988 C    -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
                0989 C         /* IceGrowthRateOpenWater */
                0990 
                0991           S_a_IGROW(I,J) = ZERO
                0992 
                0993 c         Expand ice cover if the open water growth rate is positive
                0994           IF ( IceGrowthRateOpenWater(I,J) .GT. ZERO) THEN
4dd39c50d9 Mart*0995            IF ( AREApreTH(I,J) .GT. ZERO   .OR.
                0996      &          S_h(I,J)       .GT. ZERO )  THEN
                0997 c     Determine which hemisphere for hemisphere-dependent
                0998 c     "lead closing variable", HO
                0999             IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN
                1000              S_a_IGROW(I,J) = (ONE - AREApreTH(I,J)) *
                1001      &            IceGrowthRateOpenWater(I,J)*recip_HO_south
                1002             ELSE
                1003              S_a_IGROW(I,J) = (ONE - AREApreTH(I,J)) *
                1004      &            IceGrowthRateOpenWater(I,J)*recip_HO
dc54d31829 Ian *1005             ENDIF
4dd39c50d9 Mart*1006            ENDIF
dc54d31829 Ian *1007 
                1008 C    -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
                1009 c         Contract ice cover if the open water growth rate is negative
                1010           ELSE
4dd39c50d9 Mart*1011            S_a_IGROW(I,J) = tmpscal0 *
dc54d31829 Ian *1012      &          IceGrowthRateOpenWater(I,J) * (ONE - AREApreTH(I,J))
                1013           ENDIF
                1014 
                1015           S_a(I,J) = S_a(I,J) + S_a_IGROW(I,J)
                1016 
                1017 C    -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
                1018 c         /* IceGrowthRateMixedLayer */
                1019 C         Contract ice if the IceGrowthRateMixedLayer is negative
                1020           S_a_IGRML(I,J) = ZERO
                1021 
                1022           IF ( IceGrowthRateMixedLayer(I,J) .LE. ZERO)  THEN
4dd39c50d9 Mart*1023            S_a_IGRML(I,J) = tmpscal0 * IceGrowthRateMixedLayer(I,J)
dc54d31829 Ian *1024           ENDIF
                1025 
                1026           S_a(I,J) = S_a(I,J) + S_a_IGRML(I,J)
                1027 
                1028 c
                1029 C    -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
                1030 C         Contract ice if the NetExistingIceGrowthRate is negative
                1031 C         /* NetExistingIceGrowthRate */
                1032           S_a_IGRNE(I,J) = ZERO
4dd39c50d9 Mart*1033           IF ( NetExistingIceGrowthRate(I,J) .LE. ZERO .AND.
                1034      &         HEFFpreTH(I,J)                .GT. ZERO ) THEN
dc54d31829 Ian *1035 
4dd39c50d9 Mart*1036            S_a_IGRNE(I,J) =
dc54d31829 Ian *1037      &       tmpscal0 * NetExistingIceGrowthRate(I,J) * AREApreTH(I,J)
                1038 
                1039           ENDIF
                1040 
                1041           S_a(I,J) = S_a(I,J) + S_a_IGRNE(I,J)
                1042 
2383f7d4e9 Mart*1043          ENDDO
                1044         ENDDO
dc54d31829 Ian *1045 
2383f7d4e9 Mart*1046 C     Update the area, heff, and hsnow
dc54d31829 Ian *1047         DO J=1,sNy
                1048          DO I=1,sNx
4dd39c50d9 Mart*1049           HEFF(I,J,bi,bj)  = HEFFpreTH(I,J) +
6b47d550f4 Mart*1050      &         SEAICE_deltaTtherm * S_h(I,J) * HEFFM(I,J,bi,bj)
dc54d31829 Ian *1051 
4dd39c50d9 Mart*1052           AREA(I,J,bi,bj)  = AREApreTH(I,J) +
6b47d550f4 Mart*1053      &         SEAICE_deltaTtherm * S_a(I,J) * HEFFM(I,J,bi,bj)
dc54d31829 Ian *1054 
4dd39c50d9 Mart*1055           HSNOW(I,J,bi,bj) = HSNWpreTH(I,J) +
6b47d550f4 Mart*1056      &         SEAICE_deltaTtherm * S_hsnow(I,J) * HEFFM(I,J,bi,bj)
dc54d31829 Ian *1057          ENDDO
                1058         ENDDO
                1059 
2383f7d4e9 Mart*1060 #ifdef ALLOW_AUTODIFF_TAMC
                1061 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
edb6656069 Mart*1062 CADJ STORE AREA (:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                1063 CADJ STORE HEFF (:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                1064 CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
2383f7d4e9 Mart*1065 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                1066 #endif /* ALLOW_AUTODIFF_TAMC */
dc54d31829 Ian *1067         DO J=1,sNy
4dd39c50d9 Mart*1068          DO I=1,sNx
                1069 c     Bound area, heff, and hsnow
                1070           tmpscal0 = AREA(I,J,bi,bj)
                1071           AREA(I,J,bi,bj)  = MIN(ONE,  tmpscal0)
                1072           tmpscal1 = AREA(I,J,bi,bj)
                1073           AREA(I,J,bi,bj)  = MAX(ZERO, tmpscal1)
                1074           tmpscal2 = HEFF(I,J,bi,bj)
                1075           HEFF(I,J,bi,bj)  = MAX(ZERO, tmpscal2)
                1076           tmpscal3 = HSNOW(I,J,bi,bj)
                1077           HSNOW(I,J,bi,bj) = MAX(ZERO, tmpscal3)
                1078 
                1079 c     Sanity checks
                1080           IF ( HEFF(I,J,bi,bj) .LE. ZERO .OR.
                1081      &         AREA(I,J,bi,bj) .LE. ZERO ) THEN
                1082 
                1083            AREA(I,J,bi,bj)       = 0. _d 0
                1084            HEFF(I,J,bi,bj)       = 0. _d 0
                1085            HSNOW(I,J,bi,bj)      = 0. _d 0
dc54d31829 Ian *1086 
4dd39c50d9 Mart*1087           ENDIF
dc54d31829 Ian *1088 
4dd39c50d9 Mart*1089          ENDDO
                1090         ENDDO
dc54d31829 Ian *1091 
00c7090dc0 Mart*1092 #ifdef SHORTWAVE_HEATING
                1093         IF ( selectPenetratingSW .GE. 1 ) THEN
                1094          DO J=1,sNy
                1095           DO I=1,sNx
dc54d31829 Ian *1096 
                1097 c         THE EFFECTIVE SHORTWAVE HEATING RATE
00c7090dc0 Mart*1098            QSW(I,J,bi,bj)  =
                1099      &          QSWI(I,J)  * (      AREApreTH(I,J)) +
                1100      &          QSWO(I,J)  * (ONE - AREApreTH(I,J))
                1101           ENDDO
                1102          ENDDO
                1103         ELSE
dc54d31829 Ian *1104 #endif
00c7090dc0 Mart*1105          DO J=1,sNy
                1106           DO I=1,sNx
                1107            QSW(I,J,bi,bj) = 0. _d 0
                1108           ENDDO
4dd39c50d9 Mart*1109          ENDDO
00c7090dc0 Mart*1110 #ifdef SHORTWAVE_HEATING
                1111         ENDIF
                1112 #endif
4dd39c50d9 Mart*1113 
                1114         DO J=1,sNy
                1115          DO I=1,sNx
dc54d31829 Ian *1116 
4dd39c50d9 Mart*1117 c     The actual ice volume change over the time step
                1118           ActualNewTotalVolumeChange(I,J) =
                1119      &         HEFF(I,J,bi,bj) - HEFFpreTH(I,J)
                1120 
                1121 c     The net average snow thickness melt that is actually realized. e.g.
                1122 c     hsnow_orig  = 0.25 m (e.g. 1 m of ice over a cell 1/4 covered in snow)
                1123 c     hsnow_new   = 0.20 m
                1124 c     snow accum  = 0.05 m
                1125 c     melt = 0.25 + 0.05 - 0.2 = 0.1 m
                1126 
                1127 c     since this is in mean snow thickness it might have been  0.4 of actual
                1128 c     snow thickness over the 1/4 of the cell which is ice covered.
                1129           ActualNewTotalSnowMelt(I,J) =
                1130      &         HSNWpreTH(I,J) +
                1131      &         SnowAccOverIce(I,J) -
                1132      &         HSNOW(I,J,bi,bj)
                1133 
                1134 c     The energy required to melt or form the new ice volume
                1135           EnergyInNewTotalIceVolume(I,J) =
                1136      &         ActualNewTotalVolumeChange(I,J)*recip_QI
dc54d31829 Ian *1137 c
                1138 c     This is the net energy flux out of the ice+ocean system
                1139 c     Remember -----
                1140 c     F_ia_net : Under ice/snow surface freezing conditions,
                1141 c                vertical conductive heat flux convergence (F_c < 0) balances
                1142 c                heat flux divergence to atmosphere (F_ia > 0)
                1143 c                Otherwise, F_ia_net = F_ia (pos)
                1144 c
                1145 c     F_io_net : Under ice/snow surface freezing conditions, F_c < 0.
                1146 c                Under ice surface melting conditions, F_c = 0 (no energy flux
                1147 c                from the ice to ocean)
                1148 c
                1149 c     So if we are freezing, F_io_net = the conductive flux and there
                1150 c     is energy balance at ice surface, F_ia_net =0.  If we are melting,
                1151 c     there is a convergence of energy into the ice from above
4dd39c50d9 Mart*1152           NetEnergyFluxOutOfOcean(I,J) = SEAICE_deltaTtherm *
dc54d31829 Ian *1153      &               ( AREApreTH(I,J) *
                1154      &                 (F_ia_net(I,J) + F_io_net(I,J) + QSWI(I,J))
                1155      &         +     ( ONE - AREApreTH(I,J)) *  F_ao(I,J))
                1156 
                1157 c     THE QUANTITY OF HEAT WHICH IS THE RESIDUAL TO THE QUANTITY OF
                1158 c     ML temperature.  If the net energy flux is exactly balanced by the
                1159 c     latent energy of fusion in the new ice created then we will not
                1160 c     change the ML temperature at all.
                1161 
4dd39c50d9 Mart*1162           ResidualEnergyOutOfOcean(I,J) =
dc54d31829 Ian *1163      &             NetEnergyFluxOutOfOcean(I,J) -
                1164      &             EnergyInNewTotalIceVolume(I,J)
                1165 
                1166 C     NOW FORMULATE QNET
                1167 C     THIS QNET DETERMINES THE TEMPERATURE CHANGE
                1168 C     QNET IS A DEPTH AVERAGED HEAT FLUX FOR THE OCEAN COLUMN
                1169 
4dd39c50d9 Mart*1170           QNET(I,J,bi,bj) =
                1171      &             ResidualEnergyOutOfOcean(I,J) * recip_deltaTtherm
dc54d31829 Ian *1172 
4dd39c50d9 Mart*1173 c     Like snow melt, if there is melting, this quantity is positive.
                1174 c     The change of freshwater content is per unit area over the entire
                1175 c     cell, not just over the ice covered bits.  This term is only used
                1176 c     to calculate freshwater fluxes for the purpose of changing the
                1177 c     salinity of the liquid cell.  In the case of non-zero ice salinity,
                1178 c     the amount of freshwater is reduced by the ratio of ice salinity
                1179 c     to water cell salinity.
                1180           IF  (salt(I,J,kSurface,bi,bj) .GE. SEAICE_salt0 .AND.
                1181      &         salt(I,J,kSurface,bi,bj) .GT. 0. _d 0) THEN
dc54d31829 Ian *1182 
4dd39c50d9 Mart*1183            FreshwaterContribFromIce(I,J) =
dc54d31829 Ian *1184      &            - ActualNewTotalVolumeChange(I,J) *
4dd39c50d9 Mart*1185      &              rhoIce2rhoFresh
dc54d31829 Ian *1186 C     &             * (ONE - SEAICE_salt0/salt(I,J,kSurface,bi,bj))
                1187 
4dd39c50d9 Mart*1188           ELSE
                1189 C     If the liquid cell has a lower salinity than the specified
                1190 c     salinity of sea ice then assume the sea ice is completely fresh
                1191            FreshwaterContribFromIce(I,J) =
dc54d31829 Ian *1192      &             -ActualNewTotalVolumeChange(I,J) *
4dd39c50d9 Mart*1193      &              rhoIce2rhoFresh
                1194           ENDIF
dc54d31829 Ian *1195 CAB
                1196           tmpscal3 = max( 0. _d 0,
                1197      &                    min(SEAICE_salt0,salt(I,J,kSurface,bi,bj)) )
4dd39c50d9 Mart*1198 CAB   Salt in ice + ridgid - ice from snow flood
                1199           tmpscal2 = (ActualNewTotalVolumeChange(I,J)+
dc54d31829 Ian *1200      &                   SIheffNeg(I,J,bi,bj)+
                1201      &                  0.0 )
                1202      &                   * tmpscal3
                1203      &                   * HEFFM(I,J,bi,bj)
4dd39c50d9 Mart*1204      &                   * recip_deltaTtherm * SEAICE_rhoIce
dc54d31829 Ian *1205 
4dd39c50d9 Mart*1206 CAB   test for getting the lines
                1207           saltflux(I,J,bi,bj) = tmpscal2
dc54d31829 Ian *1208 CAB
                1209 
4dd39c50d9 Mart*1210 c     The freshwater contribution from snow comes only in the form of melt
                1211 c     unlike ice, which takes freshwater upon growth and yields freshwater
                1212 c     upon melt.  This is why the the actual new average snow melt was
                1213 c     determined. In m/m^2 over the entire cell.
                1214           FreshwaterContribFromSnowMelt(I,J) =
                1215      &         ActualNewTotalSnowMelt(I,J)*rhoSnow2rhoFresh
dc54d31829 Ian *1216 
                1217 c    This seems to be in m/s, original time level 2 for area
                1218 c    Only the precip and evap need to be area weighted.  The runoff
                1219 c    and freshwater contribs from ice and snow melt are already mean
                1220 c    weighted
4dd39c50d9 Mart*1221           EmPmR(I,J,bi,bj)  = HEFFM(I,J,bi,bj)*(
                1222      &         ( EVAP(I,J,bi,bj) - PRECIP(I,J,bi,bj) )
                1223      &         * ( ONE - AREApreTH(I,J) )
                1224      &         - PrecipRateOverIceSurfaceToSea(I,J)*AREApreTH(I,J)
dc54d31829 Ian *1225 #ifdef ALLOW_RUNOFF
4dd39c50d9 Mart*1226      &         - RUNOFF(I,J,bi,bj)
dc54d31829 Ian *1227 #endif
4dd39c50d9 Mart*1228      &         - (FreshwaterContribFromIce(I,J) +
                1229      &            FreshwaterContribFromSnowMelt(I,J))
                1230      &           *recip_deltaTtherm ) * rhoConstFresh
                1231      &         +( SIheffNeg(I,J,bi,bj)*SEAICE_rhoIce
                1232      &          + SIhsnwNeg(I,J,bi,bj)*SEAICE_rhoSnow)
                1233      &          *recip_deltaTtherm * HEFFM(I,J,bi,bj)
                1234 
dc54d31829 Ian *1235          ENDDO
                1236         ENDDO
                1237 
                1238 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                1239 
d73d15cdef antn*1240 C Treat flooding after calculations of QNET and before Sea Ice loading
                1241         IF(SEAICEuseFlooding) THEN
                1242 #ifdef ALLOW_AUTODIFF_TAMC
                1243 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                1244 CADJ STORE HEFF (:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                1245 CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                1246 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                1247 #endif /* ALLOW_AUTODIFF_TAMC */
                1248 C     convert snow to ice if submerged.
                1249          DO J=1,sNy
                1250           DO I=1,sNx
                1251            tmpscal0           = ( HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
                1252      &              +HEFF(I,J,bi,bj)*SEAICE_rhoIce )*recip_rhoConst
                1253            d_HEFFbyFLOODING(I,J) = MAX(0. _d 0,tmpscal0-HEFF(I,J,bi,bj))
                1254            HEFF(I,J,bi,bj)    = HEFF(I,J,bi,bj) + d_HEFFbyFLOODING(I,J)
                1255            HSNOW(I,J,bi,bj)   = HSNOW(I,J,bi,bj)
                1256      &                           - d_HEFFbyFLOODING(I,J)*ICE2SNOW
                1257 
                1258           ENDDO
                1259          ENDDO
                1260 c     SEAICEuseFlooding
                1261         ENDIF
                1262 #ifdef ALLOW_DIAGNOSTICS
                1263         IF (useDiagnostics) THEN
                1264          tmpscal1=1. _d 0 * recip_deltaTtherm
                1265          CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyFLOODING,
                1266      &      tmpscal1,1,'SIdHbFLO',0,1,3,bi,bj,myThid)
                1267         ENDIF
                1268 #endif /* ALLOW_DIAGNOSTICS */
                1269 
4dd39c50d9 Mart*1270 C     Sea Ice Load on the sea surface.
                1271 C     =================================
dc54d31829 Ian *1272         IF ( useRealFreshWaterFlux ) THEN
                1273          DO J=1,sNy
                1274           DO I=1,sNx
                1275 #ifdef SEAICE_CAP_ICELOAD
                1276            tmpscal1 = HEFF(I,J,bi,bj)*SEAICE_rhoIce
                1277      &              + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
                1278            tmpscal2 = MIN(tmpscal1,heffTooHeavy*rhoConst)
                1279 #else
                1280            tmpscal2 = HEFF(I,J,bi,bj)*SEAICE_rhoIce
                1281      &              + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
                1282 #endif
                1283            sIceLoad(i,j,bi,bj) = tmpscal2
                1284           ENDDO
                1285          ENDDO
                1286         ENDIF
                1287 
d73d15cdef antn*1288 C Everything down here is diagnostic, and should reflect what we
                1289 C calculated above and passed out to the ocean instead of computed
                1290 C again.  Otherwise we can end up with mismatch when we modify above
                1291 C and not replicate down here.
                1292 #ifdef ALLOW_DIAGNOSTICS
                1293         IF (useDiagnostics) THEN
                1294          DO J=1,sNy
                1295           DO I=1,sNx
                1296            SItflux(I,J,bi,bj) = 0.0 +
                1297 c     &        (AREApreTH(I,J)
                1298 c     &        *(F_ia_net(I,J)
                1299 c     &        + F_io_net(I,J) + QSWI(I,J))
                1300 c     &        +( ONE - AREApreTH(I,J)) *  F_ao(I,J))
                1301      &          NetEnergyFluxOutOfOcean(I,J)*recip_deltaTtherm
                1302      &          *HEFFM(I,J,bi,bj)
                1303      &          +(SIheffNeg(I,J,bi,bj)*recip_QI
                1304      &          +( SIhsnwNeg(I,J,bi,bj) - ActualNewTotalSnowMelt(I,J)
                1305      &           + SnowAccOverIce(I,J) )*recip_QS)
                1306      &          * recip_deltaTtherm * HEFFM(I,J,bi,bj)
                1307 
                1308            SIatmFW(I,J,bi,bj) = HEFFM(I,J,bi,bj)*(
                1309      &          EVAP(I,J,bi,bj)*( ONE - AREApreTH(I,J) )
                1310      &          - PRECIP(I,J,bi,bj)
                1311 #ifdef ALLOW_RUNOFF
                1312      &          - RUNOFF(I,J,bi,bj)
                1313 #endif /* ALLOW_RUNOFF */
                1314      &          )*rhoConstFresh
                1315           ENDDO
                1316          ENDDO
                1317 
                1318          IF (useRealFreshWaterFlux.AND.(nonlinFreeSurf.NE.0)) THEN
                1319           DO J=1,sNy
                1320            DO I=1,sNx
                1321 C Heat flux term associated with EmPmR:
                1322             SIeprflx(I,J,bi,bj) =  -EmPmR(I,J,bi,bj)
                1323      &           * HeatCapacity_Cp *theta(I,J,kSurface,bi,bj)
                1324            ENDDO
                1325           ENDDO
                1326 c        ELSE
                1327 CAB I  HAVE NOT FIGURED  OUT the QNET for the ELSE case, should not matter 01/14/2021 L2275
                1328          ENDIF
                1329         ENDIF
                1330 #endif /* ALLOW_DIAGNOSTICS */
                1331 
dc54d31829 Ian *1332 #ifdef SEAICE_DEBUG
4dd39c50d9 Mart*1333         DO j=1,sNy
                1334          DO i=1,sNx
                1335 
                1336           IF ( (i .EQ. SEAICE_debugPointI)   .and.
                1337      &         (j .EQ. SEAICE_debugPointJ) ) THEN
                1338 
                1339            print *,'ifice: myTime,myIter:',myTime,myIter
                1340 
                1341            print  '(A,2i4,2(1x,1P3E15.4))',
                1342      &          'ifice i j --------------  ',i,j
                1343 
                1344            print  '(A,2i4,2(1x,1P3E15.4))',
                1345      &          'ifice i j F(mi ao), rHIA  ',
                1346      &          i,j, F_mi(i,j), F_ao(i,j),
                1347      &          recip_hiceActual(i,j)
                1348 
                1349            print  '(A,2i4,2(1x,1P3E15.4))',
                1350      &          'ifice i j Fi(a,ant2/1 ont)',
                1351      &          i,j, F_ia(i,j),
                1352      &          F_ia_net_before_snow(i,j),
                1353      &          F_ia_net(i,j),
                1354      &          F_io_net(i,j)
                1355 
                1356            print  '(A,2i4,2(1x,1P3E15.4))',
                1357      &          'ifice i j AREA2/1 HEFF2/1 ',i,j,
                1358      &          AREApreTH(I,J),
                1359      &          AREA(i,j,bi,bj),
                1360      &          HEFFpreTH(I,J),
                1361      &          HEFF(i,j,bi,bj)
                1362 
                1363            print  '(A,2i4,2(1x,1P3E15.4))',
                1364      &          'ifice i j HSNOW2/1 TMX    ',i,j,
                1365      &          HSNWpreTH(I,J),
                1366      &          HSNOW(I,J,bi,bj),
                1367      &          theta(I,J,kSurface,bi,bj)
                1368 
                1369            print  '(A,2i4,2(1x,1P3E15.4))',
                1370      &          'ifice i j TI ATP LWD      ',i,j,
                1371      &          TICES(i,j,1, bi,bj) - celsius2k,
                1372      &          ATEMP(i,j,bi,bj) - celsius2k,
                1373      &          LWDOWN(i,j,bi,bj)
                1374 
                1375            print  '(A,2i4,2(1x,1P3E15.4))',
                1376      &          'ifice i j S_a(tot,OW,ML,NE',i,j,
                1377      &          S_a(i,j),
                1378      &          S_a_IGROW(I,J),
                1379      &          S_a_IGRML(I,J),
                1380      &          S_a_IGRNE(I,J)
                1381 
                1382            print  '(A,2i4,2(1x,1P3E15.4))',
                1383      &          'ifice i j S_a S_h S_hsnow ',i,j,
                1384      &          S_a(i,j),
                1385      &          S_h(i,j),
                1386      &          S_hsnow(i,j)
                1387 
                1388            print  '(A,2i4,2(1x,1P3E15.4))',
                1389      &          'ifice i j IGR(ML OW ICE)  ',i,j,
                1390      &          IceGrowthRateMixedLayer(i,j),
                1391      &          IceGrowthRateOpenWater(i,j),
                1392      &          NetExistingIceGrowthRate(i,j)
                1393 
                1394            print  '(A,2i4,2(1x,1P3E15.4))',
                1395      &          'ifice i j IVC(A ENIN)     ',i,j,
                1396      &          ActualNewTotalVolumeChange(i,j),
                1397      &          EnergyInNewTotalIceVolume(i,j)
                1398 c     &          ExpectedIceVolumeChange(i,j),
                1399 
                1400            print  '(A,2i4,2(1x,1P3E15.4))',
                1401      &          'ifice i j EF(NOS RE) QNET ',i,j,
                1402      &          NetEnergyFluxOutOfOcean(i,j),
                1403      &          ResidualEnergyOutOfOcean(i,j),
                1404      &          QNET(I,J,bi,bj)
                1405 
                1406            print  '(A,2i4,3(1x,1P3E15.4))',
                1407      &          'ifice i j QSW QSWO QSWI   ',i,j,
                1408      &          QSW(i,j,bi,bj),
                1409      &          QSWO(i,j),
                1410      &          QSWI(i,j)
dc54d31829 Ian *1411 
                1412 c                  print  '(A,2i4,3(1x,1P3E15.4))',
                1413 c     &                 'ifice i j SW(BML IML SW)  ',i,j,
                1414 c     &                 QSW_absorb_below_first_layer(i,j),
                1415 c     &                 QSW_absorb_in_first_layer(i,j),
8e32c48b8f Mart*1416 c     &                 SEAICE_SWFrac
dc54d31829 Ian *1417 
                1418 c                  print  '(A,2i4,3(1x,1P3E15.4))',
                1419 c     &                 'ifice i j ptc(to, qsw, oa)',i,j,
                1420 c     &                 PredTempChange(i,j),
                1421 c     &                 PredTempChangeFromQSW (i,j),
                1422 c     &                 PredTempChangeFromOA_MQNET(i,j)
                1423 
                1424 c                  print  '(A,2i4,3(1x,1P3E15.4))',
                1425 c     &                 'ifice i j ptc(fion,ian,ia)',i,j,
                1426 c     &                 PredTempChangeFromF_IO_NET(i,j),
                1427 c     &                 PredTempChangeFromF_IA_NET(i,j),
                1428 c     &                 PredTempChangeFromFIA(i,j)
                1429 
                1430 c                  print  '(A,2i4,3(1x,1P3E15.4))',
                1431 c     &                 'ifice i j ptc(niv)        ',i,j,
                1432 c     &                 PredTempChangeFromNewIceVol(i,j)
                1433 
4dd39c50d9 Mart*1434            print  '(A,2i4,3(1x,1P3E15.4))',
                1435      &          'ifice i j EmPmR EVP PRE RU',i,j,
                1436      &          EmPmR(I,J,bi,bj),
                1437      &          EVAP(I,J,bi,bj),
                1438      &          PRECIP(I,J,bi,bj),
                1439      &          RUNOFF(I,J,bi,bj)
                1440 
                1441            print  '(A,2i4,3(1x,1P3E15.4))',
                1442      &          'ifice i j PRROIS,SAOI(R .)',i,j,
                1443      &          PrecipRateOverIceSurfaceToSea(I,J),
                1444      &          SnowAccRateOverIce(I,J),
                1445      &          SnowAccOverIce(I,J)
                1446 
                1447            print  '(A,2i4,4(1x,1P3E15.4))',
                1448      &          'ifice i j SM(PM PMR . .R) ',i,j,
                1449      &          PotSnowMeltFromSurf(I,J),
                1450      &          PotSnowMeltRateFromSurf(I,J),
                1451      &          SnowMeltFromSurface(I,J),
                1452      &          SnowMeltRateFromSurface(I,J)
                1453 
                1454            print  '(A,2i4,4(1x,1P3E15.4))',
                1455      &          'ifice i j TotSnwMlt       ',i,j,
                1456      &          ActualNewTotalSnowMelt(I,J)
dc54d31829 Ian *1457 c     &                 ExpectedSnowVolumeChange(I,J)
                1458 
4dd39c50d9 Mart*1459            print  '(A,2i4,4(1x,1P3E15.4))',
                1460      &          'ifice i j fw(CFICE, CFSM) ',i,j,
                1461      &          FreshwaterContribFromIce(I,J),
                1462      &          FreshwaterContribFromSnowMelt(I,J)
dc54d31829 Ian *1463 
4dd39c50d9 Mart*1464            print  '(A,2i4,2(1x,1P3E15.4))',
                1465      &          'ifice i j --------------  ',i,j
dc54d31829 Ian *1466 
4dd39c50d9 Mart*1467           ENDIF
dc54d31829 Ian *1468 
4dd39c50d9 Mart*1469          ENDDO
                1470         ENDDO
dc54d31829 Ian *1471 #endif /* SEAICE_DEBUG */
                1472 
                1473 C close bi,bj loops
                1474        ENDDO
                1475       ENDDO
                1476 
                1477 #ifdef ALLOW_DIAGNOSTICS
                1478       IF ( useDiagnostics ) THEN
                1479 C these diags need to be done outside of the bi,bj loop so that
                1480 C we may do potential global mean adjustement to them consistently.
                1481 
4dd39c50d9 Mart*1482        CALL DIAGNOSTICS_FILL(SItflux, 'SItflux ',0,1,0,1,1,myThid)
a4e168e012 antn*1483        CALL DIAGNOSTICS_FILL(SIeprflx,'SIeprflx',0,1,0,1,1,myThid)
dc54d31829 Ian *1484 
d73d15cdef antn*1485 C SIatmFW follows the same convention as empmr -- SIatmFW diag does not
4dd39c50d9 Mart*1486        tmpscal1= - 1. _d 0
                1487        CALL DIAGNOSTICS_SCALE_FILL(SIatmFW,
                1488      &      tmpscal1,1,'SIatmFW ',0,1,0,1,1,myThid)
dc54d31829 Ian *1489       ENDIF
                1490 #endif /* ALLOW_DIAGNOSTICS */
                1491 
                1492 #else  /* ALLOW_EXF and ALLOW_ATM_TEMP */
4dd39c50d9 Mart*1493       STOP 'SEAICE_GROWTH_ADX not compiled with EXF and ALLOW_ATM_TEMP'
dc54d31829 Ian *1494 #endif /* ALLOW_EXF and ALLOW_ATM_TEMP */
5337b47ce7 Jean*1495 #endif /* SEAICE_USE_GROWTH_ADX */
dc54d31829 Ian *1496 
                1497       RETURN
                1498       END