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