pkg/seaice/seaice_growth.F

33e17487ce Dimi*0001 #include "SEAICE_OPTIONS.h"
66d21a8387 Jean*0002 #ifdef ALLOW_EXF
                0003 # include "EXF_OPTIONS.h"
                0004 #endif
cafd3818b7 An T*0005 #ifdef ALLOW_SALT_PLUME
                0006 # include "SALT_PLUME_OPTIONS.h"
                0007 #endif
772b2ed80e Gael*0008 #ifdef ALLOW_AUTODIFF
                0009 # include "AUTODIFF_OPTIONS.h"
                0010 #endif
33e17487ce Dimi*0011
                0012 CBOP
                0013 C     !ROUTINE: SEAICE_GROWTH
                0014 C     !INTERFACE:
                0015       SUBROUTINE SEAICE_GROWTH( myTime, myIter, myThid )
                0016 C     !DESCRIPTION: \bv
                0017 C     *==========================================================*
                0018 C     | SUBROUTINE seaice_growth
                0019 C     | o Updata ice thickness and snow depth
                0020 C     *==========================================================*
                0021 C     \ev
                0022
                0023 C     !USES:
                0024       IMPLICIT NONE
                0025 C     === Global variables ===
                0026 #include "SIZE.h"
                0027 #include "EEPARAMS.h"
                0028 #include "PARAMS.h"
                0029 #include "DYNVARS.h"
                0030 #include "GRID.h"
                0031 #include "FFIELDS.h"
ccaa3c61f4 Patr*0032 #include "SEAICE_SIZE.h"
33e17487ce Dimi*0033 #include "SEAICE_PARAMS.h"
                0034 #include "SEAICE.h"
ccaa3c61f4 Patr*0035 #include "SEAICE_TRACER.h"
33e17487ce Dimi*0036 #ifdef ALLOW_EXF
                0037 # include "EXF_PARAM.h"
66d21a8387 Jean*0038 # include "EXF_FIELDS.h"
33e17487ce Dimi*0039 #endif
                0040 #ifdef ALLOW_SALT_PLUME
                0041 # include "SALT_PLUME.h"
                0042 #endif
                0043 #ifdef ALLOW_AUTODIFF_TAMC
                0044 # include "tamc.h"
                0045 #endif
                0046
                0047 C     !INPUT/OUTPUT PARAMETERS:
                0048 C     === Routine arguments ===
                0049 C     myTime :: Simulation time
                0050 C     myIter :: Simulation timestep number
                0051 C     myThid :: Thread no. that called this routine.
                0052       _RL myTime
                0053       INTEGER myIter, myThid
66d21a8387 Jean*0054 CEOP
33e17487ce Dimi*0055
4dd39c50d9 Mart*0056 #ifndef SEAICE_USE_GROWTH_ADX
634144d037 Jean*0057 #if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP)
2651ba3350 Jean*0058 C     !FUNCTIONS:
                0059 #ifdef ALLOW_DIAGNOSTICS
                0060       LOGICAL  DIAGNOSTICS_IS_ON
                0061       EXTERNAL DIAGNOSTICS_IS_ON
                0062 #endif
                0063
33e17487ce Dimi*0064 C     !LOCAL VARIABLES:
                0065 C     === Local variables ===
c50ad14e64 Gael*0066 C
0c0ecd4c7b Jean*0067 C unit/sign convention:
                0068 C    Within the thermodynamic computation all stocks, except HSNOW,
c50ad14e64 Gael*0069 C      are in 'effective ice meters' units, and >0 implies more ice.
0c0ecd4c7b Jean*0070 C    This holds for stocks due to ocean and atmosphere heat,
                0071 C      at the outset of 'PART 2: determine heat fluxes/stocks'
c50ad14e64 Gael*0072 C      and until 'PART 7: determine ocean model forcing'
                0073 C    This strategy minimizes the need for multiplications/divisions
0c0ecd4c7b Jean*0074 C      by ice fraction, heat capacity, etc. The only conversions that
                0075 C      occurs are for the HSNOW (in effective snow meters) and
                0076 C      PRECIP (fresh water m/s).
2651ba3350 Jean*0077 C
                0078 C HEFF is effective Hice thickness (m3/m2)
                0079 C HSNOW is Heffective snow thickness (m3/m2)
                0080 C HSALT is Heffective salt content (g/m2)
                0081 C AREA is the seaice cover fraction (0<=AREA<=1)
                0082 C Q denotes heat stocks -- converted to ice stocks (m3/m2) early on
                0083 C
                0084 C For all other stocks/increments, such as d_HEFFbyATMonOCN
                0085 C or a_QbyATM_cover, the naming convention is as follows:
                0086 C    The prefix 'a_' means available, the prefix 'd_' means delta
                0087 C       (i.e. increment), and the prefix 'r_' means residual.
                0088 C    The suffix '_cover' denotes a value for the ice covered fraction
                0089 C       of the grid cell, whereas '_open' is for the open water fraction.
                0090 C    The main part of the name states what ice/snow stock is concerned
                0091 C       (e.g. QbyATM or HEFF), and how it is affected (e.g. d_HEFFbyATMonOCN
                0092 C       is the increment of HEFF due to the ATMosphere extracting heat from the
                0093 C       OCeaN surface, or providing heat to the OCeaN surface).
c50ad14e64 Gael*0094
33e17487ce Dimi*0095 C     i,j,bi,bj :: Loop counters
                0096       INTEGER i, j, bi, bj
                0097 C     number of surface interface layer
                0098       INTEGER kSurface
8377b8ee87 Mart*0099 C     IT        :: ice thickness category index (MULTICATEGORIES and ITD code)
286983d3d2 Patr*0100       INTEGER IT
8377b8ee87 Mart*0101 C     msgBuf    :: Informational/error message buffer
75bbb16cce Jean*0102 #ifdef ALLOW_BALANCE_FLUXES
                0103       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0104 #elif (defined (SEAICE_DEBUG))
5ffca4a0e2 Patr*0105       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0106       CHARACTER*12 msgBufForm
75bbb16cce Jean*0107 #endif
33e17487ce Dimi*0108 C     constants
bb8e6379cb Mart*0109       _RL tempFrz, ICE2SNOW, SNOW2ICE
                0110       _RL QI, QS, recip_QI
1d74e34c65 Jean*0111       _RL lhSublim
                0112
                0113 C conversion factors to go from Q (W/m2) to HEFF (ice meters)
                0114       _RL convertQ2HI, convertHI2Q
                0115 C conversion factors to go from precip (m/s) unit to HEFF (ice meters)
                0116       _RL convertPRECIP2HI, convertHI2PRECIP
                0117 C     Factor by which we increase the upper ocean friction velocity (u*) when
                0118 C     ice is absent in a grid cell  (dimensionless)
                0119       _RL MixedLayerTurbulenceFactor
                0120
                0121 C     wind speed square
                0122       _RL SPEED_SQ
                0123
                0124 C     Regularization values squared
                0125       _RL area_reg_sq, hice_reg_sq
                0126 C     pathological cases thresholds
                0127       _RL heffTooHeavy
0320e25227 Mart*0128 C     local copy of surface layer thickness in meters
                0129       _RL dzSurf
1d74e34c65 Jean*0130
                0131 C     Helper variables: reciprocal of some constants
                0132       _RL recip_multDim
                0133       _RL recip_deltaTtherm
                0134       _RL recip_rhoIce
                0135 C     local value (=1/HO or 1/HO_south)
                0136       _RL recip_HO
                0137 C     local value (=1/ice thickness)
                0138       _RL recip_HH
6571f3ca98 Jean*0139 #ifndef SEAICE_ITD
74c037b5fb Mart*0140 C     facilitate multi-category snow implementation
6571f3ca98 Jean*0141       _RL pFac, pFacSnow
                0142 #endif
4b6d456764 Mart*0143 C     additional factors accounting for a non-uniform sea-ice PDF
                0144       _RL denominator, recip_denominator, areaPDFfac
1d74e34c65 Jean*0145
                0146 C     temporary variables available for the various computations
                0147       _RL tmpscal0, tmpscal1, tmpscal2, tmpscal3, tmpscal4
                0148
                0149 #ifdef ALLOW_SITRACER
                0150       INTEGER iTr
                0151 #ifdef ALLOW_DIAGNOSTICS
                0152       CHARACTER*8   diagName
                0153 #endif
f61838dfc1 Torg*0154 #ifdef SEAICE_GREASE
                0155       INTEGER iTrGrease
                0156       _RL greaseDecayTime
                0157       _RL greaseNewFrazil
                0158       _RL THIRD
                0159       PARAMETER (THIRD = 1.0 _d 0 / 3.0 _d 0)
                0160 #endif
1d74e34c65 Jean*0161 #endif /* ALLOW_SITRACER */
                0162 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0163 C     tkey :: tape key (depends on tiles)
                0164       INTEGER tkey
1d74e34c65 Jean*0165 #endif
33e17487ce Dimi*0166
1d74e34c65 Jean*0167 C==   local arrays ==
8377b8ee87 Mart*0168 C     TmixLoc   :: ocean surface/mixed-layer temperature (in K)
5b0abbe6ee Jean*0169       _RL TmixLoc       (1:sNx,1:sNy)
                0170
6571f3ca98 Jean*0171 #ifndef SEAICE_ITD
1d74e34c65 Jean*0172 C     actual ice thickness (with upper and lower limit)
                0173       _RL heffActual          (1:sNx,1:sNy)
                0174 C     actual snow thickness
                0175       _RL hsnowActual         (1:sNx,1:sNy)
6571f3ca98 Jean*0176 #endif
1d74e34c65 Jean*0177 C     actual ice thickness (with lower limit only) Reciprocal
                0178       _RL recip_heffActual    (1:sNx,1:sNy)
                0179
8377b8ee87 Mart*0180 C     AREA_PRE  :: hold sea-ice fraction field before any seaice-thermo update
1d74e34c65 Jean*0181       _RL AREApreTH           (1:sNx,1:sNy)
                0182       _RL HEFFpreTH           (1:sNx,1:sNy)
                0183       _RL HSNWpreTH           (1:sNx,1:sNy)
286983d3d2 Patr*0184 #ifdef SEAICE_ITD
                0185       _RL AREAITDpreTH        (1:sNx,1:sNy,1:nITD)
                0186       _RL HEFFITDpreTH        (1:sNx,1:sNy,1:nITD)
                0187       _RL HSNWITDpreTH        (1:sNx,1:sNy,1:nITD)
                0188       _RL areaFracFactor      (1:sNx,1:sNy,1:nITD)
                0189 #endif
1d74e34c65 Jean*0190
                0191 C     wind speed
                0192       _RL UG                  (1:sNx,1:sNy)
                0193
                0194 C     temporary variables available for the various computations
                0195       _RL tmparr1             (1:sNx,1:sNy)
                0196
f913c5a485 Mart*0197       _RL ticeInMult          (1:sNx,1:sNy,nITD)
                0198       _RL ticeOutMult         (1:sNx,1:sNy,nITD)
                0199       _RL heffActualMult      (1:sNx,1:sNy,nITD)
                0200       _RL hsnowActualMult     (1:sNx,1:sNy,nITD)
286983d3d2 Patr*0201 #ifdef SEAICE_ITD
f913c5a485 Mart*0202       _RL recip_heffActualMult(1:sNx,1:sNy,nITD)
286983d3d2 Patr*0203 #endif
f913c5a485 Mart*0204       _RL a_QbyATMmult_cover  (1:sNx,1:sNy,nITD)
                0205       _RL a_QSWbyATMmult_cover(1:sNx,1:sNy,nITD)
                0206       _RL a_FWbySublimMult    (1:sNx,1:sNy,nITD)
56d13a40ed Mart*0207 #if defined ( ALLOW_SITRACER ) && defined ( SEAICE_GREASE )
f61838dfc1 Torg*0208       _RL greaseLayerThick    (1:sNx,1:sNy)
                0209       _RL d_HEFFbyGREASE      (1:sNx,1:sNy)
56d13a40ed Mart*0210       _RL uRelW               (1:sNx,1:sNy)
                0211       _RL vRelW               (1:sNx,1:sNy)
f61838dfc1 Torg*0212 #endif
286983d3d2 Patr*0213 #ifdef SEAICE_ITD
f913c5a485 Mart*0214       _RL r_QbyATMmult_cover  (1:sNx,1:sNy,nITD)
                0215       _RL r_FWbySublimMult    (1:sNx,1:sNy,nITD)
1080be7801 Jean*0216 C for lateral melt parameterization:
f913c5a485 Mart*0217       _RL latMeltFrac         (1:sNx,1:sNy,nITD)
                0218       _RL latMeltRate         (1:sNx,1:sNy,nITD)
53b2f6dc29 Torg*0219       _RL floeAlpha
                0220       _RL floeDiameter
                0221       _RL floeDiameterMin
                0222       _RL floeDiameterMax
286983d3d2 Patr*0223 #endif
1d74e34c65 Jean*0224
0c0ecd4c7b Jean*0225 C     a_QbyATM_cover :: available heat (in W/m^2) due to the interaction of
2ae913cfea Gael*0226 C             the atmosphere and the ocean surface - for ice covered water
c50ad14e64 Gael*0227 C     a_QbyATM_open  :: same but for open water
65b7f51792 Gael*0228 C     r_QbyATM_cover :: residual of a_QbyATM_cover after freezing/melting processes
c50ad14e64 Gael*0229 C     r_QbyATM_open  :: same but for open water
292695ea58 Gael*0230       _RL a_QbyATM_cover      (1:sNx,1:sNy)
                0231       _RL a_QbyATM_open       (1:sNx,1:sNy)
                0232       _RL r_QbyATM_cover      (1:sNx,1:sNy)
b34884f5be Mart*0233       _RL r_QbyATM_open       (1:sNx,1:sNy)
2ae913cfea Gael*0234 C     a_QSWbyATM_open   - short wave heat flux over ocean in W/m^2
                0235 C     a_QSWbyATM_cover  - short wave heat flux under ice in W/m^2
292695ea58 Gael*0236       _RL a_QSWbyATM_open     (1:sNx,1:sNy)
                0237       _RL a_QSWbyATM_cover    (1:sNx,1:sNy)
286983d3d2 Patr*0238 C     a_QbyOCN :: available heat (in W/m^2) due to the
292695ea58 Gael*0239 C             interaction of the ice pack and the ocean surface
0c0ecd4c7b Jean*0240 C     r_QbyOCN :: residual of a_QbyOCN after freezing/melting
2ae913cfea Gael*0241 C             processes have been accounted for
65b7f51792 Gael*0242       _RL a_QbyOCN            (1:sNx,1:sNy)
                0243       _RL r_QbyOCN            (1:sNx,1:sNy)
292695ea58 Gael*0244
1d74e34c65 Jean*0245 C     The change of mean ice thickness due to turbulent ocean-sea ice heat fluxes
2afe30fba0 Dimi*0246       _RL d_HEFFbyOCNonICE    (1:sNx,1:sNy)
                0247
1d74e34c65 Jean*0248 C     The sum of mean ice thickness increments due to atmospheric fluxes over
                0249 C     the open water fraction and ice-covered fractions of the grid cell
2afe30fba0 Dimi*0250       _RL d_HEFFbyATMonOCN    (1:sNx,1:sNy)
1d74e34c65 Jean*0251 C     The change of mean ice thickness due to flooding by snow
2afe30fba0 Dimi*0252       _RL d_HEFFbyFLOODING    (1:sNx,1:sNy)
2651ba3350 Jean*0253
1d74e34c65 Jean*0254 C     The mean ice thickness increments due to atmospheric fluxes over the open
                0255 C     water fraction and ice-covered fractions of the grid cell, respectively
6b712295de Dimi*0256       _RL d_HEFFbyATMonOCN_open(1:sNx,1:sNy)
                0257       _RL d_HEFFbyATMonOCN_cover(1:sNx,1:sNy)
                0258
2afe30fba0 Dimi*0259       _RL d_HSNWbyATMonSNW    (1:sNx,1:sNy)
                0260       _RL d_HSNWbyOCNonSNW    (1:sNx,1:sNy)
65b7f51792 Gael*0261       _RL d_HSNWbyRAIN        (1:sNx,1:sNy)
2651ba3350 Jean*0262
65b7f51792 Gael*0263       _RL d_HFRWbyRAIN        (1:sNx,1:sNy)
1d74e34c65 Jean*0264
83ad492c2d Jean*0265 C     a_FWbySublim :: fresh water flux implied by latent heat of
b34884f5be Mart*0266 C                     sublimation to atmosphere, same sign convention
                0267 C                     as EVAP (positive upward)
                0268       _RL a_FWbySublim        (1:sNx,1:sNy)
ae36251cae Gael*0269       _RL r_FWbySublim        (1:sNx,1:sNy)
b34884f5be Mart*0270       _RL d_HEFFbySublim      (1:sNx,1:sNy)
                0271       _RL d_HSNWbySublim      (1:sNx,1:sNy)
1f07e6d037 Gael*0272
840c7fba30 Gael*0273 #ifdef SEAICE_CAP_SUBLIM
76f7eb184b Ian *0274 C     The latent heat flux which will sublimate all snow and ice
                0275 C     over one time step
                0276       _RL latentHeatFluxMax   (1:sNx,1:sNy)
f913c5a485 Mart*0277       _RL latentHeatFluxMaxMult(1:sNx,1:sNy,nITD)
2651ba3350 Jean*0278 #endif
d187c22362 Gael*0279
286983d3d2 Patr*0280 #ifdef SEAICE_ITD
8377b8ee87 Mart*0281       _RL d_HEFFbySublim_ITD        (1:sNx,1:sNy,1:nITD)
                0282       _RL d_HSNWbySublim_ITD        (1:sNx,1:sNy,1:nITD)
                0283       _RL d_HEFFbyOCNonICE_ITD      (1:sNx,1:sNy,1:nITD)
                0284       _RL d_HSNWbyATMonSNW_ITD      (1:sNx,1:sNy,1:nITD)
                0285       _RL d_HEFFbyATMonOCN_ITD      (1:sNx,1:sNy,1:nITD)
                0286       _RL d_HEFFbyATMonOCN_cover_ITD(1:sNx,1:sNy,1:nITD)
                0287       _RL d_HEFFbyATMonOCN_open_ITD (1:sNx,1:sNy,1:nITD)
                0288       _RL d_HSNWbyRAIN_ITD          (1:sNx,1:sNy,1:nITD)
                0289       _RL d_HSNWbyOCNonSNW_ITD      (1:sNx,1:sNy,1:nITD)
                0290       _RL d_HEFFbyFLOODING_ITD      (1:sNx,1:sNy,1:nITD)
286983d3d2 Patr*0291 #endif
                0292
33e17487ce Dimi*0293 #ifdef ALLOW_DIAGNOSTICS
1d74e34c65 Jean*0294 C ICE/SNOW stocks tendencies associated with the various melt/freeze processes
8377b8ee87 Mart*0295       _RL d_AREAbyATM   (1:sNx,1:sNy)
                0296       _RL d_AREAbyOCN   (1:sNx,1:sNy)
                0297       _RL d_AREAbyICE   (1:sNx,1:sNy)
1d74e34c65 Jean*0298 C     Helper variables for diagnostics
2afe30fba0 Dimi*0299       _RL DIAGarrayA    (1:sNx,1:sNy)
                0300       _RL DIAGarrayB    (1:sNx,1:sNy)
                0301       _RL DIAGarrayC    (1:sNx,1:sNy)
                0302       _RL DIAGarrayD    (1:sNx,1:sNy)
1d74e34c65 Jean*0303 #endif /* ALLOW_DIAGNOSTICS */
1ed503f8a3 Gael*0304
6509326d8c Gael*0305       _RL SItflux     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0306       _RL SIatmQnt    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0307       _RL SIatmFW     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0308 #ifdef ALLOW_BALANCE_FLUXES
                0309       _RL FWFsiTile(nSx,nSy)
                0310       _RL FWFsiGlob
                0311       _RL HFsiTile(nSx,nSy)
                0312       _RL HFsiGlob
                0313       _RL FWF2HFsiTile(nSx,nSy)
                0314       _RL FWF2HFsiGlob
                0315 #endif
                0316
cafd3818b7 An T*0317 #ifdef ALLOW_SALT_PLUME
8377b8ee87 Mart*0318       _RL localSPfrac         (1:sNx,1:sNy)
cafd3818b7 An T*0319 #ifdef SALT_PLUME_IN_LEADS
8377b8ee87 Mart*0320       _RL leadPlumeFraction   (1:sNx,1:sNy)
                0321       _RL IceGrowthRateInLeads(1:sNx,1:sNy)
cafd3818b7 An T*0322 #endif /* SALT_PLUME_IN_LEADS */
                0323 #endif /* ALLOW_SALT_PLUME */
                0324
2651ba3350 Jean*0325 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0326
                0327 C ===================================================================
                0328 C =================PART 0: constants and initializations=============
                0329 C ===================================================================
aea7db20a6 Gael*0330
0320e25227 Mart*0331       IF ( usingPCoords ) THEN
                0332        kSurface = Nr
                0333        dzSurf   = drF(kSurface)*recip_rhoConst*recip_gravity
33e17487ce Dimi*0334       ELSE
0320e25227 Mart*0335        kSurface = 1
                0336        dzSurf   = drF(kSurface)
33e17487ce Dimi*0337       ENDIF
                0338
bb8e6379cb Mart*0339 C     avoid unnecessary divisions in loops
f5282c5b03 Gael*0340       recip_multDim        = SEAICE_multDim
                0341       recip_multDim        = ONE / recip_multDim
                0342 C     above/below: double/single precision calculation of recip_multDim
                0343 c     recip_multDim        = 1./float(SEAICE_multDim)
a4bc0a0b4c Jean*0344       recip_deltaTtherm = ONE / SEAICE_deltaTtherm
                0345       recip_rhoIce      = ONE / SEAICE_rhoIce
136908bfac Ian *0346
2651ba3350 Jean*0347 C     Cutoff for iceload
0320e25227 Mart*0348       heffTooHeavy = dzSurf * 0.2 _d 0
33e17487ce Dimi*0349 C     RATIO OF SEA ICE DENSITY to SNOW DENSITY
c5a377c43d Gael*0350       ICE2SNOW     = SEAICE_rhoIce/SEAICE_rhoSnow
a4bc0a0b4c Jean*0351       SNOW2ICE     = ONE / ICE2SNOW
85b399b441 Gael*0352
33e17487ce Dimi*0353 C     HEAT OF FUSION OF ICE (J/m^3)
fff6be1885 Mart*0354       QI           = SEAICE_rhoIce*SEAICE_lhFusion
a4bc0a0b4c Jean*0355       recip_QI     = ONE / QI
33e17487ce Dimi*0356 C     HEAT OF FUSION OF SNOW (J/m^3)
fff6be1885 Mart*0357       QS           = SEAICE_rhoSnow*SEAICE_lhFusion
136908bfac Ian *0358
52ff14d141 Ian *0359 C     ICE LATENT HEAT CONSTANT
                0360       lhSublim = SEAICE_lhEvap + SEAICE_lhFusion
                0361
136908bfac Ian *0362 C     regularization constants
                0363       area_reg_sq = SEAICE_area_reg * SEAICE_area_reg
                0364       hice_reg_sq = SEAICE_hice_reg * SEAICE_hice_reg
85b399b441 Gael*0365
2651ba3350 Jean*0366 C conversion factors to go from Q (W/m2) to HEFF (ice meters)
65b7f51792 Gael*0367       convertQ2HI=SEAICE_deltaTtherm/QI
4eb4a54cba Jean*0368       convertHI2Q = ONE/convertQ2HI
2651ba3350 Jean*0369 C conversion factors to go from precip (m/s) unit to HEFF (ice meters)
292695ea58 Gael*0370       convertPRECIP2HI=SEAICE_deltaTtherm*rhoConstFresh/SEAICE_rhoIce
4eb4a54cba Jean*0371       convertHI2PRECIP = ONE/convertPRECIP2HI
4b6d456764 Mart*0372 C     compute parameters for thickness pdf (cheap enough to do it here
                0373 C     and not in seaice_readparms and store in common block)
                0374       denominator = 0. _d 0
                0375       DO IT=1,SEAICE_multDim
                0376        denominator = denominator + IT * SEAICE_pdf(IT)
                0377       ENDDO
                0378       denominator = (2.0 _d 0 * denominator) - 1.0 _d 0
                0379       recip_denominator = 1. _d 0 / denominator
                0380 #ifdef SEAICE_ITD
                0381       areaPDFfac  = 1. _d 0
                0382 #else
                0383       areaPDFfac  = denominator * recip_multDim
                0384 #endif /* SEAICE_ITD */
