pkg/seaice/seaice_growth_adx.F

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