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