53b2f6dc29 Torg*0385 #ifdef SEAICE_ITD
1080be7801 Jean*0386 C constants for lateral melt parameterization:
                0387 C following Steele (1992), Equ. 2
53b2f6dc29 Torg*0388       floeAlpha                  = 0.66 _d 0
1080be7801 Jean*0389 C typical mean diameter used in CICE 4.1:
                0390 C (this is currently computed as a function of ice concentration
                0391 C  following a suggestion by Luepkes at al. (2012))
                0392 C      floeDiameter               = 300. _d 0
                0393 C parameters needed for variable floe diameter following Luepkes et al. (2012):
53b2f6dc29 Torg*0394       floeDiameterMin            = 8. _d 0
                0395       floeDiameterMax            = 300. _d 0
                0396 #endif
292695ea58 Gael*0397
33e17487ce Dimi*0398       DO bj=myByLo(myThid),myByHi(myThid)
                0399        DO bi=myBxLo(myThid),myBxHi(myThid)
2651ba3350 Jean*0400
33e17487ce Dimi*0401 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0402         tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
33e17487ce Dimi*0403 #endif /* ALLOW_AUTODIFF_TAMC */
2ae913cfea Gael*0404
56d13a40ed Mart*0405 #if defined ( ALLOW_SITRACER ) && defined ( SEAICE_GREASE )
f61838dfc1 Torg*0406 C time scale of grease ice decline by solidification,
                0407 C with 50% grease ice becoming solid pancake ice within 1 day:
                0408         greaseDecayTime=1.44*86400. _d 0
                0409 C store position of 'grease' in tracer array
                0410         iTrGrease=-1
                0411         DO iTr = 1, SItrNumInUse
56d13a40ed Mart*0412          IF (SItrName(iTr).EQ.'grease') iTrGrease=iTr
f61838dfc1 Torg*0413         ENDDO
                0414 #endif
                0415
2ae913cfea Gael*0416 C array initializations
                0417 C =====================
                0418
8377b8ee87 Mart*0419         DO j=1,sNy
                0420          DO i=1,sNx
                0421           a_QbyATM_cover(i,j)        = 0.0 _d 0
                0422           a_QbyATM_open (i,j)        = 0.0 _d 0
                0423           r_QbyATM_cover(i,j)        = 0.0 _d 0
                0424           r_QbyATM_open (i,j)        = 0.0 _d 0
2651ba3350 Jean*0425
8377b8ee87 Mart*0426           a_QSWbyATM_open (i,j)      = 0.0 _d 0
                0427           a_QSWbyATM_cover(i,j)      = 0.0 _d 0
2651ba3350 Jean*0428
8377b8ee87 Mart*0429           a_QbyOCN (i,j)             = 0.0 _d 0
                0430           r_QbyOCN (i,j)             = 0.0 _d 0
2651ba3350 Jean*0431
381adf77df Gael*0432 #ifdef ALLOW_DIAGNOSTICS
8377b8ee87 Mart*0433           d_AREAbyATM(i,j)           = 0.0 _d 0
                0434           d_AREAbyICE(i,j)           = 0.0 _d 0
                0435           d_AREAbyOCN(i,j)           = 0.0 _d 0
381adf77df Gael*0436 #endif
2651ba3350 Jean*0437
8377b8ee87 Mart*0438           d_HEFFbyOCNonICE(i,j)      = 0.0 _d 0
                0439           d_HEFFbyATMonOCN(i,j)      = 0.0 _d 0
                0440           d_HEFFbyFLOODING(i,j)      = 0.0 _d 0
2651ba3350 Jean*0441
8377b8ee87 Mart*0442           d_HEFFbyATMonOCN_open(i,j) = 0.0 _d 0
                0443           d_HEFFbyATMonOCN_cover(i,j)= 0.0 _d 0
6b712295de Dimi*0444
8377b8ee87 Mart*0445           d_HSNWbyATMonSNW(i,j)      = 0.0 _d 0
                0446           d_HSNWbyOCNonSNW(i,j)      = 0.0 _d 0
                0447           d_HSNWbyRAIN(i,j)          = 0.0 _d 0
                0448           a_FWbySublim(i,j)          = 0.0 _d 0
                0449           r_FWbySublim(i,j)          = 0.0 _d 0
                0450           d_HEFFbySublim(i,j)        = 0.0 _d 0
                0451           d_HSNWbySublim(i,j)        = 0.0 _d 0
840c7fba30 Gael*0452 #ifdef SEAICE_CAP_SUBLIM
8377b8ee87 Mart*0453           latentHeatFluxMax(i,j)     = 0.0 _d 0
840c7fba30 Gael*0454 #endif
8377b8ee87 Mart*0455           d_HFRWbyRAIN(i,j)          = 0.0 _d 0
                0456           tmparr1(i,j)               = 0.0 _d 0
56d13a40ed Mart*0457 #if defined ( ALLOW_SITRACER ) && defined ( SEAICE_GREASE )
8377b8ee87 Mart*0458           greaseLayerThick(i,j)      = 0.0 _d 0
                0459           d_HEFFbyGREASE(i,j)        = 0.0 _d 0
f61838dfc1 Torg*0460 #endif
286983d3d2 Patr*0461           DO IT=1,SEAICE_multDim
8377b8ee87 Mart*0462             ticeInMult(i,j,IT)            = 0.0 _d 0
                0463             ticeOutMult(i,j,IT)           = 0.0 _d 0
                0464             a_QbyATMmult_cover(i,j,IT)    = 0.0 _d 0
                0465             a_QSWbyATMmult_cover(i,j,IT)  = 0.0 _d 0
                0466             a_FWbySublimMult(i,j,IT)      = 0.0 _d 0
a73db480d4 Jean*0467 #ifdef SEAICE_CAP_SUBLIM
8377b8ee87 Mart*0468             latentHeatFluxMaxMult(i,j,IT) = 0.0 _d 0
286983d3d2 Patr*0469 #endif
                0470 #ifdef SEAICE_ITD
8377b8ee87 Mart*0471             d_HEFFbySublim_ITD(i,j,IT)         = 0.0 _d 0
                0472             d_HSNWbySublim_ITD(i,j,IT)         = 0.0 _d 0
                0473             d_HEFFbyOCNonICE_ITD(i,j,IT)       = 0.0 _d 0
                0474             d_HSNWbyATMonSNW_ITD(i,j,IT)       = 0.0 _d 0
                0475             d_HEFFbyATMonOCN_ITD(i,j,IT)       = 0.0 _d 0
                0476             d_HEFFbyATMonOCN_cover_ITD(i,j,IT) = 0.0 _d 0
                0477             d_HEFFbyATMonOCN_open_ITD(i,j,IT)  = 0.0 _d 0
                0478             d_HSNWbyRAIN_ITD(i,j,IT)           = 0.0 _d 0
                0479             d_HSNWbyOCNonSNW_ITD(i,j,IT)       = 0.0 _d 0
                0480             d_HEFFbyFLOODING_ITD(i,j,IT)       = 0.0 _d 0
                0481             r_QbyATMmult_cover(i,j,IT)         = 0.0 _d 0
                0482             r_FWbySublimMult(i,j,IT)           = 0.0 _d 0
1080be7801 Jean*0483 C for lateral melt parameterization:
8377b8ee87 Mart*0484             latMeltFrac(i,j,IT)                = 0.0 _d 0
                0485             latMeltRate(i,j,IT)                = 0.0 _d 0
a73db480d4 Jean*0486 #endif
f5282c5b03 Gael*0487           ENDDO
33e17487ce Dimi*0488          ENDDO
                0489         ENDDO
                0490
2651ba3350 Jean*0491 C =====================================================================
                0492 C ===========PART 1: treat pathological cases (post advdiff)===========
                0493 C =====================================================================
aea7db20a6 Gael*0494
c6b168144e Jean*0495 C     This part has been mostly moved to S/R seaice_reg_ridge, which is
1cf549c217 Mart*0496 C     called before S/R seaice_growth
2651ba3350 Jean*0497
1cf549c217 Mart*0498 C     store regularized values of heff, hsnow, area at the onset of thermo.
8377b8ee87 Mart*0499         DO j=1,sNy
                0500          DO i=1,sNx
                0501           HEFFpreTH(i,j) = HEFF(i,j,bi,bj)
                0502           HSNWpreTH(i,j) = HSNOW(i,j,bi,bj)
                0503           AREApreTH(i,j) = AREA(i,j,bi,bj)
4213eb5769 Gael*0504 #ifdef ALLOW_DIAGNOSTICS
8377b8ee87 Mart*0505           DIAGarrayB(i,j) = AREA(i,j,bi,bj)
                0506           DIAGarrayC(i,j) = HEFF(i,j,bi,bj)
                0507           DIAGarrayD(i,j) = HSNOW(i,j,bi,bj)
4213eb5769 Gael*0508 #endif
f50f58ec54 Gael*0509 #ifdef ALLOW_SITRACER
8377b8ee87 Mart*0510           SItrHEFF(i,j,bi,bj,1)=HEFF(i,j,bi,bj)
                0511           SItrAREA(i,j,bi,bj,2)=AREA(i,j,bi,bj)
f50f58ec54 Gael*0512 #endif
581175eaf0 Gael*0513          ENDDO
                0514         ENDDO
286983d3d2 Patr*0515 #ifdef SEAICE_ITD
f913c5a485 Mart*0516         DO IT=1,SEAICE_multDim
8377b8ee87 Mart*0517          DO j=1,sNy
                0518           DO i=1,sNx
                0519            HEFFITDpreTH(i,j,IT)=HEFFITD(i,j,IT,bi,bj)
                0520            HSNWITDpreTH(i,j,IT)=HSNOWITD(i,j,IT,bi,bj)
                0521            AREAITDpreTH(i,j,IT)=AREAITD(i,j,IT,bi,bj)
286983d3d2 Patr*0522
1cf549c217 Mart*0523 C     keep track of areal and volume fraction of each ITD category
8377b8ee87 Mart*0524            IF (AREA(i,j,bi,bj) .GT. ZERO) THEN
                0525             areaFracFactor(i,j,IT)=AREAITD(i,j,IT,bi,bj)/AREA(i,j,bi,bj)
114c791332 Jean*0526            ELSE
                0527 C           if there is no ice, potential growth starts in 1st category
                0528             IF (IT .EQ. 1) THEN
8377b8ee87 Mart*0529              areaFracFactor(i,j,IT)=ONE
114c791332 Jean*0530             ELSE
8377b8ee87 Mart*0531              areaFracFactor(i,j,IT)=ZERO
114c791332 Jean*0532             ENDIF
                0533            ENDIF
286983d3d2 Patr*0534           ENDDO
                0535          ENDDO
                0536         ENDDO
                0537 #ifdef ALLOW_SITRACER
1cf549c217 Mart*0538 C     prepare SItrHEFF to be computed as cumulative sum
286983d3d2 Patr*0539         DO iTr=2,5
8377b8ee87 Mart*0540          DO j=1,sNy
                0541           DO i=1,sNx
                0542            SItrHEFF(i,j,bi,bj,iTr)=ZERO
286983d3d2 Patr*0543           ENDDO
                0544          ENDDO
                0545         ENDDO
1cf549c217 Mart*0546 C     prepare SItrAREA to be computed as cumulative sum
8377b8ee87 Mart*0547         DO j=1,sNy
                0548          DO i=1,sNx
                0549           SItrAREA(i,j,bi,bj,3)=ZERO
286983d3d2 Patr*0550          ENDDO
                0551         ENDDO
                0552 #endif
                0553 #endif /* SEAICE_ITD */
33e17487ce Dimi*0554
4213eb5769 Gael*0555 #ifdef ALLOW_DIAGNOSTICS
4bc8f4264e Mart*0556         IF ( useDiagnostics ) THEN
                0557          CALL DIAGNOSTICS_FILL(DIAGarrayB,'SIareaPT',0,1,3,bi,bj,myThid)
                0558          CALL DIAGNOSTICS_FILL(DIAGarrayC,'SIheffPT',0,1,3,bi,bj,myThid)
                0559          CALL DIAGNOSTICS_FILL(DIAGarrayD,'SIhsnoPT',0,1,3,bi,bj,myThid)
3721cfe5e4 Gael*0560 #ifdef ALLOW_SITRACER
38cfb58d85 Gael*0561          DO iTr = 1, SItrNumInUse
4bc8f4264e Mart*0562           WRITE(diagName,'(A4,I2.2,A2)') 'SItr',iTr,'PT'
4eb4a54cba Jean*0563           IF (SItrMate(iTr).EQ.'HEFF') THEN
4bc8f4264e Mart*0564            CALL DIAGNOSTICS_FRACT_FILL(
                0565      I          SItracer(1-OLx,1-OLy,bi,bj,iTr),HEFF(1-OLx,1-OLy,bi,bj),
                0566      I          ONE, 1, diagName,0,1,2,bi,bj,myThid )
4eb4a54cba Jean*0567           ELSE
4bc8f4264e Mart*0568            CALL DIAGNOSTICS_FRACT_FILL(
                0569      I          SItracer(1-OLx,1-OLy,bi,bj,iTr),AREA(1-OLx,1-OLy,bi,bj),
                0570      I          ONE, 1, diagName,0,1,2,bi,bj,myThid )
4eb4a54cba Jean*0571           ENDIF
4bc8f4264e Mart*0572          ENDDO
a73db480d4 Jean*0573 #endif /* ALLOW_SITRACER */
4bc8f4264e Mart*0574         ENDIF
a73db480d4 Jean*0575 #endif /* ALLOW_DIAGNOSTICS */
4213eb5769 Gael*0576
8377b8ee87 Mart*0577 #if (defined ALLOW_AUTODIFF && defined SEAICE_MODIFY_GROWTH_ADJ)
3a3bf6419a Gael*0578 Cgf no additional dependency of air-sea fluxes to ice
4bc8f4264e Mart*0579         IF ( SEAICEadjMODE.GE.1 ) THEN
8377b8ee87 Mart*0580          DO j=1,sNy
                0581           DO i=1,sNx
                0582            HEFFpreTH(i,j) = 0. _d 0
                0583            HSNWpreTH(i,j) = 0. _d 0
                0584            AREApreTH(i,j) = 0. _d 0
4bc8f4264e Mart*0585           ENDDO
3a3bf6419a Gael*0586          ENDDO
286983d3d2 Patr*0587 #ifdef SEAICE_ITD
f913c5a485 Mart*0588          DO IT=1,SEAICE_multDim
8377b8ee87 Mart*0589           DO j=1,sNy
                0590            DO i=1,sNx
                0591             HEFFITDpreTH(i,j,IT) = 0. _d 0
                0592             HSNWITDpreTH(i,j,IT) = 0. _d 0
                0593             AREAITDpreTH(i,j,IT) = 0. _d 0
286983d3d2 Patr*0594            ENDDO
                0595           ENDDO
                0596          ENDDO
                0597 #endif
4bc8f4264e Mart*0598         ENDIF
3a3bf6419a Gael*0599 #endif
581175eaf0 Gael*0600
1cf549c217 Mart*0601 C     COMPUTE ACTUAL ICE/SNOW THICKNESS; USE MIN/MAX VALUES
                0602 C     TO REGULARIZE SEAICE_SOLVE4TEMP/d_AREA COMPUTATIONS
2651ba3350 Jean*0603
581175eaf0 Gael*0604 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0605 CADJ STORE AREApreTH = comlev1_bibj, key = tkey, byte = isbyte
                0606 CADJ STORE HEFFpreTH = comlev1_bibj, key = tkey, byte = isbyte
                0607 CADJ STORE HSNWpreTH = comlev1_bibj, key = tkey, byte = isbyte
33e17487ce Dimi*0608 #endif /* ALLOW_AUTODIFF_TAMC */
286983d3d2 Patr*0609 #ifdef SEAICE_ITD
f913c5a485 Mart*0610         DO IT=1,SEAICE_multDim
8377b8ee87 Mart*0611          DO j=1,sNy
                0612           DO i=1,sNx
                0613            IF (HEFFITDpreTH(i,j,IT) .GT. ZERO) THEN
1cf549c217 Mart*0614 C     regularize AREA with SEAICE_area_reg
8377b8ee87 Mart*0615             tmpscal1 = SQRT( AREAITDpreTH(i,j,IT)*AREAITDpreTH(i,j,IT)
                0616      &                     + area_reg_sq )
1cf549c217 Mart*0617 C     heffActual calculated with the regularized AREA
8377b8ee87 Mart*0618             tmpscal2 = HEFFITDpreTH(i,j,IT) / tmpscal1
1cf549c217 Mart*0619 C     regularize heffActual with SEAICE_hice_reg (add lower bound)
8377b8ee87 Mart*0620             heffActualMult(i,j,IT) = SQRT(tmpscal2 * tmpscal2
1cf549c217 Mart*0621      &                                    + hice_reg_sq)
                0622 C     hsnowActual calculated with the regularized AREA
8377b8ee87 Mart*0623             hsnowActualMult(i,j,IT) = HSNWITDpreTH(i,j,IT) / tmpscal1
1cf549c217 Mart*0624 C     regularize the inverse of heffActual by hice_reg
8377b8ee87 Mart*0625             recip_heffActualMult(i,j,IT)  = AREAITDpreTH(i,j,IT) /
                0626      &           SQRT(HEFFITDpreTH(i,j,IT) * HEFFITDpreTH(i,j,IT)
1cf549c217 Mart*0627      &           + hice_reg_sq)
                0628 C     Do not regularize when HEFFpreTH = 0
                0629            ELSE
8377b8ee87 Mart*0630             heffActualMult(i,j,IT) = ZERO
                0631             hsnowActualMult(i,j,IT) = ZERO
                0632             recip_heffActualMult(i,j,IT)  = ZERO
1cf549c217 Mart*0633            ENDIF
                0634           ENDDO
                0635          ENDDO
                0636         ENDDO
                0637 #else /* ndef SEAICE_ITD */
8377b8ee87 Mart*0638         DO j=1,sNy
                0639          DO i=1,sNx
                0640           IF (HEFFpreTH(i,j) .GT. ZERO) THEN
1d74e34c65 Jean*0641 Cif        regularize AREA with SEAICE_area_reg
8377b8ee87 Mart*0642            tmpscal1 = SQRT(AREApreTH(i,j)* AREApreTH(i,j) + area_reg_sq)
1d74e34c65 Jean*0643 Cif        heffActual calculated with the regularized AREA
8377b8ee87 Mart*0644            tmpscal2 = HEFFpreTH(i,j) / tmpscal1
1d74e34c65 Jean*0645 Cif        regularize heffActual with SEAICE_hice_reg (add lower bound)
8377b8ee87 Mart*0646            heffActual(i,j) = SQRT(tmpscal2 * tmpscal2 + hice_reg_sq)
1d74e34c65 Jean*0647 Cif        hsnowActual calculated with the regularized AREA
8377b8ee87 Mart*0648            hsnowActual(i,j) = HSNWpreTH(i,j) / tmpscal1
1d74e34c65 Jean*0649 Cif        regularize the inverse of heffActual by hice_reg
8377b8ee87 Mart*0650            recip_heffActual(i,j)  = AREApreTH(i,j) /
                0651      &                 SQRT(HEFFpreTH(i,j)*HEFFpreTH(i,j) + hice_reg_sq)
1d74e34c65 Jean*0652 Cif       Do not regularize when HEFFpreTH = 0
4bc8f4264e Mart*0653           ELSE
8377b8ee87 Mart*0654            heffActual(i,j) = ZERO
                0655            hsnowActual(i,j) = ZERO
                0656            recip_heffActual(i,j)  = ZERO
4bc8f4264e Mart*0657           ENDIF
33e17487ce Dimi*0658          ENDDO
                0659         ENDDO
1cf549c217 Mart*0660 #endif /* SEAICE_ITD */
33e17487ce Dimi*0661
8377b8ee87 Mart*0662 #if (defined ALLOW_AUTODIFF && defined SEAICE_MODIFY_GROWTH_ADJ)
4bc8f4264e Mart*0663         CALL ZERO_ADJ_1D( sNx*sNy, heffActual, myThid)
                0664         CALL ZERO_ADJ_1D( sNx*sNy, hsnowActual, myThid)
                0665         CALL ZERO_ADJ_1D( sNx*sNy, recip_heffActual, myThid)
b2891f8c1a Gael*0666 #endif
581175eaf0 Gael*0667
840c7fba30 Gael*0668 #ifdef SEAICE_CAP_SUBLIM
1cf549c217 Mart*0669 C     COMPUTE MAXIMUM LATENT HEAT FLUXES FOR THE CURRENT ICE
                0670 C     AND SNOW THICKNESS
                0671 C     The latent heat flux over the sea ice which
                0672 C     will sublimate all of the snow and ice over one time
                0673 C     step (W/m^2)
286983d3d2 Patr*0674 #ifdef SEAICE_ITD
f913c5a485 Mart*0675         DO IT=1,SEAICE_multDim
8377b8ee87 Mart*0676          DO j=1,sNy
                0677           DO i=1,sNx
                0678            IF (HEFFITDpreTH(i,j,IT) .GT. ZERO) THEN
                0679             latentHeatFluxMaxMult(i,j,IT) = lhSublim*recip_deltaTtherm *
                0680      &           (HEFFITDpreTH(i,j,IT)*SEAICE_rhoIce +
                0681      &            HSNWITDpreTH(i,j,IT)*SEAICE_rhoSnow)
                0682      &           /AREAITDpreTH(i,j,IT)
1cf549c217 Mart*0683            ELSE
8377b8ee87 Mart*0684             latentHeatFluxMaxMult(i,j,IT) = ZERO
1cf549c217 Mart*0685            ENDIF
                0686           ENDDO
                0687          ENDDO
                0688         ENDDO
                0689 #else /* ndef SEAICE_ITD */
8377b8ee87 Mart*0690         DO j=1,sNy
                0691          DO i=1,sNx
                0692           IF (HEFFpreTH(i,j) .GT. ZERO) THEN
                0693            latentHeatFluxMax(i,j) = lhSublim * recip_deltaTtherm *
                0694      &          (HEFFpreTH(i,j) * SEAICE_rhoIce +
                0695      &           HSNWpreTH(i,j) * SEAICE_rhoSnow)/AREApreTH(i,j)
4bc8f4264e Mart*0696           ELSE
8377b8ee87 Mart*0697            latentHeatFluxMax(i,j) = ZERO
4bc8f4264e Mart*0698           ENDIF
52ff14d141 Ian *0699          ENDDO
                0700         ENDDO
1cf549c217 Mart*0701 #endif /* SEAICE_ITD */
a73db480d4 Jean*0702 #endif /* SEAICE_CAP_SUBLIM */
2afe30fba0 Dimi*0703
2651ba3350 Jean*0704 C ===================================================================
                0705 C ================PART 2: determine heat fluxes/stocks===============
                0706 C ===================================================================
aea7db20a6 Gael*0707
2ae913cfea Gael*0708 C determine available heat due to the atmosphere -- for open water
                0709 C ================================================================
                0710
5b0abbe6ee Jean*0711         DO j=1,sNy
                0712          DO i=1,sNx
840c7fba30 Gael*0713 C ocean surface/mixed layer temperature
5b0abbe6ee Jean*0714           TmixLoc(i,j) = theta(i,j,kSurface,bi,bj)+celsius2K
2ae913cfea Gael*0715 C wind speed from exf
8377b8ee87 Mart*0716           UG(i,j) = MAX(SEAICE_EPS,wspeed(i,j,bi,bj))
33e17487ce Dimi*0717          ENDDO
                0718         ENDDO
                0719
