C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_calc_thickn.F,v 1.30 2014/01/19 14:47:35 jmc Exp $
C $Name: $
#include "THSICE_OPTIONS.h"
CBOP
C !ROUTINE: THSICE_CALC_THICKN
C !INTERFACE:
SUBROUTINE THSICE_CALC_THICKN(
I bi, bj,
I iMin,iMax, jMin,jMax, dBugFlag,
I iceMask, tFrz, tOce, v2oc,
I snowP, prcAtm, sHeat, flxCnB,
U icFrac, hIce, hSnow1, tSrf, qIc1, qIc2,
U frwAtm, fzMlOc, flx2oc,
O frw2oc, fsalt, frzSeaWat,
I myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | S/R THSICE_CALC_THICKN
C | o Calculate ice & snow thickness changes
C *==========================================================*
C \ev
C ADAPTED FROM:
C LANL CICE.v2.0.2
C-----------------------------------------------------------------------
C.. thermodynamics (vertical physics) based on M. Winton 3-layer model
C.. See Bitz, C. M. and W. H. Lipscomb, 1999: An energy-conserving
C.. thermodynamic sea ice model for climate study.
C.. J. Geophys. Res., 104, 15669 - 15677.
C.. Winton, M., 1999: "A reformulated three-layer sea ice model."
C.. Submitted to J. Atmos. Ocean. Technol.
C.. authors Elizabeth C. Hunke and William Lipscomb
C.. Fluid Dynamics Group, Los Alamos National Laboratory
C-----------------------------------------------------------------------
Cc****subroutine thermo_winton(n,fice,fsnow,dqice,dTsfc)
C.. Compute temperature change using Winton model with 2 ice layers, of
C.. which only the top layer has a variable heat capacity.
C---------------------------------
C parameters that control the partitioning between lateral (ice area) and
C vertical (ice thickness) ice volume changes.
C a) surface melting and bottom melting (runtime parameter: fracEnMelt):
C frace is the fraction of available heat that is used for
C lateral melting (and 1-frace reduces the thickness ) when
C o hi < hThinIce & frac > lowIcFrac2 : frace=1 (lateral melting only)
C o hThinIce < hi < hThickIce & frac > lowIcFrac1 : frace=fracEnMelt
C o hi > hThickIce or frac < lowIcFrac1 : frace=0 (thinning only)
C b) ocean freezing (and ice forming):
C - conductive heat flux (below sea-ice) always increases thickness.
C - under sea-ice, freezing potential (x iceFraction) is used to increase ice
C thickness or ice fraction (lateral growth), according to:
C o hi < hThinIce : use freezing potential to grow ice vertically;
C o hThinIce < hi < hThickIce : use partition factor fracEnFreez for lateral growth
c and (1-fracEnFreez) to increase thickness.
C o hi > hThickIce : use all freezing potential to grow ice laterally
C (up to areaMax)
C - over open ocean, use freezing potential [x(1-iceFraction)] to grow ice laterally
C - lateral growth forms ice of the same or =hNewIceMax thickness, the less of the 2.
C - starts to form sea-ice over fraction iceMaskMin, as minimum ice-volume is reached
C---------------------------------
C !USES:
IMPLICIT NONE
C == Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "THSICE_SIZE.h"
#include "THSICE_PARAMS.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
# include "tamc_keys.h"
#endif
C !INPUT/OUTPUT PARAMETERS:
C == Routine Arguments ==
C bi,bj :: tile indices
C iMin,iMax :: computation domain: 1rst index range
C jMin,jMax :: computation domain: 2nd index range
C dBugFlag :: allow to print debugging stuff (e.g. on 1 grid point).
C--- Input:
C iceMask :: sea-ice fractional mask [0-1]
C tFrz :: sea-water freezing temperature [oC] (function of S)
C tOce :: surface level oceanic temperature [oC]
C v2oc :: square of ocean surface-level velocity [m2/s2]
C snowP :: snow precipitation [kg/m2/s]
C prcAtm :: total precip from the atmosphere [kg/m2/s]
C sHeat :: surf heating flux left to melt snow or ice (= Atmos-conduction)
C flxCnB :: heat flux conducted through the ice to bottom surface
C--- Modified (input&output):
C icFrac :: fraction of grid area covered in ice
C hIce :: ice height [m]
C hSnow1 :: snow height [m]
C tSrf :: surface (ice or snow) temperature
C qIc1 (qicen) :: ice enthalpy (J/kg), 1rst level
C qIc2 (qicen) :: ice enthalpy (J/kg), 2nd level
C frwAtm (evpAtm):: evaporation to the atmosphere [kg/m2/s] (>0 if evaporate)
C fzMlOc :: ocean mixed-layer freezing/melting potential [W/m2]
C flx2oc :: net heat flux to ocean [W/m2] (> 0 downward)
C--- Output
C frw2oc :: Total fresh water flux to ocean [kg/m2/s] (> 0 downward)
C fsalt :: salt flux to ocean [g/m2/s] (> 0 downward)
C frzSeaWat :: seawater freezing rate (expressed as mass flux) [kg/m^2/s]
C--- Input:
C myTime :: current Time of simulation [s]
C myIter :: current Iteration number in simulation
C myThid :: my Thread Id number
INTEGER bi,bj
INTEGER iMin, iMax
INTEGER jMin, jMax
LOGICAL dBugFlag
_RL iceMask(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL tFrz (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL tOce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL v2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL snowP (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL prcAtm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL sHeat (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL flxCnB (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL icFrac (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL hIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL hSnow1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL tSrf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL qIc1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL qIc2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL frwAtm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL fzMlOc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL flx2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL frw2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL fsalt (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL frzSeaWat(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL myTime
INTEGER myIter
INTEGER myThid
CEOP
#ifdef ALLOW_THSICE
C !LOCAL VARIABLES:
C--- local copy of input/output argument list variables (see description above)
_RL qicen(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nlyr)
C == Local Variables ==
C i,j,k :: loop indices
C rec_nlyr :: reciprocal of number of ice layers (real value)
C evapLoc :: evaporation over snow/ice [kg/m2/s] (>0 if evaporate)
C Fbot :: oceanic heat flux used to melt/form ice [W/m2]
C etop :: energy for top melting (J m-2)
C ebot :: energy for bottom melting (J m-2)
C etope :: energy (from top) for lateral melting (J m-2)
C ebote :: energy (from bottom) for lateral melting (J m-2)
C extend :: total energy for lateral melting (J m-2)
C hnew(nlyr) :: new ice layer thickness (m)
C hlyr :: individual ice layer thickness (m)
C dhi :: change in ice thickness
C dhs :: change in snow thickness
C rq :: rho * q for a layer
C rqh :: rho * q * h for a layer
C qbot :: enthalpy for new ice at bottom surf (J/kg)
C dt :: timestep
C esurp :: surplus energy from melting (J m-2)
C mwater0 :: fresh water mass gained/lost (kg/m^2)
C msalt0 :: salt gained/lost (kg/m^2)
C freshe :: fresh water gain from extension melting
C salte :: salt gained from extension melting
C lowIcFrac1 :: ice-fraction lower limit above which partial (lowIcFrac1)
C lowIcFrac2 :: or full (lowIcFrac2) lateral melting is allowed.
C from THSICE_RESHAPE_LAYERS
C f1 :: Fraction of upper layer ice in new layer
C qh1, qh2 :: qice*h for layers 1 and 2
C qhtot :: qh1 + qh2
C q2tmp :: Temporary value of qice for layer 2
INTEGER i,j,k
_RL rec_nlyr
_RL evapLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL Fbot (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL etop (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL ebot (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL etope (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL ebote (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL esurp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL extend
_RL hnew (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nlyr)
_RL hlyr
_RL dhi
_RL dhs
_RL rq
_RL rqh
_RL qbot
_RL dt
_RL mwater0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL msalt0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL freshe
_RL salte
_RL lowIcFrac1, lowIcFrac2
_RL f1
_RL qh1, qh2
_RL qhtot
_RL q2tmp
#ifdef CHECK_ENERGY_CONSERV
_RL qaux(nlyr)
#endif /* CHECK_ENERGY_CONSERV */
_RL ustar, cpchr
_RL chi
_RL frace, rs, hq
#ifdef THSICE_FRACEN_POWERLAW
INTEGER powerLaw
_RL rec_pLaw
_RL c1Mlt, c2Mlt, aMlt, hMlt
_RL c1Frz, c2Frz, aFrz, hFrz
_RL enFrcMlt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL xxMlt, tmpMlt
_RL enFrcFrz(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL xxFrz, tmpFrz
#endif
C- define grid-point location where to print debugging values
#include "THSICE_DEBUG.h"
1020 FORMAT(A,1P4E11.3)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_AUTODIFF_TAMC
act1 = bi - myBxLo(myThid)
max1 = myBxHi(myThid) - myBxLo(myThid) + 1
act2 = bj - myByLo(myThid)
max2 = myByHi(myThid) - myByLo(myThid) + 1
act3 = myThid - 1
max3 = nTx*nTy
act4 = ikey_dynamics - 1
ticekey = (act1 + 1) + act2*max1
& + act3*max1*max2
& + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */
rec_nlyr = nlyr
rec_nlyr = 1. _d 0 / rec_nlyr
dt = thSIce_deltaT
C for now, use hard coded threshold (iceMaskMin +1.% and +10.%)
lowIcFrac1 = iceMaskMin*1.01 _d 0
lowIcFrac2 = iceMaskMin*1.10 _d 0
#ifdef THSICE_FRACEN_POWERLAW
IF ( powerLawExp2 .GE. 1 ) THEN
powerLaw = 1 + 2**powerLawExp2
rec_pLaw = powerLaw
rec_pLaw = 1. _d 0 / rec_pLaw
C- Coef for melting:
C lateral-melting energy fraction = fracEnMelt - [ aMlt*(hi-hMlt) ]^powerLaw
c1Mlt = fracEnMelt**rec_pLaw
c2Mlt = (1. _d 0 - fracEnMelt)**rec_pLaw
aMlt = (c1Mlt+c2Mlt)/(hThickIce-hThinIce)
hMlt = hThinIce+c2Mlt/aMlt
C- Coef for freezing:
C thickening energy fraction = fracEnFreez - [ aFrz*(hi-hFrz) ]^powerLaw
c1Frz = fracEnFreez**rec_pLaw
c2Frz = (1. _d 0 -fracEnFreez)**rec_pLaw
aFrz = (c1Frz+c2Frz)/(hThickIce-hThinIce)
hFrz = hThinIce+c2Frz/aFrz
ELSE
C- Linear relation
powerLaw = 1
aMlt = -1. _d 0 /(hThickIce-hThinIce)
hMlt = hThickIce
aFrz = -1. _d 0 /(hThickIce-hThinIce)
hFrz = hThickIce
ENDIF
#endif /* THSICE_FRACEN_POWERLAW */
C initialise local arrays
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
evapLoc(i,j) = 0. _d 0
Fbot (i,j) = 0. _d 0
etop (i,j) = 0. _d 0
ebot (i,j) = 0. _d 0
etope (i,j) = 0. _d 0
ebote (i,j) = 0. _d 0
esurp (i,j) = 0. _d 0
mwater0(i,j) = 0. _d 0
msalt0 (i,j) = 0. _d 0
#ifdef THSICE_FRACEN_POWERLAW
enFrcMlt(i,j)= 0. _d 0
enFrcFrz(i,j)= 0. _d 0
#endif
ENDDO
ENDDO
DO k = 1,nlyr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
qicen(i,j,k) = 0. _d 0
hnew (i,j,k) = 0. _d 0
ENDDO
ENDDO
ENDDO
DO j = jMin, jMax
DO i = iMin, iMax
CML#ifdef ALLOW_AUTODIFF_TAMC
CML ikey_1 = i
CML & + sNx*(j-1)
CML & + sNx*sNy*act1
CML & + sNx*sNy*max1*act2
CML & + sNx*sNy*max1*max2*act3
CML & + sNx*sNy*max1*max2*max3*act4
CML#endif /* ALLOW_AUTODIFF_TAMC */
CML#ifdef ALLOW_AUTODIFF_TAMC
CMLCADJ STORE frwatm(i,j) = comlev1_thsice_1, key=ikey_1
CMLCADJ STORE fzmloc(i,j) = comlev1_thsice_1, key=ikey_1
CMLCADJ STORE hice(i,j) = comlev1_thsice_1, key=ikey_1
CMLCADJ STORE hSnow1(i,j) = comlev1_thsice_1, key=ikey_1
CMLCADJ STORE icfrac(i,j) = comlev1_thsice_1, key=ikey_1
CMLCADJ STORE qic1(i,j) = comlev1_thsice_1, key=ikey_1
CMLCADJ STORE qic2(i,j) = comlev1_thsice_1, key=ikey_1
CML#endif
IF (iceMask(i,j).GT.0. _d 0) THEN
qicen(i,j,1)= qIc1(i,j)
qicen(i,j,2)= qIc2(i,j)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C initialize energies
esurp(i,j) = 0. _d 0
c make a local copy of evaporation
evapLoc(i,j) = frwAtm(i,j)
C------------------------------------------------------------------------
C-- Compute growth and/or melting at the top and bottom surfaces
C------------------------------------------------------------------------
#ifdef THSICE_FRACEN_POWERLAW
xxMlt = aMlt*(hIce(i,j)-hMlt)
xxFrz = aFrz*(hIce(i,j)-hFrz)
c--
IF ( powerLawExp2 .GE. 1 ) THEN
#ifdef TARGET_NEC_SX
C avoid the short inner loop that cannot be vectorized
xxMlt = xxMlt**powerLaw
xxFrz = xxFrz**powerLaw
#else
tmpMlt = xxMlt
tmpFrz = xxFrz
DO k=1,powerLawExp2
tmpMlt = tmpMlt*tmpMlt
tmpFrz = tmpFrz*tmpFrz
ENDDO
xxMlt = xxMlt*tmpMlt
xxFrz = xxFrz*tmpFrz
#endif /* TARGET_NEC_SX */
xxMlt = fracEnMelt -xxMlt
xxFrz = fracEnFreez-xxFrz
ENDIF
enFrcMlt(i,j) = MAX( 0. _d 0, MIN( xxMlt, 1. _d 0 ) )
enFrcFrz(i,j) = MAX( 0. _d 0, MIN( xxFrz, 1. _d 0 ) )
#endif /* THSICE_FRACEN_POWERLAW */
IF (fzMlOc(i,j).GE. 0. _d 0) THEN
C !-----------------------------------------------------------------
C ! freezing conditions
C !-----------------------------------------------------------------
Fbot(i,j) = fzMlOc(i,j)
IF ( icFrac(i,j).LT.iceMaskMax ) THEN
#ifdef THSICE_FRACEN_POWERLAW
Fbot(i,j) = enFrcFrz(i,j)*fzMlOc(i,j)
#else /* THSICE_FRACEN_POWERLAW */
IF (hIce(i,j).GT.hThickIce) THEN
C if higher than hThickIce, use all fzMlOc energy to grow extra ice
Fbot(i,j) = 0. _d 0
ELSEIF (hIce(i,j).GE.hThinIce) THEN
C between hThinIce & hThickIce, use partition factor fracEnFreez
Fbot(i,j) = (1. _d 0 - fracEnFreez)*fzMlOc(i,j)
ENDIF
#endif /* THSICE_FRACEN_POWERLAW */
ENDIF
ELSE
C !-----------------------------------------------------------------
C ! melting conditions
C !-----------------------------------------------------------------
C for no currents:
ustar = 5. _d -2
C frictional velocity between ice and water
IF (v2oc(i,j) .NE. 0.)
& ustar = SQRT(0.00536 _d 0*v2oc(i,j))
ustar=max(5. _d -3,ustar)
cpchr =cpWater*rhosw*bMeltCoef
Fbot(i,j) = cpchr*(tFrz(i,j)-tOce(i,j))*ustar
C fzMlOc < Fbot < 0
Fbot(i,j) = max(Fbot(i,j),fzMlOc(i,j))
Fbot(i,j) = min(Fbot(i,j),0. _d 0)
ENDIF
C mass of fresh water and salt initially present in ice
mwater0(i,j) = rhos*hSnow1(i,j) + rhoi*hIce(i,j)
msalt0 (i,j) = rhoi*hIce(i,j)*saltIce
#ifdef ALLOW_DBUG_THSICE
IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
& 'ThSI_CALC_TH: evpAtm, fzMlOc, Fbot =',
& frwAtm(i,j),fzMlOc(i,j),Fbot(i,j)
#endif
C endif iceMask > 0
ENDIF
C end i/j-loops
ENDDO
ENDDO
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
DO j = jMin, jMax
DO i = iMin, iMax
IF (iceMask(i,j).GT.0. _d 0) THEN
C Compute energy available for melting/growth.
#ifdef THSICE_FRACEN_POWERLAW
IF ( fracEnMelt.EQ.0. _d 0 ) THEN
frace = 0. _d 0
ELSE
frace = (icFrac(i,j) - lowIcFrac1)/(lowIcFrac2-iceMaskMin)
frace = MIN( enFrcMlt(i,j), MAX( 0. _d 0, frace ) )
ENDIF
#else /* THSICE_FRACEN_POWERLAW */
IF ( hIce(i,j).GT.hThickIce .OR. fracEnMelt.EQ.0. _d 0 ) THEN
C above certain height (or when no ice fractionation), only melt from top
frace = 0. _d 0
ELSEIF (hIce(i,j).LT.hThinIce) THEN
C below a certain height, all energy goes to changing ice extent
frace = 1. _d 0
ELSE
frace = fracEnMelt
ENDIF
C Reduce lateral melting when ice fraction is low : the purpose is to avoid
C disappearing of (up to hThinIce thick) sea-ice by over doing lateral melting
C (which would bring icFrac below iceMaskMin).
IF ( icFrac(i,j).LE.lowIcFrac1 ) THEN
frace = 0. _d 0
ELSEIF (icFrac(i,j).LE.lowIcFrac2 ) THEN
frace = MIN( frace, fracEnMelt )
ENDIF
#endif /* THSICE_FRACEN_POWERLAW */
c IF (tSrf(i,j) .EQ. 0. _d 0 .AND. sHeat(i,j).GT.0. _d 0) THEN
IF ( sHeat(i,j).GT.0. _d 0 ) THEN
etop(i,j) = (1. _d 0-frace)*sHeat(i,j) * dt
etope(i,j) = frace*sHeat(i,j) * dt
ELSE
etop(i,j) = 0. _d 0
etope(i,j) = 0. _d 0
C jmc: found few cases where tSrf=0 & sHeat < 0 : add this line to conserv energy:
esurp(i,j) = sHeat(i,j) * dt
ENDIF
C-- flux at the base of sea-ice:
C conduction H.flx= flxCnB (+ =down); oceanic turbulent H.flx= Fbot (+ =down).
C- ==> energy available(+ => melt)= (flxCnB-Fbot)*dt
c IF (fzMlOc(i,j).LT.0. _d 0) THEN
c ebot(i,j) = (1. _d 0-frace)*(flxCnB-Fbot(i,j)) * dt
c ebote(i,j) = frace*(flxCnB-Fbot(i,j)) * dt
c ELSE
c ebot(i,j) = (flxCnB-Fbot(i,j)) * dt
c ebote(i,j) = 0. _d 0
c ENDIF
C- original formulation(above): Loose energy when flxCnB < Fbot < 0
ebot(i,j) = (flxCnB(i,j)-Fbot(i,j)) * dt
IF (ebot(i,j).GT.0. _d 0) THEN
ebote(i,j) = frace*ebot(i,j)
ebot(i,j) = ebot(i,j)-ebote(i,j)
ELSE
ebote(i,j) = 0. _d 0
ENDIF
#ifdef ALLOW_DBUG_THSICE
IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
& 'ThSI_CALC_TH: etop,etope,ebot,ebote=',
& etop(i,j),etope(i,j),ebot(i,j),ebote(i,j)
#endif
C endif iceMask > 0
ENDIF
C end i/j-loops
ENDDO
ENDDO
C Initialize layer thicknesses. Divide total thickness equally between
C layers
DO k = 1, nlyr
DO j = jMin, jMax
DO i = iMin, iMax
hnew(i,j,k) = hIce(i,j) * rec_nlyr
ENDDO
ENDDO
ENDDO
DO j = jMin, jMax
DO i = iMin, iMax
CML#ifdef ALLOW_AUTODIFF_TAMC
CMLCADJ STORE etop(i,j) = comlev1_thsice_1, key=ikey_1
CML#endif
IF (iceMask(i,j) .GT. 0. _d 0 .AND.
& etop(i,j) .GT. 0. _d 0 .AND.
& hSnow1(i,j) .GT. 0. _d 0) THEN
C Make sure internal ice temperatures do not exceed Tmlt.
C If they do, then eliminate the layer. (Dont think this will happen
C for reasonable values of i0.)
C Top melt: snow, then ice.
rq = rhos * qsnow
rqh = rq * hSnow1(i,j)
IF (etop(i,j) .LT. rqh) THEN
hSnow1(i,j) = hSnow1(i,j) - etop(i,j)/rq
etop(i,j) = 0. _d 0
ELSE
etop(i,j) = etop(i,j) - rqh
hSnow1(i,j) = 0. _d 0
ENDIF
C endif iceMask > 0, etc.
ENDIF
C end i/j-loops
ENDDO
ENDDO
C two layers of ice
DO k = 1, nlyr
DO j = jMin, jMax
DO i = iMin, iMax
IF (iceMask(i,j).GT.0. _d 0) THEN
CML#ifdef ALLOW_AUTODIFF_TAMC
CML ikey_2 = k
CML & + nlyr*(i-1)
CML & + nlyr*sNx*(j-1)
CML & + nlyr*sNx*sNy*act1
CML & + nlyr*sNx*sNy*max1*act2
CML & + nlyr*sNx*sNy*max1*max2*act3
CML & + nlyr*sNx*sNy*max1*max2*max3*act4
CML#endif
CML#ifdef ALLOW_AUTODIFF_TAMC
CMLCADJ STORE etop(i,j) = comlev1_thsice_2, key=ikey_2
CMLCADJ STORE hnew(i,j,k) = comlev1_thsice_2, key=ikey_2
CML#endif
IF (etop(i,j) .GT. 0. _d 0) THEN
rq = rhoi * qicen(i,j,k)
rqh = rq * hnew(i,j,k)
IF (etop(i,j) .LT. rqh) THEN
hnew(i,j,k) = hnew(i,j,k) - etop(i,j) / rq
etop(i,j) = 0. _d 0
ELSE
etop(i,j) = etop(i,j) - rqh
hnew(i,j,k) = 0. _d 0
ENDIF
ELSE
etop(i,j)=0. _d 0
ENDIF
C If ice is gone and melting energy remains
c IF (etop(i,j) .GT. 0. _d 0) THEN
c WRITE (6,*) 'QQ All ice melts from top ', i,j
c hIce(i,j)=0. _d 0
c go to 200
c ENDIF
C endif iceMask > 0
ENDIF
C end i/j-loops
ENDDO
ENDDO
C end k-loop
ENDDO
C Bottom growth:
DO j = jMin, jMax
DO i = iMin, iMax
IF (iceMask(i,j).GT.0. _d 0 .AND. ebot(i,j) .LT. 0. _d 0) THEN
C Compute enthalpy of new ice growing at bottom surface.
qbot = -cpIce *tFrz(i,j) + Lfresh
dhi = -ebot(i,j) / (qbot * rhoi)
ebot(i,j) = 0. _d 0
cph k = nlyr
CML#ifdef ALLOW_AUTODIFF_TAMC
CMLCADJ STORE hnew(i,j,:) = comlev1_thsice_1, key=ikey_1
CMLCADJ STORE qicen(i,j,:) = comlev1_thsice_1, key=ikey_1
CML#endif
qicen(i,j,nlyr) =
& (hnew(i,j,nlyr)*qicen(i,j,nlyr)+dhi*qbot) /
& (hnew(i,j,nlyr)+dhi)
hnew(i,j,nlyr) = hnew(i,j,nlyr) + dhi
frzSeaWat(i,j) = rhoi*dhi/dt
C endif iceMask > 0 and ebot < 0
ENDIF
C end i/j-loops
ENDDO
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE etop(:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE ebot(:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE hnew(:,:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE qicen(:,:,:) = comlev1_bibj, key=ticekey, byte=isbyte
#endif
C Bottom melt:
DO k = nlyr, 1, -1
CML#ifdef ALLOW_AUTODIFF_TAMC
CML ikey_2 = (nlyr-k+1)
CML & + nlyr*(i-1)
CML & + nlyr*sNx*(j-1)
CML & + nlyr*sNx*sNy*act1
CML & + nlyr*sNx*sNy*max1*act2
CML & + nlyr*sNx*sNy*max1*max2*act3
CML & + nlyr*sNx*sNy*max1*max2*max3*act4
CML#endif
CML#ifdef ALLOW_AUTODIFF_TAMC
CMLCADJ STORE ebot(i,j) = comlev1_thsice_2, key=ikey_2
CMLCADJ STORE hnew(i,j,k) = comlev1_thsice_2, key=ikey_2
CMLCADJ STORE qicen(i,j,k) = comlev1_thsice_2, key=ikey_2
CML#endif
DO j = jMin, jMax
DO i = iMin, iMax
IF (iceMask(i,j) .GT. 0. _d 0 .AND.
& ebot(i,j) .GT. 0. _d 0 .AND.
& hnew(i,j,k) .GT. 0. _d 0) THEN
rq = rhoi * qicen(i,j,k)
rqh = rq * hnew(i,j,k)
IF (ebot(i,j) .LT. rqh) THEN
hnew(i,j,k) = hnew(i,j,k) - ebot(i,j) / rq
ebot(i,j) = 0. _d 0
ELSE
ebot(i,j) = ebot(i,j) - rqh
hnew(i,j,k) = 0. _d 0
ENDIF
C endif iceMask > 0 etc.
ENDIF
C end i/j-loops
ENDDO
ENDDO
C end k-loop
ENDDO
C If ice melts completely and snow is left, remove the snow with
C energy from the mixed layer
DO j = jMin, jMax
DO i = iMin, iMax
IF (iceMask(i,j) .GT. 0. _d 0 .AND.
& ebot(i,j) .GT. 0. _d 0 .AND.
& hSnow1(i,j) .GT. 0. _d 0) THEN
rq = rhos * qsnow
rqh = rq * hSnow1(i,j)
IF (ebot(i,j) .LT. rqh) THEN
hSnow1(i,j) = hSnow1(i,j) - ebot(i,j) / rq
ebot(i,j) = 0. _d 0
ELSE
ebot(i,j) = ebot(i,j) - rqh
hSnow1(i,j) = 0. _d 0
ENDIF
c IF (ebot(i,j) .GT. 0. _d 0) THEN
c IF (dBug) WRITE(6,*) 'All ice (& snow) melts from bottom'
c hIce(i,j)=0. _d 0
c go to 200
c ENDIF
C endif iceMask > 0, etc.
ENDIF
C end i/j-loops
ENDDO
ENDDO
DO j = jMin, jMax
DO i = iMin, iMax
IF (iceMask(i,j).GT.0. _d 0) THEN
C Compute new total ice thickness.
hIce(i,j) = hnew(i,j,1) + hnew(i,j,2)
#ifdef ALLOW_DBUG_THSICE
IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
& 'ThSI_CALC_TH: etop, ebot, hIce, hSnow1 =',
& etop(i,j), ebot(i,j), hIce(i,j), hSnow1(i,j)
#endif
C If hIce < hIceMin, melt the ice.
IF ( hIce(i,j).LT.hIceMin
& .AND. (hIce(i,j)+hSnow1(i,j)).GT.0. _d 0 ) THEN
esurp(i,j) = esurp(i,j) - rhos*qsnow*hSnow1(i,j)
& - rhoi*qicen(i,j,1)*hnew(i,j,1)
& - rhoi*qicen(i,j,2)*hnew(i,j,2)
hIce(i,j) = 0. _d 0
hSnow1(i,j) = 0. _d 0
tSrf(i,j) = 0. _d 0
icFrac(i,j) = 0. _d 0
qicen(i,j,1) = 0. _d 0
qicen(i,j,2) = 0. _d 0
#ifdef ALLOW_DBUG_THSICE
IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
& 'ThSI_CALC_TH: -1 : esurp=',esurp(i,j)
#endif
ENDIF
C endif iceMask > 0
ENDIF
C end i/j-loops
ENDDO
ENDDO
DO j = jMin, jMax
DO i = iMin, iMax
IF (iceMask(i,j).GT.0. _d 0) THEN
C-- do a mass-budget of sea-ice to compute "fresh" = the fresh-water flux
C that is returned to the ocean ; needs to be done before snow/evap
frw2oc(i,j) = (mwater0(i,j)
& - (rhos*hSnow1(i,j)+rhoi*hIce(i,j)))/dt
IF ( hIce(i,j) .LE. 0. _d 0 ) THEN
C- return snow to the ocean (account for Latent heat of freezing)
frw2oc(i,j) = frw2oc(i,j) + snowP(i,j)
flx2oc(i,j) = flx2oc(i,j) - snowP(i,j)*Lfresh
ENDIF
C endif iceMask > 0
ENDIF
C end i/j-loops
ENDDO
ENDDO
C- else: hIce > 0
DO j = jMin, jMax
DO i = iMin, iMax
IF (iceMask(i,j).GT.0. _d 0) THEN
IF ( hIce(i,j) .GT. 0. _d 0 ) THEN
C Let it snow
hSnow1(i,j) = hSnow1(i,j) + dt*snowP(i,j)/rhos
C If ice evap is used to sublimate surface snow/ice or
C if no ice pass on to ocean
CML#ifdef ALLOW_AUTODIFF_TAMC
CMLCADJ STORE evapLoc(i,j) = comlev1_thsice_1, key=ikey_1
CMLCADJ STORE hSnow1(i,j) = comlev1_thsice_1, key=ikey_1
CML#endif
IF (hSnow1(i,j).GT.0. _d 0) THEN
IF (evapLoc(i,j)/rhos *dt.GT.hSnow1(i,j)) THEN
evapLoc(i,j)=evapLoc(i,j)-hSnow1(i,j)*rhos/dt
hSnow1(i,j)=0. _d 0
ELSE
hSnow1(i,j) = hSnow1(i,j) - evapLoc(i,j)/rhos *dt
evapLoc(i,j)=0. _d 0
ENDIF
ENDIF
C endif hice > 0
ENDIF
C endif iceMask > 0
ENDIF
C end i/j-loops
ENDDO
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE evaploc(:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE hnew(:,:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE hSnow1(:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE icfrac(:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE hice(:,:) = comlev1_bibj, key=ticekey, byte=isbyte
#endif
C- else: hIce > 0
DO k = 1, nlyr
DO j = jMin, jMax
DO i = iMin, iMax
IF (iceMask(i,j).GT.0. _d 0 ) THEN
CML#ifdef ALLOW_AUTODIFF_TAMC
CMLCADJ STORE evapLoc(i,j) = comlev1_thsice_1, key=ikey_1
CMLCADJ STORE hIce(i,j) = comlev1_thsice_1, key=ikey_1
CMLCADJ STORE hnew(i,j,:) = comlev1_thsice_1, key=ikey_1
CMLCADJ STORE qicen(i,j,:) = comlev1_thsice_1, key=ikey_1
CML#endif
IF (hIce(i,j).GT.0. _d 0.AND.evapLoc(i,j).GT.0. _d 0) THEN
CML#ifdef ALLOW_AUTODIFF_TAMC
CML ikey_2 = k
CML & + nlyr*(i-1)
CML & + nlyr*sNx*(j-1)
CML & + nlyr*sNx*sNy*act1
CML & + nlyr*sNx*sNy*max1*act2
CML & + nlyr*sNx*sNy*max1*max2*act3
CML & + nlyr*sNx*sNy*max1*max2*max3*act4
CML#endif
CMLC--
CML#ifdef ALLOW_AUTODIFF_TAMC
CMLCADJ STORE evapLoc(i,j) = comlev1_thsice_2, key=ikey_2
CMLCADJ STORE hnew(i,j,k) = comlev1_thsice_2, key=ikey_2
CMLCADJ STORE qicen(i,j,k) = comlev1_thsice_2, key=ikey_2
CML#endif
C IF (evapLoc(i,j) .GT. 0. _d 0) THEN
C-- original scheme, does not care about ice temp.
C- this can produce small error (< 1.W/m2) in the Energy budget
c IF (evapLoc(i,j)/rhoi *dt.GT.hnew(i,j,k)) THEN
c evapLoc(i,j)=evapLoc(i,j)-hnew(i,j,k)*rhoi/dt
c hnew(i,j,k)=0. _d 0
c ELSE
c hnew(i,j,k) = hnew(i,j,k) - evapLoc(i,j)/rhoi *dt
c evapLoc(i,j)=0. _d 0
c ENDIF
C-- modified scheme. taking into account Ice enthalpy
dhi = evapLoc(i,j)/rhoi*dt
IF (dhi.GE.hnew(i,j,k)) THEN
evapLoc(i,j)=evapLoc(i,j)-hnew(i,j,k)*rhoi/dt
esurp(i,j) = esurp(i,j)
& - hnew(i,j,k)*rhoi*(qicen(i,j,k)-Lfresh)
hnew(i,j,k)=0. _d 0
ELSE
CML#ifdef ALLOW_AUTODIFF_TAMC
CMLCADJ STORE hnew(i,j,k) = comlev1_thsice_2, key=ikey_2
CML#endif
hq = hnew(i,j,k)*qicen(i,j,k)-dhi*Lfresh
hnew(i,j,k) = hnew(i,j,k) - dhi
qicen(i,j,k)=hq/hnew(i,j,k)
evapLoc(i,j)=0. _d 0
ENDIF
C-------
c IF (evapLoc(i,j) .GT. 0. _d 0) THEN
c WRITE (6,*) 'BB All ice sublimates', i,j
c hIce(i,j)=0. _d 0
c go to 200
c ENDIF
C endif hice > 0 and evaploc > 0
ENDIF
C endif iceMask > 0
ENDIF
C end i/j-loops
ENDDO
ENDDO
C end k-loop
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE etop(:,:) = comlev1_bibj, key=ticekey, byte=isbyte
cphCADJ STORE icemask(:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE hice(:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE hnew(:,:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE qicen(:,:,:) = comlev1_bibj, key=ticekey, byte=isbyte
#endif
C still else: hice > 0
DO j = jMin, jMax
DO i = iMin, iMax
IF (iceMask(i,j).GT.0. _d 0) THEN
IF (hIce(i,j) .GT. 0. _d 0) THEN
C Compute new total ice thickness.
hIce(i,j) = hnew(i,j,1) + hnew(i,j,2)
C If hIce < hIceMin, melt the ice.
IF ( hIce(i,j).GT.0. _d 0 .AND. hIce(i,j).LT.hIceMin ) THEN
frw2oc(i,j) = frw2oc(i,j)
& + (rhos*hSnow1(i,j) + rhoi*hIce(i,j))/dt
esurp(i,j) = esurp(i,j) - rhos*qsnow*hSnow1(i,j)
& - rhoi*qicen(i,j,1)*hnew(i,j,1)
& - rhoi*qicen(i,j,2)*hnew(i,j,2)
hIce(i,j) = 0. _d 0
hSnow1(i,j) = 0. _d 0
tSrf(i,j) = 0. _d 0
icFrac(i,j) = 0. _d 0
qicen(i,j,1) = 0. _d 0
qicen(i,j,2) = 0. _d 0
#ifdef ALLOW_DBUG_THSICE
IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
& 'ThSI_CALC_TH: -2 : esurp,frw2oc=',
& esurp(i,j), frw2oc(i,j)
#endif
ENDIF
C-- else hIce > 0: end
ENDIF
C endif iceMask > 0
ENDIF
C end i/j-loops
ENDDO
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
cphCADJ STORE icemask(:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE hice(:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE hnew(:,:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE hSnow1(:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE qicen(:,:,:) = comlev1_bibj, key=ticekey, byte=isbyte
#endif
DO j = jMin, jMax
DO i = iMin, iMax
IF (iceMask(i,j).GT.0. _d 0) THEN
IF ( hIce(i,j) .GT. 0. _d 0 ) THEN
C If there is enough snow to lower the ice/snow interface below
C freeboard, convert enough snow to ice to bring the interface back
C to sea-level OR if snow height is larger than hsMax, snow is
C converted to ice to bring hSnow1 down to hsMax. Largest change is
C applied and enthalpy of top ice layer adjusted accordingly.
#ifdef ALLOW_AUTODIFF_TAMC
ikey_1 = i
& + sNx*(j-1)
& + sNx*sNy*act1
& + sNx*sNy*max1*act2
& + sNx*sNy*max1*max2*act3
& + sNx*sNy*max1*max2*max3*act4
#endif /* ALLOW_AUTODIFF_TAMC */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hIce(i,j) = comlev1_thsice_1, key=ikey_1
CADJ STORE hSnow1(i,j) = comlev1_thsice_1, key=ikey_1
CADJ STORE hnew(i,j,:) = comlev1_thsice_1, key=ikey_1
CADJ STORE qicen(i,j,:) = comlev1_thsice_1, key=ikey_1
#endif
IF ( hSnow1(i,j) .GT. hIce(i,j)*floodFac
& .OR. hSnow1(i,j) .GT. hsMax ) THEN
cBB WRITE (6,*) 'Freeboard adjusts'
c dhi = (hSnow1(i,j) * rhos - hIce(i,j) * rhoiw) / rhosw
c dhs = dhi * rhoi / rhos
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hnew(i,j,:) = comlev1_thsice_1, key=ikey_1
#endif
dhs = (hSnow1(i,j) - hIce(i,j)*floodFac) * rhoi / rhosw
dhs = MAX( hSnow1(i,j) - hsMax, dhs )
dhi = dhs * rhos / rhoi
rqh = rhoi*qicen(i,j,1)*hnew(i,j,1) + rhos*qsnow*dhs
hnew(i,j,1) = hnew(i,j,1) + dhi
qicen(i,j,1) = rqh / (rhoi*hnew(i,j,1))
hIce(i,j) = hIce(i,j) + dhi
hSnow1(i,j) = hSnow1(i,j) - dhs
ENDIF
C limit ice height
C- NOTE: this part does not conserve Energy ;
C but surplus of fresh water and salt are taken into account.
IF (hIce(i,j).GT.hiMax) THEN
cBB print*,'BBerr, hIce>hiMax',i,j,hIce(i,j)
chi=hIce(i,j)-hiMax
hnew(i,j,1)=hnew(i,j,1)-chi/2. _d 0
hnew(i,j,2)=hnew(i,j,2)-chi/2. _d 0
frw2oc(i,j) = frw2oc(i,j) + chi*rhoi/dt
ENDIF
c IF (hSnow1(i,j).GT.hsMax) THEN
cc print*,'BBerr, hSnow1>hsMax',i,j,hSnow1(i,j)
c chs=hSnow1(i,j)-hsMax
c hSnow1(i,j)=hsMax
c frw2oc(i,j) = frw2oc(i,j) + chs*rhos/dt
c ENDIF
C Compute new total ice thickness.
hIce(i,j) = hnew(i,j,1) + hnew(i,j,2)
#ifdef ALLOW_DBUG_THSICE
IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
& 'ThSI_CALC_TH: b-Winton: hnew, qice =',
& hnew(i,j,1), hnew(i,j,2),
& qicen(i,j,1), qicen(i,j,2)
#endif
hlyr = hIce(i,j) * rec_nlyr
CML CALL THSICE_RESHAPE_LAYERS(
CML U qicen(i,j,:),
CML I hlyr, hnew(i,j,:), myThid )
C inlined version of S/R THSICE_RESHAPE_LAYERS
C | Repartition into equal-thickness layers, conserving energy.
C *==========================================================*
C | This is the 2-layer version (formerly "NEW_LAYERS_WINTON")
C | from M. Winton 1999, JAOT, sea-ice model.
if (hnew(i,j,1).gt.hnew(i,j,2)) then
C-- Layer 1 gives ice to layer 2
f1 = (hnew(i,j,1)-hlyr)/hlyr
q2tmp = f1*qicen(i,j,1) + (1. _d 0-f1)*qicen(i,j,2)
if (q2tmp.gt.Lfresh) then
qicen(i,j,2) = q2tmp
else
C- Keep q2 fixed to avoid q20
qh2 = hlyr*qicen(i,j,2)
qhtot = hnew(i,j,1)*qicen(i,j,1) + hnew(i,j,2)*qicen(i,j,2)
qh1 = qhtot - qh2
qicen(i,j,1) = qh1/hlyr
endif
else
C- Layer 2 gives ice to layer 1
f1 = hnew(i,j,1)/hlyr
qicen(i,j,1) = f1*qicen(i,j,1) + (1. _d 0-f1)*qicen(i,j,2)
endif
C end of inlined S/R THSICE_RESHAPE_LAYERS
#ifdef ALLOW_DBUG_THSICE
IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
& 'ThSI_CALC_TH: icFrac,hIce, qtot, hSnow1 =',
& icFrac(i,j),hIce(i,j), (qicen(i,j,1)+qicen(i,j,2))*0.5,
& hSnow1(i,j)
#endif
C- if hIce > 0 : end
ENDIF
c200 CONTINUE
C endif iceMask > 0
ENDIF
C end i/j-loops
ENDDO
ENDDO
DO j = jMin, jMax
DO i = iMin, iMax
IF (iceMask(i,j).GT.0. _d 0) THEN
C- Compute surplus energy left over from melting.
IF (hIce(i,j).LE.0. _d 0) icFrac(i,j)=0. _d 0
C.. heat fluxes left over for ocean
flx2oc(i,j) = flx2oc(i,j)
& + (Fbot(i,j)+(esurp(i,j)+etop(i,j)+ebot(i,j))/dt)
#ifdef ALLOW_DBUG_THSICE
IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
& 'ThSI_CALC_TH: [esurp,etop+ebot]/dt =',
& esurp(i,j)/dt,etop(i,j)/dt,ebot(i,j)/dt
#endif
C-- Evaporation left to the ocean :
frw2oc(i,j) = frw2oc(i,j) - evapLoc(i,j)
C-- Correct Atmos. fluxes for this different latent heat:
C evap was computed over freezing surf.(tSrf<0), latent heat = Lvap+Lfresh
C but should be Lvap only for the fraction "evap" that is left to the ocean.
flx2oc(i,j) = flx2oc(i,j) + evapLoc(i,j)*Lfresh
C fresh and salt fluxes
c frw2oc(i,j) = (mwater0(i,j) - (rhos*(hSnow1(i,j))
c & + rhoi*(hIce(i,j))))/dt-evapLoc(i,j)
c fsalt = (msalt0(i,j) - rhoi*hIce(i,j)*saltIce)/35. _d 0/dt ! for same units as frw2oc
C note (jmc): frw2oc is computed from a sea-ice mass budget that already
C contains, at this point, snow & evaporation (of snow & ice)
C but are not meant to be part of ice/ocean fresh-water flux.
C fix: a) like below or b) by making the budget before snow/evap is added
c frw2oc(i,j) = (mwater0(i,j) - (rhos*(hSnow1(i,j)) + rhoi*(hIce(i,j))))/dt
c & + snow(i,j,bi,bj)*rhos - frwAtm
fsalt(i,j) = (msalt0(i,j) - rhoi*hIce(i,j)*saltIce)/dt
#ifdef ALLOW_DBUG_THSICE
IF (dBug(i,j,bi,bj) ) THEN
WRITE(6,1020)'ThSI_CALC_TH:dH2O,Ev[kg],frw2oc,fsalt',
& (mwater0(i,j)-(rhos*hSnow1(i,j)+rhoi*hIce(i,j)))/dt,
& evapLoc(i,j),frw2oc(i,j),fsalt(i,j)
WRITE(6,1020)'ThSI_CALC_TH: flx2oc,Fbot,extend/dt =',
& flx2oc(I,J),Fbot(i,j),(etope(i,j)+ebote(i,j))/dt
ENDIF
#endif
C-- add remaining liquid Precip (rain+RunOff) directly to ocean:
frw2oc(i,j) = frw2oc(i,j) + (prcAtm(i,j)-snowP(i,j))
C endif iceMask > 0
ENDIF
C end i/j-loops
ENDDO
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
cphCADJ STORE icemask(:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE icfrac(:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE hnew(:,:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE hSnow1(:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE qicen(:,:,:) = comlev1_bibj, key=ticekey, byte=isbyte
#endif
DO j = jMin, jMax
DO i = iMin, iMax
IF (iceMask(i,j).GT.0. _d 0) THEN
C-- note: at this point, icFrac has not been changed (unless reset to zero)
C and it can only be reduced by lateral melting in the following part:
C calculate extent changes
extend=etope(i,j)+ebote(i,j)
IF (icFrac(i,j).GT.0. _d 0.AND.extend.GT.0. _d 0) THEN
rq = rhoi * 0.5 _d 0*(qicen(i,j,1)+qicen(i,j,2))
rs = rhos * qsnow
rqh = rq * hIce(i,j) + rs * hSnow1(i,j)
freshe=(rhos*hSnow1(i,j)+rhoi*hIce(i,j))/dt
salte=(rhoi*hIce(i,j)*saltIce)/dt
IF ( extend.LT.rqh ) THEN
icFrac(i,j)=(1. _d 0-extend/rqh)*icFrac(i,j)
ENDIF
IF ( extend.LT.rqh .AND. icFrac(i,j).GE.iceMaskMin ) THEN
frw2oc(i,j)=frw2oc(i,j)+extend/rqh*freshe
fsalt(i,j)=fsalt(i,j)+extend/rqh*salte
ELSE
icFrac(i,j)=0. _d 0
hIce(i,j) =0. _d 0
hSnow1(i,j) =0. _d 0
flx2oc(i,j)=flx2oc(i,j)+(extend-rqh)/dt
frw2oc(i,j)=frw2oc(i,j)+freshe
fsalt(i,j)=fsalt(i,j)+salte
ENDIF
ELSEIF (extend.GT.0. _d 0) THEN
flx2oc(i,j)=flx2oc(i,j)+extend/dt
ENDIF
C endif iceMask > 0
ENDIF
C end i/j-loops
ENDDO
ENDDO
DO j = jMin, jMax
DO i = iMin, iMax
IF (iceMask(i,j).GT.0. _d 0) THEN
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-- Update output variables :
C-- Diagnostic of Atmos. fresh water flux (E-P) over sea-ice :
C substract precip from Evap (<- stored in frwAtm array)
frwAtm(i,j) = frwAtm(i,j) - prcAtm(i,j)
C-- update Mixed-Layer Freezing potential heat flux by substracting the
C part which has already been accounted for (Fbot):
fzMlOc(i,j) = fzMlOc(i,j) - Fbot(i,j)*iceMask(i,j)
C-- Update Sea-Ice state output:
qIc1(i,j) = qicen(i,j,1)
qIc2(i,j) = qicen(i,j,2)
#ifdef ALLOW_DBUG_THSICE
IF (dBug(i,j,bi,bj) ) WRITE(6,1020)
& 'ThSI_CALC_TH: icFrac,flx2oc,fsalt,frw2oc=',
& icFrac(i,j), flx2oc(i,j), fsalt(i,j), frw2oc(i,j)
#endif
C endif iceMask > 0
ENDIF
ENDDO
ENDDO
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef CHECK_ENERGY_CONSERV
DO j = jMin, jMax
DO i = iMin, iMax
IF (iceMask(i,j).GT.0. _d 0) THEN
qaux(1)=qIc1(i,j)
qaux(2)=qIc2(i,j)
CALL THSICE_CHECK_CONSERV( dBugFlag, i, j, bi, bj, 0,
I iceMask(i,j), icFrac(i,j), hIce(i,j), hSnow1(i,j),
I qaux,
I flx2oc(i,j), frw2oc(i,j), fsalt,
I myTime, myIter, myThid )
C endif iceMask > 0
ENDIF
C end i/j-loops
ENDDO
ENDDO
#endif /* CHECK_ENERGY_CONSERV */
#endif /* ALLOW_THSICE */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
RETURN
END