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