3a3bf6419a Gael*0720 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0721 CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj, key = tkey,byte=isbyte
                0722 CADJ STORE qsw(:,:,bi,bj)  = comlev1_bibj, key = tkey,byte=isbyte
                0723 cCADJ STORE UG = comlev1_bibj, key = tkey,byte=isbyte
                0724 cCADJ STORE TmixLoc = comlev1_bibj, key = tkey,byte=isbyte
3a3bf6419a Gael*0725 #endif /* ALLOW_AUTODIFF_TAMC */
                0726
2ae913cfea Gael*0727         CALL SEAICE_BUDGET_OCEAN(
                0728      I       UG,
5b0abbe6ee Jean*0729      I       TmixLoc,
2ae913cfea Gael*0730      O       a_QbyATM_open, a_QSWbyATM_open,
                0731      I       bi, bj, myTime, myIter, myThid )
                0732
                0733 C determine available heat due to the atmosphere -- for ice covered water
                0734 C =======================================================================
33e17487ce Dimi*0735
358649780a Gael*0736         IF (useRelativeWind.AND.useAtmWind) THEN
33e17487ce Dimi*0737 C     Compute relative wind speed over sea ice.
8377b8ee87 Mart*0738          DO j=1,sNy
                0739           DO i=1,sNx
33e17487ce Dimi*0740            SPEED_SQ =
8377b8ee87 Mart*0741      &          (uWind(i,j,bi,bj)
33e17487ce Dimi*0742      &          -0.5 _d 0*(uice(i,j,bi,bj)+uice(i+1,j,bi,bj)))**2
8377b8ee87 Mart*0743      &          +(vWind(i,j,bi,bj)
33e17487ce Dimi*0744      &          -0.5 _d 0*(vice(i,j,bi,bj)+vice(i,j+1,bi,bj)))**2
                0745            IF ( SPEED_SQ .LE. SEAICE_EPS_SQ ) THEN
8377b8ee87 Mart*0746              UG(i,j)=SEAICE_EPS
33e17487ce Dimi*0747            ELSE
8377b8ee87 Mart*0748              UG(i,j)=SQRT(SPEED_SQ)
33e17487ce Dimi*0749            ENDIF
                0750           ENDDO
                0751          ENDDO
                0752         ENDIF
2ae913cfea Gael*0753
33e17487ce Dimi*0754 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0755 CADJ STORE hsnowActual = comlev1_bibj, key = tkey, byte = isbyte
                0756 CADJ STORE heffActual  = comlev1_bibj, key = tkey, byte = isbyte
                0757 CADJ STORE UG          = comlev1_bibj, key = tkey, byte = isbyte
66d21a8387 Jean*0758 CADJ STORE tices(:,:,:,bi,bj)
edb6656069 Mart*0759 CADJ &     = comlev1_bibj, key = tkey, byte = isbyte
4fd8e94be0 Gael*0760 CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
edb6656069 Mart*0761 CADJ &                       key = tkey, byte = isbyte
057ebb1030 Mart*0762 #endif /* ALLOW_AUTODIFF_TAMC */
                0763
f5282c5b03 Gael*0764 C--   Start loop over multi-categories
cbd0ee24a8 Mart*0765         DO IT=1,SEAICE_multDim
8377b8ee87 Mart*0766          DO j=1,sNy
                0767           DO i=1,sNx
                0768            ticeInMult(i,j,IT)  = TICES(i,j,IT,bi,bj)
                0769            ticeOutMult(i,j,IT) = TICES(i,j,IT,bi,bj)
                0770            TICES(i,j,IT,bi,bj) = ZERO
286983d3d2 Patr*0771           ENDDO
                0772          ENDDO
cbd0ee24a8 Mart*0773 #ifndef SEAICE_ITD
                0774 C     for SEAICE_ITD heffActualMult and latentHeatFluxMaxMult have been
4b6d456764 Mart*0775 C     calculated above (instead of heffActual and latentHeatFluxMax)
cbd0ee24a8 Mart*0776 C--   assume homogeneous distribution between 0 and 2 x heffActual
4b6d456764 Mart*0777          pFac = (2.0 _d 0*IT - 1.0 _d 0)*recip_denominator
74c037b5fb Mart*0778          pFacSnow = 1. _d 0
                0779          IF ( SEAICE_useMultDimSnow ) pFacSnow=pFac
8377b8ee87 Mart*0780          DO j=1,sNy
                0781           DO i=1,sNx
                0782            heffActualMult(i,j,IT)        = heffActual(i,j)*pFac
                0783            hsnowActualMult(i,j,IT)       = hsnowActual(i,j)*pFacSnow
840c7fba30 Gael*0784 #ifdef SEAICE_CAP_SUBLIM
8377b8ee87 Mart*0785            latentHeatFluxMaxMult(i,j,IT) = latentHeatFluxMax(i,j)*pFac
52ff14d141 Ian *0786 #endif
33e17487ce Dimi*0787           ENDDO
                0788          ENDDO
cbd0ee24a8 Mart*0789 #endif /* ndef SEAICE_ITD */
f5282c5b03 Gael*0790         ENDDO
                0791
                0792 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0793 CADJ STORE heffActualMult = comlev1_bibj, key = tkey, byte = isbyte
                0794 CADJ STORE hsnowActualMult= comlev1_bibj, key = tkey, byte = isbyte
                0795 CADJ STORE ticeInMult     = comlev1_bibj, key = tkey, byte = isbyte
f5282c5b03 Gael*0796 # ifdef SEAICE_CAP_SUBLIM
4eb4a54cba Jean*0797 CADJ STORE latentHeatFluxMaxMult
edb6656069 Mart*0798 CADJ &     = comlev1_bibj, key = tkey, byte = isbyte
f5282c5b03 Gael*0799 # endif
4eb4a54cba Jean*0800 CADJ STORE a_QbyATMmult_cover   =
edb6656069 Mart*0801 CADJ &     comlev1_bibj, key = tkey, byte = isbyte
4eb4a54cba Jean*0802 CADJ STORE a_QSWbyATMmult_cover =
edb6656069 Mart*0803 CADJ &     comlev1_bibj, key = tkey, byte = isbyte
4eb4a54cba Jean*0804 CADJ STORE a_FWbySublimMult     =
edb6656069 Mart*0805 CADJ &     comlev1_bibj, key = tkey, byte = isbyte
f5282c5b03 Gael*0806 #endif /* ALLOW_AUTODIFF_TAMC */
                0807
286983d3d2 Patr*0808         DO IT=1,SEAICE_multDim
9637aec598 Jean*0809          CALL SEAICE_SOLVE4TEMP(
286983d3d2 Patr*0810      I        UG, heffActualMult(1,1,IT), hsnowActualMult(1,1,IT),
840c7fba30 Gael*0811 #ifdef SEAICE_CAP_SUBLIM
286983d3d2 Patr*0812      I        latentHeatFluxMaxMult(1,1,IT),
52ff14d141 Ian *0813 #endif
4dd39c50d9 Mart*0814      I        ticeInMult(1,1,IT),
                0815      O        ticeOutMult(1,1,IT),
286983d3d2 Patr*0816      O        a_QbyATMmult_cover(1,1,IT),
                0817      O        a_QSWbyATMmult_cover(1,1,IT),
                0818      O        a_FWbySublimMult(1,1,IT),
33e17487ce Dimi*0819      I        bi, bj, myTime, myIter, myThid )
f5282c5b03 Gael*0820         ENDDO
                0821
                0822 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0823 CADJ STORE heffActualMult = comlev1_bibj, key = tkey, byte = isbyte
                0824 CADJ STORE hsnowActualMult= comlev1_bibj, key = tkey, byte = isbyte
                0825 CADJ STORE ticeOutMult    = comlev1_bibj, key = tkey, byte = isbyte
f5282c5b03 Gael*0826 # ifdef SEAICE_CAP_SUBLIM
4eb4a54cba Jean*0827 CADJ STORE latentHeatFluxMaxMult
edb6656069 Mart*0828 CADJ &     = comlev1_bibj, key = tkey, byte = isbyte
f5282c5b03 Gael*0829 # endif
4eb4a54cba Jean*0830 CADJ STORE a_QbyATMmult_cover   =
edb6656069 Mart*0831 CADJ &     comlev1_bibj, key = tkey, byte = isbyte
4eb4a54cba Jean*0832 CADJ STORE a_QSWbyATMmult_cover =
edb6656069 Mart*0833 CADJ &     comlev1_bibj, key = tkey, byte = isbyte
4eb4a54cba Jean*0834 CADJ STORE a_FWbySublimMult     =
edb6656069 Mart*0835 CADJ &     comlev1_bibj, key = tkey, byte = isbyte
f5282c5b03 Gael*0836 #endif /* ALLOW_AUTODIFF_TAMC */
4eb4a54cba Jean*0837
286983d3d2 Patr*0838         DO IT=1,SEAICE_multDim
8377b8ee87 Mart*0839          DO j=1,sNy
                0840           DO i=1,sNx
b69fbfd195 Mart*0841 C     update TICES
                0842 CMLC     and tIce, if required later on (currently not the case)
                0843 CML#ifdef SEAICE_ITD
                0844 CMLC     calculate area weighted mean
                0845 CMLC     (although the ice temperature relates to its energy content
                0846 CMLC      and hence should be averaged weighted by ice volume,
                0847 CMLC      the temperature here is a result of the fluxes through the ice surface
                0848 CMLC      computed individually for each single category in SEAICE_SOLVE4TEMP
                0849 CMLC      and hence is averaged area weighted [areaFracFactor])
                0850 CML           tIce(I,J,bi,bj) = tIce(I,J,bi,bj)
                0851 CML     &          +  ticeOutMult(I,J,IT)*areaFracFactor(I,J,IT)
                0852 CML#else
                0853 CML           tIce(I,J,bi,bj) = tIce(I,J,bi,bj)
4b6d456764 Mart*0854 CML     &          +  ticeOutMult(I,J,IT)*SEAICE_PDF(IT)
b69fbfd195 Mart*0855 CML#endif
8377b8ee87 Mart*0856            TICES(i,j,IT,bi,bj) = ticeOutMult(i,j,IT)
2ae913cfea Gael*0857 C     average over categories
286983d3d2 Patr*0858 #ifdef SEAICE_ITD
                0859 C     calculate area weighted mean
                0860 C     (fluxes are per unit (ice surface) area and are thus area weighted)
8377b8ee87 Mart*0861            a_QbyATM_cover   (i,j) = a_QbyATM_cover(i,j)
                0862      &          + a_QbyATMmult_cover(i,j,IT)*areaFracFactor(i,j,IT)
                0863            a_QSWbyATM_cover (i,j) = a_QSWbyATM_cover(i,j)
                0864      &          + a_QSWbyATMmult_cover(i,j,IT)*areaFracFactor(i,j,IT)
                0865            a_FWbySublim     (i,j) = a_FWbySublim(i,j)
                0866      &          + a_FWbySublimMult(i,j,IT)*areaFracFactor(i,j,IT)
286983d3d2 Patr*0867 #else
8377b8ee87 Mart*0868            a_QbyATM_cover   (i,j) = a_QbyATM_cover(i,j)
                0869      &          + a_QbyATMmult_cover(i,j,IT)*SEAICE_PDF(IT)
                0870            a_QSWbyATM_cover (i,j) = a_QSWbyATM_cover(i,j)
                0871      &          + a_QSWbyATMmult_cover(i,j,IT)*SEAICE_PDF(IT)
                0872            a_FWbySublim     (i,j) = a_FWbySublim(i,j)
                0873      &          + a_FWbySublimMult(i,j,IT)*SEAICE_PDF(IT)
286983d3d2 Patr*0874 #endif
33e17487ce Dimi*0875           ENDDO
                0876          ENDDO
                0877         ENDDO
                0878
840c7fba30 Gael*0879 #ifdef SEAICE_CAP_SUBLIM
                0880 # ifdef ALLOW_DIAGNOSTICS
8377b8ee87 Mart*0881         DO j=1,sNy
                0882          DO i=1,sNx
1d74e34c65 Jean*0883 C          The actual latent heat flux realized by SOLVE4TEMP
8377b8ee87 Mart*0884            DIAGarrayA(i,j) = a_FWbySublim(i,j) * lhSublim
52ff14d141 Ian *0885          ENDDO
                0886         ENDDO
1d74e34c65 Jean*0887 Cif     The actual vs. maximum latent heat flux
52ff14d141 Ian *0888         IF ( useDiagnostics ) THEN
                0889           CALL DIAGNOSTICS_FILL(DIAGarrayA,
                0890      &     'SIactLHF',0,1,3,bi,bj,myThid)
                0891           CALL DIAGNOSTICS_FILL(latentHeatFluxMax,
                0892      &     'SImaxLHF',0,1,3,bi,bj,myThid)
                0893         ENDIF
a73db480d4 Jean*0894 # endif /* ALLOW_DIAGNOSTICS */
                0895 #endif /* SEAICE_CAP_SUBLIM */
52ff14d141 Ian *0896
4fd8e94be0 Gael*0897 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0898 CADJ STORE AREApreTH       = comlev1_bibj, key = tkey, byte = isbyte
                0899 CADJ STORE a_QbyATM_cover  = comlev1_bibj, key = tkey, byte = isbyte
                0900 CADJ STORE a_QSWbyATM_cover= comlev1_bibj, key = tkey, byte = isbyte
                0901 CADJ STORE a_QbyATM_open   = comlev1_bibj, key = tkey, byte = isbyte
                0902 CADJ STORE a_QSWbyATM_open = comlev1_bibj, key = tkey, byte = isbyte
                0903 CADJ STORE a_FWbySublim    = comlev1_bibj, key = tkey, byte = isbyte
4fd8e94be0 Gael*0904 #endif /* ALLOW_AUTODIFF_TAMC */
                0905
2651ba3350 Jean*0906 C switch heat fluxes from W/m2 to 'effective' ice meters
286983d3d2 Patr*0907 #ifdef SEAICE_ITD
f913c5a485 Mart*0908         DO IT=1,SEAICE_multDim
8377b8ee87 Mart*0909          DO j=1,sNy
                0910           DO i=1,sNx
                0911            a_QbyATMmult_cover(i,j,IT)   = a_QbyATMmult_cover(i,j,IT)
                0912      &          * convertQ2HI * AREAITDpreTH(i,j,IT)
                0913            a_QSWbyATMmult_cover(i,j,IT) = a_QSWbyATMmult_cover(i,j,IT)
                0914      &          * convertQ2HI * AREAITDpreTH(i,j,IT)
286983d3d2 Patr*0915 C and initialize r_QbyATMmult_cover
8377b8ee87 Mart*0916            r_QbyATMmult_cover(i,j,IT)=a_QbyATMmult_cover(i,j,IT)
286983d3d2 Patr*0917 C     Convert fresh water flux by sublimation to 'effective' ice meters.
                0918 C     Negative sublimation is resublimation and will be added as snow.
                0919 #ifdef SEAICE_DISABLE_SUBLIM
8377b8ee87 Mart*0920            a_FWbySublimMult(i,j,IT) = ZERO
286983d3d2 Patr*0921 #endif
8377b8ee87 Mart*0922            a_FWbySublimMult(i,j,IT) = SEAICE_deltaTtherm*recip_rhoIce
                0923      &            * a_FWbySublimMult(i,j,IT)*AREAITDpreTH(i,j,IT)
                0924            r_FWbySublimMult(i,j,IT)=a_FWbySublimMult(i,j,IT)
114c791332 Jean*0925           ENDDO
286983d3d2 Patr*0926          ENDDO
                0927         ENDDO
8377b8ee87 Mart*0928         DO j=1,sNy
                0929          DO i=1,sNx
                0930           a_QbyATM_open(i,j)    = a_QbyATM_open(i,j)
                0931      &         * convertQ2HI * ( ONE - AREApreTH(i,j) )
                0932           a_QSWbyATM_open(i,j)  = a_QSWbyATM_open(i,j)
                0933      &         * convertQ2HI * ( ONE - AREApreTH(i,j) )
286983d3d2 Patr*0934 C and initialize r_QbyATM_open
8377b8ee87 Mart*0935           r_QbyATM_open(i,j)=a_QbyATM_open(i,j)
286983d3d2 Patr*0936          ENDDO
                0937         ENDDO
                0938 #else /* SEAICE_ITD */
8377b8ee87 Mart*0939         DO j=1,sNy
                0940          DO i=1,sNx
                0941           a_QbyATM_cover(i,j)   = a_QbyATM_cover(i,j)
                0942      &         * convertQ2HI * AREApreTH(i,j)
                0943           a_QSWbyATM_cover(i,j) = a_QSWbyATM_cover(i,j)
                0944      &         * convertQ2HI * AREApreTH(i,j)
                0945           a_QbyATM_open(i,j)    = a_QbyATM_open(i,j)
                0946      &         * convertQ2HI * ( ONE - AREApreTH(i,j) )
                0947           a_QSWbyATM_open(i,j)  = a_QSWbyATM_open(i,j)
                0948      &         * convertQ2HI * ( ONE - AREApreTH(i,j) )
2651ba3350 Jean*0949 C and initialize r_QbyATM_cover/r_QbyATM_open
8377b8ee87 Mart*0950           r_QbyATM_cover(i,j)=a_QbyATM_cover(i,j)
                0951           r_QbyATM_open(i,j)=a_QbyATM_open(i,j)
83ad492c2d Jean*0952 C     Convert fresh water flux by sublimation to 'effective' ice meters.
b34884f5be Mart*0953 C     Negative sublimation is resublimation and will be added as snow.
840c7fba30 Gael*0954 #ifdef SEAICE_DISABLE_SUBLIM
1d74e34c65 Jean*0955 Cgf just for those who may need to omit this term to reproduce old results
8377b8ee87 Mart*0956           a_FWbySublim(i,j) = ZERO
933b1d4757 Jean*0957 #endif /* SEAICE_DISABLE_SUBLIM */
8377b8ee87 Mart*0958           a_FWbySublim(i,j) = SEAICE_deltaTtherm*recip_rhoIce
                0959      &           * a_FWbySublim(i,j)*AREApreTH(i,j)
                0960           r_FWbySublim(i,j)=a_FWbySublim(i,j)
bea7d9d588 Gael*0961          ENDDO
                0962         ENDDO
286983d3d2 Patr*0963 #endif /* SEAICE_ITD */
2ae913cfea Gael*0964
4fd8e94be0 Gael*0965 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0966 CADJ STORE AREApreTH       = comlev1_bibj, key = tkey, byte = isbyte
                0967 CADJ STORE a_QbyATM_cover  = comlev1_bibj, key = tkey, byte = isbyte
                0968 CADJ STORE a_QSWbyATM_cover= comlev1_bibj, key = tkey, byte = isbyte
                0969 CADJ STORE a_QbyATM_open   = comlev1_bibj, key = tkey, byte = isbyte
                0970 CADJ STORE a_QSWbyATM_open = comlev1_bibj, key = tkey, byte = isbyte
                0971 CADJ STORE a_FWbySublim    = comlev1_bibj, key = tkey, byte = isbyte
                0972 CADJ STORE r_QbyATM_cover  = comlev1_bibj, key = tkey, byte = isbyte
                0973 CADJ STORE r_QbyATM_open   = comlev1_bibj, key = tkey, byte = isbyte
                0974 CADJ STORE r_FWbySublim    = comlev1_bibj, key = tkey, byte = isbyte
4fd8e94be0 Gael*0975 #endif /* ALLOW_AUTODIFF_TAMC */
52ff14d141 Ian *0976
8377b8ee87 Mart*0977 #if (defined ALLOW_AUTODIFF && defined SEAICE_MODIFY_GROWTH_ADJ)
83ad492c2d Jean*0978 Cgf no additional dependency through ice cover
4bc8f4264e Mart*0979         IF ( SEAICEadjMODE.GE.3 ) THEN
286983d3d2 Patr*0980 #ifdef SEAICE_ITD
f913c5a485 Mart*0981          DO IT=1,SEAICE_multDim
8377b8ee87 Mart*0982           DO j=1,sNy
                0983            DO i=1,sNx
                0984             a_QbyATMmult_cover(i,j,IT)   = 0. _d 0
                0985             r_QbyATMmult_cover(i,j,IT)   = 0. _d 0
                0986             a_QSWbyATMmult_cover(i,j,IT) = 0. _d 0
286983d3d2 Patr*0987            ENDDO
                0988           ENDDO
114c791332 Jean*0989          ENDDO
cbd0ee24a8 Mart*0990 #else /* ndef SEAICE_ITD */
8377b8ee87 Mart*0991          DO j=1,sNy
                0992           DO i=1,sNx
                0993            a_QbyATM_cover(i,j)   = 0. _d 0
                0994            r_QbyATM_cover(i,j)   = 0. _d 0
                0995            a_QSWbyATM_cover(i,j) = 0. _d 0
4bc8f4264e Mart*0996           ENDDO
3a3bf6419a Gael*0997          ENDDO
cbd0ee24a8 Mart*0998 #endif /* SEAICE_ITD */
4bc8f4264e Mart*0999         ENDIF
3a3bf6419a Gael*1000 #endif
                1001
0c0ecd4c7b Jean*1002 C determine available heat due to the ice pack tying the
2ae913cfea Gael*1003 C underlying surface water temperature to freezing point
                1004 C ======================================================
                1005
33e17487ce Dimi*1006 #ifdef ALLOW_AUTODIFF_TAMC
4fd8e94be0 Gael*1007 CADJ STORE theta(:,:,kSurface,bi,bj) = comlev1_bibj,
edb6656069 Mart*1008 CADJ &                       key = tkey, byte = isbyte
4fd8e94be0 Gael*1009 CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
edb6656069 Mart*1010 CADJ &                       key = tkey, byte = isbyte
33e17487ce Dimi*1011 #endif
01e3cf59a2 Gael*1012
8377b8ee87 Mart*1013         DO j=1,sNy
                1014          DO i=1,sNx
1d74e34c65 Jean*1015 C         FREEZING TEMP. OF SEA WATER (deg C)
a4bc0a0b4c Jean*1016           tempFrz = SEAICE_tempFrz0 +
8377b8ee87 Mart*1017      &              SEAICE_dTempFrz_dS *salt(i,j,kSurface,bi,bj)
1d74e34c65 Jean*1018 C efficiency of turbulent fluxes : dependency to sign of THETA-TBC
8377b8ee87 Mart*1019           IF ( theta(i,j,kSurface,bi,bj) .GE. tempFrz ) THEN
ceae9498ad Gael*1020            tmpscal1 = SEAICE_mcPheePiston
bb8e6379cb Mart*1021           ELSE
0320e25227 Mart*1022            tmpscal1 =SEAICE_frazilFrac*dzSurf/SEAICE_deltaTtherm
bb8e6379cb Mart*1023           ENDIF
1d74e34c65 Jean*1024 C efficiency of turbulent fluxes : dependency to AREA (McPhee cases)
8377b8ee87 Mart*1025           IF ( (AREApreTH(i,j) .GT. 0. _d 0).AND.
ceae9498ad Gael*1026      &         (.NOT.SEAICE_mcPheeStepFunc) ) THEN
a4bc0a0b4c Jean*1027            MixedLayerTurbulenceFactor = ONE -
8377b8ee87 Mart*1028      &          SEAICE_mcPheeTaper * AREApreTH(i,j)
                1029           ELSEIF ( (AREApreTH(i,j) .GT. 0. _d 0).AND.
ceae9498ad Gael*1030      &             (SEAICE_mcPheeStepFunc) ) THEN
                1031            MixedLayerTurbulenceFactor = ONE - SEAICE_mcPheeTaper
bb8e6379cb Mart*1032           ELSE
                1033            MixedLayerTurbulenceFactor = ONE
                1034           ENDIF
