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