pkg/thsice/thsice_calc_thickn.F

fc7306ba7d Jean*0001 #include "THSICE_OPTIONS.h"
6b47d550f4 Mart*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
a85293d087 Mart*0004 # define ALLOW_AUTODIFF_TAMC_MORE
6b47d550f4 Mart*0005 #endif
fc7306ba7d Jean*0006
87ea84cac6 Jean*0007 CBOP
                0008 C     !ROUTINE: THSICE_CALC_THICKN
                0009 C     !INTERFACE:
                0010       SUBROUTINE THSICE_CALC_THICKN(
6dc8890c80 Patr*0011      I                  bi, bj,
7269783f6f Jean*0012      I                  iMin,iMax, jMin,jMax, dBugFlag,
                0013      I                  iceMask, tFrz, tOce, v2oc,
                0014      I                  snowP, prcAtm, sHeat, flxCnB,
c1c3d0f9d7 Patr*0015      U                  icFrac, hIce, hSnow1, tSrf, qIc1, qIc2,
7269783f6f Jean*0016      U                  frwAtm, fzMlOc, flx2oc,
281cce82f4 Jean*0017      O                  frw2oc, fsalt, frzSeaWat,
7269783f6f Jean*0018      I                  myTime, myIter, myThid )
87ea84cac6 Jean*0019 C     !DESCRIPTION: \bv
                0020 C     *==========================================================*
                0021 C     | S/R  THSICE_CALC_THICKN
                0022 C     | o Calculate ice & snow thickness changes
                0023 C     *==========================================================*
                0024 C     \ev
7269783f6f Jean*0025 C ADAPTED FROM:
                0026 C LANL CICE.v2.0.2
                0027 C-----------------------------------------------------------------------
                0028 C.. thermodynamics (vertical physics) based on M. Winton 3-layer model
88f72205aa Jean*0029 C.. See Bitz, C. M. and W. H. Lipscomb, 1999:  An energy-conserving
                0030 C..       thermodynamic sea ice model for climate study.
                0031 C..       J. Geophys. Res., 104, 15669 - 15677.
7269783f6f Jean*0032 C..     Winton, M., 1999:  "A reformulated three-layer sea ice model."
                0033 C..       Submitted to J. Atmos. Ocean. Technol.
                0034 C.. authors Elizabeth C. Hunke and William Lipscomb
                0035 C..         Fluid Dynamics Group, Los Alamos National Laboratory
                0036 C-----------------------------------------------------------------------
                0037 Cc****subroutine thermo_winton(n,fice,fsnow,dqice,dTsfc)
                0038 C.. Compute temperature change using Winton model with 2 ice layers, of
                0039 C.. which only the top layer has a variable heat capacity.
87ea84cac6 Jean*0040
c53e1eb37c Jean*0041 C---------------------------------
                0042 C  parameters that control the partitioning between lateral (ice area) and
                0043 C    vertical (ice thickness) ice volume changes.
6053aec1b8 Jean*0044 C a) surface melting and bottom melting (runtime parameter: fracEnMelt):
c53e1eb37c Jean*0045 C  frace is the fraction of available heat that is used for
                0046 C  lateral melting (and 1-frace reduces the thickness ) when
6053aec1b8 Jean*0047 C o       hi < hThinIce        & frac > lowIcFrac2 : frace=1 (lateral melting only)
                0048 C o hThinIce < hi < hThickIce  & frac > lowIcFrac1 : frace=fracEnMelt
                0049 C o            hi > hThickIce or frac < lowIcFrac1 : frace=0 (thinning only)
c53e1eb37c Jean*0050 C b) ocean freezing (and ice forming):
6053aec1b8 Jean*0051 C - conductive heat flux (below sea-ice) always increases thickness.
                0052 C - under sea-ice, freezing potential (x iceFraction) is used to increase ice
                0053 C                  thickness or ice fraction (lateral growth), according to:
                0054 C o       hi < hThinIce       : use freezing potential to grow ice vertically;
                0055 C o hThinIce < hi < hThickIce : use partition factor fracEnFreez for lateral growth
                0056 c                               and (1-fracEnFreez) to increase thickness.
                0057 C o            hi > hThickIce : use all freezing potential to grow ice laterally
                0058 C                                (up to areaMax)
                0059 C - over open ocean, use freezing potential [x(1-iceFraction)] to grow ice laterally
                0060 C - lateral growth forms ice of the same or =hNewIceMax thickness, the less of the 2.
                0061 C - starts to form sea-ice over fraction iceMaskMin, as minimum ice-volume is reached
c53e1eb37c Jean*0062 C---------------------------------
87ea84cac6 Jean*0063 C     !USES:
fc7306ba7d Jean*0064       IMPLICIT NONE
                0065
87ea84cac6 Jean*0066 C     == Global variables ===
d09af74739 Mart*0067 #include "SIZE.h"
dbce8fc2d4 Jean*0068 #include "EEPARAMS.h"
fc7306ba7d Jean*0069 #include "THSICE_SIZE.h"
                0070 #include "THSICE_PARAMS.h"
d6f06800ae Patr*0071 #ifdef ALLOW_AUTODIFF_TAMC
                0072 # include "tamc.h"
                0073 #endif
fc7306ba7d Jean*0074
87ea84cac6 Jean*0075 C     !INPUT/OUTPUT PARAMETERS:
                0076 C     == Routine Arguments ==
7269783f6f Jean*0077 C     bi,bj       :: tile indices
                0078 C     iMin,iMax   :: computation domain: 1rst index range
                0079 C     jMin,jMax   :: computation domain: 2nd  index range
                0080 C     dBugFlag    :: allow to print debugging stuff (e.g. on 1 grid point).
                0081 C---  Input:
                0082 C         iceMask :: sea-ice fractional mask [0-1]
d09af74739 Mart*0083 C  tFrz           :: sea-water freezing temperature [oC] (function of S)
                0084 C  tOce           :: surface level oceanic temperature [oC]
                0085 C  v2oc           :: square of ocean surface-level velocity [m2/s2]
                0086 C  snowP          :: snow precipitation                [kg/m2/s]
                0087 C  prcAtm         :: total precip from the atmosphere [kg/m2/s]
                0088 C  sHeat          :: surf heating flux left to melt snow or ice (= Atmos-conduction)
                0089 C  flxCnB         :: heat flux conducted through the ice to bottom surface
7269783f6f Jean*0090 C---  Modified (input&output):
d09af74739 Mart*0091 C  icFrac         :: fraction of grid area covered in ice
                0092 C  hIce           :: ice height [m]
c1c3d0f9d7 Patr*0093 C  hSnow1          :: snow height [m]
d09af74739 Mart*0094 C  tSrf           :: surface (ice or snow) temperature
7269783f6f Jean*0095 C  qIc1   (qicen) :: ice enthalpy (J/kg), 1rst level
                0096 C  qIc2   (qicen) :: ice enthalpy (J/kg), 2nd level
                0097 C  frwAtm (evpAtm):: evaporation to the atmosphere [kg/m2/s] (>0 if evaporate)
d09af74739 Mart*0098 C  fzMlOc         :: ocean mixed-layer freezing/melting potential [W/m2]
                0099 C  flx2oc         :: net heat flux to ocean    [W/m2]          (> 0 downward)
7269783f6f Jean*0100 C---  Output
d09af74739 Mart*0101 C  frw2oc         :: Total fresh water flux to ocean [kg/m2/s] (> 0 downward)
                0102 C  fsalt          :: salt flux to ocean        [g/m2/s]        (> 0 downward)
281cce82f4 Jean*0103 C  frzSeaWat      :: seawater freezing rate (expressed as mass flux) [kg/m^2/s]
7269783f6f Jean*0104 C---  Input:
                0105 C     myTime      :: current Time of simulation [s]
                0106 C     myIter      :: current Iteration number in simulation
                0107 C     myThid      :: my Thread Id number
                0108       INTEGER bi,bj
                0109       INTEGER iMin, iMax
                0110       INTEGER jMin, jMax
                0111       LOGICAL dBugFlag