1d74e34c65 Jean*1035 C maximum turbulent flux, in ice meters
bb8e6379cb Mart*1036           tmpscal2= - (HeatCapacity_Cp*rhoConst * recip_QI)
8377b8ee87 Mart*1037      &         * (theta(i,j,kSurface,bi,bj)-tempFrz)
ec0d7df165 Mart*1038      &         * SEAICE_deltaTtherm * HEFFM(i,j,bi,bj)
1d74e34c65 Jean*1039 C available turbulent flux
bb8e6379cb Mart*1040           a_QbyOCN(i,j) =
                1041      &         tmpscal1 * tmpscal2 * MixedLayerTurbulenceFactor
                1042           r_QbyOCN(i,j) = a_QbyOCN(i,j)
01e3cf59a2 Gael*1043          ENDDO
                1044         ENDDO
83ad492c2d Jean*1045
53b2f6dc29 Torg*1046 #ifdef SEAICE_ITD
                1047 C determine lateral melt rate at floe edges based on an
6571f3ca98 Jean*1048 C average floe diameter or a floe size distribution
53b2f6dc29 Torg*1049 C following Steele (1992, Tab. 2)
                1050 C ======================================================
f913c5a485 Mart*1051         DO IT=1,SEAICE_multDim
8377b8ee87 Mart*1052          DO j=1,sNy
                1053           DO i=1,sNx
9a1fd902e7 Mart*1054            tempFrz = SEAICE_tempFrz0 +
8377b8ee87 Mart*1055      &          SEAICE_dTempFrz_dS *salt(i,j,kSurface,bi,bj)
                1056            tmpscal1=(theta(i,j,kSurface,bi,bj)-tempFrz)
                1057            tmpscal2=SQRT(0.87 + 0.067*UG(i,j)) * UG(i,j)
9a1fd902e7 Mart*1058
c6b168144e Jean*1059 C     variable floe diameter following Luepkes et al. (2012, JGR, Equ. 26)
9a1fd902e7 Mart*1060 C     with beta=1
                1061 CML           tmpscal3=ONE/(ONE-(floeDiameterMin/floeDiameterMax))
                1062 CML           floeDiameter = floeDiameterMin
                1063 CML     &          * (tmpscal3 / (tmpscal3-AREApreTH(I,J)))
                1064 C     this form involves fewer divisions but gives the same result
                1065            floeDiameter = floeDiameterMin * floeDiameterMax
8377b8ee87 Mart*1066      &          / ( floeDiameterMax*( 1. _d 0 - AREApreTH(i,j) )
                1067      &            + floeDiameterMin*AREApreTH(i,j) )
9a1fd902e7 Mart*1068 C     following the idea of SEAICE_areaLossFormula == 1:
8377b8ee87 Mart*1069            IF (a_QbyATMmult_cover(i,j,IT).LT.ZERO .OR.
53b2f6dc29 Torg*1070      &         a_QbyATM_open(i,j) .LT.ZERO .OR.
                1071      &         a_QbyOCN(i,j)      .LT.ZERO) THEN
9a1fd902e7 Mart*1072 C     lateral melt rate as suggested by Perovich, 1983 (PhD thesis)
8377b8ee87 Mart*1073 c           latMeltRate(i,j,IT) = 1.6 _d -6 * tmpscal1**1.36
47771cb6d0 Mart*1074 C     The following for does the same, but is faster
8377b8ee87 Mart*1075             latMeltRate(i,j,IT) = ZERO
a24915ab1a Jean*1076             IF (tmpscal1 .GT. ZERO)
8377b8ee87 Mart*1077      &         latMeltRate(i,j,IT) = 1.6 _d -6 * exp(1.36*log(tmpscal1))
c6b168144e Jean*1078 C     lateral melt rate as suggested by Maykut and Perovich, 1987
9a1fd902e7 Mart*1079 C     (JGR 92(C7)), Equ. 24
8377b8ee87 Mart*1080 c           latMeltRate(i,j,IT) = 13.5 _d -6 * tmpscal2 * tmpscal1**1.3
c6b168144e Jean*1081 C     further suggestion by Maykut and Perovich to avoid
9a1fd902e7 Mart*1082 C     latMeltRate -> 0 for UG -> 0
8377b8ee87 Mart*1083 c           latMeltRate(i,j,IT) = (1.6 _d -6 + 13.5 _d -6 * tmpscal2)
9a1fd902e7 Mart*1084 c    &                          * tmpscal1**1.3
                1085 C     factor determining fraction of area and ice volume reduction
                1086 C     due to lateral melt
8377b8ee87 Mart*1087             latMeltFrac(i,j,IT) =
                1088      &       latMeltRate(i,j,IT)*SEAICE_deltaTtherm*PI /
6571f3ca98 Jean*1089      &       (floeAlpha * floeDiameter)
8377b8ee87 Mart*1090             latMeltFrac(i,j,IT)=max(ZERO, min(latMeltFrac(i,j,IT),ONE))
6571f3ca98 Jean*1091            ELSE
8377b8ee87 Mart*1092             latMeltRate(i,j,IT)=0.0 _d 0
                1093             latMeltFrac(i,j,IT)=0.0 _d 0
6571f3ca98 Jean*1094            ENDIF
                1095           ENDDO
                1096          ENDDO
53b2f6dc29 Torg*1097         ENDDO
cbd0ee24a8 Mart*1098 #endif /* SEAICE_ITD */
53b2f6dc29 Torg*1099
8377b8ee87 Mart*1100 #if (defined ALLOW_AUTODIFF && defined SEAICE_MODIFY_GROWTH_ADJ)
4bc8f4264e Mart*1101         CALL ZERO_ADJ_1D( sNx*sNy, r_QbyOCN, myThid)
d187c22362 Gael*1102 #endif
aea7db20a6 Gael*1103
2651ba3350 Jean*1104 C ===================================================================
                1105 C =========PART 3: determine effective thicknesses increments========
                1106 C ===================================================================
aea7db20a6 Gael*1107
56d13a40ed Mart*1108 #if defined ( ALLOW_SITRACER ) && defined ( SEAICE_GREASE )
f61838dfc1 Torg*1109 C convert SItracer 'grease' from ratio to grease ice volume:
                1110 C ==========================================================
56d13a40ed Mart*1111         DO j=1,sNy
                1112          DO i=1,sNx
                1113           SItracer(i,j,bi,bj,iTrGrease) =
                1114      &     SItracer(i,j,bi,bj,iTrGrease) * HEFF(i,j,bi,bj)
                1115           uRelW(i,j) = uWind(i,j,bi,bj)
                1116           vRelW(i,j) = vWind(i,j,bi,bj)
f61838dfc1 Torg*1117          ENDDO
                1118         ENDDO
56d13a40ed Mart*1119         IF ( useRelativeWind ) THEN
                1120          DO j=1,sNy
                1121           DO i=1,sNx
                1122            uRelW(i,j) = uWind(i,j,bi,bj)
                1123      &          - 0.5 _d 0*( uIce(i,j,bi,bj) + uIce(i+1,j,bi,bj) )
                1124            vRelW(i,j) = vWind(i,j,bi,bj)
                1125      &          - 0.5 _d 0*( vIce(i,j,bi,bj) + vIce(i,j+1,bi,bj) )
                1126           ENDDO
                1127          ENDDO
                1128         ENDIF
f61838dfc1 Torg*1129 C compute actual grease ice layer thickness
                1130 C as a function of wind stress
                1131 C following Smedsrud [2011, Ann.Glac.]
                1132 C =========================================
56d13a40ed Mart*1133          DO j=1,sNy
                1134           DO i=1,sNx
f61838dfc1 Torg*1135 C
c6b168144e Jean*1136 C          computing compaction force acting on grease
f61838dfc1 Torg*1137 C           (= air and water stress acting on ice)
                1138 C          wind stress (using calculation of SPEED_SQ & UG as template)
                1139 C          u_a**2 + v_a**2:
56d13a40ed Mart*1140            tmpscal1 = uRelW(i,j)*uRelW(i,j) + vRelW(i,j)*vRelW(i,j)
f61838dfc1 Torg*1141            IF ( tmpscal1 .LE. SEAICE_EPS_SQ ) THEN
                1142              tmpscal1=SEAICE_EPS
                1143            ELSE
                1144              tmpscal1=SQRT(tmpscal2)
                1145            ENDIF
                1146            tmpscal1 = 1.4 _d 0 * 1.3 _d -3 * tmpscal1
                1147 C          water stress
                1148 C          u_w - u_i:
c6b168144e Jean*1149            tmpscal2 =
f61838dfc1 Torg*1150      &        0.5 _d 0*(uVel(i,j,kSurface,bi,bj)
                1151      &                 +uVel(i+1,j,kSurface,bi,bj))
56d13a40ed Mart*1152      &       -0.5 _d 0*(uIce(i,j,bi,bj)+uIce(i+1,j,bi,bj))
f61838dfc1 Torg*1153 C          v_w - v_i:
c6b168144e Jean*1154            tmpscal3 =
f61838dfc1 Torg*1155      &        0.5 _d 0*(vVel(i,j,kSurface,bi,bj)
                1156      &                 +vVel(i,j+1,kSurface,bi,bj))
56d13a40ed Mart*1157      &       -0.5 _d 0*(vIce(i,j,bi,bj)+vIce(i,j+1,bi,bj))
f61838dfc1 Torg*1158 C          (u_w - u_i)**2 + (v_w - v_i)**2:
                1159            tmpscal4 = (tmpscal2*tmpscal2 + tmpscal3*tmpscal3)
                1160            IF ( tmpscal4 .LE. SEAICE_EPS_SQ ) THEN
                1161              tmpscal4=SEAICE_EPS
                1162            ELSE
                1163              tmpscal4=SQRT(tmpscal4)
                1164            ENDIF
                1165            tmpscal4 = 1027.0 _d 0 * 6.0 _d -3 * tmpscal4
                1166 C          magnitude of compined stresses:
c6b168144e Jean*1167            tmpscal0 =
56d13a40ed Mart*1168      &         ( tmpscal1 * uRelW(i,j) + tmpscal4 * tmpscal2 )**2
                1169      &       + ( tmpscal1 * vRelW(i,j) + tmpscal4 * tmpscal3 )**2
f61838dfc1 Torg*1170            IF ( tmpscal0 .LE. SEAICE_EPS_SQ ) THEN
                1171              tmpscal0=SEAICE_EPS
                1172            ELSE
                1173              tmpscal0=SQRT(tmpscal0)
                1174            ENDIF
                1175 C
                1176 C          mean grid cell width between tracer points
56d13a40ed Mart*1177            tmpscal3 = 0.5 _d 0 * (dxC(i,j,bi,bj)+dyC(i,j,bi,bj))
f61838dfc1 Torg*1178 C          grease ice volume Vg [m^3/m] as in Smedsrud [2011]
                1179 C           is SItracer * gridcell_area / lead_width
c6b168144e Jean*1180 C           where lead width is lead extent along lead,
f61838dfc1 Torg*1181 C           i.e. perpendicular to L_lead in Smedsrud (2011),
                1182 C           which is in the absence of lead length statistics
                1183 C           simply the grid cell length
56d13a40ed Mart*1184            tmpscal4 = 4. _d 0 * SItracer(i,j,bi,bj,iTrGrease) * tmpscal3
c6b168144e Jean*1185 C
f61838dfc1 Torg*1186 C          mean grease ice layer thickness <h_g>, Eq.10 in Smedsrud [2011] but incl. water stress
56d13a40ed Mart*1187            greaseLayerThick(i,j) = 0.763 _d 0
f61838dfc1 Torg*1188 C          grease ice volume
                1189      &           * ( tmpscal4
                1190 C          times magnitude of vector sum of air and water stresses
                1191 C          (in fact, only the component perpendicular to the lead edge, i.e. along L_lead counts)
cbeffc6b9a Jean*1192 C          ATTENTION: since we do not have lead orientation with respect to wind
f61838dfc1 Torg*1193 C                     we use 80% of the total force assuming angles 45-90deg between wind and lead edge
                1194      &           * 0.8 _d 0 * tmpscal0
                1195 C          devided by K_r = 100 N/m^3 (resistance of grease against compression)
                1196      &           * 0.01 _d 0 )**THIRD
                1197 C
                1198 C          assure a minimum thickness of 4 cm (equals HO=0.01):
56d13a40ed Mart*1199            greaseLayerThick(i,j)=max(4. _d -2, greaseLayerThick(i,j))
f61838dfc1 Torg*1200 C          ... and a maximum thickness of 4 m (equals HO=1.0):
56d13a40ed Mart*1201            greaseLayerThick(i,j)=min(4. _d 0 , greaseLayerThick(i,j))
f61838dfc1 Torg*1202 C
                1203           ENDDO
                1204          ENDDO
                1205 #endif /* SEAICE_GREASE */
                1206
ae36251cae Gael*1207 C compute snow/ice tendency due to sublimation
                1208 C ============================================
2ae913cfea Gael*1209
b34884f5be Mart*1210 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1211 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                1212 CADJ STORE r_FWbySublim     = comlev1_bibj,key=tkey,byte=isbyte
b34884f5be Mart*1213 #endif /* ALLOW_AUTODIFF_TAMC */
286983d3d2 Patr*1214 #ifdef SEAICE_ITD
f913c5a485 Mart*1215         DO IT=1,SEAICE_multDim
cbd0ee24a8 Mart*1216 #endif /* SEAICE_ITD */
8377b8ee87 Mart*1217         DO j=1,sNy
                1218          DO i=1,sNx
ae36251cae Gael*1219 C     First sublimate/deposite snow
a4bc0a0b4c Jean*1220           tmpscal2 =
286983d3d2 Patr*1221 #ifdef SEAICE_ITD
8377b8ee87 Mart*1222      &     MAX(MIN(r_FWbySublimMult(i,j,IT),HSNOWITD(i,j,IT,bi,bj)
286983d3d2 Patr*1223      &             *SNOW2ICE),ZERO)
8377b8ee87 Mart*1224           d_HSNWbySublim_ITD(i,j,IT) = - tmpscal2 * ICE2SNOW
286983d3d2 Patr*1225 C         accumulate change over ITD categories
8377b8ee87 Mart*1226           d_HSNWbySublim(i,j)     = d_HSNWbySublim(i,j)      - tmpscal2
286983d3d2 Patr*1227      &                                                       *ICE2SNOW
8377b8ee87 Mart*1228           r_FWbySublimMult(i,j,IT)= r_FWbySublimMult(i,j,IT) - tmpscal2
cbd0ee24a8 Mart*1229 #else /* ndef SEAICE_ITD */
8377b8ee87 Mart*1230      &     MAX(MIN(r_FWbySublim(i,j),HSNOW(i,j,bi,bj)*SNOW2ICE),ZERO)
                1231           d_HSNWbySublim(i,j) = - tmpscal2 * ICE2SNOW
                1232           HSNOW(i,j,bi,bj)    = HSNOW(i,j,bi,bj)  - tmpscal2*ICE2SNOW
                1233           r_FWbySublim(i,j)   = r_FWbySublim(i,j) - tmpscal2
cbd0ee24a8 Mart*1234 #endif /* SEAICE_ITD */
b34884f5be Mart*1235          ENDDO
                1236         ENDDO
ae36251cae Gael*1237 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1238 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                1239 CADJ STORE r_FWbySublim    = comlev1_bibj,key=tkey,byte=isbyte
ae36251cae Gael*1240 #endif /* ALLOW_AUTODIFF_TAMC */
8377b8ee87 Mart*1241         DO j=1,sNy
                1242          DO i=1,sNx
ae36251cae Gael*1243 C     If anything is left, sublimate ice
a4bc0a0b4c Jean*1244           tmpscal2 =
286983d3d2 Patr*1245 #ifdef SEAICE_ITD
8377b8ee87 Mart*1246      &     MAX(MIN(r_FWbySublimMult(i,j,IT),HEFFITD(i,j,IT,bi,bj)),ZERO)
                1247           d_HEFFbySublim_ITD(i,j,IT) = - tmpscal2
286983d3d2 Patr*1248 C         accumulate change over ITD categories
8377b8ee87 Mart*1249           d_HEFFbySublim(i,j)      = d_HEFFbySublim(i,j)      - tmpscal2
                1250           r_FWbySublimMult(i,j,IT) = r_FWbySublimMult(i,j,IT) - tmpscal2
cbd0ee24a8 Mart*1251 #else /* ndef SEAICE_ITD */
8377b8ee87 Mart*1252      &     MAX(MIN(r_FWbySublim(i,j),HEFF(i,j,bi,bj)),ZERO)
                1253           d_HEFFbySublim(i,j) = - tmpscal2
                1254           HEFF(i,j,bi,bj)     = HEFF(i,j,bi,bj)   - tmpscal2
                1255           r_FWbySublim(i,j)   = r_FWbySublim(i,j) - tmpscal2
cbd0ee24a8 Mart*1256 #endif /* SEAICE_ITD */
ae36251cae Gael*1257          ENDDO
                1258         ENDDO
8377b8ee87 Mart*1259         DO j=1,sNy
                1260          DO i=1,sNx
eed5d4f5a4 Gael*1261 C     If anything is left, it will be evaporated from the ocean rather than sublimated.
286983d3d2 Patr*1262 C     Since a_QbyATM_cover was computed for sublimation, not simple evaporation, we need to
eed5d4f5a4 Gael*1263 C     remove the fusion part for the residual (that happens to be precisely r_FWbySublim).
286983d3d2 Patr*1264 #ifdef SEAICE_ITD
8377b8ee87 Mart*1265           a_QbyATMmult_cover(i,j,IT) = a_QbyATMmult_cover(i,j,IT)
                1266      &                               - r_FWbySublimMult(i,j,IT)
                1267           r_QbyATMmult_cover(i,j,IT) = r_QbyATMmult_cover(i,j,IT)
                1268      &                               - r_FWbySublimMult(i,j,IT)
cbd0ee24a8 Mart*1269 #else /* ndef SEAICE_ITD */
8377b8ee87 Mart*1270           a_QbyATM_cover(i,j) = a_QbyATM_cover(i,j)-r_FWbySublim(i,j)
                1271           r_QbyATM_cover(i,j) = r_QbyATM_cover(i,j)-r_FWbySublim(i,j)
cbd0ee24a8 Mart*1272 #endif /* SEAICE_ITD */
eed5d4f5a4 Gael*1273          ENDDO
                1274         ENDDO
286983d3d2 Patr*1275 #ifdef SEAICE_ITD
                1276 C       end IT loop
114c791332 Jean*1277         ENDDO
cbd0ee24a8 Mart*1278 #endif /* SEAICE_ITD */
b34884f5be Mart*1279
52ff14d141 Ian *1280 C compute ice thickness tendency due to ice-ocean interaction
                1281 C ===========================================================
                1282
                1283 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1284 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                1285 CADJ STORE r_QbyOCN = comlev1_bibj,key=tkey,byte=isbyte
52ff14d141 Ian *1286 #endif /* ALLOW_AUTODIFF_TAMC */
                1287
62cc8945c8 Gael*1288        IF (.NOT.SEAICE_growMeltByConv) THEN
6571f3ca98 Jean*1289
286983d3d2 Patr*1290 #ifdef SEAICE_ITD
f913c5a485 Mart*1291         DO IT=1,SEAICE_multDim
8377b8ee87 Mart*1292          DO j=1,sNy
                1293           DO i=1,sNx
114c791332 Jean*1294 C          ice growth/melt due to ocean heat r_QbyOCN (W/m^2) is
                1295 C          equally distributed under the ice and hence weighted by
286983d3d2 Patr*1296 C          fractional area of each thickness category
8377b8ee87 Mart*1297            tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(i,j,IT),
                1298      &                               -HEFFITD(i,j,IT,bi,bj))
                1299            d_HEFFbyOCNonICE_ITD(i,j,IT)=tmpscal1
                1300            d_HEFFbyOCNonICE(i,j) = d_HEFFbyOCNonICE(i,j) + tmpscal1
286983d3d2 Patr*1301           ENDDO
                1302          ENDDO
                1303         ENDDO
                1304 #ifdef ALLOW_SITRACER
8377b8ee87 Mart*1305         DO j=1,sNy
                1306          DO i=1,sNx
                1307           SItrHEFF(i,j,bi,bj,2) = HEFFpreTH(i,j)
                1308      &                          + d_HEFFbySublim(i,j)
                1309      &                          + d_HEFFbyOCNonICE(i,j)
286983d3d2 Patr*1310          ENDDO
                1311         ENDDO
                1312 #endif
8377b8ee87 Mart*1313         DO j=1,sNy
                1314          DO i=1,sNx
                1315           r_QbyOCN(i,j)=r_QbyOCN(i,j)-d_HEFFbyOCNonICE(i,j)
286983d3d2 Patr*1316          ENDDO
                1317         ENDDO
                1318 #else /* SEAICE_ITD */
8377b8ee87 Mart*1319         DO j=1,sNy
                1320          DO i=1,sNx
                1321           d_HEFFbyOCNonICE(i,j)=MAX(r_QbyOCN(i,j), -HEFF(i,j,bi,bj))
                1322           r_QbyOCN(i,j)=r_QbyOCN(i,j)-d_HEFFbyOCNonICE(i,j)
                1323           HEFF(i,j,bi,bj)=HEFF(i,j,bi,bj) + d_HEFFbyOCNonICE(i,j)
52ff14d141 Ian *1324 #ifdef ALLOW_SITRACER
8377b8ee87 Mart*1325           SItrHEFF(i,j,bi,bj,2)=HEFF(i,j,bi,bj)
52ff14d141 Ian *1326 #endif
                1327          ENDDO
                1328         ENDDO
286983d3d2 Patr*1329 #endif /* SEAICE_ITD */
52ff14d141 Ian *1330
634144d037 Jean*1331       ENDIF !SEAICE_growMeltByConv
62cc8945c8 Gael*1332
ae36251cae Gael*1333 C compute snow melt tendency due to snow-atmosphere interaction
                1334 C ==================================================================
                1335
33e17487ce Dimi*1336 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1337 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                1338 CADJ STORE r_QbyATM_cover = comlev1_bibj,key=tkey,byte=isbyte
33e17487ce Dimi*1339 #endif /* ALLOW_AUTODIFF_TAMC */
                1340
286983d3d2 Patr*1341 #ifdef SEAICE_ITD
f913c5a485 Mart*1342         DO IT=1,SEAICE_multDim
8377b8ee87 Mart*1343          DO j=1,sNy
                1344           DO i=1,sNx
286983d3d2 Patr*1345 C     Convert to standard units (meters of ice) rather than to meters
                1346 C     of snow. This appears to be more robust.
8377b8ee87 Mart*1347            tmpscal1=MAX(r_QbyATMmult_cover(i,j,IT),
                1348      &                  -HSNOWITD(i,j,IT,bi,bj)*SNOW2ICE)
286983d3d2 Patr*1349            tmpscal2=MIN(tmpscal1,0. _d 0)
                1350 #ifdef SEAICE_MODIFY_GROWTH_ADJ
                1351 Cgf no additional dependency through snow
                1352            IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
                1353 #endif
8377b8ee87 Mart*1354            d_HSNWbyATMonSNW_ITD(i,j,IT) = tmpscal2*ICE2SNOW
                1355            d_HSNWbyATMonSNW(i,j) = d_HSNWbyATMonSNW(i,j)
286983d3d2 Patr*1356      &                           + tmpscal2*ICE2SNOW
8377b8ee87 Mart*1357            r_QbyATMmult_cover(i,j,IT)=r_QbyATMmult_cover(i,j,IT)
286983d3d2 Patr*1358      &                           - tmpscal2
114c791332 Jean*1359           ENDDO
                1360          ENDDO
                1361         ENDDO
286983d3d2 Patr*1362 #else /* SEAICE_ITD */
8377b8ee87 Mart*1363         DO j=1,sNy
                1364          DO i=1,sNx
bb8e6379cb Mart*1365 C     Convert to standard units (meters of ice) rather than to meters
                1366 C     of snow. This appears to be more robust.
8377b8ee87 Mart*1367           tmpscal1=MAX(r_QbyATM_cover(i,j),-HSNOW(i,j,bi,bj)*SNOW2ICE)
3a3bf6419a Gael*1368           tmpscal2=MIN(tmpscal1,0. _d 0)
                1369 #ifdef SEAICE_MODIFY_GROWTH_ADJ
                1370 Cgf no additional dependency through snow
4bc8f4264e Mart*1371           IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
3a3bf6419a Gael*1372 #endif
8377b8ee87 Mart*1373           d_HSNWbyATMonSNW(i,j)= tmpscal2*ICE2SNOW
                1374           HSNOW(i,j,bi,bj) = HSNOW(i,j,bi,bj) + tmpscal2*ICE2SNOW
                1375           r_QbyATM_cover(i,j)=r_QbyATM_cover(i,j) - tmpscal2
33e17487ce Dimi*1376          ENDDO
                1377         ENDDO
286983d3d2 Patr*1378 #endif /* SEAICE_ITD */
33e17487ce Dimi*1379
c50ad14e64 Gael*1380 C compute ice thickness tendency due to the atmosphere
                1381 C ====================================================
2ae913cfea Gael*1382
b34884f5be Mart*1383 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1384 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                1385 CADJ STORE r_QbyATM_cover  = comlev1_bibj,key=tkey,byte=isbyte
3a3bf6419a Gael*1386 #endif /* ALLOW_AUTODIFF_TAMC */
2ae913cfea Gael*1387
2651ba3350 Jean*1388 Cgf note: this block is not actually tested by lab_sea
                1389 Cgf where all experiments start in January. So even though
                1390 Cgf the v1.81=>v1.82 revision would change results in
                1391 Cgf warming conditions, the lab_sea results were not changed.
e27b57218b Gael*1392
286983d3d2 Patr*1393 #ifdef SEAICE_ITD
f913c5a485 Mart*1394         DO IT=1,SEAICE_multDim
8377b8ee87 Mart*1395          DO j=1,sNy
                1396           DO i=1,sNx
                1397            tmpscal1 = HEFFITDpreTH(i,j,IT)
                1398      &              + d_HEFFbySublim_ITD(i,j,IT)
                1399      &              + d_HEFFbyOCNonICE_ITD(i,j,IT)
286983d3d2 Patr*1400            tmpscal2 = MAX(-tmpscal1,
8377b8ee87 Mart*1401      &                     r_QbyATMmult_cover(i,j,IT)
1080be7801 Jean*1402 C         Limit ice growth by potential melt by ocean
8377b8ee87 Mart*1403      &              + AREAITDpreTH(i,j,IT) * r_QbyOCN(i,j))
                1404            d_HEFFbyATMonOCN_cover_ITD(i,j,IT) = tmpscal2
                1405            d_HEFFbyATMonOCN_cover(i,j) = d_HEFFbyATMonOCN_cover(i,j)
286983d3d2 Patr*1406      &                                 + tmpscal2
8377b8ee87 Mart*1407            d_HEFFbyATMonOCN_ITD(i,j,IT) = d_HEFFbyATMonOCN_ITD(i,j,IT)
286983d3d2 Patr*1408      &                                 + tmpscal2
8377b8ee87 Mart*1409            d_HEFFbyATMonOCN(i,j)       = d_HEFFbyATMonOCN(i,j)
286983d3d2 Patr*1410      &                                 + tmpscal2
8377b8ee87 Mart*1411            r_QbyATMmult_cover(i,j,IT)  = r_QbyATMmult_cover(i,j,IT)
286983d3d2 Patr*1412      &                                 - tmpscal2
114c791332 Jean*1413           ENDDO
                1414          ENDDO
                1415         ENDDO
286983d3d2 Patr*1416 #ifdef ALLOW_SITRACER
8377b8ee87 Mart*1417         DO j=1,sNy
                1418          DO i=1,sNx
                1419           SItrHEFF(i,j,bi,bj,3) = SItrHEFF(i,j,bi,bj,2)
                1420      &                          + d_HEFFbyATMonOCN_cover(i,j)
114c791332 Jean*1421          ENDDO
                1422         ENDDO
286983d3d2 Patr*1423 #endif
cbd0ee24a8 Mart*1424 #else /* ndef SEAICE_ITD */
8377b8ee87 Mart*1425         DO j=1,sNy
                1426          DO i=1,sNx
2afe30fba0 Dimi*1427
8377b8ee87 Mart*1428           tmpscal2 = MAX(-HEFF(i,j,bi,bj),r_QbyATM_cover(i,j)+
1d74e34c65 Jean*1429 C         Limit ice growth by potential melt by ocean
8377b8ee87 Mart*1430      &        AREApreTH(i,j) * r_QbyOCN(i,j))
2afe30fba0 Dimi*1431
8377b8ee87 Mart*1432           d_HEFFbyATMonOCN_cover(i,j)=tmpscal2
                1433           d_HEFFbyATMonOCN(i,j)=d_HEFFbyATMonOCN(i,j)+tmpscal2
                1434           r_QbyATM_cover(i,j)=r_QbyATM_cover(i,j)-tmpscal2
                1435           HEFF(i,j,bi,bj) = HEFF(i,j,bi,bj) + tmpscal2
2afe30fba0 Dimi*1436
f50f58ec54 Gael*1437 #ifdef ALLOW_SITRACER
8377b8ee87 Mart*1438           SItrHEFF(i,j,bi,bj,3)=HEFF(i,j,bi,bj)
f50f58ec54 Gael*1439 #endif
ee38904d3a Gael*1440          ENDDO
                1441         ENDDO
286983d3d2 Patr*1442 #endif /* SEAICE_ITD */
33e17487ce Dimi*1443
ed104d6f21 Dimi*1444 C add snow precipitation to HSNOW.
2ae913cfea Gael*1445 C =================================================
634144d037 Jean*1446 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1447 CADJ STORE a_QbyATM_cover = comlev1_bibj,key=tkey,byte=isbyte
                1448 CADJ STORE PRECIP(:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                1449 CADJ STORE AREApreTH = comlev1_bibj,key=tkey,byte=isbyte
634144d037 Jean*1450 #endif /* ALLOW_AUTODIFF_TAMC */
ed104d6f21 Dimi*1451         IF ( snowPrecipFile .NE. ' ' ) THEN
                1452 C add snowPrecip to HSNOW
8377b8ee87 Mart*1453          DO j=1,sNy
                1454           DO i=1,sNx
                1455            d_HSNWbyRAIN(i,j) = convertPRECIP2HI * ICE2SNOW *
                1456      &          snowPrecip(i,j,bi,bj) * AREApreTH(i,j)
                1457            d_HFRWbyRAIN(i,j) = -convertPRECIP2HI *
                1458      &          ( PRECIP(i,j,bi,bj) - snowPrecip(i,j,bi,bj) ) *
                1459      &          AREApreTH(i,j)
                1460            HSNOW(i,j,bi,bj) = HSNOW(i,j,bi,bj) + d_HSNWbyRAIN(i,j)
ed104d6f21 Dimi*1461           ENDDO
                1462          ENDDO
                1463         ELSE
                1464 C attribute precip to fresh water or snow stock,
                1465 C depending on atmospheric conditions.
8377b8ee87 Mart*1466          DO j=1,sNy
                1467           DO i=1,sNx
2651ba3350 Jean*1468 C possible alternatives to the a_QbyATM_cover criterium
b69fbfd195 Mart*1469 c          IF (tIce(I,J,bi,bj) .LT. TMIX) THEN ! would require computing tIce
c50ad14e64 Gael*1470 c          IF (atemp(I,J,bi,bj) .LT. celsius2K) THEN
8377b8ee87 Mart*1471            IF ( a_QbyATM_cover(i,j).GE. 0. _d 0 ) THEN
3b4bb9d1ee Gael*1472 C           add precip as snow
8377b8ee87 Mart*1473             d_HFRWbyRAIN(i,j)=0. _d 0
                1474             d_HSNWbyRAIN(i,j)=convertPRECIP2HI*ICE2SNOW*
                1475      &            PRECIP(i,j,bi,bj)*AREApreTH(i,j)
ed104d6f21 Dimi*1476            ELSE
2651ba3350 Jean*1477 C           add precip to the fresh water bucket
8377b8ee87 Mart*1478             d_HFRWbyRAIN(i,j)=-convertPRECIP2HI*
                1479      &            PRECIP(i,j,bi,bj)*AREApreTH(i,j)
                1480             d_HSNWbyRAIN(i,j)=0. _d 0
ed104d6f21 Dimi*1481            ENDIF
286983d3d2 Patr*1482          ENDDO
                1483         ENDDO
                1484 #ifdef SEAICE_ITD
f913c5a485 Mart*1485         DO IT=1,SEAICE_multDim
8377b8ee87 Mart*1486          DO j=1,sNy
                1487           DO i=1,sNx
                1488            d_HSNWbyRAIN_ITD(i,j,IT)
                1489      &     = d_HSNWbyRAIN(i,j)*areaFracFactor(i,j,IT)
ed104d6f21 Dimi*1490           ENDDO
33e17487ce Dimi*1491          ENDDO
114c791332 Jean*1492         ENDDO
cbd0ee24a8 Mart*1493 #else /* ndef SEAICE_ITD */
8377b8ee87 Mart*1494         DO j=1,sNy
                1495          DO i=1,sNx
                1496           HSNOW(i,j,bi,bj) = HSNOW(i,j,bi,bj) + d_HSNWbyRAIN(i,j)
286983d3d2 Patr*1497          ENDDO
                1498         ENDDO
cbd0ee24a8 Mart*1499 #endif /* SEAICE_ITD */
2651ba3350 Jean*1500 Cgf note: this does not affect air-sea heat flux,
                1501 Cgf since the implied air heat gain to turn
                1502 Cgf rain to snow is not a surface process.
286983d3d2 Patr*1503 C end of IF statement snowPrecipFile:
ed104d6f21 Dimi*1504         ENDIF
33e17487ce Dimi*1505
c50ad14e64 Gael*1506 C compute snow melt due to heat available from ocean.
2ae913cfea Gael*1507 C =================================================================
33e17487ce Dimi*1508
2651ba3350 Jean*1509 Cgf do we need to keep this comment and cpp bracket?
                1510 Cph( very sensitive bit here by JZ
33e17487ce Dimi*1511 #ifndef SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING
3a3bf6419a Gael*1512 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1513 CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                1514 CADJ STORE r_QbyOCN = comlev1_bibj,key=tkey,byte=isbyte
3a3bf6419a Gael*1515 #endif /* ALLOW_AUTODIFF_TAMC */
286983d3d2 Patr*1516
62cc8945c8 Gael*1517       IF (.NOT.SEAICE_growMeltByConv) THEN
                1518
286983d3d2 Patr*1519 #ifdef SEAICE_ITD
f913c5a485 Mart*1520         DO IT=1,SEAICE_multDim
8377b8ee87 Mart*1521          DO j=1,sNy
                1522           DO i=1,sNx
                1523            tmpscal4 = HSNWITDpreTH(i,j,IT)
                1524      &              + d_HSNWbySublim_ITD(i,j,IT)
                1525      &              + d_HSNWbyATMonSNW_ITD(i,j,IT)
                1526      &              + d_HSNWbyRAIN_ITD(i,j,IT)
                1527            tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(i,j,IT),
286983d3d2 Patr*1528      &                  -tmpscal4)
                1529            tmpscal2=MIN(tmpscal1,0. _d 0)
                1530 #ifdef SEAICE_MODIFY_GROWTH_ADJ
                1531 Cgf no additional dependency through snow
8377b8ee87 Mart*1532            IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
286983d3d2 Patr*1533 #endif
8377b8ee87 Mart*1534            d_HSNWbyOCNonSNW_ITD(i,j,IT) = tmpscal2
                1535            d_HSNWbyOCNonSNW(i,j) = d_HSNWbyOCNonSNW(i,j) + tmpscal2
                1536            r_QbyOCN(i,j)=r_QbyOCN(i,j) - tmpscal2*SNOW2ICE
286983d3d2 Patr*1537           ENDDO
                1538          ENDDO
                1539         ENDDO
cbd0ee24a8 Mart*1540 #else /* ndef SEAICE_ITD */
8377b8ee87 Mart*1541         DO j=1,sNy
                1542          DO i=1,sNx
                1543           tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW, -HSNOW(i,j,bi,bj))
3a3bf6419a Gael*1544           tmpscal2=MIN(tmpscal1,0. _d 0)
                1545 #ifdef SEAICE_MODIFY_GROWTH_ADJ
                1546 Cgf no additional dependency through snow
8377b8ee87 Mart*1547           IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
3a3bf6419a Gael*1548 #endif
8377b8ee87 Mart*1549           d_HSNWbyOCNonSNW(i,j) = tmpscal2
                1550           r_QbyOCN(i,j)=r_QbyOCN(i,j)
                1551      &                               -d_HSNWbyOCNonSNW(i,j)*SNOW2ICE
                1552           HSNOW(i,j,bi,bj) = HSNOW(i,j,bi,bj)+d_HSNWbyOCNonSNW(i,j)
33e17487ce Dimi*1553          ENDDO
                1554         ENDDO
286983d3d2 Patr*1555 #endif /* SEAICE_ITD */
62cc8945c8 Gael*1556
634144d037 Jean*1557       ENDIF !SEAICE_growMeltByConv
62cc8945c8 Gael*1558
2ae913cfea Gael*1559 #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */
2651ba3350 Jean*1560 Cph)
85586adda4 Gael*1561
                1562 C gain of new ice over open water
                1563 C ===============================
                1564 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1565 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                1566 CADJ STORE r_QbyATM_open = comlev1_bibj,key=tkey,byte=isbyte
                1567 CADJ STORE r_QbyOCN = comlev1_bibj,key=tkey,byte=isbyte
                1568 CADJ STORE a_QSWbyATM_cover = comlev1_bibj,key=tkey,byte=isbyte
                1569 CADJ STORE a_QSWbyATM_open = comlev1_bibj,key=tkey,byte=isbyte
3a3bf6419a Gael*1570 #endif /* ALLOW_AUTODIFF_TAMC */
6ec4646d60 Gael*1571
8377b8ee87 Mart*1572         DO j=1,sNy
                1573          DO i=1,sNx
286983d3d2 Patr*1574 #ifdef SEAICE_ITD
                1575 C         HEFF will be updated at the end of PART 3,
                1576 C         hence sum of tendencies so far is needed
8377b8ee87 Mart*1577           tmpscal4 = HEFFpreTH(i,j)
                1578      &             + d_HEFFbySublim(i,j)
                1579      &             + d_HEFFbyOCNonICE(i,j)
                1580      &             + d_HEFFbyATMonOCN(i,j)
cbd0ee24a8 Mart*1581 #else /* ndef SEAICE_ITD */
286983d3d2 Patr*1582 C         HEFF is updated step by step throughout seaice_growth
8377b8ee87 Mart*1583           tmpscal4 = HEFF(i,j,bi,bj)
cbd0ee24a8 Mart*1584 #endif /* SEAICE_ITD */
1d74e34c65 Jean*1585 C           Initial ice growth is triggered by open water
                1586 C           heat flux overcoming potential melt by ocean
8377b8ee87 Mart*1587             tmpscal1=r_QbyATM_open(i,j)+r_QbyOCN(i,j) *
                1588      &            (1.0 _d 0 - AREApreTH(i,j))
1d74e34c65 Jean*1589 C           Penetrative shortwave flux beyond first layer
                1590 C           that is therefore not available to ice growth/melt
8e32c48b8f Mart*1591             tmpscal2=SEAICE_SWFrac * a_QSWbyATM_open(i,j)
6ec4646d60 Gael*1592 C           impose -HEFF as the maxmum melting if SEAICE_doOpenWaterMelt
                1593 C           or 0. otherwise (no melting if not SEAICE_doOpenWaterMelt)
                1594             tmpscal3=facOpenGrow*MAX(tmpscal1-tmpscal2,
8377b8ee87 Mart*1595      &           -tmpscal4*facOpenMelt)*HEFFM(i,j,bi,bj)
56d13a40ed Mart*1596 #if defined ( ALLOW_SITRACER ) && defined ( SEAICE_GREASE )
f61838dfc1 Torg*1597 C Grease ice is a tracer or "bucket" for newly formed frazil ice
                1598 C  that instead of becoming solid sea ice instantly has a half-time
c6b168144e Jean*1599 C  of 1 day (see greaseDecayTime) before solidifying.
                1600 C The most important effect is that for fluxes the grease ice area/volume
f61838dfc1 Torg*1601 C  acts like open water.
                1602 C
                1603 C store freezing/melting condition:
                1604           greaseNewFrazil = max(0.0 _d 0, tmpscal3)
                1605 C
                1606 C 1) mechanical removal of grease by ridging:
                1607 C    if there is no open water left after advection, there cannot be any grease ice
8377b8ee87 Mart*1608           IF ((1.0 _d 0 - AREApreTH(i,j)).LE.siEps) THEN
                1609            tmpscal3 = tmpscal3 + SItracer(i,j,bi,bj,iTrGrease)
                1610            SItracer(i,j,bi,bj,iTrGrease) = 0. _d 0
c6b168144e Jean*1611
8377b8ee87 Mart*1612           ELSEIF (greaseNewFrazil .GT. 0. _d 0) THEN
f61838dfc1 Torg*1613 C    new ice growth goes into grease tracer not HEFF:
c6b168144e Jean*1614            tmpscal3=0. _d 0
f61838dfc1 Torg*1615 C
                1616 C 2) solidification of "old" grease ice
                1617 C    (only when cold enough for ice growth)
                1618 C
c6b168144e Jean*1619 C    time scale dependent solidification
f61838dfc1 Torg*1620 C    (50%  of grease ice area become solid ice within 24h):
                1621            tmpscal1=exp(-SEAICE_deltaTtherm/greaseDecayTime)
                1622 C    gain in solid sea ice volume due to solidified grease:
8377b8ee87 Mart*1623            d_HEFFbyGREASE(i,j) =
                1624      &             SItracer(i,j,bi,bj,iTrGrease)
f61838dfc1 Torg*1625      &           * (1.0 _d 0 - tmpscal1)
                1626 C    ... and related loss of grease ice (tracer) volume
8377b8ee87 Mart*1627            SItracer(i,j,bi,bj,iTrGrease) =
                1628      &             SItracer(i,j,bi,bj,iTrGrease) * tmpscal1
f61838dfc1 Torg*1629 C    the solidified grease ice volume needs to be added to HEFF:
8377b8ee87 Mart*1630            SItrBucket(i,j,bi,bj,iTrGrease) =
                1631      &             SItrBucket(i,j,bi,bj,iTrGrease)
                1632      &           + d_HEFFbyGREASE(i,j)
f61838dfc1 Torg*1633 C
                1634 C 3) grease ice growth from new frazil ice:
                1635 C
c6b168144e Jean*1636            SItracer(i,j,bi,bj,iTrGrease) =
f61838dfc1 Torg*1637      &      SItracer(i,j,bi,bj,iTrGrease) + greaseNewFrazil
8377b8ee87 Mart*1638           ENDIF
f61838dfc1 Torg*1639 C 4) mapping SItrBucket to external variable, in this case HEFF, ...
8377b8ee87 Mart*1640           tmpscal3=tmpscal3+SItrBucket(i,j,bi,bj,iTrGrease)
f61838dfc1 Torg*1641 C    ... and empty SItrBucket for tracer 'grease'
8377b8ee87 Mart*1642           SItrBucket(i,j,bi,bj,iTrGrease)=0. _d 0
f61838dfc1 Torg*1643 #endif /* SEAICE_GREASE */
286983d3d2 Patr*1644 #ifdef SEAICE_ITD
                1645 C         ice growth in open water adds to first category
8377b8ee87 Mart*1646           d_HEFFbyATMonOCN_open_ITD(i,j,1)=tmpscal3
                1647           d_HEFFbyATMonOCN_ITD(i,j,1)     =d_HEFFbyATMonOCN_ITD(i,j,1)
286983d3d2 Patr*1648      &                                    +tmpscal3
cbd0ee24a8 Mart*1649 #endif /* SEAICE_ITD */
8377b8ee87 Mart*1650           d_HEFFbyATMonOCN_open(i,j)=tmpscal3
                1651           d_HEFFbyATMonOCN(i,j)=d_HEFFbyATMonOCN(i,j)+tmpscal3
                1652           r_QbyATM_open(i,j)=r_QbyATM_open(i,j)-tmpscal3
                1653           HEFF(i,j,bi,bj) = HEFF(i,j,bi,bj) + tmpscal3
85586adda4 Gael*1654          ENDDO
                1655         ENDDO
                1656
f50f58ec54 Gael*1657 #ifdef ALLOW_SITRACER
8377b8ee87 Mart*1658         DO j=1,sNy
                1659          DO i=1,sNx
1d74e34c65 Jean*1660 C needs to be here to allow use also with LEGACY branch
286983d3d2 Patr*1661 #ifdef SEAICE_ITD
8377b8ee87 Mart*1662           SItrHEFF(i,j,bi,bj,4)=SItrHEFF(i,j,bi,bj,3)
                1663      &                         +d_HEFFbyATMonOCN_open(i,j)
cbd0ee24a8 Mart*1664 #else /* ndef SEAICE_ITD */
8377b8ee87 Mart*1665           SItrHEFF(i,j,bi,bj,4)=HEFF(i,j,bi,bj)
cbd0ee24a8 Mart*1666 #endif /* SEAICE_ITD */
f50f58ec54 Gael*1667          ENDDO
                1668         ENDDO
a73db480d4 Jean*1669 #endif /* ALLOW_SITRACER */
f50f58ec54 Gael*1670
581175eaf0 Gael*1671 C convert snow to ice if submerged.
                1672 C =================================
                1673
2651ba3350 Jean*1674 C note: in legacy, this process is done at the end
3a3bf6419a Gael*1675 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1676 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                1677 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
3a3bf6419a Gael*1678 #endif /* ALLOW_AUTODIFF_TAMC */
581175eaf0 Gael*1679         IF ( SEAICEuseFlooding ) THEN
286983d3d2 Patr*1680 #ifdef SEAICE_ITD
f913c5a485 Mart*1681          DO IT=1,SEAICE_multDim
8377b8ee87 Mart*1682           DO j=1,sNy
                1683            DO i=1,sNx
                1684             tmpscal3 = HEFFITDpreTH(i,j,IT)
                1685      &               + d_HEFFbySublim_ITD(i,j,IT)
                1686      &               + d_HEFFbyOCNonICE_ITD(i,j,IT)
                1687      &               + d_HEFFbyATMonOCN_ITD(i,j,IT)
                1688             tmpscal4 = HSNWITDpreTH(i,j,IT)
                1689      &               + d_HSNWbySublim_ITD(i,j,IT)
                1690      &               + d_HSNWbyATMonSNW_ITD(i,j,IT)
                1691      &               + d_HSNWbyRAIN_ITD(i,j,IT)
286983d3d2 Patr*1692             tmpscal0 = (tmpscal4*SEAICE_rhoSnow
                1693      &               +  tmpscal3*SEAICE_rhoIce)
                1694      &               * recip_rhoConst
                1695             tmpscal1 = MAX( 0. _d 0, tmpscal0 - tmpscal3)
8377b8ee87 Mart*1696             d_HEFFbyFLOODING_ITD(i,j,IT) = tmpscal1
                1697             d_HEFFbyFLOODING(i,j) = d_HEFFbyFLOODING(i,j)  + tmpscal1
114c791332 Jean*1698            ENDDO
                1699           ENDDO
                1700          ENDDO
cbd0ee24a8 Mart*1701 #else /* ndef SEAICE_ITD */
8377b8ee87 Mart*1702          DO j=1,sNy
                1703           DO i=1,sNx
                1704            tmpscal0 = (HSNOW(i,j,bi,bj)*SEAICE_rhoSnow
                1705      &              +HEFF(i,j,bi,bj)*SEAICE_rhoIce)*recip_rhoConst
                1706            tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFF(i,j,bi,bj))
                1707            d_HEFFbyFLOODING(i,j)=tmpscal1
                1708            HEFF(i,j,bi,bj) = HEFF(i,j,bi,bj)+d_HEFFbyFLOODING(i,j)
                1709            HSNOW(i,j,bi,bj) = HSNOW(i,j,bi,bj)-
                1710      &                           d_HEFFbyFLOODING(i,j)*ICE2SNOW
581175eaf0 Gael*1711           ENDDO
                1712          ENDDO