6dc8890c80 Patr*0112       _RL iceMask(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0113       _RL tFrz   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0114       _RL tOce   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0115       _RL v2oc   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0116       _RL snowP  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0117       _RL prcAtm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0118       _RL sHeat  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0119       _RL flxCnB (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0120       _RL icFrac (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0121       _RL hIce   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0707ba3b3d Jean*0122       _RL hSnow1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
6dc8890c80 Patr*0123       _RL tSrf   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0124       _RL qIc1   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0125       _RL qIc2   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0126       _RL frwAtm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0127       _RL fzMlOc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0128       _RL flx2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0129       _RL frw2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0130       _RL fsalt  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
281cce82f4 Jean*0131       _RL frzSeaWat(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
7269783f6f Jean*0132       _RL  myTime
                0133       INTEGER myIter
                0134       INTEGER myThid
                0135 CEOP
                0136
                0137 #ifdef ALLOW_THSICE
                0138
                0139 C     !LOCAL VARIABLES:
                0140 C---  local copy of input/output argument list variables (see description above)
281cce82f4 Jean*0141       _RL qicen(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nlyr)
98b4e0ca2d Jean*0142
87ea84cac6 Jean*0143 C     == Local Variables ==
98b4e0ca2d Jean*0144 C     i,j,k      :: loop indices
                0145 C     rec_nlyr   :: reciprocal of number of ice layers (real value)
d09af74739 Mart*0146 C     evapLoc    :: evaporation over snow/ice [kg/m2/s] (>0 if evaporate)
98b4e0ca2d Jean*0147 C     Fbot       :: oceanic heat flux used to melt/form ice [W/m2]
                0148 C     etop       :: energy for top melting    (J m-2)
                0149 C     ebot       :: energy for bottom melting (J m-2)
                0150 C     etope      :: energy (from top)    for lateral melting (J m-2)
                0151 C     ebote      :: energy (from bottom) for lateral melting (J m-2)
                0152 C     extend     :: total energy for lateral melting (J m-2)
                0153 C     hnew(nlyr) :: new ice layer thickness (m)
                0154 C     hlyr       :: individual ice layer thickness (m)
                0155 C     dhi        :: change in ice thickness
                0156 C     dhs        :: change in snow thickness
                0157 C     rq         :: rho * q for a layer
                0158 C     rqh        :: rho * q * h for a layer
                0159 C     qbot       :: enthalpy for new ice at bottom surf (J/kg)
                0160 C     dt         :: timestep
                0161 C     esurp      :: surplus energy from melting (J m-2)
                0162 C     mwater0    :: fresh water mass gained/lost (kg/m^2)
                0163 C     msalt0     :: salt gained/lost  (kg/m^2)
                0164 C     freshe     :: fresh water gain from extension melting
                0165 C     salte      :: salt gained from extension melting
                0166 C     lowIcFrac1 :: ice-fraction lower limit above which partial (lowIcFrac1)
                0167 C     lowIcFrac2 :: or full (lowIcFrac2) lateral melting is allowed.
d09af74739 Mart*0168 C     from THSICE_RESHAPE_LAYERS
                0169 C     f1         :: Fraction of upper layer ice in new layer
                0170 C     qh1, qh2   :: qice*h for layers 1 and 2
                0171 C     qhtot      :: qh1 + qh2
                0172 C     q2tmp      :: Temporary value of qice for layer 2
98b4e0ca2d Jean*0173       INTEGER  i,j,k
d09af74739 Mart*0174       _RL rec_nlyr
281cce82f4 Jean*0175       _RL evapLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0176       _RL Fbot   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0177       _RL etop   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0178       _RL ebot   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0179       _RL etope  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0180       _RL ebote  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0181       _RL esurp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
d09af74739 Mart*0182       _RL extend
281cce82f4 Jean*0183       _RL hnew   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nlyr)
a85293d087 Mart*0184 C     help TAF with these variables
                0185       _RL hnewTmp, icFracTmp, hIceTmp, hSnwTmp
d09af74739 Mart*0186       _RL hlyr
                0187       _RL dhi
                0188       _RL dhs
                0189       _RL rq
                0190       _RL rqh
                0191       _RL qbot
                0192       _RL dt
281cce82f4 Jean*0193       _RL mwater0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0194       _RL msalt0  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
d09af74739 Mart*0195       _RL freshe
                0196       _RL salte
98b4e0ca2d Jean*0197       _RL lowIcFrac1, lowIcFrac2
d09af74739 Mart*0198       _RL  f1
                0199       _RL  qh1, qh2
                0200       _RL  qhtot
                0201       _RL  q2tmp
                0202 #ifdef CHECK_ENERGY_CONSERV
                0203       _RL  qaux(nlyr)
                0204 #endif /* CHECK_ENERGY_CONSERV */
fc7306ba7d Jean*0205
87ea84cac6 Jean*0206       _RL  ustar, cpchr
98b4e0ca2d Jean*0207       _RL  chi
fc7306ba7d Jean*0208       _RL  frace, rs, hq
4a5f035778 Jean*0209 #ifdef THSICE_FRACEN_POWERLAW
                0210       INTEGER powerLaw
                0211       _RL rec_pLaw
                0212       _RL c1Mlt, c2Mlt, aMlt, hMlt
                0213       _RL c1Frz, c2Frz, aFrz, hFrz
281cce82f4 Jean*0214       _RL enFrcMlt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
a85293d087 Mart*0215       _RL xxMlt
281cce82f4 Jean*0216       _RL enFrcFrz(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
a85293d087 Mart*0217       _RL xxFrz
                0218 # if ( !defined TARGET_NEC_SX && !defined ALLOW_AUTODIFF )
                0219       _RL tmpMlt, tmpFrz
                0220 # endif
4a5f035778 Jean*0221 #endif
7269783f6f Jean*0222
7c50f07931 Mart*0223 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0224 C     tkey :: tape key (depends on tiles)
                0225 C     kkey :: tape key (depends on levels and tiles)
                0226       INTEGER tkey, kkey
7c50f07931 Mart*0227 #endif
f6de6620bc Mart*0228 #ifdef THSICE_REGULARIZE_CALC_THICKN
a85293d087 Mart*0229       _RL kScal
f6de6620bc Mart*0230       PARAMETER ( kScal = 1.D0 )
a85293d087 Mart*0231 #endif
                0232 #ifdef ALLOW_DBUG_THSICE
7269783f6f Jean*0233 C-    define grid-point location where to print debugging values
                0234 #include "THSICE_DEBUG.h"
fc7306ba7d Jean*0235
87ea84cac6 Jean*0236  1020 FORMAT(A,1P4E11.3)
a85293d087 Mart*0237 #endif
87ea84cac6 Jean*0238
7269783f6f Jean*0239 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0240
1818702d6f Patr*0241 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0242       tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
1818702d6f Patr*0243 #endif /* ALLOW_AUTODIFF_TAMC */
d6f06800ae Patr*0244
7269783f6f Jean*0245       rec_nlyr = nlyr
                0246       rec_nlyr = 1. _d 0 / rec_nlyr
fc7306ba7d Jean*0247       dt  = thSIce_deltaT
a85293d087 Mart*0248       cpchr = cpWater*rhosw*bMeltCoef
fc7306ba7d Jean*0249
6053aec1b8 Jean*0250 C     for now, use hard coded threshold (iceMaskMin +1.% and +10.%)
                0251       lowIcFrac1 = iceMaskMin*1.01 _d 0
                0252       lowIcFrac2 = iceMaskMin*1.10 _d 0
4a5f035778 Jean*0253 #ifdef THSICE_FRACEN_POWERLAW
                0254       IF ( powerLawExp2 .GE. 1 ) THEN
                0255         powerLaw = 1 + 2**powerLawExp2
                0256         rec_pLaw = powerLaw
                0257         rec_pLaw = 1. _d 0 / rec_pLaw
                0258 C-    Coef for melting:
                0259 C     lateral-melting energy fraction = fracEnMelt - [ aMlt*(hi-hMlt) ]^powerLaw
                0260         c1Mlt = fracEnMelt**rec_pLaw
                0261         c2Mlt = (1. _d 0 - fracEnMelt)**rec_pLaw
                0262         aMlt = (c1Mlt+c2Mlt)/(hThickIce-hThinIce)
                0263         hMlt = hThinIce+c2Mlt/aMlt
                0264 C-    Coef for freezing:
                0265 C     thickening energy fraction     = fracEnFreez - [ aFrz*(hi-hFrz) ]^powerLaw
                0266         c1Frz = fracEnFreez**rec_pLaw
                0267         c2Frz = (1. _d 0 -fracEnFreez)**rec_pLaw
                0268         aFrz = (c1Frz+c2Frz)/(hThickIce-hThinIce)
                0269         hFrz = hThinIce+c2Frz/aFrz
                0270       ELSE
                0271 C-    Linear relation
                0272         powerLaw = 1
                0273         aMlt = -1. _d 0 /(hThickIce-hThinIce)
                0274         hMlt = hThickIce
                0275         aFrz = -1. _d 0 /(hThickIce-hThinIce)
                0276         hFrz = hThickIce
                0277       ENDIF
                0278 #endif /* THSICE_FRACEN_POWERLAW */
                0279
d09af74739 Mart*0280 C     initialise local arrays
281cce82f4 Jean*0281       DO j=1-OLy,sNy+OLy
                0282        DO i=1-OLx,sNx+OLx
d09af74739 Mart*0283         evapLoc(i,j) = 0. _d 0
                0284         Fbot   (i,j) = 0. _d 0
                0285         etop   (i,j) = 0. _d 0
                0286         ebot   (i,j) = 0. _d 0
                0287         etope  (i,j) = 0. _d 0
                0288         ebote  (i,j) = 0. _d 0
                0289         esurp  (i,j) = 0. _d 0
                0290         mwater0(i,j) = 0. _d 0
                0291         msalt0 (i,j) = 0. _d 0
                0292 #ifdef THSICE_FRACEN_POWERLAW
                0293         enFrcMlt(i,j)= 0. _d 0
                0294         enFrcFrz(i,j)= 0. _d 0
                0295 #endif
                0296        ENDDO
                0297       ENDDO
                0298       DO k = 1,nlyr
281cce82f4 Jean*0299        DO j=1-OLy,sNy+OLy
                0300         DO i=1-OLx,sNx+OLx
0707ba3b3d Jean*0301          qicen(i,j,k) = 0. _d 0
                0302          hnew (i,j,k) = 0. _d 0
d09af74739 Mart*0303         ENDDO
                0304        ENDDO
                0305       ENDDO
                0306
a85293d087 Mart*0307 #ifdef ALLOW_AUTODIFF_TAMC
                0308 C     The commented store directives ("cCADJ") are not picked up by TAF,
edb6656069 Mart*0309 C     that does not mean that under different circumstances they might be,
a85293d087 Mart*0310 C     so I leave them here as a reminder.
edb6656069 Mart*0311 CADJ STORE hice   = comlev1_bibj, key=tkey, kind=isbyte
                0312 CADJ STORE icfrac = comlev1_bibj, key=tkey, kind=isbyte
                0313 cCADJ STORE frwatm = comlev1_bibj, key=tkey, kind=isbyte
                0314 cCADJ STORE fzmloc = comlev1_bibj, key=tkey, kind=isbyte
                0315 cCADJ STORE hSnow1 = comlev1_bibj, key=tkey, kind=isbyte
                0316 cCADJ STORE qic1   = comlev1_bibj, key=tkey, kind=isbyte
                0317 cCADJ STORE qic2   = comlev1_bibj, key=tkey, kind=isbyte
a85293d087 Mart*0318 #endif
7269783f6f Jean*0319       DO j = jMin, jMax
                0320        DO i = iMin, iMax
d6f06800ae Patr*0321
7269783f6f Jean*0322         IF (iceMask(i,j).GT.0. _d 0) THEN
d09af74739 Mart*0323          qicen(i,j,1)= qIc1(i,j)
                0324          qicen(i,j,2)= qIc2(i,j)
7269783f6f Jean*0325 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
87ea84cac6 Jean*0326 C     initialize energies
d09af74739 Mart*0327          esurp(i,j) = 0. _d 0
fc7306ba7d Jean*0328
d09af74739 Mart*0329 c     make a local copy of evaporation
                0330          evapLoc(i,j) = frwAtm(i,j)
fc7306ba7d Jean*0331
d09af74739 Mart*0332 C------------------------------------------------------------------------
                0333 C--   Compute growth and/or melting at the top and bottom surfaces
                0334 C------------------------------------------------------------------------
fc7306ba7d Jean*0335
4a5f035778 Jean*0336 #ifdef THSICE_FRACEN_POWERLAW
d09af74739 Mart*0337          xxMlt = aMlt*(hIce(i,j)-hMlt)
                0338          xxFrz = aFrz*(hIce(i,j)-hFrz)
4a5f035778 Jean*0339 c--
                0340          IF ( powerLawExp2 .GE. 1 ) THEN
a85293d087 Mart*0341 #if ( defined TARGET_NEC_SX || defined ALLOW_AUTODIFF )
d09af74739 Mart*0342 C     avoid the short inner loop that cannot be vectorized
                0343           xxMlt = xxMlt**powerLaw
                0344           xxFrz = xxFrz**powerLaw
                0345 #else
                0346           tmpMlt = xxMlt
                0347           tmpFrz = xxFrz
                0348           DO k=1,powerLawExp2
                0349            tmpMlt = tmpMlt*tmpMlt
                0350            tmpFrz = tmpFrz*tmpFrz
                0351           ENDDO
                0352           xxMlt = xxMlt*tmpMlt
                0353           xxFrz = xxFrz*tmpFrz
                0354 #endif /* TARGET_NEC_SX */
                0355           xxMlt = fracEnMelt -xxMlt
                0356           xxFrz = fracEnFreez-xxFrz
4a5f035778 Jean*0357          ENDIF
d09af74739 Mart*0358          enFrcMlt(i,j) = MAX( 0. _d 0, MIN( xxMlt, 1. _d 0 ) )
                0359          enFrcFrz(i,j) = MAX( 0. _d 0, MIN( xxFrz, 1. _d 0 ) )
4a5f035778 Jean*0360 #endif /* THSICE_FRACEN_POWERLAW */
                0361
d09af74739 Mart*0362          IF (fzMlOc(i,j).GE. 0. _d 0) THEN
fc7306ba7d Jean*0363 C     !-----------------------------------------------------------------
                0364 C     ! freezing conditions
                0365 C     !-----------------------------------------------------------------
d09af74739 Mart*0366           Fbot(i,j) = fzMlOc(i,j)
                0367           IF ( icFrac(i,j).LT.iceMaskMax ) THEN
4a5f035778 Jean*0368 #ifdef THSICE_FRACEN_POWERLAW
d09af74739 Mart*0369            Fbot(i,j) = enFrcFrz(i,j)*fzMlOc(i,j)
4a5f035778 Jean*0370 #else /* THSICE_FRACEN_POWERLAW */
d09af74739 Mart*0371            IF (hIce(i,j).GT.hThickIce) THEN
                0372 C if higher than hThickIce, use all fzMlOc energy to grow extra ice
                0373             Fbot(i,j) = 0. _d 0
                0374            ELSEIF (hIce(i,j).GE.hThinIce) THEN
5884801cd0 Jean*0375 C between hThinIce & hThickIce, use partition factor fracEnFreez
d09af74739 Mart*0376             Fbot(i,j) = (1. _d 0 - fracEnFreez)*fzMlOc(i,j)
                0377            ENDIF
4a5f035778 Jean*0378 #endif /* THSICE_FRACEN_POWERLAW */
d09af74739 Mart*0379           ENDIF
                0380          ELSE
fc7306ba7d Jean*0381 C     !-----------------------------------------------------------------
                0382 C     ! melting conditions
                0383 C     !-----------------------------------------------------------------
f6de6620bc Mart*0384 #ifdef THSICE_REGULARIZE_CALC_THICKN
a85293d087 Mart*0385 C     smooth regulatization
f6de6620bc Mart*0386           ustar = SQRT(0.00536 _d 0 *v2oc(i,j) + 25. _d -6)
a85293d087 Mart*0387 #else
98b4e0ca2d Jean*0388 C     for no currents:
f6de6620bc Mart*0389           ustar = 5. _d -3
fc7306ba7d Jean*0390 C frictional velocity between ice and water
d4b7a0d2bc Patr*0391           IF (v2oc(i,j) .NE. 0.)
                0392      &     ustar = SQRT(0.00536 _d 0*v2oc(i,j))
f6de6620bc Mart*0393           ustar = MAX(5. _d -3,ustar)
a85293d087 Mart*0394 #endif
                0395 c         cpchr =cpWater*rhosw*bMeltCoef
d09af74739 Mart*0396           Fbot(i,j) = cpchr*(tFrz(i,j)-tOce(i,j))*ustar
                0397 C     fzMlOc < Fbot < 0
f6de6620bc Mart*0398 #ifdef THSICE_REGULARIZE_CALC_THICKN
a85293d087 Mart*0399 C     smooth regulatization
                0400 C     we want to have some soft maximum for better differentiability:
                0401 C     max(a,b;k) = (a*exp(k*a)+b*exp(k*b))/(exp(k*a)+exp(k*b))
                0402           Fbot(i,j) = (
f6de6620bc Mart*0403      &             Fbot(i,j) * EXP(kScal*  Fbot(i,j))
                0404      &         + fzMlOc(i,j) * EXP(kScal*fzMlOc(i,j))
                0405      &         )/( EXP(kScal*Fbot(i,j)) + EXP(kScal*fzMlOc(i,j)) )
a85293d087 Mart*0406 #else
f6de6620bc Mart*0407           Fbot(i,j) = MAX(Fbot(i,j),fzMlOc(i,j))
a85293d087 Mart*0408 #endif
f6de6620bc Mart*0409           Fbot(i,j) = MIN(Fbot(i,j),0. _d 0)
d09af74739 Mart*0410          ENDIF
fc7306ba7d Jean*0411
                0412 C  mass of fresh water and salt initially present in ice
c1c3d0f9d7 Patr*0413          mwater0(i,j) = rhos*hSnow1(i,j) + rhoi*hIce(i,j)
d09af74739 Mart*0414          msalt0 (i,j) = rhoi*hIce(i,j)*saltIce
fc7306ba7d Jean*0415
7269783f6f Jean*0416 #ifdef ALLOW_DBUG_THSICE
d09af74739 Mart*0417          IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
6fc136ac68 Jean*0418      &        'ThSI_CALC_TH: evpAtm, fzMlOc, Fbot =',
d09af74739 Mart*0419      &        frwAtm(i,j),fzMlOc(i,j),Fbot(i,j)
7269783f6f Jean*0420 #endif
6fc136ac68 Jean*0421 C     endif iceMask > 0
d09af74739 Mart*0422         ENDIF
                0423 C     end i/j-loops
                0424        ENDDO
                0425       ENDDO
a85293d087 Mart*0426 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0427 CADJ STORE icFrac = comlev1_bibj, key = tkey, kind = isbyte
a85293d087 Mart*0428 #endif
87ea84cac6 Jean*0429 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
d09af74739 Mart*0430       DO j = jMin, jMax
                0431        DO i = iMin, iMax
                0432         IF (iceMask(i,j).GT.0. _d 0) THEN
fc7306ba7d Jean*0433
d09af74739 Mart*0434 C     Compute energy available for melting/growth.
fc7306ba7d Jean*0435
4a5f035778 Jean*0436 #ifdef THSICE_FRACEN_POWERLAW
d09af74739 Mart*0437          IF ( fracEnMelt.EQ.0. _d 0 ) THEN
                0438           frace = 0. _d 0
                0439          ELSE
                0440           frace = (icFrac(i,j) - lowIcFrac1)/(lowIcFrac2-iceMaskMin)
                0441           frace = MIN( enFrcMlt(i,j), MAX( 0. _d 0, frace ) )
                0442          ENDIF
4a5f035778 Jean*0443 #else /* THSICE_FRACEN_POWERLAW */
d09af74739 Mart*0444          IF ( hIce(i,j).GT.hThickIce .OR. fracEnMelt.EQ.0. _d 0 ) THEN
5884801cd0 Jean*0445 C above certain height (or when no ice fractionation), only melt from top
d09af74739 Mart*0446           frace = 0. _d 0
                0447          ELSEIF (hIce(i,j).LT.hThinIce) THEN
fc7306ba7d Jean*0448 C below a certain height, all energy goes to changing ice extent
d09af74739 Mart*0449           frace = 1. _d 0
                0450          ELSE
                0451           frace = fracEnMelt
                0452          ENDIF
6053aec1b8 Jean*0453 C     Reduce lateral melting when ice fraction is low : the purpose is to avoid
                0454 C     disappearing of (up to hThinIce thick) sea-ice by over doing lateral melting
d09af74739 Mart*0455 C     (which would bring icFrac below iceMaskMin).
                0456          IF ( icFrac(i,j).LE.lowIcFrac1 ) THEN
                0457           frace = 0. _d 0
                0458          ELSEIF (icFrac(i,j).LE.lowIcFrac2 ) THEN
                0459           frace = MIN( frace, fracEnMelt )
                0460          ENDIF
4a5f035778 Jean*0461 #endif /* THSICE_FRACEN_POWERLAW */
fc7306ba7d Jean*0462
d09af74739 Mart*0463 c     IF (tSrf(i,j) .EQ. 0. _d 0 .AND. sHeat(i,j).GT.0. _d 0) THEN
                0464          IF ( sHeat(i,j).GT.0. _d 0 ) THEN
                0465           etop(i,j) = (1. _d 0-frace)*sHeat(i,j) * dt
                0466           etope(i,j) = frace*sHeat(i,j) * dt
                0467          ELSE
                0468           etop(i,j) =  0. _d 0
                0469           etope(i,j) = 0. _d 0
                0470 C jmc: found few cases where tSrf=0 & sHeat < 0 : add this line to conserv energy:
                0471           esurp(i,j) = sHeat(i,j) * dt
                0472          ENDIF
7269783f6f Jean*0473 C--   flux at the base of sea-ice:
87ea84cac6 Jean*0474 C     conduction H.flx= flxCnB (+ =down); oceanic turbulent H.flx= Fbot (+ =down).
                0475 C-    ==> energy available(+ => melt)= (flxCnB-Fbot)*dt
d09af74739 Mart*0476 c     IF (fzMlOc(i,j).LT.0. _d 0) THEN
                0477 c         ebot(i,j) = (1. _d 0-frace)*(flxCnB-Fbot(i,j)) * dt
                0478 c         ebote(i,j) = frace*(flxCnB-Fbot(i,j)) * dt
9a68c0a761 Jean*0479 c     ELSE
d09af74739 Mart*0480 c         ebot(i,j) = (flxCnB-Fbot(i,j)) * dt
                0481 c         ebote(i,j) = 0. _d 0
9a68c0a761 Jean*0482 c     ENDIF
87ea84cac6 Jean*0483 C- original formulation(above): Loose energy when flxCnB < Fbot < 0
d09af74739 Mart*0484          ebot(i,j) = (flxCnB(i,j)-Fbot(i,j)) * dt
                0485          IF (ebot(i,j).GT.0. _d 0) THEN
                0486           ebote(i,j) = frace*ebot(i,j)
                0487           ebot(i,j)  = ebot(i,j)-ebote(i,j)
                0488          ELSE
                0489           ebote(i,j) = 0. _d 0
                0490          ENDIF
7269783f6f Jean*0491 #ifdef ALLOW_DBUG_THSICE
d09af74739 Mart*0492          IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
                0493      &        'ThSI_CALC_TH: etop,etope,ebot,ebote=',
                0494      &        etop(i,j),etope(i,j),ebot(i,j),ebote(i,j)
7269783f6f Jean*0495 #endif
6fc136ac68 Jean*0496 C     endif iceMask > 0
d09af74739 Mart*0497         ENDIF
                0498 C     end i/j-loops
                0499        ENDDO
                0500       ENDDO
fc7306ba7d Jean*0501
d09af74739 Mart*0502 C     Initialize layer thicknesses. Divide total thickness equally between
                0503 C     layers
9a68c0a761 Jean*0504       DO k = 1, nlyr
d09af74739 Mart*0505        DO j = jMin, jMax
                0506         DO i = iMin, iMax
                0507          hnew(i,j,k) = hIce(i,j) * rec_nlyr
                0508         ENDDO
                0509        ENDDO
9a68c0a761 Jean*0510       ENDDO
fc7306ba7d Jean*0511
a85293d087 Mart*0512 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0513 CADJ STORE hSnow1 = comlev1_bibj, key = tkey, kind = isbyte
                0514 CADJ STORE etop   = comlev1_bibj, key = tkey, kind = isbyte
a85293d087 Mart*0515 #endif
d09af74739 Mart*0516       DO j = jMin, jMax
                0517        DO i = iMin, iMax
                0518         IF (iceMask(i,j) .GT. 0. _d 0 .AND.
                0519      &         etop(i,j) .GT. 0. _d 0 .AND.
f6de6620bc Mart*0520      &       hSnow1(i,j) .GT. 0. _d 0) THEN
d09af74739 Mart*0521
                0522 C     Make sure internal ice temperatures do not exceed Tmlt.
                0523 C     If they do, then eliminate the layer.  (Dont think this will happen
                0524 C     for reasonable values of i0.)
                0525 C     Top melt: snow, then ice.
                0526          rq =  rhos * qsnow
c1c3d0f9d7 Patr*0527          rqh = rq * hSnow1(i,j)
d09af74739 Mart*0528          IF (etop(i,j) .LT. rqh) THEN
c1c3d0f9d7 Patr*0529           hSnow1(i,j) = hSnow1(i,j) - etop(i,j)/rq
d09af74739 Mart*0530           etop(i,j) = 0. _d 0
                0531          ELSE
                0532           etop(i,j) = etop(i,j) - rqh
c1c3d0f9d7 Patr*0533           hSnow1(i,j) = 0. _d 0
9a68c0a761 Jean*0534          ENDIF
d09af74739 Mart*0535 C     endif iceMask > 0, etc.
                0536         ENDIF
                0537 C     end i/j-loops
                0538        ENDDO
                0539       ENDDO
a85293d087 Mart*0540 #ifdef ALLOW_AUTODIFF_TAMC_MORE
edb6656069 Mart*0541 CADJ STORE hSnow1 = comlev1_bibj, key = tkey, kind = isbyte
a85293d087 Mart*0542 #endif
d09af74739 Mart*0543 C     two layers of ice
                0544       DO k = 1, nlyr
a85293d087 Mart*0545 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0546        kkey = k + (tkey-1)*nlyr
                0547 CADJ STORE etop         = comlev1_thsice_nlyr, key=kkey, kind=isbyte
                0548 CADJ STORE hnew (:,:,k) = comlev1_thsice_nlyr, key=kkey, kind=isbyte
                0549 CADJ STORE qicen(:,:,k) = comlev1_thsice_nlyr, key=kkey, kind=isbyte
a85293d087 Mart*0550 #endif
d09af74739 Mart*0551        DO j = jMin, jMax
                0552         DO i = iMin, iMax
                0553          IF (iceMask(i,j).GT.0. _d 0) THEN
                0554           IF (etop(i,j) .GT. 0. _d 0) THEN
                0555            rq =  rhoi * qicen(i,j,k)
                0556            rqh = rq * hnew(i,j,k)
                0557            IF (etop(i,j) .LT. rqh) THEN
                0558             hnew(i,j,k) = hnew(i,j,k) - etop(i,j) / rq
                0559             etop(i,j) = 0. _d 0
                0560            ELSE
                0561             etop(i,j) = etop(i,j) - rqh
                0562             hnew(i,j,k) = 0. _d 0
                0563            ENDIF
                0564           ELSE
                0565            etop(i,j)=0. _d 0
                0566           ENDIF
fc7306ba7d Jean*0567 C If ice is gone and melting energy remains
d09af74739 Mart*0568 c     IF (etop(i,j) .GT. 0. _d 0) THEN
9a68c0a761 Jean*0569 c        WRITE (6,*)  'QQ All ice melts from top  ', i,j
d09af74739 Mart*0570 c        hIce(i,j)=0. _d 0
fc7306ba7d Jean*0571 c        go to 200
9a68c0a761 Jean*0572 c     ENDIF
fc7306ba7d Jean*0573
6fc136ac68 Jean*0574 C     endif iceMask > 0
d09af74739 Mart*0575          ENDIF
                0576 C     end i/j-loops
                0577         ENDDO
                0578        ENDDO
                0579 C     end k-loop
                0580       ENDDO
a85293d087 Mart*0581 #ifdef ALLOW_AUTODIFF_TAMC_MORE
                0582 C     if isbyte=4, these change the AD-monitor output a little (11-14
                0583 C     digits of agreement)
edb6656069 Mart*0584 CADJ STORE ebot  = comlev1_bibj, key = tkey, kind = isbyte
                0585 CADJ STORE hnew  = comlev1_bibj, key = tkey, kind = isbyte
                0586 CADJ STORE qicen = comlev1_bibj, key = tkey, kind = isbyte
a85293d087 Mart*0587 #endif
fc7306ba7d Jean*0588
281cce82f4 Jean*0589 C Bottom growth:
d09af74739 Mart*0590       DO j = jMin, jMax
                0591        DO i = iMin, iMax
                0592         IF (iceMask(i,j).GT.0. _d 0 .AND. ebot(i,j) .LT. 0. _d 0) THEN
fc7306ba7d Jean*0593 C Compute enthalpy of new ice growing at bottom surface.
d09af74739 Mart*0594          qbot =  -cpIce *tFrz(i,j) + Lfresh
                0595          dhi = -ebot(i,j) / (qbot * rhoi)
                0596          ebot(i,j) = 0. _d 0
6fc136ac68 Jean*0597          qicen(i,j,nlyr) =
d09af74739 Mart*0598      &        (hnew(i,j,nlyr)*qicen(i,j,nlyr)+dhi*qbot) /
                0599      &        (hnew(i,j,nlyr)+dhi)
                0600          hnew(i,j,nlyr) = hnew(i,j,nlyr) + dhi
281cce82f4 Jean*0601          frzSeaWat(i,j) = rhoi*dhi/dt
d09af74739 Mart*0602
                0603 C     endif iceMask > 0 and ebot < 0
                0604         ENDIF
                0605 C     end i/j-loops
                0606        ENDDO
                0607       ENDDO
1818702d6f Patr*0608
a85293d087 Mart*0609 #ifdef ALLOW_AUTODIFF_TAMC_MORE
                0610 C     these change the results (increase gradient error sligthly);
                0611 C     they do not improve performance too much, so we leave commented
edb6656069 Mart*0612 cCADJ STORE hnew (:,:,nlyr) = comlev1_bibj, key = tkey, kind = isbyte
                0613 cCADJ STORE qicen(:,:,nlyr) = comlev1_bibj, key = tkey, kind = isbyte
1818702d6f Patr*0614 #endif
                0615
281cce82f4 Jean*0616 C Bottom melt:
d09af74739 Mart*0617       DO k = nlyr, 1, -1
a85293d087 Mart*0618 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0619        kkey = k + (tkey-1)*nlyr
                0620 CADJ STORE ebot = comlev1_thsice_nlyr, key=kkey, kind = isbyte
a85293d087 Mart*0621 #endif
d09af74739 Mart*0622        DO j = jMin, jMax
                0623         DO i = iMin, iMax
                0624          IF (iceMask(i,j) .GT. 0. _d 0 .AND.
6fc136ac68 Jean*0625      &        ebot(i,j)   .GT. 0. _d 0 .AND.
d09af74739 Mart*0626      &        hnew(i,j,k) .GT. 0. _d 0) THEN
                0627           rq =  rhoi * qicen(i,j,k)
                0628           rqh = rq * hnew(i,j,k)
                0629           IF (ebot(i,j) .LT. rqh) THEN
                0630            hnew(i,j,k) = hnew(i,j,k) - ebot(i,j) / rq
                0631            ebot(i,j) = 0. _d 0
                0632           ELSE
                0633            ebot(i,j) = ebot(i,j) - rqh
                0634            hnew(i,j,k) = 0. _d 0
                0635           ENDIF
                0636 C     endif iceMask > 0 etc.
                0637          ENDIF
                0638 C     end i/j-loops
                0639         ENDDO
                0640        ENDDO
                0641 C     end k-loop
                0642       ENDDO
a85293d087 Mart*0643 #ifdef ALLOW_AUTODIFF_TAMC_MORE
                0644 C     if isbyte=4, these change the AD-monitor output a little (10-11
                0645 C     digits of agreement)
edb6656069 Mart*0646 CADJ STORE hnew = comlev1_bibj, key = tkey, kind = isbyte
                0647 CADJ STORE ebot = comlev1_bibj, key = tkey, kind = isbyte
a85293d087 Mart*0648 #endif
d09af74739 Mart*0649 C     If ice melts completely and snow is left, remove the snow with
                0650 C     energy from the mixed layer
                0651       DO j = jMin, jMax
                0652        DO i = iMin, iMax
6fc136ac68 Jean*0653         IF (iceMask(i,j) .GT. 0. _d 0 .AND.
                0654      &       ebot(i,j)   .GT. 0. _d 0 .AND.
f6de6620bc Mart*0655      &       hSnow1(i,j) .GT. 0. _d 0) THEN
d09af74739 Mart*0656          rq =  rhos * qsnow
c1c3d0f9d7 Patr*0657          rqh = rq * hSnow1(i,j)
d09af74739 Mart*0658          IF (ebot(i,j) .LT. rqh) THEN
c1c3d0f9d7 Patr*0659           hSnow1(i,j) = hSnow1(i,j) - ebot(i,j) / rq
d09af74739 Mart*0660           ebot(i,j) = 0. _d 0
                0661          ELSE
                0662           ebot(i,j) = ebot(i,j) - rqh
c1c3d0f9d7 Patr*0663           hSnow1(i,j) = 0. _d 0
9a68c0a761 Jean*0664          ENDIF
d09af74739 Mart*0665 c        IF (ebot(i,j) .GT. 0. _d 0) THEN
fc7306ba7d Jean*0666 c           IF (dBug) WRITE(6,*) 'All ice (& snow) melts from bottom'
d09af74739 Mart*0667 c           hIce(i,j)=0. _d 0
fc7306ba7d Jean*0668 c           go to 200
9a68c0a761 Jean*0669 c        ENDIF
fc7306ba7d Jean*0670
d09af74739 Mart*0671 C     endif iceMask > 0, etc.
                0672         ENDIF
                0673 C     end i/j-loops
                0674        ENDDO
9a68c0a761 Jean*0675       ENDDO
a85293d087 Mart*0676 #ifdef ALLOW_AUTODIFF_TAMC_MORE
edb6656069 Mart*0677 CADJ STORE hSnow1 = comlev1_bibj, key = tkey, kind = isbyte
a85293d087 Mart*0678 #endif
d09af74739 Mart*0679       DO j = jMin, jMax
                0680        DO i = iMin, iMax
                0681         IF (iceMask(i,j).GT.0. _d 0) THEN
                0682 C Compute new total ice thickness.
                0683          hIce(i,j) = hnew(i,j,1) + hnew(i,j,2)
7269783f6f Jean*0684 #ifdef ALLOW_DBUG_THSICE
d09af74739 Mart*0685          IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
c1c3d0f9d7 Patr*0686      &        'ThSI_CALC_TH:   etop, ebot, hIce, hSnow1 =',
                0687      &        etop(i,j), ebot(i,j), hIce(i,j), hSnow1(i,j)
7269783f6f Jean*0688 #endif
fc7306ba7d Jean*0689
d09af74739 Mart*0690 C If hIce < hIceMin, melt the ice.
6fc136ac68 Jean*0691          IF ( hIce(i,j).LT.hIceMin
c1c3d0f9d7 Patr*0692      &        .AND. (hIce(i,j)+hSnow1(i,j)).GT.0. _d 0 ) THEN
                0693           esurp(i,j) = esurp(i,j) - rhos*qsnow*hSnow1(i,j)
d09af74739 Mart*0694      &         - rhoi*qicen(i,j,1)*hnew(i,j,1)
                0695      &         - rhoi*qicen(i,j,2)*hnew(i,j,2)
                0696           hIce(i,j)   = 0. _d 0
c1c3d0f9d7 Patr*0697           hSnow1(i,j)  = 0. _d 0
d09af74739 Mart*0698           tSrf(i,j)   = 0. _d 0
                0699           icFrac(i,j) = 0. _d 0
                0700           qicen(i,j,1) = 0. _d 0
                0701           qicen(i,j,2) = 0. _d 0
7269783f6f Jean*0702 #ifdef ALLOW_DBUG_THSICE
d09af74739 Mart*0703           IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
                0704      &         'ThSI_CALC_TH: -1 : esurp=',esurp(i,j)
7269783f6f Jean*0705 #endif
d09af74739 Mart*0706          ENDIF
                0707
6fc136ac68 Jean*0708 C     endif iceMask > 0
d09af74739 Mart*0709         ENDIF
                0710 C     end i/j-loops
                0711        ENDDO
                0712       ENDDO
1818702d6f Patr*0713
d09af74739 Mart*0714       DO j = jMin, jMax
                0715        DO i = iMin, iMax
                0716         IF (iceMask(i,j).GT.0. _d 0) THEN
fc7306ba7d Jean*0717
                0718 C--   do a mass-budget of sea-ice to compute "fresh" = the fresh-water flux
                0719 C     that is returned to the ocean ; needs to be done before snow/evap
6fc136ac68 Jean*0720          frw2oc(i,j) = (mwater0(i,j)
c1c3d0f9d7 Patr*0721      &        - (rhos*hSnow1(i,j)+rhoi*hIce(i,j)))/dt
fc7306ba7d Jean*0722
d09af74739 Mart*0723          IF ( hIce(i,j) .LE. 0. _d 0 ) THEN
9a68c0a761 Jean*0724 C-    return  snow to the ocean (account for Latent heat of freezing)
d09af74739 Mart*0725           frw2oc(i,j) = frw2oc(i,j) + snowP(i,j)
                0726           flx2oc(i,j) = flx2oc(i,j) - snowP(i,j)*Lfresh
                0727          ENDIF
6fc136ac68 Jean*0728
                0729 C     endif iceMask > 0
d09af74739 Mart*0730         ENDIF
                0731 C     end i/j-loops
                0732        ENDDO
                0733       ENDDO
                0734 C-    else: hIce > 0
                0735       DO j = jMin, jMax
                0736        DO i = iMin, iMax
                0737         IF (iceMask(i,j).GT.0. _d 0) THEN
fc7306ba7d Jean*0738
d09af74739 Mart*0739          IF ( hIce(i,j) .GT. 0. _d 0 ) THEN
fc7306ba7d Jean*0740 C Let it snow
c1c3d0f9d7 Patr*0741           hSnow1(i,j) = hSnow1(i,j) + dt*snowP(i,j)/rhos
fc7306ba7d Jean*0742 C If ice evap is used to sublimate surface snow/ice or
                0743 C if no ice pass on to ocean
c1c3d0f9d7 Patr*0744           IF (hSnow1(i,j).GT.0. _d 0) THEN
                0745            IF (evapLoc(i,j)/rhos *dt.GT.hSnow1(i,j)) THEN
                0746             evapLoc(i,j)=evapLoc(i,j)-hSnow1(i,j)*rhos/dt
                0747             hSnow1(i,j)=0. _d 0
d09af74739 Mart*0748            ELSE
c1c3d0f9d7 Patr*0749             hSnow1(i,j) = hSnow1(i,j) - evapLoc(i,j)/rhos *dt
d09af74739 Mart*0750             evapLoc(i,j)=0. _d 0
                0751            ENDIF
                0752           ENDIF
                0753 C     endif hice > 0
9a68c0a761 Jean*0754          ENDIF
6fc136ac68 Jean*0755 C     endif iceMask > 0
d09af74739 Mart*0756         ENDIF
                0757 C     end i/j-loops
                0758        ENDDO
                0759       ENDDO
1818702d6f Patr*0760
d09af74739 Mart*0761 C-    else: hIce > 0
                0762       DO k = 1, nlyr
a85293d087 Mart*0763 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0764        kkey = k + (tkey-1)*nlyr
                0765 CADJ STORE evaploc      = comlev1_thsice_nlyr, key=kkey, kind=isbyte
a85293d087 Mart*0766 #endif
d09af74739 Mart*0767        DO j = jMin, jMax
                0768         DO i = iMin, iMax
                0769          IF (iceMask(i,j).GT.0. _d 0 ) THEN
                0770
a85293d087 Mart*0771            hnewTmp = hnew(i,j,k)
d09af74739 Mart*0772            IF (hIce(i,j).GT.0. _d 0.AND.evapLoc(i,j).GT.0. _d 0) THEN
                0773 C            IF (evapLoc(i,j) .GT. 0. _d 0) THEN
fc7306ba7d Jean*0774 C-- original scheme, does not care about ice temp.
                0775 C-  this can produce small error (< 1.W/m2) in the Energy budget
a85293d087 Mart*0776 c              IF (evapLoc(i,j)/rhoi *dt.GT.hnewTmp) THEN
                0777 c                evapLoc(i,j)=evapLoc(i,j)-hnewTmp*rhoi/dt
d09af74739 Mart*0778 c                hnew(i,j,k)=0. _d 0
9a68c0a761 Jean*0779 c              ELSE
a85293d087 Mart*0780 c                hnew(i,j,k) = hnewTmp - evapLoc(i,j)/rhoi *dt
d09af74739 Mart*0781 c                evapLoc(i,j)=0. _d 0
9a68c0a761 Jean*0782 c              ENDIF
fc7306ba7d Jean*0783 C-- modified scheme. taking into account Ice enthalpy
d09af74739 Mart*0784              dhi = evapLoc(i,j)/rhoi*dt
a85293d087 Mart*0785              IF (dhi.GE.hnewTmp) THEN
                0786               evapLoc(i,j)=evapLoc(i,j)-hnewTmp*rhoi/dt
6fc136ac68 Jean*0787               esurp(i,j) = esurp(i,j)
a85293d087 Mart*0788      &             - hnewTmp*rhoi*(qicen(i,j,k)-Lfresh)
d09af74739 Mart*0789               hnew(i,j,k)=0. _d 0
                0790              ELSE
a85293d087 Mart*0791               hq = hnewTmp*qicen(i,j,k)-dhi*Lfresh
                0792               hnew(i,j,k) = hnewTmp - dhi
d09af74739 Mart*0793               qicen(i,j,k)=hq/hnew(i,j,k)
                0794               evapLoc(i,j)=0. _d 0
                0795              ENDIF
fc7306ba7d Jean*0796 C-------
d09af74739 Mart*0797 c     IF (evapLoc(i,j) .GT. 0. _d 0) THEN
9a68c0a761 Jean*0798 c           WRITE (6,*)  'BB All ice sublimates', i,j
d09af74739 Mart*0799 c           hIce(i,j)=0. _d 0
fc7306ba7d Jean*0800 c           go to 200
9a68c0a761 Jean*0801 c     ENDIF
d09af74739 Mart*0802 C     endif hice > 0 and evaploc > 0
                0803           ENDIF
6fc136ac68 Jean*0804 C     endif iceMask > 0
d09af74739 Mart*0805          ENDIF
                0806 C     end i/j-loops
                0807         ENDDO
9a68c0a761 Jean*0808        ENDDO
d09af74739 Mart*0809 C     end k-loop
                0810       ENDDO
1818702d6f Patr*0811
a85293d087 Mart*0812 #ifdef ALLOW_AUTODIFF_TAMC_MORE
                0813 C     if isbyte=4, these change the AD-monitor output a little (8-10
                0814 C     digits of agreement)
edb6656069 Mart*0815 CADJ STORE evapLoc = comlev1_bibj, key = tkey, kind = isbyte
                0816 CADJ STORE hSnow1  = comlev1_bibj, key = tkey, kind = isbyte
                0817 CADJ STORE hIce    = comlev1_bibj, key = tkey, kind = isbyte
                0818 CADJ STORE hnew    = comlev1_bibj, key = tkey, kind = isbyte
                0819 CADJ STORE qicen   = comlev1_bibj, key = tkey, kind = isbyte
1818702d6f Patr*0820 #endif
                0821
6fc136ac68 Jean*0822 C     still else: hice > 0
d09af74739 Mart*0823       DO j = jMin, jMax
                0824        DO i = iMin, iMax
                0825         IF (iceMask(i,j).GT.0. _d 0) THEN
                0826          IF (hIce(i,j) .GT. 0. _d 0) THEN
                0827 C Compute new total ice thickness.
6fc136ac68 Jean*0828           hIce(i,j) = hnew(i,j,1) + hnew(i,j,2)
d09af74739 Mart*0829 C If hIce < hIceMin, melt the ice.
                0830           IF ( hIce(i,j).GT.0. _d 0 .AND. hIce(i,j).LT.hIceMin ) THEN
6fc136ac68 Jean*0831            frw2oc(i,j) = frw2oc(i,j)
c1c3d0f9d7 Patr*0832      &          + (rhos*hSnow1(i,j) + rhoi*hIce(i,j))/dt
                0833            esurp(i,j) = esurp(i,j) - rhos*qsnow*hSnow1(i,j)
d09af74739 Mart*0834      &          - rhoi*qicen(i,j,1)*hnew(i,j,1)
                0835      &          - rhoi*qicen(i,j,2)*hnew(i,j,2)
                0836            hIce(i,j)   = 0. _d 0
c1c3d0f9d7 Patr*0837            hSnow1(i,j)  = 0. _d 0
d09af74739 Mart*0838            tSrf(i,j)   = 0. _d 0
                0839            icFrac(i,j) = 0. _d 0
                0840            qicen(i,j,1) = 0. _d 0
                0841            qicen(i,j,2) = 0. _d 0
7269783f6f Jean*0842 #ifdef ALLOW_DBUG_THSICE
d09af74739 Mart*0843            IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
6fc136ac68 Jean*0844      &          'ThSI_CALC_TH: -2 : esurp,frw2oc=',
d09af74739 Mart*0845      &          esurp(i,j), frw2oc(i,j)
7269783f6f Jean*0846 #endif
d09af74739 Mart*0847           ENDIF
9a68c0a761 Jean*0848
d09af74739 Mart*0849 C--   else hIce > 0: end
                0850          ENDIF
                0851
6fc136ac68 Jean*0852 C     endif iceMask > 0
d09af74739 Mart*0853         ENDIF
                0854 C     end i/j-loops
                0855        ENDDO
                0856       ENDDO
1818702d6f Patr*0857
d09af74739 Mart*0858       DO j = jMin, jMax
                0859        DO i = iMin, iMax
a85293d087 Mart*0860 C     These helper variables avoid TAF recomputation warnings
                0861         hIceTmp = hIce(i,j)
                0862         hSnwTmp = hSnow1(i,j)
                0863         hnewTmp = hnew(i,j,1)
d09af74739 Mart*0864         IF (iceMask(i,j).GT.0. _d 0) THEN
9a68c0a761 Jean*0865
a85293d087 Mart*0866          IF ( hIceTmp .GT. 0. _d 0 ) THEN
fc7306ba7d Jean*0867
7269783f6f Jean*0868 C If there is enough snow to lower the ice/snow interface below
                0869 C freeboard, convert enough snow to ice to bring the interface back
e6b6bab319 Davi*0870 C to sea-level OR if snow height is larger than hsMax, snow is
c1c3d0f9d7 Patr*0871 C converted to ice to bring hSnow1 down to hsMax. Largest change is
e6b6bab319 Davi*0872 C applied and enthalpy of top ice layer adjusted accordingly.
fc7306ba7d Jean*0873
a85293d087 Mart*0874           IF ( hSnwTmp .GT. hIceTmp*floodFac
                0875      &         .OR. hSnwTmp .GT. hsMax ) THEN
9a68c0a761 Jean*0876 cBB               WRITE (6,*)  'Freeboard adjusts'
a85293d087 Mart*0877 c          dhi = (hSnwTmp * rhos - hIceTmp * rhoiw) / rhosw
d09af74739 Mart*0878 c          dhs = dhi * rhoi / rhos
a85293d087 Mart*0879            dhs = (hSnwTmp - hIceTmp*floodFac) * rhoi / rhosw
                0880            dhs = MAX( hSnwTmp - hsMax, dhs )
d09af74739 Mart*0881            dhi = dhs * rhos / rhoi
a85293d087 Mart*0882            rqh = rhoi*qicen(i,j,1)*hnewTmp + rhos*qsnow*dhs
                0883            hnew(i,j,1)  = hnewTmp + dhi
                0884            qicen(i,j,1) = rqh / (rhoi*hnew(i,j,1))
                0885            hIce(i,j)    = hIceTmp + dhi
                0886            hSnow1(i,j)  = hSnwTmp - dhs
d09af74739 Mart*0887           ENDIF
a85293d087 Mart*0888 #ifdef ALLOW_AUTODIFF
                0889 C     This is awkward, but makes storing a few variables (hIce, hnew,
                0890 C     qicen) unnecessary. Splitting the loop also changes the condition
                0891 C     hIce>0 because hIce has been modified above. This could be handled
                0892 C     more cleanly, but because hIce>0 is required for the previous
                0893 C     block to be executed where snow is turned into ice hIce>0 is still
                0894 C     unchanged, so splitting this loop with the if conditions should be
                0895 C     safe. If this turns out to be robust, we can use this code also
                0896 C     for non-AD simulation.
fc7306ba7d Jean*0897
a85293d087 Mart*0898 C-    if hIce > 0 : end
                0899          ENDIF
                0900 C     endif iceMask > 0
                0901         ENDIF
                0902 C     end i/j-loops
                0903        ENDDO
                0904       ENDDO
                0905 #ifdef ALLOW_AUTODIFF_TAMC_MORE
                0906 C     if isbyte=4, these change the AD-monitor output a little (8-10
                0907 C     digits of agreement), but avoid quite a few recomputations
edb6656069 Mart*0908 CADJ STORE qicen       = comlev1_bibj, key = tkey, kind = isbyte
                0909 CADJ STORE hice        = comlev1_bibj, key = tkey, kind = isbyte
                0910 CADJ STORE hnew(:,:,1) = comlev1_bibj, key = tkey, kind = isbyte
a85293d087 Mart*0911 #endif
                0912       DO j = jMin, jMax
                0913        DO i = iMin, iMax
                0914 C     This helper variable avoids TAF recomputation warnings
                0915         hIceTmp = hIce(i,j)
                0916         IF (iceMask(i,j).GT.0. _d 0) THEN
                0917          IF ( hIceTmp .GT. 0. _d 0 ) THEN
                0918 C limit ice height
                0919 C- NOTE: this part does not conserve Energy ;
                0920 C        but surplus of fresh water and salt are taken into account.
                0921           IF (hIceTmp.GT.hiMax) THEN
                0922 cBB      print*,'BBerr, hIce>hiMax',i,j,hIce(i,j)
                0923            chi=hIceTmp-hiMax
                0924 #else /* ndef ALLOW_AUTODIFF */
fc7306ba7d Jean*0925 C limit ice height
                0926 C- NOTE: this part does not conserve Energy ;
                0927 C        but surplus of fresh water and salt are taken into account.
d09af74739 Mart*0928           IF (hIce(i,j).GT.hiMax) THEN
                0929 cBB      print*,'BBerr, hIce>hiMax',i,j,hIce(i,j)
                0930            chi=hIce(i,j)-hiMax
a85293d087 Mart*0931 #endif /* ALLOW_AUTODIFF */
d09af74739 Mart*0932            hnew(i,j,1)=hnew(i,j,1)-chi/2. _d 0
                0933            hnew(i,j,2)=hnew(i,j,2)-chi/2. _d 0
                0934            frw2oc(i,j) = frw2oc(i,j) + chi*rhoi/dt
                0935           ENDIF
c1c3d0f9d7 Patr*0936 c       IF (hSnow1(i,j).GT.hsMax) THEN
                0937 cc        print*,'BBerr, hSnow1>hsMax',i,j,hSnow1(i,j)
                0938 c         chs=hSnow1(i,j)-hsMax
                0939 c         hSnow1(i,j)=hsMax
d09af74739 Mart*0940 c         frw2oc(i,j) = frw2oc(i,j) + chs*rhos/dt
e6b6bab319 Davi*0941 c       ENDIF
fc7306ba7d Jean*0942
                0943 C Compute new total ice thickness.
d09af74739 Mart*0944           hIce(i,j) = hnew(i,j,1) + hnew(i,j,2)
fc7306ba7d Jean*0945
7269783f6f Jean*0946 #ifdef ALLOW_DBUG_THSICE
d09af74739 Mart*0947           IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
6fc136ac68 Jean*0948      &         'ThSI_CALC_TH: b-Winton: hnew, qice =',
                0949      &         hnew(i,j,1), hnew(i,j,2),
                0950      &         qicen(i,j,1), qicen(i,j,2)
7269783f6f Jean*0951 #endif
                0952
d09af74739 Mart*0953           hlyr = hIce(i,j) * rec_nlyr
                0954 CML          CALL THSICE_RESHAPE_LAYERS(
                0955 CML     U         qicen(i,j,:),
                0956 CML     I         hlyr, hnew(i,j,:), myThid )
                0957 C     inlined version of S/R THSICE_RESHAPE_LAYERS
                0958 C     | Repartition into equal-thickness layers, conserving energy.
                0959 C     *==========================================================*
6fc136ac68 Jean*0960 C     | This is the 2-layer version (formerly "NEW_LAYERS_WINTON")
d09af74739 Mart*0961 C     |  from M. Winton 1999, JAOT, sea-ice model.
                0962           if (hnew(i,j,1).gt.hnew(i,j,2)) then
                0963 C--   Layer 1 gives ice to layer 2
                0964            f1 = (hnew(i,j,1)-hlyr)/hlyr
                0965            q2tmp = f1*qicen(i,j,1) + (1. _d 0-f1)*qicen(i,j,2)
                0966            if (q2tmp.gt.Lfresh) then
                0967             qicen(i,j,2) = q2tmp
                0968            else
                0969 C-    Keep q2 fixed to avoid q2<Lfresh and T2>0
                0970             qh2 = hlyr*qicen(i,j,2)
                0971             qhtot = hnew(i,j,1)*qicen(i,j,1) + hnew(i,j,2)*qicen(i,j,2)
                0972             qh1 = qhtot - qh2
                0973             qicen(i,j,1) = qh1/hlyr
                0974            endif
                0975           else
                0976 C-    Layer 2 gives ice to layer 1
                0977            f1 = hnew(i,j,1)/hlyr
                0978            qicen(i,j,1) = f1*qicen(i,j,1) + (1. _d 0-f1)*qicen(i,j,2)
                0979           endif
                0980 C     end of inlined S/R THSICE_RESHAPE_LAYERS
fc7306ba7d Jean*0981
7269783f6f Jean*0982 #ifdef ALLOW_DBUG_THSICE
d09af74739 Mart*0983           IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
c1c3d0f9d7 Patr*0984      &         'ThSI_CALC_TH: icFrac,hIce, qtot, hSnow1 =',
6fc136ac68 Jean*0985      &         icFrac(i,j),hIce(i,j), (qicen(i,j,1)+qicen(i,j,2))*0.5,
c1c3d0f9d7 Patr*0986      &         hSnow1(i,j)
7269783f6f Jean*0987 #endif
fc7306ba7d Jean*0988
d09af74739 Mart*0989 C-    if hIce > 0 : end
                0990          ENDIF
29188f9691 Jean*0991 c200     CONTINUE
fc7306ba7d Jean*0992
6fc136ac68 Jean*0993 C     endif iceMask > 0
d09af74739 Mart*0994         ENDIF
                0995 C     end i/j-loops
                0996        ENDDO
                0997       ENDDO
a85293d087 Mart*0998 #ifdef ALLOW_AUTODIFF_TAMC_MORE
edb6656069 Mart*0999 CADJ STORE qicen = comlev1_bibj, key = tkey, kind = isbyte
                1000 CADJ STORE hnew  = comlev1_bibj, key = tkey, kind = isbyte
a85293d087 Mart*1001 #endif
1818702d6f Patr*1002
d09af74739 Mart*1003       DO j = jMin, jMax
                1004        DO i = iMin, iMax
                1005         IF (iceMask(i,j).GT.0. _d 0) THEN
fc7306ba7d Jean*1006 C-  Compute surplus energy left over from melting.
                1007
d09af74739 Mart*1008          IF (hIce(i,j).LE.0. _d 0) icFrac(i,j)=0. _d 0
fc7306ba7d Jean*1009
                1010 C.. heat fluxes left over for ocean
6fc136ac68 Jean*1011          flx2oc(i,j) = flx2oc(i,j)
d09af74739 Mart*1012      &        + (Fbot(i,j)+(esurp(i,j)+etop(i,j)+ebot(i,j))/dt)
7269783f6f Jean*1013 #ifdef ALLOW_DBUG_THSICE
d09af74739 Mart*1014          IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
                1015      &        'ThSI_CALC_TH: [esurp,etop+ebot]/dt =',
                1016      &        esurp(i,j)/dt,etop(i,j)/dt,ebot(i,j)/dt
7269783f6f Jean*1017 #endif
fc7306ba7d Jean*1018
d09af74739 Mart*1019 C--   Evaporation left to the ocean :
                1020          frw2oc(i,j) = frw2oc(i,j) - evapLoc(i,j)
                1021 C--   Correct Atmos. fluxes for this different latent heat:
                1022 C     evap was computed over freezing surf.(tSrf<0), latent heat = Lvap+Lfresh
                1023 C     but should be Lvap only for the fraction "evap" that is left to the ocean.
                1024          flx2oc(i,j) = flx2oc(i,j) + evapLoc(i,j)*Lfresh
fc7306ba7d Jean*1025
                1026 C fresh and salt fluxes
c1c3d0f9d7 Patr*1027 c     frw2oc(i,j) = (mwater0(i,j) - (rhos*(hSnow1(i,j))
d09af74739 Mart*1028 c    &              + rhoi*(hIce(i,j))))/dt-evapLoc(i,j)
                1029 c     fsalt = (msalt0(i,j) - rhoi*hIce(i,j)*saltIce)/35. _d 0/dt  ! for same units as frw2oc
                1030 C     note (jmc): frw2oc is computed from a sea-ice mass budget that already
                1031 C     contains, at this point, snow & evaporation (of snow & ice)
                1032 C     but are not meant to be part of ice/ocean fresh-water flux.
                1033 C     fix: a) like below or b) by making the budget before snow/evap is added
c1c3d0f9d7 Patr*1034 c     frw2oc(i,j) = (mwater0(i,j) - (rhos*(hSnow1(i,j)) + rhoi*(hIce(i,j))))/dt
7269783f6f Jean*1035 c    &      + snow(i,j,bi,bj)*rhos - frwAtm
d09af74739 Mart*1036          fsalt(i,j) = (msalt0(i,j) - rhoi*hIce(i,j)*saltIce)/dt
7269783f6f Jean*1037
                1038 #ifdef ALLOW_DBUG_THSICE
d09af74739 Mart*1039          IF (dBug(i,j,bi,bj) ) THEN
                1040           WRITE(6,1020)'ThSI_CALC_TH:dH2O,Ev[kg],frw2oc,fsalt',
c1c3d0f9d7 Patr*1041      &         (mwater0(i,j)-(rhos*hSnow1(i,j)+rhoi*hIce(i,j)))/dt,
d09af74739 Mart*1042      &         evapLoc(i,j),frw2oc(i,j),fsalt(i,j)
                1043           WRITE(6,1020)'ThSI_CALC_TH: flx2oc,Fbot,extend/dt =',
f6de6620bc Mart*1044      &         flx2oc(i,j),Fbot(i,j),(etope(i,j)+ebote(i,j))/dt
d09af74739 Mart*1045          ENDIF
7269783f6f Jean*1046 #endif
fc7306ba7d Jean*1047
d09af74739 Mart*1048 C--   add remaining liquid Precip (rain+RunOff) directly to ocean:
                1049          frw2oc(i,j) = frw2oc(i,j) + (prcAtm(i,j)-snowP(i,j))
                1050
6fc136ac68 Jean*1051 C     endif iceMask > 0
d09af74739 Mart*1052         ENDIF
                1053 C     end i/j-loops
                1054        ENDDO
                1055       ENDDO
1818702d6f Patr*1056
d09af74739 Mart*1057       DO j = jMin, jMax
                1058        DO i = iMin, iMax
a85293d087 Mart*1059 C     this copy really only helps TAF
                1060         icFracTmp = icFrac(i,j)
d09af74739 Mart*1061         IF (iceMask(i,j).GT.0. _d 0) THEN
                1062 C--   note: at this point, icFrac has not been changed (unless reset to zero)
                1063 C     and it can only be reduced by lateral melting in the following part:
                1064
                1065 C     calculate extent changes
                1066          extend=etope(i,j)+ebote(i,j)
a85293d087 Mart*1067          IF (icFracTmp.GT.0. _d 0.AND.extend.GT.0. _d 0) THEN
d09af74739 Mart*1068           rq =  rhoi * 0.5 _d 0*(qicen(i,j,1)+qicen(i,j,2))
                1069           rs =  rhos * qsnow
c1c3d0f9d7 Patr*1070           rqh = rq * hIce(i,j) + rs * hSnow1(i,j)
                1071           freshe=(rhos*hSnow1(i,j)+rhoi*hIce(i,j))/dt
d09af74739 Mart*1072           salte=(rhoi*hIce(i,j)*saltIce)/dt
                1073           IF ( extend.LT.rqh ) THEN
a85293d087 Mart*1074            icFrac(i,j)=(1. _d 0-extend/rqh)*icFracTmp
d09af74739 Mart*1075           ENDIF
                1076           IF ( extend.LT.rqh .AND. icFrac(i,j).GE.iceMaskMin ) THEN
                1077            frw2oc(i,j)=frw2oc(i,j)+extend/rqh*freshe
7269783f6f Jean*1078            fsalt(i,j)=fsalt(i,j)+extend/rqh*salte
d09af74739 Mart*1079           ELSE
                1080            icFrac(i,j)=0. _d 0
                1081            hIce(i,j)  =0. _d 0
c1c3d0f9d7 Patr*1082            hSnow1(i,j) =0. _d 0
d09af74739 Mart*1083            flx2oc(i,j)=flx2oc(i,j)+(extend-rqh)/dt
                1084            frw2oc(i,j)=frw2oc(i,j)+freshe
7269783f6f Jean*1085            fsalt(i,j)=fsalt(i,j)+salte
d09af74739 Mart*1086           ENDIF
                1087          ELSEIF (extend.GT.0. _d 0) THEN
                1088           flx2oc(i,j)=flx2oc(i,j)+extend/dt
9a68c0a761 Jean*1089          ENDIF
6fc136ac68 Jean*1090 C     endif iceMask > 0
d09af74739 Mart*1091         ENDIF
                1092 C     end i/j-loops
                1093        ENDDO
                1094       ENDDO
                1095       DO j = jMin, jMax
                1096        DO i = iMin, iMax
                1097         IF (iceMask(i,j).GT.0. _d 0) THEN
7269783f6f Jean*1098 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
d09af74739 Mart*1099 C--   Update output variables :
7269783f6f Jean*1100
d09af74739 Mart*1101 C--   Diagnostic of Atmos. fresh water flux (E-P) over sea-ice :
7269783f6f Jean*1102 C     substract precip from Evap (<- stored in frwAtm array)
d09af74739 Mart*1103          frwAtm(i,j) = frwAtm(i,j) - prcAtm(i,j)
7269783f6f Jean*1104
d09af74739 Mart*1105 C--   update Mixed-Layer Freezing potential heat flux by substracting the
                1106 C     part which has already been accounted for (Fbot):
                1107          fzMlOc(i,j) = fzMlOc(i,j) - Fbot(i,j)*iceMask(i,j)
7269783f6f Jean*1108
d09af74739 Mart*1109 C-- Update Sea-Ice state output:
                1110          qIc1(i,j)   = qicen(i,j,1)
                1111          qIc2(i,j)   = qicen(i,j,2)
7269783f6f Jean*1112 #ifdef ALLOW_DBUG_THSICE
d09af74739 Mart*1113          IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
                1114      &        'ThSI_CALC_TH: icFrac,flx2oc,fsalt,frw2oc=',
                1115      &        icFrac(i,j), flx2oc(i,j), fsalt(i,j), frw2oc(i,j)
7269783f6f Jean*1116 #endif
d09af74739 Mart*1117 C     endif iceMask > 0
                1118         ENDIF
                1119        ENDDO
                1120       ENDDO
6fc136ac68 Jean*1121
7269783f6f Jean*1122 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
6fc136ac68 Jean*1123 #ifdef CHECK_ENERGY_CONSERV
d09af74739 Mart*1124       DO j = jMin, jMax
                1125        DO i = iMin, iMax
                1126         IF (iceMask(i,j).GT.0. _d 0) THEN
                1127          qaux(1)=qIc1(i,j)
                1128          qaux(2)=qIc2(i,j)
                1129          CALL THSICE_CHECK_CONSERV( dBugFlag, i, j, bi, bj, 0,
c1c3d0f9d7 Patr*1130      I        iceMask(i,j), icFrac(i,j), hIce(i,j), hSnow1(i,j),
d09af74739 Mart*1131      I        qaux,
                1132      I        flx2oc(i,j), frw2oc(i,j), fsalt,
                1133      I        myTime, myIter, myThid )
6fc136ac68 Jean*1134 C     endif iceMask > 0
7269783f6f Jean*1135         ENDIF
d09af74739 Mart*1136 C     end i/j-loops
7269783f6f Jean*1137        ENDDO
                1138       ENDDO
d09af74739 Mart*1139 #endif /* CHECK_ENERGY_CONSERV */
                1140
fc7306ba7d Jean*1141 #endif  /* ALLOW_THSICE */
                1142
                1143 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                1144
                1145       RETURN
                1146       END