cbd0ee24a8 Mart*1713 #endif /* SEAICE_ITD */
581175eaf0 Gael*1714         ENDIF
050eb90cc6 Gael*1715
286983d3d2 Patr*1716 #ifdef SEAICE_ITD
                1717 C apply ice and snow thickness changes
                1718 C =================================================================
f913c5a485 Mart*1719          DO IT=1,SEAICE_multDim
8377b8ee87 Mart*1720           DO j=1,sNy
                1721            DO i=1,sNx
                1722             HEFFITD(i,j,IT,bi,bj) = HEFFITD(i,j,IT,bi,bj)
                1723      &                            + d_HEFFbySublim_ITD(i,j,IT)
                1724      &                            + d_HEFFbyOCNonICE_ITD(i,j,IT)
                1725      &                            + d_HEFFbyATMonOCN_ITD(i,j,IT)
                1726      &                            + d_HEFFbyFLOODING_ITD(i,j,IT)
                1727             HSNOWITD(i,j,IT,bi,bj) = HSNOWITD(i,j,IT,bi,bj)
                1728      &                            + d_HSNWbySublim_ITD(i,j,IT)
                1729      &                            + d_HSNWbyATMonSNW_ITD(i,j,IT)
                1730      &                            + d_HSNWbyRAIN_ITD(i,j,IT)
                1731      &                            + d_HSNWbyOCNonSNW_ITD(i,j,IT)
                1732      &                            - d_HEFFbyFLOODING_ITD(i,j,IT)
286983d3d2 Patr*1733      &                            * ICE2SNOW
114c791332 Jean*1734            ENDDO
                1735           ENDDO
                1736          ENDDO
cbd0ee24a8 Mart*1737 #endif /* SEAICE_ITD */
581175eaf0 Gael*1738
2651ba3350 Jean*1739 C ===================================================================
                1740 C ==========PART 4: determine ice cover fraction increments=========-
                1741 C ===================================================================
2ae913cfea Gael*1742
33e17487ce Dimi*1743 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1744 CADJ STORE d_HEFFbyATMonOCN = comlev1_bibj,key=tkey,byte=isbyte
                1745 CADJ STORE d_HEFFbyATMonOCN_cover = comlev1_bibj,key=tkey,byte=isbyte
                1746 CADJ STORE d_HEFFbyATMonOCN_open = comlev1_bibj,key=tkey,byte=isbyte
                1747 CADJ STORE d_HEFFbyOCNonICE = comlev1_bibj,key=tkey,byte=isbyte
                1748 CADJ STORE recip_heffActual = comlev1_bibj,key=tkey,byte=isbyte
                1749 CADJ STORE d_hsnwbyatmonsnw = comlev1_bibj,key=tkey,byte=isbyte
4fcdd931bd Patr*1750 cph(
edb6656069 Mart*1751 cphCADJ STORE d_AREAbyATM  = comlev1_bibj,key=tkey,byte=isbyte
                1752 cphCADJ STORE d_AREAbyICE  = comlev1_bibj,key=tkey,byte=isbyte
                1753 cphCADJ STORE d_AREAbyOCN  = comlev1_bibj,key=tkey,byte=isbyte
4fcdd931bd Patr*1754 cph)
edb6656069 Mart*1755 CADJ STORE a_QbyATM_open = comlev1_bibj,key=tkey,byte=isbyte
                1756 CADJ STORE heffActual = comlev1_bibj,key=tkey,byte=isbyte
                1757 CADJ STORE AREApreTH = comlev1_bibj,key=tkey,byte=isbyte
                1758 CADJ STORE HEFF(:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                1759 CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                1760 CADJ STORE AREA(:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
33e17487ce Dimi*1761 #endif /* ALLOW_AUTODIFF_TAMC */
                1762
286983d3d2 Patr*1763 #ifdef SEAICE_ITD
c6b168144e Jean*1764 C--   in thinnest category account for lateral ice growth and melt the
f913c5a485 Mart*1765 C--   "non-ITD" way, so that the ITD simulation with SEAICE_multDim=1
                1766 C--   is identical with the non-ITD simulation;
cbd0ee24a8 Mart*1767 C--   use HEFF, ARE, HSNOW, etc. as temporal storage for 1st category
8377b8ee87 Mart*1768         DO j=1,sNy
                1769          DO i=1,sNx
                1770           HEFF(i,j,bi,bj)=HEFFITD(i,j,1,bi,bj)
                1771           AREA(i,j,bi,bj)=AREAITD(i,j,1,bi,bj)
                1772           HSNOW(i,j,bi,bj)=HSNOWITD(i,j,1,bi,bj)
                1773           HEFFpreTH(i,j)=HEFFITDpreTH(i,j,1)
                1774           AREApreTH(i,j)=AREAITDpreTH(i,j,1)
                1775           recip_heffActual(i,j)=recip_heffActualMult(i,j,1)
114c791332 Jean*1776          ENDDO
                1777         ENDDO
cbd0ee24a8 Mart*1778 #endif /* SEAICE_ITD */
8377b8ee87 Mart*1779         DO j=1,sNy
                1780          DO i=1,sNx
2651ba3350 Jean*1781
56d13a40ed Mart*1782 #if defined ( ALLOW_SITRACER ) && defined ( SEAICE_GREASE )
c6b168144e Jean*1783 C         grease ice layer thickness (Hg) includes
f61838dfc1 Torg*1784 C         1 part frazil ice and 3 parts sea water
                1785 C         i.e. HO = 0.25 * Hg
8377b8ee87 Mart*1786           recip_HO=4. _d 0 / greaseLayerThick(i,j)
f61838dfc1 Torg*1787 #else /* SEAICE_GREASE */
8377b8ee87 Mart*1788           IF ( YC(i,j,bi,bj) .LT. ZERO ) THEN
6ec4646d60 Gael*1789            recip_HO=1. _d 0 / HO_south
c7b18bd1fe Gael*1790           ELSE
6ec4646d60 Gael*1791            recip_HO=1. _d 0 / HO
2afe30fba0 Dimi*1792           ENDIF
f61838dfc1 Torg*1793 #endif /* SEAICE_GREASE */
8377b8ee87 Mart*1794           recip_HH = recip_heffActual(i,j)
2afe30fba0 Dimi*1795
a4bc0a0b4c Jean*1796 C gain of ice over open water : computed from
56d13a40ed Mart*1797 #if defined ( ALLOW_SITRACER ) && defined ( SEAICE_GREASE )
f61838dfc1 Torg*1798 C   from growth by ATM with grease ice time delay
8377b8ee87 Mart*1799           tmpscal4 = MAX(ZERO,d_HEFFbyGREASE(i,j))
f61838dfc1 Torg*1800 #else /* SEAICE_GREASE */
6ec4646d60 Gael*1801 C   (SEAICE_areaGainFormula.EQ.1) from growth by ATM
                1802 C   (SEAICE_areaGainFormula.EQ.2) from predicted growth by ATM
4eb4a54cba Jean*1803           IF (SEAICE_areaGainFormula.EQ.1) THEN
8377b8ee87 Mart*1804             tmpscal4 = MAX(ZERO,d_HEFFbyATMonOCN_open(i,j))
4eb4a54cba Jean*1805           ELSE
8377b8ee87 Mart*1806             tmpscal4=MAX(ZERO,a_QbyATM_open(i,j))
4eb4a54cba Jean*1807           ENDIF
f61838dfc1 Torg*1808 #endif /* SEAICE_GREASE */
65f34462d4 Gael*1809
a4bc0a0b4c Jean*1810 C loss of ice cover by melting : computed from
6ec4646d60 Gael*1811 C   (SEAICE_areaLossFormula.EQ.1) from all but only melt conributions by ATM and OCN
                1812 C   (SEAICE_areaLossFormula.EQ.2) from net melt-growth>0 by ATM and OCN
                1813 C   (SEAICE_areaLossFormula.EQ.3) from predicted melt by ATM
4eb4a54cba Jean*1814           IF (SEAICE_areaLossFormula.EQ.1) THEN
8377b8ee87 Mart*1815             tmpscal3 = MIN( 0. _d 0 , d_HEFFbyATMonOCN_cover(i,j) )
                1816      &        + MIN( 0. _d 0 , d_HEFFbyATMonOCN_open(i,j) )
                1817      &        + MIN( 0. _d 0 , d_HEFFbyOCNonICE(i,j) )
4eb4a54cba Jean*1818           ELSEIF (SEAICE_areaLossFormula.EQ.2) THEN
8377b8ee87 Mart*1819             tmpscal3 = MIN( 0. _d 0 , d_HEFFbyATMonOCN_cover(i,j)
                1820      &        + d_HEFFbyATMonOCN_open(i,j) + d_HEFFbyOCNonICE(i,j) )
4eb4a54cba Jean*1821           ELSE
6ec4646d60 Gael*1822 C           compute heff after ice melt by ocn:
8377b8ee87 Mart*1823             tmpscal0=HEFF(i,j,bi,bj) - d_HEFFbyATMonOCN(i,j)
6ec4646d60 Gael*1824 C           compute available heat left after snow melt by atm:
8377b8ee87 Mart*1825             tmpscal1= a_QbyATM_open(i,j)+a_QbyATM_cover(i,j)
                1826      &            - d_HSNWbyATMonSNW(i,j)*SNOW2ICE
6ec4646d60 Gael*1827 C           could not melt more than all the ice
                1828             tmpscal2 = MAX(-tmpscal0,tmpscal1)
                1829             tmpscal3 = MIN(ZERO,tmpscal2)
4eb4a54cba Jean*1830           ENDIF
a4bc0a0b4c Jean*1831
2651ba3350 Jean*1832 C apply tendency
85586adda4 Gael*1833           IF ( (HEFF(i,j,bi,bj).GT.0. _d 0).OR.
                1834      &        (HSNOW(i,j,bi,bj).GT.0. _d 0) ) THEN
8377b8ee87 Mart*1835            AREA(i,j,bi,bj)=MAX(0. _d 0,
                1836      &      MIN( SEAICE_area_max, AREA(i,j,bi,bj)
4b6d456764 Mart*1837      &       + recip_HO*tmpscal4+HALF*recip_HH*tmpscal3
                1838      &       * areaPDFfac ))
85586adda4 Gael*1839           ELSE
8377b8ee87 Mart*1840            AREA(i,j,bi,bj)=0. _d 0
85586adda4 Gael*1841           ENDIF
bb24b8a3e6 Gael*1842 #ifdef ALLOW_SITRACER
8377b8ee87 Mart*1843           SItrAREA(i,j,bi,bj,3)=AREA(i,j,bi,bj)
a73db480d4 Jean*1844 #endif /* ALLOW_SITRACER */
381adf77df Gael*1845 #ifdef ALLOW_DIAGNOSTICS
8377b8ee87 Mart*1846           d_AREAbyATM(i,j)=
                1847      &       recip_HO*MAX(ZERO,d_HEFFbyATMonOCN_open(i,j))
                1848      &       +HALF*recip_HH*MIN(0. _d 0,d_HEFFbyATMonOCN_open(i,j))
4b6d456764 Mart*1849      &       *areaPDFfac
8377b8ee87 Mart*1850           d_AREAbyICE(i,j)=
                1851      &        HALF*recip_HH*MIN(0. _d 0,d_HEFFbyATMonOCN_cover(i,j))
4b6d456764 Mart*1852      &       *areaPDFfac
8377b8ee87 Mart*1853           d_AREAbyOCN(i,j)=
                1854      &        HALF*recip_HH*MIN( 0. _d 0,d_HEFFbyOCNonICE(i,j) )
4b6d456764 Mart*1855      &       *areaPDFfac
a73db480d4 Jean*1856 #endif /* ALLOW_DIAGNOSTICS */
33e17487ce Dimi*1857          ENDDO
                1858         ENDDO
286983d3d2 Patr*1859 #ifdef SEAICE_ITD
                1860 C       transfer 1st category values back into ITD variables
8377b8ee87 Mart*1861         DO j=1,sNy
                1862          DO i=1,sNx
                1863           HEFFITD(i,j,1,bi,bj)=HEFF(i,j,bi,bj)
                1864           AREAITD(i,j,1,bi,bj)=AREA(i,j,bi,bj)
                1865           HSNOWITD(i,j,1,bi,bj)=HSNOW(i,j,bi,bj)
286983d3d2 Patr*1866          ENDDO
                1867         ENDDO
f322f85e9b Torg*1868 C       now melt ice laterally in all other thickness categories
                1869 C       (areal growth, i.e. new ice formation, only occurrs in 1st category)
8377b8ee87 Mart*1870         IF (SEAICE_multDim .GT. 1) THEN
f913c5a485 Mart*1871          DO IT=2,SEAICE_multDim
8377b8ee87 Mart*1872           DO j=1,sNy
                1873            DO i=1,sNx
                1874             IF (HEFFITD(i,j,IT,bi,bj).LE.ZERO) THEN
f322f85e9b Torg*1875 C       when thickness is zero, area should be zero, too:
8377b8ee87 Mart*1876              AREAITD(i,j,IT,bi,bj)=ZERO
f322f85e9b Torg*1877             ELSE
3da5adc21e Mart*1878 C     tmpscal1 is the minimal ice concentration after lateral melt that will
                1879 C     not lead to an unphysical increase of ice thickness by lateral melt;
a24915ab1a Jean*1880 C     estimated as the concentration before thermodynamics scaled by the
3da5adc21e Mart*1881 C     ratio of new ice thickness and ice thickness before thermodynamics
8377b8ee87 Mart*1882              IF ( HEFFITDpreTH(i,j,IT).LE.ZERO ) THEN
3da5adc21e Mart*1883               tmpscal1=0. _d 0
                1884              ELSE
8377b8ee87 Mart*1885               tmpscal1=AREAITDpreTH(i,j,IT)*
                1886      &             HEFFITD(i,j,IT,bi,bj)/HEFFITDpreTH(i,j,IT)
3da5adc21e Mart*1887              ENDIF
f322f85e9b Torg*1888 C       melt ice laterally based on an average floe sice
                1889 C       following Steele (1992)
8377b8ee87 Mart*1890              AREAITD(i,j,IT,bi,bj) = AREAITD(i,j,IT,bi,bj)
                1891      &                             * (ONE - latMeltFrac(i,j,IT))
3da5adc21e Mart*1892 CML   not necessary:
                1893 CML          AREAITD(I,J,IT,bi,bj) = max(ZERO,AREAITD(I,J,IT,bi,bj))
f322f85e9b Torg*1894 C       limit area reduction so that actual ice thickness does not increase
8377b8ee87 Mart*1895              AREAITD(i,j,IT,bi,bj) = max(AREAITD(i,j,IT,bi,bj),
3da5adc21e Mart*1896      &                                   tmpscal1)
f322f85e9b Torg*1897             ENDIF
                1898 #ifdef ALLOW_SITRACER
8377b8ee87 Mart*1899             SItrAREA(i,j,bi,bj,3)=SItrAREA(i,j,bi,bj,3)
                1900      &                           +AREAITD(i,j,IT,bi,bj)
f322f85e9b Torg*1901 #endif /* ALLOW_SITRACER */
                1902            ENDDO
                1903           ENDDO
                1904          ENDDO
                1905         ENDIF
cbd0ee24a8 Mart*1906 #endif /* SEAICE_ITD */
33e17487ce Dimi*1907
8377b8ee87 Mart*1908 #if (defined ALLOW_AUTODIFF && defined SEAICE_MODIFY_GROWTH_ADJ)
83ad492c2d Jean*1909 Cgf 'bulk' linearization of area=f(HEFF)
4bc8f4264e Mart*1910         IF ( SEAICEadjMODE.GE.1 ) THEN
286983d3d2 Patr*1911 #ifdef SEAICE_ITD
f913c5a485 Mart*1912          DO IT=1,SEAICE_multDim
8377b8ee87 Mart*1913           DO j=1,sNy
                1914            DO i=1,sNx
                1915             AREAITD(i,j,IT,bi,bj) = AREAITDpreTH(i,j,IT) + 0.1 _d 0 *
                1916      &               ( HEFFITD(i,j,IT,bi,bj) - HEFFITDpreTH(i,j,IT) )
286983d3d2 Patr*1917            ENDDO
                1918           ENDDO
                1919          ENDDO
cbd0ee24a8 Mart*1920 #else /* ndef SEAICE_ITD */
8377b8ee87 Mart*1921          DO j=1,sNy
                1922           DO i=1,sNx
3a3bf6419a Gael*1923 C            AREA(I,J,bi,bj) = 0.1 _d 0 * HEFF(I,J,bi,bj)
8377b8ee87 Mart*1924            AREA(i,j,bi,bj) = AREApreTH(i,j) + 0.1 _d 0 *
                1925      &               ( HEFF(i,j,bi,bj) - HEFFpreTH(i,j) )
4bc8f4264e Mart*1926           ENDDO
3a3bf6419a Gael*1927          ENDDO
cbd0ee24a8 Mart*1928 #endif /* SEAICE_ITD */
4bc8f4264e Mart*1929         ENDIF
3a3bf6419a Gael*1930 #endif
286983d3d2 Patr*1931 #ifdef SEAICE_ITD
5df73465ef Torg*1932 C check categories for consistency with limits after growth/melt ...
ed2f6fecc4 Mart*1933         IF ( SEAICEuseLinRemapITD ) CALL SEAICE_ITD_REMAP(
                1934      I     heffitdPreTH, areaitdPreTH,
                1935      I     bi, bj, myTime, myIter, myThid )
5df73465ef Torg*1936         CALL SEAICE_ITD_REDIST(bi, bj, myTime, myIter, myThid)
                1937 C ... and update total AREA, HEFF, HSNOW
                1938 C     (the updated HEFF is used below for ice salinity increments)
                1939         CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)
cbd0ee24a8 Mart*1940 #endif /* SEAICE_ITD */
56d13a40ed Mart*1941 #if defined ( ALLOW_SITRACER ) && defined ( SEAICE_GREASE )
f61838dfc1 Torg*1942 C convert SItracer 'grease' from grease ice volume back to ratio:
                1943 C ===============================================================
8377b8ee87 Mart*1944         DO j=1,sNy
                1945          DO i=1,sNx
                1946           IF (HEFF(i,j,bi,bj).GT.siEps) THEN
                1947            SItracer(i,j,bi,bj,iTrGrease) =
                1948      &      SItracer(i,j,bi,bj,iTrGrease) / HEFF(i,j,bi,bj)
                1949           ELSE
                1950            SItracer(i,j,bi,bj,iTrGrease) = 0. _d 0
                1951           ENDIF
f61838dfc1 Torg*1952          ENDDO
                1953         ENDDO
                1954 #endif /* SEAICE_GREASE */
3a3bf6419a Gael*1955
2651ba3350 Jean*1956 C ===================================================================
                1957 C =============PART 5: determine ice salinity increments=============
                1958 C ===================================================================
aea7db20a6 Gael*1959
a98c4b8072 Ian *1960 #ifndef SEAICE_VARIABLE_SALINITY
8120fff0c1 Mart*1961 # ifdef ALLOW_AUTODIFF_TAMC
1cf549c217 Mart*1962 CADJ STORE d_HEFFbyNEG(:,:,bi,bj) = comlev1_bibj,
edb6656069 Mart*1963 CADJ &                              key = tkey, byte = isbyte
8120fff0c1 Mart*1964 #  ifdef ALLOW_SALT_PLUME
edb6656069 Mart*1965 CADJ STORE d_HEFFbyOCNonICE = comlev1_bibj,key=tkey,byte=isbyte
                1966 CADJ STORE d_HEFFbyATMonOCN = comlev1_bibj,key=tkey,byte=isbyte
                1967 CADJ STORE d_HEFFbyATMonOCN_open = comlev1_bibj,key=tkey,byte=isbyte
                1968 CADJ STORE d_HEFFbyATMonOCN_cover = comlev1_bibj,key=tkey,byte=isbyte
                1969 CADJ STORE d_HEFFbyFLOODING = comlev1_bibj,key=tkey,byte=isbyte
                1970 CADJ STORE d_HEFFbySublim = comlev1_bibj,key=tkey,byte=isbyte
6b89e1b973 Gael*1971 CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
edb6656069 Mart*1972 CADJ &                       key = tkey, byte = isbyte
8120fff0c1 Mart*1973 #  endif /* ALLOW_SALT_PLUME */
                1974 # endif /* ALLOW_AUTODIFF_TAMC */
8377b8ee87 Mart*1975         DO j=1,sNy
                1976          DO i=1,sNx
                1977           tmpscal1 = d_HEFFbyNEG(i,j,bi,bj) + d_HEFFbyOCNonICE(i,j) +
                1978      &               d_HEFFbyATMonOCN(i,j) + d_HEFFbyFLOODING(i,j)
                1979      &             + d_HEFFbySublim(i,j)
6510a54854 Jean*1980 #ifdef EXF_SEAICE_FRACTION
8377b8ee87 Mart*1981      &             + d_HEFFbyRLX(i,j,bi,bj)
d32fe07ad8 Patr*1982 #endif
cafd3818b7 An T*1983 Catn: can not take out more that surface salinity when SSS<SEAICE_salt0
                1984           tmpscal3 = max( 0. _d 0,
8377b8ee87 Mart*1985      &                    min(SEAICE_salt0,salt(i,j,kSurface,bi,bj)) )
                1986           tmpscal2 = tmpscal1 * tmpscal3 * HEFFM(i,j,bi,bj)
bb8e6379cb Mart*1987      &            * recip_deltaTtherm * SEAICE_rhoIce
8377b8ee87 Mart*1988           saltFlux(i,j,bi,bj) = tmpscal2
8a423a90df Gael*1989 #ifdef ALLOW_SALT_PLUME
ef6001d9aa An T*1990 #ifdef SALT_PLUME_SPLIT_BASIN
                1991 catn attempt to split East/West basins in Arctic
8377b8ee87 Mart*1992           localSPfrac(i,j) = SPsalFRAC(1)
ef6001d9aa An T*1993           IF ( SaltPlumeSplitBasin ) THEN
8377b8ee87 Mart*1994             localSPfrac(i,j) = SPsalFRAC(2)
                1995             IF(YC(i,j,bi,bj).LT. 85.0 .AND. YC(i,j,bi,bj).GT. 71.0
                1996      &       .AND. XC(i,j,bi,bj) .LT. -90.0) THEN
                1997              localSPfrac(i,j) = SPsalFRAC(1)
ef6001d9aa An T*1998             ENDIF
                1999           ENDIF
                2000 #else
8377b8ee87 Mart*2001           localSPfrac(i,j) = SPsalFRAC
ef6001d9aa An T*2002 #endif /* SALT_PLUME_SPLIT_BASIN */
cafd3818b7 An T*2003 #ifdef SALT_PLUME_IN_LEADS
                2004 Catn: Only d_HEFFbyATMonOCN should contribute to plume.
                2005 C     By redefining tmpscal1 here, saltPlumeFlux is smaller in case
                2006 C     define inLeads than case undef inLeads.  Physical interpretation
                2007 C     is that when d_HEFF is formed from below via ocean freezing, it
                2008 C     occurs more uniform over grid cell and not inLeads, thus not
                2009 C     participating in pkg/salt_plume.
                2010 C     Note: tmpscal1 is defined only after saltFlux is calculated.
8377b8ee87 Mart*2011           IceGrowthRateInLeads(i,j)=max(0. _d 0,d_HEFFbyATMonOCN(i,j))
                2012           tmpscal1 = IceGrowthRateInLeads(i,j)
                2013           leadPlumeFraction(i,j) =
                2014      &      (ONE + EXP( (SPinflectionPoint - AREApreTH(i,j))*5.0
cafd3818b7 An T*2015      &                 /(ONE - SPinflectionPoint) ))**(-ONE)
8377b8ee87 Mart*2016           localSPfrac(i,j)=localSPfrac(i,j)*leadPlumeFraction(i,j)
cafd3818b7 An T*2017 #endif /* SALT_PLUME_IN_LEADS */
8377b8ee87 Mart*2018           tmpscal3 = tmpscal1*salt(i,j,kSurface,bi,bj)*HEFFM(i,j,bi,bj)
bb8e6379cb Mart*2019      &            * recip_deltaTtherm * SEAICE_rhoIce
8377b8ee87 Mart*2020           saltPlumeFlux(i,j,bi,bj) = MAX( tmpscal3-tmpscal2 , 0. _d 0)
                2021      &            *localSPfrac(i,j)
cafd3818b7 An T*2022 C if SaltPlumeSouthernOcean=.FALSE. turn off salt plume in Southern Ocean
                2023           IF ( .NOT. SaltPlumeSouthernOcean ) THEN
8377b8ee87 Mart*2024            IF ( YC(i,j,bi,bj) .LT. 0.0 _d 0 )
cafd3818b7 An T*2025      &          saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
                2026           ENDIF
8a423a90df Gael*2027 #endif /* ALLOW_SALT_PLUME */
                2028          ENDDO
                2029         ENDDO
a73db480d4 Jean*2030 #endif /* ndef SEAICE_VARIABLE_SALINITY */
8a423a90df Gael*2031
a98c4b8072 Ian *2032 #ifdef SEAICE_VARIABLE_SALINITY
33e17487ce Dimi*2033
                2034 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*2035 CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
33e17487ce Dimi*2036 #endif /* ALLOW_AUTODIFF_TAMC */
                2037
8377b8ee87 Mart*2038         DO j=1,sNy
                2039          DO i=1,sNx
2651ba3350 Jean*2040 C sum up the terms that affect the salt content of the ice pack
8377b8ee87 Mart*2041          tmpscal1=d_HEFFbyOCNonICE(i,j)+d_HEFFbyATMonOCN(i,j)
2afe30fba0 Dimi*2042
4eb4a54cba Jean*2043 C recompute HEFF before thermodynamic updates (which is not AREApreTH in legacy code)
8377b8ee87 Mart*2044          tmpscal2=HEFF(i,j,bi,bj)-tmpscal1-d_HEFFbyFLOODING(i,j)
65b7f51792 Gael*2045 C tmpscal1 > 0 : m of sea ice that is created
                2046           IF ( tmpscal1 .GE. 0.0 ) THEN
8377b8ee87 Mart*2047              saltFlux(i,j,bi,bj) =
                2048      &            HEFFM(i,j,bi,bj)*recip_deltaTtherm
                2049      &            *SEAICE_saltFrac*salt(i,j,kSurface,bi,bj)
3d15111920 Ian *2050      &            *tmpscal1*SEAICE_rhoIce
33e17487ce Dimi*2051 #ifdef ALLOW_SALT_PLUME
ef6001d9aa An T*2052 #ifdef SALT_PLUME_SPLIT_BASIN
                2053 catn attempt to split East/West basins in Arctic
8377b8ee87 Mart*2054              localSPfrac(i,j) = SPsalFRAC(1)
ef6001d9aa An T*2055              IF ( SaltPlumeSplitBasin ) THEN
8377b8ee87 Mart*2056                localSPfrac(i,j) = SPsalFRAC(2)
                2057                IF(YC(i,j,bi,bj).LT. 85.0 .AND. YC(i,j,bi,bj).GT. 71.0
                2058      &                  .AND. XC(i,j,bi,bj) .LT. -90.0) THEN
                2059                  localSPfrac(i,j) = SPsalFRAC(1)
ef6001d9aa An T*2060                ENDIF
                2061              ENDIF
                2062 #else
8377b8ee87 Mart*2063              localSPfrac(i,j) = SPsalFRAC
ef6001d9aa An T*2064 #endif /* SALT_PLUME_SPLIT_BASIN */
cafd3818b7 An T*2065 #ifndef SALT_PLUME_IN_LEADS
33e17487ce Dimi*2066 C saltPlumeFlux is defined only during freezing:
8377b8ee87 Mart*2067              saltPlumeFlux(i,j,bi,bj)=
                2068      &            HEFFM(i,j,bi,bj)*recip_deltaTtherm
                2069      &            *(ONE-SEAICE_saltFrac)*salt(i,j,kSurface,bi,bj)
3d15111920 Ian *2070      &            *tmpscal1*SEAICE_rhoIce
8377b8ee87 Mart*2071      &            *localSPfrac(i,j)
cafd3818b7 An T*2072 #endif /* ndef SALT_PLUME_IN_LEADS */
33e17487ce Dimi*2073 #endif /* ALLOW_SALT_PLUME */
2651ba3350 Jean*2074
65b7f51792 Gael*2075 C tmpscal1 < 0 : m of sea ice that is melted
33e17487ce Dimi*2076           ELSE
8377b8ee87 Mart*2077              saltFlux(i,j,bi,bj) =
                2078      &         HEFFM(i,j,bi,bj)*recip_deltaTtherm
                2079      &         *HSALT(i,j,bi,bj)
581175eaf0 Gael*2080      &         *tmpscal1/tmpscal2
33e17487ce Dimi*2081 #ifdef ALLOW_SALT_PLUME
cafd3818b7 An T*2082 #ifndef SALT_PLUME_IN_LEADS
33e17487ce Dimi*2083              saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
cafd3818b7 An T*2084 #endif /* ndef SALT_PLUME_IN_LEADS */
33e17487ce Dimi*2085 #endif /* ALLOW_SALT_PLUME */
                2086           ENDIF
cafd3818b7 An T*2087
                2088 #ifdef ALLOW_SALT_PLUME
                2089 #ifdef SALT_PLUME_IN_LEADS
                2090 Catn: only d_HEFFbyATMonOCN should contribute to plume
                2091 C     By redefining tmpscal1 here, saltPlumeFlux is smaller in case
                2092 C     define inLeads than case undef inLeads.  Physical interpretation
                2093 C     is that when d_HEFF is formed from below via ocean freezing, it
                2094 C     occurs more uniform over grid cell and not inLeads, thus not
                2095 C     participating in pkg/salt_plume.
                2096 C     Note: tmpscal1 is defined only after saltFlux is calculated.
8377b8ee87 Mart*2097           IceGrowthRateInLeads(i,j)=max(0. _d 0,d_HEFFbyATMonOCN(i,j))
                2098           tmpscal1 = IceGrowthRateInLeads(i,j)
                2099           leadPlumeFraction(i,j) =
                2100      &      (ONE + EXP( (SPinflectionPoint - AREApreTH(i,j))*5.0
cafd3818b7 An T*2101      &                 /(ONE - SPinflectionPoint) ))**(-ONE)
8377b8ee87 Mart*2102           localSPfrac(i,j)=localSPfrac(i,j)*leadPlumeFraction(i,j)
cafd3818b7 An T*2103           IF ( tmpscal1 .GE. 0.0) THEN
                2104 C saltPlumeFlux is defined only during freezing:
8377b8ee87 Mart*2105              saltPlumeFlux(i,j,bi,bj)=
                2106      &            HEFFM(i,j,bi,bj)*recip_deltaTtherm
                2107      &            *(ONE-SEAICE_saltFrac)*salt(i,j,kSurface,bi,bj)
cafd3818b7 An T*2108      &            *tmpscal1*SEAICE_rhoIce
8377b8ee87 Mart*2109      &            *localSPfrac(i,j)
cafd3818b7 An T*2110           ELSE
8377b8ee87 Mart*2111              saltPlumeFlux(i,j,bi,bj) = 0. _d 0
cafd3818b7 An T*2112           ENDIF
                2113 #endif /* SALT_PLUME_IN_LEADS */
                2114 C if SaltPlumeSouthernOcean=.FALSE. turn off salt plume in Southern Ocean
                2115           IF ( .NOT. SaltPlumeSouthernOcean ) THEN
8377b8ee87 Mart*2116            IF ( YC(i,j,bi,bj) .LT. 0.0 _d 0 )
cafd3818b7 An T*2117      &          saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
                2118           ENDIF
                2119 #endif /* ALLOW_SALT_PLUME */
33e17487ce Dimi*2120 C update HSALT based on surface saltFlux
8377b8ee87 Mart*2121           HSALT(i,j,bi,bj) = HSALT(i,j,bi,bj) +
                2122      &         saltFlux(i,j,bi,bj) * SEAICE_deltaTtherm
                2123           saltFlux(i,j,bi,bj) =
                2124      &         saltFlux(i,j,bi,bj) + saltFluxAdjust(i,j,bi,bj)
33e17487ce Dimi*2125          ENDDO
                2126         ENDDO
a98c4b8072 Ian *2127 #endif /* SEAICE_VARIABLE_SALINITY */
33e17487ce Dimi*2128
f50f58ec54 Gael*2129 #ifdef ALLOW_SITRACER
8377b8ee87 Mart*2130         DO j=1,sNy
                2131          DO i=1,sNx
1d74e34c65 Jean*2132 C needs to be here to allow use also with LEGACY branch
8377b8ee87 Mart*2133           SItrHEFF(i,j,bi,bj,5)=HEFF(i,j,bi,bj)
f50f58ec54 Gael*2134          ENDDO
                2135         ENDDO
a73db480d4 Jean*2136 #endif /* ALLOW_SITRACER */
2ae913cfea Gael*2137
2651ba3350 Jean*2138 C ===================================================================
                2139 C ==============PART 7: determine ocean model forcing================
                2140 C ===================================================================
aea7db20a6 Gael*2141
0c0ecd4c7b Jean*2142 C compute net heat flux leaving/entering the ocean,
aea7db20a6 Gael*2143 C accounting for the part used in melt/freeze processes
                2144 C =====================================================
                2145
286983d3d2 Patr*2146 #ifdef SEAICE_ITD
                2147 C compute total of "mult" fluxes for ocean forcing
8377b8ee87 Mart*2148         DO j=1,sNy
                2149          DO i=1,sNx
                2150           a_QbyATM_cover(i,j)   = 0.0 _d 0
                2151           r_QbyATM_cover(i,j)   = 0.0 _d 0
                2152           a_QSWbyATM_cover(i,j) = 0.0 _d 0
                2153           r_FWbySublim(i,j)     = 0.0 _d 0
286983d3d2 Patr*2154          ENDDO
                2155         ENDDO
f913c5a485 Mart*2156         DO IT=1,SEAICE_multDim
8377b8ee87 Mart*2157          DO j=1,sNy
                2158           DO i=1,sNx
53b2f6dc29 Torg*2159 C if fluxes in W/m^2 then use:
114c791332 Jean*2160 c           a_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)
286983d3d2 Patr*2161 c     &      + a_QbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
114c791332 Jean*2162 c           r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)
286983d3d2 Patr*2163 c     &      + r_QbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
114c791332 Jean*2164 c           a_QSWbyATM_cover(I,J)=a_QSWbyATM_cover(I,J)
286983d3d2 Patr*2165 c     &      + a_QSWbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
114c791332 Jean*2166 c           r_FWbySublim(I,J)=r_FWbySublim(I,J)
286983d3d2 Patr*2167 c     &      + r_FWbySublimMult(I,J,IT) * areaFracFactor(I,J,IT)
53b2f6dc29 Torg*2168 C if fluxes in effective ice meters, i.e. ice volume per area, then use:
8377b8ee87 Mart*2169            a_QbyATM_cover(i,j)=a_QbyATM_cover(i,j)
                2170      &      + a_QbyATMmult_cover(i,j,IT)
                2171            r_QbyATM_cover(i,j)=r_QbyATM_cover(i,j)
                2172      &      + r_QbyATMmult_cover(i,j,IT)
                2173            a_QSWbyATM_cover(i,j)=a_QSWbyATM_cover(i,j)
                2174      &      + a_QSWbyATMmult_cover(i,j,IT)
                2175            r_FWbySublim(i,j)=r_FWbySublim(i,j)
                2176      &      + r_FWbySublimMult(i,j,IT)
286983d3d2 Patr*2177           ENDDO
                2178          ENDDO
                2179         ENDDO
cbd0ee24a8 Mart*2180 #endif /* SEAICE_ITD */
286983d3d2 Patr*2181
4fd8e94be0 Gael*2182 #ifdef ALLOW_AUTODIFF_TAMC
1cf549c217 Mart*2183 CADJ STORE d_hsnwbyneg(:,:,bi,bj) = comlev1_bibj,
edb6656069 Mart*2184 CADJ &                              key = tkey, byte = isbyte
                2185 CADJ STORE d_hsnwbyocnonsnw = comlev1_bibj,key=tkey,byte=isbyte
4fd8e94be0 Gael*2186 #endif /* ALLOW_AUTODIFF_TAMC */
                2187
8377b8ee87 Mart*2188         DO j=1,sNy
                2189          DO i=1,sNx
                2190           QNET(i,j,bi,bj) = r_QbyATM_cover(i,j) + r_QbyATM_open(i,j)
                2191      &         +   a_QSWbyATM_cover(i,j)
                2192      &         - ( d_HEFFbyOCNonICE(i,j)
                2193      &           + d_HSNWbyOCNonSNW(i,j)*SNOW2ICE
                2194      &           + d_HEFFbyNEG(i,j,bi,bj)
6510a54854 Jean*2195 #ifdef EXF_SEAICE_FRACTION
8377b8ee87 Mart*2196      &           + d_HEFFbyRLX(i,j,bi,bj)
d32fe07ad8 Patr*2197 #endif
8377b8ee87 Mart*2198      &           + d_HSNWbyNEG(i,j,bi,bj)*SNOW2ICE
7e924a7b23 Jean*2199      &           - convertPRECIP2HI *
8377b8ee87 Mart*2200      &             snowPrecip(i,j,bi,bj) * (ONE-AREApreTH(i,j))
                2201      &           ) * HEFFM(i,j,bi,bj)
7e924a7b23 Jean*2202          ENDDO
                2203         ENDDO
8377b8ee87 Mart*2204         DO j=1,sNy
                2205          DO i=1,sNx
                2206           QSW(i,j,bi,bj)  = a_QSWbyATM_cover(i,j) + a_QSWbyATM_open(i,j)
aea7db20a6 Gael*2207          ENDDO
                2208         ENDDO
                2209
2651ba3350 Jean*2210 C switch heat fluxes from 'effective' ice meters to W/m2
aea7db20a6 Gael*2211 C ======================================================
bea7d9d588 Gael*2212
8377b8ee87 Mart*2213         DO j=1,sNy
                2214          DO i=1,sNx
                2215           QNET(i,j,bi,bj) = QNET(i,j,bi,bj)*convertHI2Q
                2216           QSW(i,j,bi,bj)  = QSW(i,j,bi,bj)*convertHI2Q
bea7d9d588 Gael*2217          ENDDO
                2218         ENDDO
0c0ecd4c7b Jean*2219
381adf77df Gael*2220 #ifndef SEAICE_DISABLE_HEATCONSFIX
                2221 C treat advective heat flux by ocean to ice water exchange (at 0decC)
                2222 C ===================================================================
4fd8e94be0 Gael*2223 # ifdef ALLOW_AUTODIFF_TAMC
1cf549c217 Mart*2224 CADJ STORE d_HEFFbyNEG(:,:,bi,bj) = comlev1_bibj,
edb6656069 Mart*2225 CADJ &                              key = tkey, byte = isbyte
                2226 CADJ STORE d_HEFFbyOCNonICE = comlev1_bibj,key=tkey,byte=isbyte
                2227 CADJ STORE d_HEFFbyATMonOCN = comlev1_bibj,key=tkey,byte=isbyte
1cf549c217 Mart*2228 CADJ STORE d_HSNWbyNEG(:,:,bi,bj) = comlev1_bibj,
edb6656069 Mart*2229 CADJ &                              key = tkey, byte = isbyte
                2230 CADJ STORE d_HSNWbyOCNonSNW = comlev1_bibj,key=tkey,byte=isbyte
                2231 CADJ STORE d_HSNWbyATMonSNW = comlev1_bibj,key=tkey,byte=isbyte
4fd8e94be0 Gael*2232 CADJ STORE theta(:,:,kSurface,bi,bj) = comlev1_bibj,
edb6656069 Mart*2233 CADJ &                       key = tkey, byte = isbyte
4fd8e94be0 Gael*2234 # endif /* ALLOW_AUTODIFF_TAMC */
1080be7801 Jean*2235 Cgf Unlike for evap and precip, the temperature of gained/lost
1d74e34c65 Jean*2236 C ocean liquid water due to melt/freeze of solid water cannot be chosen
114c791332 Jean*2237 C arbitrarily to be e.g. the ocean SST. Indeed the present seaice model
                2238 C implies a constant ice temperature of 0degC. If melt/freeze water is exchanged
                2239 C at a different temperature, it leads to a loss of conservation in the
                2240 C ocean+ice system. While this is mostly a serious issue in the
ebf1adc97e Gael*2241 C real fresh water + non linear free surface framework, a mismatch
                2242 C between ice and ocean boundary condition can result in all cases.
114c791332 Jean*2243 C Below we therefore anticipate on external_forcing_surf.F
ebf1adc97e Gael*2244 C to diagnoze and/or apply the correction to QNET.
8377b8ee87 Mart*2245         DO j=1,sNy
                2246          DO i=1,sNx
cafd3818b7 An T*2247 catn: initialize tmpscal1
                2248            tmpscal1 = ZERO
ebf1adc97e Gael*2249 C ocean water going to ice/snow, in precip units
8377b8ee87 Mart*2250            tmpscal3=rhoConstFresh*HEFFM(i,j,bi,bj)*(
                2251      &       ( d_HSNWbyATMonSNW(i,j)*SNOW2ICE
                2252      &       + d_HSNWbyOCNonSNW(i,j)*SNOW2ICE
                2253      &       + d_HEFFbyOCNonICE(i,j) + d_HEFFbyATMonOCN(i,j)
                2254      &       + d_HEFFbyNEG(i,j,bi,bj)+ d_HSNWbyNEG(i,j,bi,bj)*SNOW2ICE )
c0fe1e7a2d Gael*2255      &       * convertHI2PRECIP
8377b8ee87 Mart*2256      &       - snowPrecip(i,j,bi,bj) * (ONE-AREApreTH(i,j)) )
ebf1adc97e Gael*2257 C factor in the heat content as done in external_forcing_surf.F
634144d037 Jean*2258            IF ( (temp_EvPrRn.NE.UNSET_RL).AND.
                2259      &         useRealFreshWaterFlux.AND.(nonlinFreeSurf.NE.0) ) THEN
ebf1adc97e Gael*2260              tmpscal1 = - tmpscal3*
c0fe1e7a2d Gael*2261      &         HeatCapacity_Cp * temp_EvPrRn
634144d037 Jean*2262            ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL).AND.
                2263      &         useRealFreshWaterFlux.AND.(nonlinFreeSurf.NE.0) ) THEN
ebf1adc97e Gael*2264              tmpscal1 = - tmpscal3*
8377b8ee87 Mart*2265      &         HeatCapacity_Cp * theta(i,j,kSurface,bi,bj)
634144d037 Jean*2266            ELSEIF ( (temp_EvPrRn.NE.UNSET_RL) ) THEN
                2267              tmpscal1 = - tmpscal3*HeatCapacity_Cp*
8377b8ee87 Mart*2268      &       ( temp_EvPrRn - theta(i,j,kSurface,bi,bj) )
634144d037 Jean*2269            ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL) ) THEN
                2270              tmpscal1 = ZERO
                2271            ENDIF
381adf77df Gael*2272 #ifdef ALLOW_DIAGNOSTICS
ebf1adc97e Gael*2273 C in all cases, diagnoze the boundary condition mismatch to SIaaflux
8377b8ee87 Mart*2274            DIAGarrayA(i,j)=tmpscal1
381adf77df Gael*2275 #endif
ebf1adc97e Gael*2276 C remove the mismatch when real fresh water is exchanged (at 0degC here)
634144d037 Jean*2277            IF ( useRealFreshWaterFlux.AND.(nonlinFreeSurf.GT.0)
                2278      &          .AND.SEAICEheatConsFix )
8377b8ee87 Mart*2279      &       QNET(i,j,bi,bj)=QNET(i,j,bi,bj)+tmpscal1
c0fe1e7a2d Gael*2280          ENDDO
                2281         ENDDO
                2282 #ifdef ALLOW_DIAGNOSTICS
1080be7801 Jean*2283         IF ( useDiagnostics ) THEN
                2284           CALL DIAGNOSTICS_FILL(DIAGarrayA,
                2285      &        'SIaaflux',0,1,3,bi,bj,myThid)
                2286         ENDIF
c0fe1e7a2d Gael*2287 #endif
a73db480d4 Jean*2288 #endif /* ndef SEAICE_DISABLE_HEATCONSFIX */
c0fe1e7a2d Gael*2289
6509326d8c Gael*2290 C compute the net heat flux, incl. adv. by water, entering ocean+ice
                2291 C ===================================================================
8377b8ee87 Mart*2292         DO j=1,sNy
                2293          DO i=1,sNx
1080be7801 Jean*2294 Cgf 1) SIatmQnt (analogous to qnet; excl. adv. by water exch.)
6509326d8c Gael*2295 CML If I consider the atmosphere above the ice, the surface flux
                2296 CML which is relevant for the air temperature dT/dt Eq
                2297 CML accounts for sensible and radiation (with different treatment
                2298 CML according to wave-length) fluxes but not for "latent heat flux",
                2299 CML since it does not contribute to heating the air.
                2300 CML So this diagnostic is only good for heat budget calculations within
                2301 CML the ice-ocean system.
8377b8ee87 Mart*2302            SIatmQnt(i,j,bi,bj) =
                2303      &            HEFFM(i,j,bi,bj)*convertHI2Q*(
                2304      &            a_QSWbyATM_cover(i,j) +
                2305      &            a_QbyATM_cover(i,j) + a_QbyATM_open(i,j) )
1080be7801 Jean*2306 Cgf 2) SItflux (analogous to tflux; includes advection by water
ebf1adc97e Gael*2307 C             exchanged between atmosphere and ocean+ice)
                2308 C solid water going to atm, in precip units
8377b8ee87 Mart*2309            tmpscal1 = rhoConstFresh*HEFFM(i,j,bi,bj)
                2310      &       * convertHI2PRECIP * ( - d_HSNWbyRAIN(i,j)*SNOW2ICE
                2311      &       + a_FWbySublim(i,j) - r_FWbySublim(i,j) )
ebf1adc97e Gael*2312 C liquid water going to atm, in precip units
8377b8ee87 Mart*2313            tmpscal2=rhoConstFresh*HEFFM(i,j,bi,bj)*
                2314      &       ( ( EVAP(i,j,bi,bj)-PRECIP(i,j,bi,bj) )
                2315      &         * ( ONE - AREApreTH(i,j) )
6509326d8c Gael*2316 #ifdef ALLOW_RUNOFF
8377b8ee87 Mart*2317      &         - RUNOFF(i,j,bi,bj)
6509326d8c Gael*2318 #endif /* ALLOW_RUNOFF */
8377b8ee87 Mart*2319      &         + ( d_HFRWbyRAIN(i,j) + r_FWbySublim(i,j) )
6509326d8c Gael*2320      &         *convertHI2PRECIP )
ebf1adc97e Gael*2321 C In real fresh water flux + nonlinFS, we factor in the advected specific
                2322 C energy (referenced to 0 for 0deC liquid water). In virtual salt flux or
                2323 C linFS, rain/evap get a special treatment (see external_forcing_surf.F).
6509326d8c Gael*2324            tmpscal1= - tmpscal1*
                2325      &       ( -SEAICE_lhFusion + HeatCapacity_Cp * ZERO )
634144d037 Jean*2326            IF ( (temp_EvPrRn.NE.UNSET_RL).AND.
                2327      &          useRealFreshWaterFlux.AND.(nonlinFreeSurf.NE.0) ) THEN
                2328              tmpscal2= - tmpscal2*
                2329      &        ( ZERO + HeatCapacity_Cp * temp_EvPrRn )
                2330            ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL).AND.
                2331      &          useRealFreshWaterFlux.AND.(nonlinFreeSurf.NE.0) ) THEN
                2332              tmpscal2= - tmpscal2*
8377b8ee87 Mart*2333      &        ( ZERO + HeatCapacity_Cp * theta(i,j,kSurface,bi,bj) )
634144d037 Jean*2334            ELSEIF ( (temp_EvPrRn.NE.UNSET_RL) ) THEN
                2335              tmpscal2= - tmpscal2*HeatCapacity_Cp*
8377b8ee87 Mart*2336      &        ( temp_EvPrRn - theta(i,j,kSurface,bi,bj) )
634144d037 Jean*2337            ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL) ) THEN
                2338              tmpscal2= ZERO
                2339            ENDIF
8377b8ee87 Mart*2340            SItflux(i,j,bi,bj)= SIatmQnt(i,j,bi,bj)-tmpscal1-tmpscal2
6509326d8c Gael*2341          ENDDO
634144d037 Jean*2342         ENDDO
6509326d8c Gael*2343
0c0ecd4c7b Jean*2344 C compute net fresh water flux leaving/entering
aea7db20a6 Gael*2345 C the ocean, accounting for fresh/salt water stocks.
                2346 C ==================================================
                2347
8377b8ee87 Mart*2348         DO j=1,sNy
                2349          DO i=1,sNx
                2350           tmpscal1= d_HSNWbyATMonSNW(i,j)*SNOW2ICE
                2351      &             +d_HFRWbyRAIN(i,j)
                2352      &             +d_HSNWbyOCNonSNW(i,j)*SNOW2ICE
                2353      &             +d_HEFFbyOCNonICE(i,j)
                2354      &             +d_HEFFbyATMonOCN(i,j)
                2355      &             +d_HEFFbyNEG(i,j,bi,bj)
6510a54854 Jean*2356 #ifdef EXF_SEAICE_FRACTION
8377b8ee87 Mart*2357      &             +d_HEFFbyRLX(i,j,bi,bj)
6ec718a0fc Dimi*2358 #endif
8377b8ee87 Mart*2359      &             +d_HSNWbyNEG(i,j,bi,bj)*SNOW2ICE
6ec718a0fc Dimi*2360 C     If r_FWbySublim>0, then it is evaporated from ocean.
8377b8ee87 Mart*2361      &             +r_FWbySublim(i,j)
                2362           EmPmR(i,j,bi,bj)  = HEFFM(i,j,bi,bj)*(
                2363      &         ( EVAP(i,j,bi,bj)-PRECIP(i,j,bi,bj) )
                2364      &         * ( ONE - AREApreTH(i,j) )
6ec718a0fc Dimi*2365 #ifdef ALLOW_RUNOFF
8377b8ee87 Mart*2366      &         - RUNOFF(i,j,bi,bj)
6ec718a0fc Dimi*2367 #endif /* ALLOW_RUNOFF */
b377c9ea62 Dimi*2368      &         + tmpscal1*convertHI2PRECIP
                2369      &         )*rhoConstFresh
cc528f098c Mart*2370 #ifdef SEAICE_ITD
                2371 C     beware of the sign: fw2ObyRidge is snow mass moved into the ocean
                2372 C     by ridging, so requires a minus sign
8377b8ee87 Mart*2373      &         - fw2ObyRidge(i,j,bi,bj)*recip_deltaTtherm
                2374      &           * HEFFM(i,j,bi,bj)
cc528f098c Mart*2375 #endif /* SEAICE_ITD */
1080be7801 Jean*2376 C and the flux leaving/entering the ocean+ice
8377b8ee87 Mart*2377            SIatmFW(i,j,bi,bj) = HEFFM(i,j,bi,bj)*(
                2378      &          EVAP(i,j,bi,bj)*( ONE - AREApreTH(i,j) )
                2379      &          - PRECIP(i,j,bi,bj)
6509326d8c Gael*2380 #ifdef ALLOW_RUNOFF
8377b8ee87 Mart*2381      &          - RUNOFF(i,j,bi,bj)
6509326d8c Gael*2382 #endif /* ALLOW_RUNOFF */
                2383      &           )*rhoConstFresh
8377b8ee87 Mart*2384      &     + a_FWbySublim(i,j) * SEAICE_rhoIce * recip_deltaTtherm
6509326d8c Gael*2385
6ec718a0fc Dimi*2386          ENDDO
114c791332 Jean*2387         ENDDO
aea7db20a6 Gael*2388
d187c22362 Gael*2389 #ifdef SEAICE_DEBUG
a24915ab1a Jean*2390        IF ( plotLevel.GE.debLevC ) THEN
4bc8f4264e Mart*2391         CALL PLOT_FIELD_XYRL( QSW,'Current QSW ', myIter, myThid )
                2392         CALL PLOT_FIELD_XYRL( QNET,'Current QNET ', myIter, myThid )
                2393         CALL PLOT_FIELD_XYRL( EmPmR,'Current EmPmR ', myIter, myThid )
a24915ab1a Jean*2394        ENDIF
d187c22362 Gael*2395 #endif /* SEAICE_DEBUG */
aea7db20a6 Gael*2396
                2397 C Sea Ice Load on the sea surface.
                2398 C =================================
                2399
3a3bf6419a Gael*2400 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*2401 CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj,key=tkey,byte=isbyte
                2402 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
3a3bf6419a Gael*2403 #endif /* ALLOW_AUTODIFF_TAMC */
                2404
aea7db20a6 Gael*2405         IF ( useRealFreshWaterFlux ) THEN
8377b8ee87 Mart*2406          DO j=1,sNy
                2407           DO i=1,sNx
0f709d7d9d Gael*2408 #ifdef SEAICE_CAP_ICELOAD
8377b8ee87 Mart*2409            tmpscal1 = HEFF(i,j,bi,bj)*SEAICE_rhoIce
                2410      &              + HSNOW(i,j,bi,bj)*SEAICE_rhoSnow
4eb4a54cba Jean*2411            tmpscal2 = MIN(tmpscal1,heffTooHeavy*rhoConst)
4fecc0d4b4 Jean*2412 #else
8377b8ee87 Mart*2413            tmpscal2 = HEFF(i,j,bi,bj)*SEAICE_rhoIce
                2414      &              + HSNOW(i,j,bi,bj)*SEAICE_rhoSnow
0f709d7d9d Gael*2415 #endif
                2416            sIceLoad(i,j,bi,bj) = tmpscal2
aea7db20a6 Gael*2417           ENDDO
                2418          ENDDO
                2419         ENDIF
                2420
6509326d8c Gael*2421 #ifdef ALLOW_BALANCE_FLUXES
                2422 C Compute tile integrals of heat/fresh water fluxes to/from atm.
                2423 C ==============================================================
634144d037 Jean*2424         FWFsiTile(bi,bj) = 0. _d 0
7e00d7e8f9 Jean*2425         IF ( selectBalanceEmPmR.EQ.1 ) THEN
634144d037 Jean*2426          DO j=1,sNy
                2427           DO i=1,sNx
                2428            FWFsiTile(bi,bj) =
                2429      &       FWFsiTile(bi,bj) + SIatmFW(i,j,bi,bj)
                2430      &       * rA(i,j,bi,bj) * maskInC(i,j,bi,bj)
                2431           ENDDO
6509326d8c Gael*2432          ENDDO
634144d037 Jean*2433         ENDIF
1080be7801 Jean*2434 C to translate global mean FWF adjustements (see below) we may need :
634144d037 Jean*2435         FWF2HFsiTile(bi,bj) = 0. _d 0
7e00d7e8f9 Jean*2436         IF ( selectBalanceEmPmR.EQ.1 .AND.
                2437      &       temp_EvPrRn.EQ.UNSET_RL ) THEN
634144d037 Jean*2438          DO j=1,sNy
                2439           DO i=1,sNx
                2440            FWF2HFsiTile(bi,bj) = FWF2HFsiTile(bi,bj) +
8377b8ee87 Mart*2441      &       HeatCapacity_Cp * theta(i,j,kSurface,bi,bj)
634144d037 Jean*2442      &       * rA(i,j,bi,bj) * maskInC(i,j,bi,bj)
                2443           ENDDO
6509326d8c Gael*2444          ENDDO
634144d037 Jean*2445         ENDIF
                2446         HFsiTile(bi,bj) = 0. _d 0
                2447         IF ( balanceQnet ) THEN
                2448          DO j=1,sNy
                2449           DO i=1,sNx
                2450            HFsiTile(bi,bj) =
                2451      &       HFsiTile(bi,bj) + SItflux(i,j,bi,bj)
                2452      &       * rA(i,j,bi,bj) * maskInC(i,j,bi,bj)
                2453           ENDDO
6509326d8c Gael*2454          ENDDO
634144d037 Jean*2455         ENDIF
                2456 #endif /* ALLOW_BALANCE_FLUXES */
6509326d8c Gael*2457
ae36251cae Gael*2458 C ===================================================================
                2459 C ======================PART 8: diagnostics==========================
                2460 C ===================================================================
                2461
                2462 #ifdef ALLOW_DIAGNOSTICS
                2463         IF ( useDiagnostics ) THEN
bb8e6379cb Mart*2464          tmpscal1=1. _d 0 * recip_deltaTtherm
82e7b1b526 Gael*2465          CALL DIAGNOSTICS_SCALE_FILL(a_QbyATM_cover,
                2466      &      tmpscal1,1,'SIaQbATC',0,1,3,bi,bj,myThid)
                2467          CALL DIAGNOSTICS_SCALE_FILL(a_QbyATM_open,
                2468      &      tmpscal1,1,'SIaQbATO',0,1,3,bi,bj,myThid)
                2469          CALL DIAGNOSTICS_SCALE_FILL(a_QbyOCN,
                2470      &      tmpscal1,1,'SIaQbOCN',0,1,3,bi,bj,myThid)
                2471          CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyOCNonICE,
                2472      &      tmpscal1,1,'SIdHbOCN',0,1,3,bi,bj,myThid)
                2473          CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyATMonOCN_cover,
                2474      &      tmpscal1,1,'SIdHbATC',0,1,3,bi,bj,myThid)
                2475          CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyATMonOCN_open,
                2476      &      tmpscal1,1,'SIdHbATO',0,1,3,bi,bj,myThid)
                2477          CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyFLOODING,
                2478      &      tmpscal1,1,'SIdHbFLO',0,1,3,bi,bj,myThid)
                2479          CALL DIAGNOSTICS_SCALE_FILL(d_HSNWbyOCNonSNW,
                2480      &      tmpscal1,1,'SIdSbOCN',0,1,3,bi,bj,myThid)
                2481          CALL DIAGNOSTICS_SCALE_FILL(d_HSNWbyATMonSNW,
                2482      &      tmpscal1,1,'SIdSbATC',0,1,3,bi,bj,myThid)
                2483          CALL DIAGNOSTICS_SCALE_FILL(d_AREAbyATM,
                2484      &      tmpscal1,1,'SIdAbATO',0,1,3,bi,bj,myThid)
                2485          CALL DIAGNOSTICS_SCALE_FILL(d_AREAbyICE,
                2486      &      tmpscal1,1,'SIdAbATC',0,1,3,bi,bj,myThid)
                2487          CALL DIAGNOSTICS_SCALE_FILL(d_AREAbyOCN,
83ad492c2d Jean*2488      &      tmpscal1,1,'SIdAbOCN',0,1,3,bi,bj,myThid)
ae36251cae Gael*2489          CALL DIAGNOSTICS_SCALE_FILL(r_QbyATM_open,
82e7b1b526 Gael*2490      &      convertHI2Q,1, 'SIqneto ',0,1,3,bi,bj,myThid)
ae36251cae Gael*2491          CALL DIAGNOSTICS_SCALE_FILL(r_QbyATM_cover,
82e7b1b526 Gael*2492      &      convertHI2Q,1, 'SIqneti ',0,1,3,bi,bj,myThid)
ae36251cae Gael*2493 C three that actually need intermediate storage
8377b8ee87 Mart*2494          DO j=1,sNy
                2495           DO i=1,sNx
                2496             DIAGarrayA(i,j) = HEFFM(i,j,bi,bj)
                2497      &        * d_HSNWbyRAIN(i,j)*SEAICE_rhoSnow*recip_deltaTtherm
                2498             DIAGarrayB(i,j) =  AREA(i,j,bi,bj)-AREApreTH(i,j)
ae36251cae Gael*2499           ENDDO
                2500          ENDDO
82e7b1b526 Gael*2501          CALL DIAGNOSTICS_FILL(DIAGarrayA,
                2502      &      'SIsnPrcp',0,1,3,bi,bj,myThid)
                2503          CALL DIAGNOSTICS_SCALE_FILL(DIAGarrayB,
                2504      &      tmpscal1,1,'SIdA    ',0,1,3,bi,bj,myThid)
8377b8ee87 Mart*2505          DO j=1,sNy
                2506           DO i=1,sNx
                2507            DIAGarrayB(i,j) = HEFFM(i,j,bi,bj) *
                2508      &       a_FWbySublim(i,j) * SEAICE_rhoIce * recip_deltaTtherm
e770811cc5 Gael*2509           ENDDO
                2510          ENDDO
                2511          CALL DIAGNOSTICS_FILL(DIAGarrayB,
                2512      &      'SIfwSubl',0,1,3,bi,bj,myThid)
                2513 C
8377b8ee87 Mart*2514          DO j=1,sNy
                2515           DO i=1,sNx
83ad492c2d Jean*2516 C the actual Freshwater flux of sublimated ice, >0 decreases ice
8377b8ee87 Mart*2517            DIAGarrayA(i,j) = HEFFM(i,j,bi,bj)
                2518      &       * (a_FWbySublim(i,j)-r_FWbySublim(i,j))
bb8e6379cb Mart*2519      &       * SEAICE_rhoIce * recip_deltaTtherm
1d74e34c65 Jean*2520 C the residual Freshwater flux of sublimated ice
8377b8ee87 Mart*2521            DIAGarrayC(i,j) = HEFFM(i,j,bi,bj)
                2522      &       * r_FWbySublim(i,j)
bb8e6379cb Mart*2523      &       * SEAICE_rhoIce * recip_deltaTtherm
e770811cc5 Gael*2524 C the latent heat flux
8377b8ee87 Mart*2525            tmpscal1= EVAP(i,j,bi,bj)*( ONE - AREApreTH(i,j) )
                2526      &             + r_FWbySublim(i,j)*convertHI2PRECIP
                2527            tmpscal2= ( a_FWbySublim(i,j)-r_FWbySublim(i,j) )
4bc8f4264e Mart*2528      &             * convertHI2PRECIP
                2529            tmpscal3= SEAICE_lhEvap+SEAICE_lhFusion
8377b8ee87 Mart*2530            DIAGarrayB(i,j) = -HEFFM(i,j,bi,bj)*rhoConstFresh
4bc8f4264e Mart*2531      &             * ( tmpscal1*SEAICE_lhEvap + tmpscal2*tmpscal3 )
                2532           ENDDO
e770811cc5 Gael*2533          ENDDO
4bc8f4264e Mart*2534          CALL DIAGNOSTICS_FILL(DIAGarrayA,'SIacSubl',0,1,3,bi,bj,myThid)
                2535          CALL DIAGNOSTICS_FILL(DIAGarrayC,'SIrsSubl',0,1,3,bi,bj,myThid)
                2536          CALL DIAGNOSTICS_FILL(DIAGarrayB,'SIhl    ',0,1,3,bi,bj,myThid)
56d13a40ed Mart*2537 # if defined ( ALLOW_SITRACER ) && defined ( SEAICE_GREASE )
f61838dfc1 Torg*2538 C actual grease ice layer thickness
                2539          CALL DIAGNOSTICS_FILL(greaseLayerThick,
                2540      &        'SIgrsLT ',0,1,3,bi,bj,myThid)
                2541 # endif /* SEAICE_GREASE */
e770811cc5 Gael*2542         ENDIF
                2543 #endif /* ALLOW_DIAGNOSTICS */
c0fe1e7a2d Gael*2544
aea7db20a6 Gael*2545 C close bi,bj loops
33e17487ce Dimi*2546        ENDDO
                2547       ENDDO
66d21a8387 Jean*2548
6509326d8c Gael*2549 C ===================================================================
                2550 C =========PART 9: HF/FWF global integrals and balancing=============
                2551 C ===================================================================
                2552
                2553 #ifdef ALLOW_BALANCE_FLUXES
                2554
1080be7801 Jean*2555 C 1)  global sums
6509326d8c Gael*2556 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*2557 CADJ STORE FWFsiTile    = comlev1, key=ikey_dynamics, kind=isbyte
                2558 CADJ STORE HFsiTile     = comlev1, key=ikey_dynamics, kind=isbyte
6509326d8c Gael*2559 CADJ STORE FWF2HFsiTile = comlev1, key=ikey_dynamics, kind=isbyte
                2560 # endif /* ALLOW_AUTODIFF_TAMC */
7e00d7e8f9 Jean*2561       FWFsiGlob = 0. _d 0
                2562       FWF2HFsiGlob = 0. _d 0
                2563       IF ( selectBalanceEmPmR.EQ.1 ) THEN
                2564         CALL GLOBAL_SUM_TILE_RL( FWFsiTile, FWFsiGlob, myThid )
                2565         IF ( temp_EvPrRn.EQ.UNSET_RL ) THEN
                2566           CALL GLOBAL_SUM_TILE_RL( FWF2HFsiTile, FWF2HFsiGlob, myThid )
                2567         ELSE
                2568           FWF2HFsiGlob = HeatCapacity_Cp * temp_EvPrRn * globalArea
                2569         ENDIF
6509326d8c Gael*2570       ENDIF
                2571       HFsiGlob=0. _d 0
                2572       IF ( balanceQnet )
                2573      &   CALL GLOBAL_SUM_TILE_RL( HFsiTile, HFsiGlob, myThid )
                2574
1080be7801 Jean*2575 C 2) global means
                2576 C mean SIatmFW
6509326d8c Gael*2577       tmpscal0=FWFsiGlob / globalArea
1080be7801 Jean*2578 C corresponding mean advection by atm to ocean+ice water exchange
                2579 C        (if mean SIatmFW was removed uniformely from ocean)
6509326d8c Gael*2580       tmpscal1=FWFsiGlob / globalArea * FWF2HFsiGlob / globalArea
1080be7801 Jean*2581 C mean SItflux (before potential adjustement due to SIatmFW)
6509326d8c Gael*2582       tmpscal2=HFsiGlob / globalArea
1080be7801 Jean*2583 C mean SItflux (after potential adjustement due to SIatmFW)
7e00d7e8f9 Jean*2584       IF ( selectBalanceEmPmR.EQ.1 ) tmpscal2=tmpscal2-tmpscal1
6509326d8c Gael*2585
1080be7801 Jean*2586 C 3) balancing adjustments
7e00d7e8f9 Jean*2587       IF ( selectBalanceEmPmR.EQ.1 ) THEN
634144d037 Jean*2588        DO bj=myByLo(myThid),myByHi(myThid)
                2589         DO bi=myBxLo(myThid),myBxHi(myThid)
                2590          DO j=1-OLy,sNy+OLy
                2591           DO i=1-OLx,sNx+OLx
6509326d8c Gael*2592             empmr(i,j,bi,bj) = empmr(i,j,bi,bj) - tmpscal0
                2593             SIatmFW(i,j,bi,bj) = SIatmFW(i,j,bi,bj) - tmpscal0
1080be7801 Jean*2594 C           adjust SItflux consistently
2f464b8a56 Gael*2595             IF ( (temp_EvPrRn.NE.UNSET_RL).AND.
                2596      &        useRealFreshWaterFlux.AND.(nonlinFreeSurf.NE.0) ) THEN
6509326d8c Gael*2597             tmpscal1=
                2598      &       ( ZERO + HeatCapacity_Cp * temp_EvPrRn )
2f464b8a56 Gael*2599             ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL).AND.
                2600      &        useRealFreshWaterFlux.AND.(nonlinFreeSurf.NE.0) ) THEN
6509326d8c Gael*2601             tmpscal1=
8377b8ee87 Mart*2602      &       ( ZERO + HeatCapacity_Cp * theta(i,j,kSurface,bi,bj) )
2f464b8a56 Gael*2603             ELSEIF ( (temp_EvPrRn.NE.UNSET_RL) ) THEN
                2604             tmpscal1=
8377b8ee87 Mart*2605      &       HeatCapacity_Cp*(temp_EvPrRn - theta(i,j,kSurface,bi,bj))
2f464b8a56 Gael*2606             ELSE
                2607             tmpscal1=ZERO
6509326d8c Gael*2608             ENDIF
                2609             SItflux(i,j,bi,bj) = SItflux(i,j,bi,bj) - tmpscal0*tmpscal1
1080be7801 Jean*2610 C           no qnet or tflux adjustement is needed
634144d037 Jean*2611           ENDDO
6509326d8c Gael*2612          ENDDO
                2613         ENDDO
                2614        ENDDO
634144d037 Jean*2615        IF ( balancePrintMean ) THEN
                2616         _BEGIN_MASTER( myThid )
0320e25227 Mart*2617         WRITE(msgBuf,'(2A,1PE21.14,A,I10)') 'rm Global mean of',
                2618      &               ' SIatmFW = ', tmpscal0, '  @ it=', myIter
634144d037 Jean*2619         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                2620      &                      SQUEEZE_RIGHT, myThid )
                2621         _END_MASTER( myThid )
                2622        ENDIF
6509326d8c Gael*2623       ENDIF
                2624       IF ( balanceQnet ) THEN
634144d037 Jean*2625        DO bj=myByLo(myThid),myByHi(myThid)
                2626         DO bi=myBxLo(myThid),myBxHi(myThid)
                2627          DO j=1-OLy,sNy+OLy
                2628           DO i=1-OLx,sNx+OLx
6509326d8c Gael*2629             SItflux(i,j,bi,bj) = SItflux(i,j,bi,bj) - tmpscal2
                2630             qnet(i,j,bi,bj) = qnet(i,j,bi,bj) - tmpscal2
                2631             SIatmQnt(i,j,bi,bj) = SIatmQnt(i,j,bi,bj) - tmpscal2
634144d037 Jean*2632           ENDDO
6509326d8c Gael*2633          ENDDO
                2634         ENDDO
                2635        ENDDO
634144d037 Jean*2636        IF ( balancePrintMean ) THEN
                2637         _BEGIN_MASTER( myThid )
0320e25227 Mart*2638         WRITE(msgBuf,'(2A,1PE21.14,A,I10)') 'rm Global mean of',
                2639      &               ' SItflux = ', tmpscal2, '  @ it=', myIter
634144d037 Jean*2640         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                2641      &                      SQUEEZE_RIGHT, myThid )
                2642         _END_MASTER( myThid )
                2643        ENDIF
6509326d8c Gael*2644       ENDIF
75bbb16cce Jean*2645 #endif /* ALLOW_BALANCE_FLUXES */
6509326d8c Gael*2646
                2647 #ifdef ALLOW_DIAGNOSTICS
1080be7801 Jean*2648       IF ( useDiagnostics ) THEN
                2649 C these diags need to be done outside of the bi,bj loop so that
                2650 C we may do potential global mean adjustement to them consistently.
                2651         CALL DIAGNOSTICS_FILL(SItflux,
                2652      &        'SItflux ',0,1,0,1,1,myThid)
                2653         CALL DIAGNOSTICS_FILL(SIatmQnt,
                2654      &        'SIatmQnt',0,1,0,1,1,myThid)
                2655 C SIatmFW follows the same convention as empmr -- SIatmFW diag does not
                2656         tmpscal1= - 1. _d 0
                2657         CALL DIAGNOSTICS_SCALE_FILL(SIatmFW,
                2658      &        tmpscal1,1,'SIatmFW ',0,1,0,1,1,myThid)
                2659       ENDIF
6509326d8c Gael*2660 #endif /* ALLOW_DIAGNOSTICS */
                2661
634144d037 Jean*2662 #else  /* ALLOW_EXF and ALLOW_ATM_TEMP */
                2663       STOP 'SEAICE_GROWTH not compiled without EXF and ALLOW_ATM_TEMP'
                2664 #endif /* ALLOW_EXF and ALLOW_ATM_TEMP */
5337b47ce7 Jean*2665 #endif /* ndef SEAICE_USE_GROWTH_ADX */
2a4b53eed1 Jean*2666
33e17487ce Dimi*2667       RETURN
                2668       END