pkg/shelfice/shelfice_thermodynamics.F

ffe464dc7d Mart*0001 #include "SHELFICE_OPTIONS.h"
26cd020e03 Jean*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
                0005 #ifdef ALLOW_CTRL
                0006 # include "CTRL_OPTIONS.h"
                0007 #endif
ead4e7560e Jean*0008
ffe464dc7d Mart*0009 CBOP
                0010 C     !ROUTINE: SHELFICE_THERMODYNAMICS
                0011 C     !INTERFACE:
ead4e7560e Jean*0012       SUBROUTINE SHELFICE_THERMODYNAMICS(
ffe464dc7d Mart*0013      I                        myTime, myIter, myThid )
                0014 C     !DESCRIPTION: \bv
                0015 C     *=============================================================*
ead4e7560e Jean*0016 C     | S/R  SHELFICE_THERMODYNAMICS
                0017 C     | o shelf-ice main routine.
                0018 C     |   compute temperature and (virtual) salt flux at the
ffe464dc7d Mart*0019 C     |   shelf-ice ocean interface
                0020 C     |
                0021 C     | stresses at the ice/water interface are computed in separate
                0022 C     | routines that are called from mom_fluxform/mom_vecinv
                0023 C     *=============================================================*
5744811f23 Mart*0024 C     \ev
ffe464dc7d Mart*0025
                0026 C     !USES:
                0027       IMPLICIT NONE
                0028
                0029 C     === Global variables ===
                0030 #include "SIZE.h"
                0031 #include "EEPARAMS.h"
                0032 #include "PARAMS.h"
8eaa195dfe Jean*0033 #include "EOS.h"
ffe464dc7d Mart*0034 #include "GRID.h"
                0035 #include "DYNVARS.h"
                0036 #include "FFIELDS.h"
                0037 #include "SHELFICE.h"
6b47d550f4 Mart*0038 #ifdef ALLOW_CTRL
83a7138cb5 Jean*0039 # include "CTRL_SIZE.h"
4d72283393 Mart*0040 # include "CTRL.h"
edcd27be69 Mart*0041 # include "CTRL_DUMMY.h"
6b47d550f4 Mart*0042 # ifdef ALLOW_GENTIM2D_CONTROL
                0043 #  include "CTRL_GENARR.h"
                0044 # endif
                0045 #endif /* ALLOW_CTRL */
7ce1881910 Mart*0046 #ifdef ALLOW_AUTODIFF_TAMC
c96d63ff5a Jean*0047 # include "tamc.h"
7ce1881910 Mart*0048 #endif /* ALLOW_AUTODIFF_TAMC */
ead4e7560e Jean*0049
ffe464dc7d Mart*0050 C     !INPUT/OUTPUT PARAMETERS:
                0051 C     === Routine arguments ===
                0052 C     myIter :: iteration counter for this thread
                0053 C     myTime :: time counter for this thread
                0054 C     myThid :: thread number for this instance of the routine.
                0055       _RL  myTime
                0056       INTEGER myIter
                0057       INTEGER myThid
                0058
                0059 #ifdef ALLOW_SHELFICE
40b188dddd Mart*0060 C     !LOCAL VARIABLES :
ffe464dc7d Mart*0061 C     === Local variables ===
f39e9f371c Mart*0062 C     i,j,k,kp1,bi,bj  :: loop counters
71125c711d Mart*0063 C     tLoc, sLoc, pLoc :: local in-situ temperature, salinity, pressure
ead4e7560e Jean*0064 C     theta/saltFreeze :: temperature and salinity of water at the
71125c711d Mart*0065 C                         ice-ocean interface (at the freezing point)
e4305b0f18 Patr*0066 C     freshWaterFlux   :: local variable for fresh water melt flux due
823dd639fd Jean*0067 C                         to melting in kg/m^2/s
e4305b0f18 Patr*0068 C                         (negative density x melt rate)
7108c05fdf Mart*0069 C     convertFW2SaltLoc:: local copy of convertFW2Salt
a8eec31802 Mart*0070 C     cFac             :: 1 for conservative form, 0, otherwise
823dd639fd Jean*0071 C     rFac             :: realFreshWaterFlux factor
5744811f23 Mart*0072 C     dFac             :: 0 for diffusive heat flux (Holland and Jenkins, 1999,
                0073 C                           eq21)
37e32ac0e0 Mart*0074 C                         1 for advective and diffusive heat flux (eq22, 26, 31)
5744811f23 Mart*0075 C     fwflxFac         :: only effective for dFac=1, 1 if we expect a melting
                0076 C                         fresh water flux, 0 otherwise
9952f046d7 dngo*0077 C     rFWinBL          :: = 1 when realFreshWaterFlux is used with BoundaryLayer
71125c711d Mart*0078 C     auxiliary variables and abbreviations:
f39e9f371c Mart*0079 C     a0, a1, a2, b0, c0
                0080 C     eps1, eps2, eps3, eps3a, eps4, eps6, eps7, eps8, epsq
7108c05fdf Mart*0081 C     aqe, bqe, cqe, discrim, recip_aqe
71125c711d Mart*0082 C     drKp1, recip_drLoc
f39e9f371c Mart*0083       INTEGER i,j,k,kp1
ffe464dc7d Mart*0084       INTEGER bi,bj
04ba99fb18 Mart*0085       _RL tLoc(1:sNx,1:sNy)
                0086       _RL sLoc(1:sNx,1:sNy)
                0087       _RL pLoc(1:sNx,1:sNy)
9952f046d7 dngo*0088       _RL uLoc(1:sNx+1,1:sNy+1)
                0089       _RL vLoc(1:sNx+1,1:sNy+1)
43bfccd244 Jean*0090       _RL velSq(1:sNx,1:sNy)
7b8b86ab99 Timo*0091       _RL tTMP, sTMP, uTMP, vTMP, vSqTmp
52d51b46b8 Jean*0092       _RL thetaFreeze, saltFreeze, recip_Cp
9952f046d7 dngo*0093       _RL freshWaterFlux
                0094 #ifdef ALLOW_ISOMIP_TD
                0095       _RL convertFW2SaltLoc
                0096 #endif
f39e9f371c Mart*0097       _RL a0, a1, a2, b0, c0
                0098       _RL eps1, eps2, eps3, eps3a, eps4, eps6, eps7, eps8, epsq
9952f046d7 dngo*0099       _RL cFac, rFac, dFac, fwflxFac, rFWinBL
71125c711d Mart*0100       _RL aqe, bqe, cqe, discrim, recip_aqe
40b188dddd Mart*0101       _RL drKp1, recip_drLoc
6c8ffc8ddb Mart*0102       _RL recip_latentHeat
a86b6152f4 Dimi*0103       _RL tmpFac
ffe464dc7d Mart*0104
e4305b0f18 Patr*0105 #ifdef SHI_ALLOW_GAMMAFRICT
784f4a51b3 Mart*0106       _RL shiPr, shiSc, shiLo, recip_shiKarman, shiTwoThirds
                0107       _RL gammaTmoleT, gammaTmoleS, gammaTurb, gammaTurbConst
e4305b0f18 Patr*0108       _RL ustar, ustarSq, etastar
9952f046d7 dngo*0109       _RL u_tmp, v_tmp
784f4a51b3 Mart*0110       PARAMETER ( shiTwoThirds = 0.66666666666666666666666666667D0 )
5744811f23 Mart*0111 #ifdef ALLOW_DIAGNOSTICS
                0112       _RL uStarDiag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0113 #endif /* ALLOW_DIAGNOSTICS */
9952f046d7 dngo*0114 #endif /*SHI_ALLOW_GAMMAFRICT */
e4305b0f18 Patr*0115
5357e1c86c Jean*0116 #ifndef ALLOW_OPENAD
e3465203db Mart*0117       _RL SW_TEMP
                0118       EXTERNAL SW_TEMP
bb49380557 Patr*0119 #endif
e3465203db Mart*0120
6b47d550f4 Mart*0121 #ifdef ALLOW_GENTIM2D_CONTROL
                0122       INTEGER iarr
9952f046d7 dngo*0123       _RL xx_shifwflx_loc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
9e75f57067 Jean*0124 #endif
dcd3df3e6a Ou W*0125 #ifdef ALLOW_AUTODIFF_TAMC
                0126 # ifdef SHI_ALLOW_GAMMAFRICT
                0127       INTEGER ikey
                0128 # endif
                0129 #endif
5744811f23 Mart*0130 CEOP
ffe464dc7d Mart*0131 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0132
e4305b0f18 Patr*0133 #ifdef SHI_ALLOW_GAMMAFRICT
26cd020e03 Jean*0134 #ifdef ALLOW_AUTODIFF
7ce1881910 Mart*0135 C     re-initialize here again, curtesy to TAF
                0136       DO bj = myByLo(myThid), myByHi(myThid)
                0137        DO bi = myBxLo(myThid), myBxHi(myThid)
00f81e6785 Ou W*0138         DO j = 1-OLy,sNy+OLy
                0139          DO i = 1-OLx,sNx+OLx
7ce1881910 Mart*0140           shiTransCoeffT(i,j,bi,bj) = SHELFICEheatTransCoeff
                0141           shiTransCoeffS(i,j,bi,bj) = SHELFICEsaltTransCoeff
                0142          ENDDO
                0143         ENDDO
                0144        ENDDO
823dd639fd Jean*0145       ENDDO
26cd020e03 Jean*0146 #endif /* ALLOW_AUTODIFF */
784f4a51b3 Mart*0147       IF ( SHELFICEuseGammaFrict ) THEN
                0148 C     Implement friction velocity-dependent transfer coefficient
                0149 C     of Holland and Jenkins, JPO, 1999
                0150        recip_shiKarman= 1. _d 0 / 0.4 _d 0
                0151        shiLo = 0. _d 0
                0152        shiPr = shiPrandtl**shiTwoThirds
                0153        shiSc = shiSchmidt**shiTwoThirds
                0154 cph      shiPr = (viscArNr(1)/diffKrNrT(1))**shiTwoThirds
                0155 cph      shiSc = (viscArNr(1)/diffKrNrS(1))**shiTwoThirds
                0156        gammaTmoleT = 12.5 _d 0 * shiPr - 6. _d 0
                0157        gammaTmoleS = 12.5 _d 0 * shiSc - 6. _d 0
7462d20066 Mart*0158 C     instead of etastar = sqrt(1+zetaN*ustar./(f*Lo*Rc))
                0159        etastar = 1. _d 0
823dd639fd Jean*0160        gammaTurbConst  = 1. _d 0 / (2. _d 0 * shiZetaN*etastar)
784f4a51b3 Mart*0161      &      - recip_shiKarman
                0162       ENDIF
                0163 #endif /* SHI_ALLOW_GAMMAFRICT */
e4305b0f18 Patr*0164
6c8ffc8ddb Mart*0165       recip_latentHeat = 0. _d 0
                0166       IF ( SHELFICElatentHeat .NE. 0. _d 0 )
                0167      &     recip_latentHeat = 1. _d 0/SHELFICElatentHeat
7108c05fdf Mart*0168 C     are we doing the conservative form of Jenkins et al. (2001)?
52d51b46b8 Jean*0169       recip_Cp = 1. _d 0 / HeatCapacity_Cp
7108c05fdf Mart*0170       cFac = 0. _d 0
                0171       IF ( SHELFICEconserve ) cFac = 1. _d 0
9952f046d7 dngo*0172 C     with "real fresh water flux" (affecting ETAN), there is more to modify
7108c05fdf Mart*0173       rFac = 1. _d 0
                0174       IF ( SHELFICEconserve .AND. useRealFreshWaterFlux ) rFac = 0. _d 0
9952f046d7 dngo*0175       rFWinBL = 0. _d 0
                0176       IF ( SHI_withBL_realFWflux ) rFWinBL = 1. _d 0
823dd639fd Jean*0177 C     heat flux into the ice shelf, default is diffusive flux
37e32ac0e0 Mart*0178 C     (Holland and Jenkins, 1999, eq.21)
                0179       dFac = 0. _d 0
                0180       IF ( SHELFICEadvDiffHeatFlux ) dFac = 1. _d 0
                0181       fwflxFac = 0. _d 0
17292dde13 Mart*0182 C     linear dependence of freezing point on salinity
ffe464dc7d Mart*0183       a0 = -0.0575   _d  0
17292dde13 Mart*0184       a1 =  0.0      _d -0
                0185       a2 =  0.0      _d -0
                0186       c0 =  0.0901   _d  0
f39e9f371c Mart*0187       b0 =  -7.61    _d -4
17292dde13 Mart*0188 #ifdef ALLOW_ISOMIP_TD
                0189       IF ( useISOMIPTD ) THEN
                0190 C     non-linear dependence of freezing point on salinity
                0191        a0 = -0.0575   _d  0
                0192        a1 = 1.710523  _d -3
                0193        a2 = -2.154996 _d -4
f39e9f371c Mart*0194        b0 = -7.53     _d -4
17292dde13 Mart*0195        c0 = 0. _d 0
                0196       ENDIF
7108c05fdf Mart*0197       convertFW2SaltLoc = convertFW2Salt
                0198 C     hardcoding this value here is OK because it only applies to ISOMIP
                0199 C     where this value is part of the protocol
                0200       IF ( convertFW2SaltLoc .EQ. -1. ) convertFW2SaltLoc = 33.4 _d 0
f34b2e0f36 Jean*0201 #endif /* ALLOW_ISOMIP_TD */
a8eec31802 Mart*0202
ffe464dc7d Mart*0203       DO bj = myByLo(myThid), myByHi(myThid)
                0204        DO bi = myBxLo(myThid), myBxHi(myThid)
f39e9f371c Mart*0205         DO j = 1-OLy,sNy+OLy
                0206          DO i = 1-OLx,sNx+OLx
                0207           shelfIceHeatFlux      (i,j,bi,bj) = 0. _d 0
                0208           shelfIceFreshWaterFlux(i,j,bi,bj) = 0. _d 0
                0209           shelficeForcingT      (i,j,bi,bj) = 0. _d 0
                0210           shelficeForcingS      (i,j,bi,bj) = 0. _d 0
74a4ab463e Mart*0211 #if (defined SHI_ALLOW_GAMMAFRICT && defined ALLOW_DIAGNOSTICS)
f39e9f371c Mart*0212           uStarDiag             (i,j,bi,bj) = 0. _d 0
74a4ab463e Mart*0213 #endif /* SHI_ALLOW_GAMMAFRICT and ALLOW_DIAGNOSTICS */
3bafcf6020 Timo*0214 #ifdef ALLOW_GENTIM2D_CONTROL
f39e9f371c Mart*0215           xx_shifwflx_loc       (i,j,bi,bj) = 0. _d 0
3bafcf6020 Timo*0216 #endif /* ALLOW_GENTIM2D_CONTROL */
40b188dddd Mart*0217          ENDDO
                0218         ENDDO
f47cd69e8b Mart*0219        ENDDO
                0220       ENDDO
cf705a6c8e Mart*0221 #if (defined ALLOW_CTRL && defined ALLOW_GENTIM2D_CONTROL)
                0222       IF ( useCTRL ) THEN
6b47d550f4 Mart*0223        DO iarr = 1, maxCtrlTim2D
                0224         IF (xx_gentim2d_file(iarr)(1:11).EQ.'xx_shifwflx') THEN
                0225          DO bj = myByLo(myThid),myByHi(myThid)
                0226           DO bi = myBxLo(myThid),myBxHi(myThid)
f39e9f371c Mart*0227            DO j = 1,sNy
                0228             DO i = 1,sNx
                0229              xx_shifwflx_loc(i,j,bi,bj)=xx_gentim2d(i,j,bi,bj,iarr)
6b47d550f4 Mart*0230             ENDDO
                0231            ENDDO
                0232           ENDDO
f47cd69e8b Mart*0233          ENDDO
6b47d550f4 Mart*0234         ENDIF
f47cd69e8b Mart*0235        ENDDO
6b47d550f4 Mart*0236       ENDIF
cf705a6c8e Mart*0237 #endif
c96d63ff5a Jean*0238
005af54e38 Jean*0239 #ifdef ALLOW_SHELFICE_REMESHING
                0240       IF ( SHI_update_kTopC ) THEN
c96d63ff5a Jean*0241 C--   Deal with ice-shelf edge advance or retreat when allowing ice-shelf
                0242 C     mass to change. Since current ice-shelf representation is missing
                0243 C     a fractional ice-shelf cover representation, this turns into making
                0244 C     or removing a "thin ice-shelf" grid-cell. Here we update "kTopC"
                0245 C     as ice-shelf mass appears or vanishes in any surface grid-cell.
                0246         DO bj = myByLo(myThid), myByHi(myThid)
                0247          DO bi = myBxLo(myThid), myBxHi(myThid)
                0248           DO j = 1-OLy,sNy+OLy
                0249            DO i = 1-OLx,sNx+OLx
                0250              IF ( kSurfC(i,j,bi,bj).LE.Nr .AND.
                0251      &            shelficeMass(i,j,bi,bj).GT.zeroRL ) THEN
                0252                kTopC(i,j,bi,bj) = kSurfC(i,j,bi,bj)
                0253              ELSE
                0254                kTopC(i,j,bi,bj) = 0
                0255              ENDIF
                0256            ENDDO
                0257           ENDDO
                0258          ENDDO
                0259         ENDDO
                0260       ENDIF
005af54e38 Jean*0261 # ifdef ALLOW_AUTODIFF_TAMC
c96d63ff5a Jean*0262 CADJ STORE kTopC        = comlev1, key=ikey_dynamics
005af54e38 Jean*0263 # endif
                0264 #endif /* ALLOW_SHELFICE_REMESHING */
                0265 #ifdef ALLOW_AUTODIFF_TAMC
c96d63ff5a Jean*0266 CADJ STORE shelficeMass = comlev1, key=ikey_dynamics, kind=isbyte
                0267 #endif
                0268
f47cd69e8b Mart*0269       DO bj = myByLo(myThid), myByHi(myThid)
                0270        DO bi = myBxLo(myThid), myBxHi(myThid)
9952f046d7 dngo*0271
7ce1881910 Mart*0272 #ifdef ALLOW_AUTODIFF_TAMC
                0273 # ifdef SHI_ALLOW_GAMMAFRICT
20dee61641 Mart*0274         ikey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
7ce1881910 Mart*0275 # endif /* SHI_ALLOW_GAMMAFRICT */
                0276 #endif /* ALLOW_AUTODIFF_TAMC */
c96d63ff5a Jean*0277
823dd639fd Jean*0278 C--   make local copies of temperature, salinity and depth (pressure in deci-bar)
04ba99fb18 Mart*0279 C--   underneath the ice
f39e9f371c Mart*0280         DO j = 1, sNy
                0281          DO i = 1, sNx
                0282           k         = MAX(1,kTopC(i,j,bi,bj))
8eaa195dfe Jean*0283 C-    original (inaccurate) pLoc expression below:
                0284 C     assumes rhoConst*gravity*SItodBar=1 + missing both "shelficeLoadAnomaly"
                0285 C      static contribution and etaN dynamical contribution:
f39e9f371c Mart*0286 c         pLoc(i,j) = ABS(R_shelfIce(i,j,bi,bj))
8eaa195dfe Jean*0287 C-    new and accurate pLoc expression (just the weight of the ice):
f39e9f371c Mart*0288           pLoc(i,j) = shelficeMass(i,j,bi,bj)*gravity*SItodBar
                0289           tLoc(i,j) = theta(i,j,k,bi,bj)
                0290           sLoc(i,j) = MAX(salt(i,j,k,bi,bj), zeroRL)
                0291           velSq(i,j)= 0.
b06dffee6b Jean*0292          ENDDO
                0293         ENDDO
f39e9f371c Mart*0294         DO j = 1, sNy+1
                0295          DO i = 1, sNx+1
                0296           uLoc(i,j) = 0.
                0297           vLoc(i,j) = 0.
9952f046d7 dngo*0298          ENDDO
                0299         ENDDO
                0300 #ifdef SHI_ALLOW_GAMMAFRICT
                0301         IF ( SHELFICEuseGammaFrict .AND. SHELFICE_oldCalcUStar ) THEN
                0302 C-    Original averaging expression for uStar:
f39e9f371c Mart*0303          DO j = 1, sNy
                0304           DO i = 1, sNx
                0305            k = MAX(1,kTopC(i,j,bi,bj))
                0306            uLoc(i,j) = recip_hFacC(i,j,k,bi,bj) * halfRL *
                0307      &         ( uVel(i,  j,k,bi,bj) * _hFacW(i,  j,k,bi,bj)
                0308      &         + uVel(i+1,j,k,bi,bj) * _hFacW(i+1,j,k,bi,bj) )
                0309            vLoc(i,j) = recip_hFacC(i,j,k,bi,bj) * halfRL *
                0310      &         ( vVel(i,j,  k,bi,bj) * _hFacS(i,j,  k,bi,bj)
                0311      &         + vVel(i,j+1,k,bi,bj) * _hFacS(i,j+1,k,bi,bj) )
                0312            velSq(i,j) = uLoc(i,j)*uLoc(i,j)+vLoc(i,j)*vLoc(i,j)
9952f046d7 dngo*0313           ENDDO
                0314          ENDDO
                0315         ELSEIF ( SHELFICEuseGammaFrict ) THEN
b06dffee6b Jean*0316 C-    New (more accurate) averaging expression for uStar:
f39e9f371c Mart*0317          DO j = 1, sNy
                0318           DO i = 1, sNx
                0319            uLoc(i,j) = 0.
                0320            vLoc(i,j) = 0.
                0321            velSq(i,j) = 0.
                0322            k = MAX(1,kTopC(i,j,bi,bj))
                0323            tmpFac = _hFacW(i,  j,k,bi,bj) + _hFacW(i+1,j,k,bi,bj)
b06dffee6b Jean*0324            IF ( tmpFac.GT.0. _d 0 )
f39e9f371c Mart*0325      &     velSq(i,j) = (
                0326      &     uVel( i, j,k,bi,bj)*uVel( i, j,k,bi,bj)*_hFacW( i, j,k,bi,bj)
                0327      &   + uVel(i+1,j,k,bi,bj)*uVel(i+1,j,k,bi,bj)*_hFacW(i+1,j,k,bi,bj)
43bfccd244 Jean*0328      &                  )/tmpFac
f39e9f371c Mart*0329            tmpFac = _hFacS(i,j,  k,bi,bj) + _hFacS(i,j+1,k,bi,bj)
b06dffee6b Jean*0330            IF ( tmpFac.GT.0. _d 0 )
f39e9f371c Mart*0331      &     velSq(i,j) = velSq(i,j) + (
                0332      &     vVel(i, j, k,bi,bj)*vVel(i, j, k,bi,bj)*_hFacS(i, j, k,bi,bj)
                0333      &   + vVel(i,j+1,k,bi,bj)*vVel(i,j+1,k,bi,bj)*_hFacS(i,j+1,k,bi,bj)
43bfccd244 Jean*0334      &                               )/tmpFac
b06dffee6b Jean*0335           ENDDO
                0336          ENDDO
                0337         ENDIF
9952f046d7 dngo*0338 #endif /* SHI_ALLOW_GAMMAFRICT */
                0339
04ba99fb18 Mart*0340         IF ( SHELFICEBoundaryLayer ) THEN
                0341 C--   average over boundary layer width
f39e9f371c Mart*0342          DO j = 1, sNy
                0343           DO i = 1, sNx
                0344            k   = kTopC(i,j,bi,bj)
                0345            IF ( k .NE. 0 .AND. k .LT. Nr ) THEN
7b8b86ab99 Timo*0346 C           This local variable assignment/reassignment is simply to
                0347 C           avoid TAF recomputations... Oh the games we play...
f39e9f371c Mart*0348             tTMP   = tLoc (i,j)
                0349             sTMP   = sLoc (i,j)
                0350             uTMP   = uLoc (i,j)
                0351             vTMP   = vLoc (i,j)
                0352             vSqTMP = velSq(i,j)
7b8b86ab99 Timo*0353
f39e9f371c Mart*0354             kp1 = MIN(Nr,k+1)
40b188dddd Mart*0355 C--   overlap into lower cell
f39e9f371c Mart*0356             drKp1 = drF(k)*( 1. _d 0 - _hFacC(i,j,k,bi,bj) )
40b188dddd Mart*0357 C--   lower cell may not be as thick as required
f39e9f371c Mart*0358             drKp1 = MIN( drKp1, drF(kp1) * _hFacC(i,j,kp1,bi,bj) )
09df6ae39f Dani*0359             drKp1 = MAX( drKp1, 0. _d 0 )
ead4e7560e Jean*0360             recip_drLoc = 1. _d 0 /
f39e9f371c Mart*0361      &           ( drF(k)*_hFacC(i,j,k,bi,bj) + drKp1 )
                0362             tTMP = ( tTMP * drF(k)*_hFacC(i,j,k,bi,bj)
                0363      &           + theta(i,j,kp1,bi,bj) *drKp1 )
04ba99fb18 Mart*0364      &           * recip_drLoc
f39e9f371c Mart*0365             sTMP = ( sTMP * drF(k)*_hFacC(i,j,k,bi,bj)
                0366      &           + MAX(salt(i,j,kp1,bi,bj), zeroRL) * drKp1 )
04ba99fb18 Mart*0367      &           * recip_drLoc
f39e9f371c Mart*0368             uTMP = ( uTMP * drF(k)*_hFacC(i,j,k,bi,bj)
                0369      &           + drKp1 * recip_hFacC(i,j,kp1,bi,bj) * halfRL *
                0370      &           ( uVel(i,  j,kp1,bi,bj) * _hFacW(i,  j,kp1,bi,bj)
                0371      &           + uVel(i+1,j,kp1,bi,bj) * _hFacW(i+1,j,kp1,bi,bj) )
7ce1881910 Mart*0372      &           ) * recip_drLoc
f39e9f371c Mart*0373             vTMP = ( vTMP * drF(k)*_hFacC(i,j,k,bi,bj)
                0374      &           + drKp1 * recip_hFacC(i,j,kp1,bi,bj) * halfRL *
                0375      &           ( vVel(i,j,  kp1,bi,bj) * _hFacS(i,j,  kp1,bi,bj)
                0376      &           + vVel(i,j+1,kp1,bi,bj) * _hFacS(i,j+1,kp1,bi,bj) )
7ce1881910 Mart*0377      &           ) * recip_drLoc
7b8b86ab99 Timo*0378             vSqTMP = uTMP*uTMP + vTMP*vTMP
                0379
f39e9f371c Mart*0380             tLoc (i,j) = tTMP
                0381             sLoc (i,j) = sTMP
                0382             uLoc (i,j) = uTMP
                0383             vLoc (i,j) = vTMP
                0384             velSq(i,j) = vSqTMP
04ba99fb18 Mart*0385            ENDIF
                0386           ENDDO
                0387          ENDDO
                0388         ENDIF
                0389
9952f046d7 dngo*0390 #ifdef SHI_ALLOW_GAMMAFRICT
                0391         IF ( SHI_withBL_uStarTopDz ) THEN
                0392 C--  TOPDR averages U/V over boundary layer at U/V points, then averages
                0393 C    (as opposed to averaging horizontally then vertically)
                0394 C    Average at u- and v- points over deltaR. use later to override uLoc/vLoc
f39e9f371c Mart*0395          DO j = 1, sNy+1
                0396           DO i = 1, sNx+1
                0397            k = kSurfW(i,j,bi,bj)
                0398            IF (k.LT.Nr) THEN
                0399             kp1 = k+1
                0400             drKp1 = drF(k)*( 1. _d 0 - _hFacW(i,j,k,bi,bj) )
                0401             drKp1 = MIN( drKp1, drF(kp1) * _hFacW(i,j,kp1,bi,bj) )
9952f046d7 dngo*0402             drKp1 = MAX( drKp1, 0. _d 0 )
                0403             recip_drLoc = 1. _d 0
f39e9f371c Mart*0404      &                  / ( drF(k)*_hFacW(i,j,k,bi,bj) + drKp1 )
                0405             uLoc(i,j) = ( drF(k)*_hFacW(i,j,k,bi,bj)*uVel(i,j,k,bi,bj)
                0406      &                  + drKp1*uVel(i,j,kp1,bi,bj)
9952f046d7 dngo*0407      &                  )*recip_drLoc
f39e9f371c Mart*0408 c           u_topDr(i,j) =
                0409 c    &           ( uVel(i,j,k,bi,bj)*drF(k)*_hFacW(i,j,k,bi,bj)
                0410 c    &           + uVel(i,j,kp1,bi,bj)*drKp1
9952f046d7 dngo*0411 c    &           )*recip_drLoc
f39e9f371c Mart*0412            ELSEIF (k.EQ.Nr) THEN
                0413             uLoc(i,j) = uVel(i,j,k,bi,bj)
9952f046d7 dngo*0414            ELSE
f39e9f371c Mart*0415             uLoc(i,j) = 0. _d 0
9952f046d7 dngo*0416            ENDIF
                0417
f39e9f371c Mart*0418            k = kSurfS(i,j,bi,bj)
                0419            IF (k.LT.Nr) THEN
                0420             kp1 = k+1
                0421             drKp1 = drF(k)*( 1. _d 0 - _hFacS(i,j,k,bi,bj) )
                0422             drKp1 = MIN( drKp1, drF(kp1) * _hFacS(i,j,kp1,bi,bj) )
9952f046d7 dngo*0423             drKp1 = MAX( drKp1, 0. _d 0 )
                0424             recip_drLoc = 1. _d 0
f39e9f371c Mart*0425      &                  / ( drF(k)*_hFacS(i,j,k,bi,bj) + drKp1 )
                0426             vLoc(i,j) = ( drF(k)*_hFacS(i,j,k,bi,bj)*vVel(i,j,k,bi,bj)
                0427      &                  + drKp1*vVel(i,j,kp1,bi,bj)
9952f046d7 dngo*0428      &                  )*recip_drLoc
f39e9f371c Mart*0429 c           v_topDr(i,j) =
                0430 c    &           ( vVel(i,j,k,bi,bj)*drF(k)*_hFacS(i,j,k,bi,bj)
                0431 c    &           + vVel(i,j,kp1,bi,bj)*drKp1
9952f046d7 dngo*0432 c    &           )*recip_drLoc
f39e9f371c Mart*0433            ELSEIF (k.EQ.Nr) THEN
                0434             vLoc(i,j) = vVel(i,j,k,bi,bj)
9952f046d7 dngo*0435            ELSE
f39e9f371c Mart*0436             vLoc(i,j) = 0. _d 0
9952f046d7 dngo*0437            ENDIF
                0438
                0439           ENDDO
                0440          ENDDO
f39e9f371c Mart*0441          DO j = 1, sNy
                0442           DO i = 1, sNx
                0443            u_tmp = halfRL*( uLoc(i,j) + uLoc(i+1,j) )
                0444            v_tmp = halfRL*( vLoc(i,j) + vLoc(i,j+1) )
                0445            velSq(i,j) = u_tmp*u_tmp + v_tmp*v_tmp
9952f046d7 dngo*0446           ENDDO
                0447          ENDDO
                0448         ENDIF
                0449 #endif /* SHI_ALLOW_GAMMAFRICT */
                0450
ead4e7560e Jean*0451 C--   turn potential temperature into in-situ temperature relative
e3465203db Mart*0452 C--   to the surface
f39e9f371c Mart*0453         DO j = 1, sNy
                0454          DO i = 1, sNx
9952f046d7 dngo*0455 #ifdef ALLOW_OPENAD
f39e9f371c Mart*0456           CALL SW_TEMP(sLoc(i,j),tLoc(i,j),pLoc(i,j),zeroRL,tLoc(i,j))
9952f046d7 dngo*0457 #else
f39e9f371c Mart*0458           tLoc(i,j) = SW_TEMP(sLoc(i,j),tLoc(i,j),pLoc(i,j),zeroRL)
bb49380557 Patr*0459 #endif
e3465203db Mart*0460          ENDDO
                0461         ENDDO
                0462
7462d20066 Mart*0463 #ifdef SHI_ALLOW_GAMMAFRICT
                0464         IF ( SHELFICEuseGammaFrict ) THEN
f39e9f371c Mart*0465          DO j = 1, sNy
                0466           DO i = 1, sNx
                0467            k = kTopC(i,j,bi,bj)
                0468            IF ( k .NE. 0 .AND. pLoc(i,j) .GT. 0. _d 0 ) THEN
                0469             ustarSq = shiCdragfld(i,j,bi,bj) * MAX( 1.D-6, velSq(i,j) )
7462d20066 Mart*0470             ustar   = SQRT(ustarSq)
5744811f23 Mart*0471 #ifdef ALLOW_DIAGNOSTICS
f39e9f371c Mart*0472             uStarDiag(i,j,bi,bj) = ustar
5744811f23 Mart*0473 #endif /* ALLOW_DIAGNOSTICS */
7462d20066 Mart*0474 C     instead of etastar = sqrt(1+zetaN*ustar./(f*Lo*Rc))
                0475 C           etastar = 1. _d 0
823dd639fd Jean*0476 C           gammaTurbConst  = 1. _d 0 / (2. _d 0 * shiZetaN*etastar)
7462d20066 Mart*0477 C    &           - recip_shiKarman
f39e9f371c Mart*0478             IF ( fCori(i,j,bi,bj) .NE. 0. _d 0 ) THEN
7462d20066 Mart*0479              gammaTurb = LOG( ustarSq * shiZetaN * etastar**2
f39e9f371c Mart*0480      &            / ABS(fCori(i,j,bi,bj) * 5.0 _d 0 * shiKinVisc))
823dd639fd Jean*0481      &            * recip_shiKarman
7462d20066 Mart*0482      &            + gammaTurbConst
823dd639fd Jean*0483 C     Do we need to catch the unlikely case of very small ustar
7462d20066 Mart*0484 C     that can lead to negative gammaTurb?
                0485 C            gammaTurb = MAX(0.D0, gammaTurb)
                0486             ELSE
                0487              gammaTurb = gammaTurbConst
                0488             ENDIF
823dd639fd Jean*0489             shiTransCoeffT(i,j,bi,bj) = MAX( zeroRL,
7462d20066 Mart*0490      &           ustar/(gammaTurb + gammaTmoleT) )
823dd639fd Jean*0491             shiTransCoeffS(i,j,bi,bj) = MAX( zeroRL,
7462d20066 Mart*0492      &           ustar/(gammaTurb + gammaTmoleS) )
                0493            ENDIF
                0494           ENDDO
                0495          ENDDO
                0496         ENDIF
                0497 #endif /* SHI_ALLOW_GAMMAFRICT */
                0498
7ce1881910 Mart*0499 #ifdef ALLOW_AUTODIFF_TAMC
                0500 # ifdef SHI_ALLOW_GAMMAFRICT
                0501 CADJ STORE shiTransCoeffS(:,:,bi,bj) = comlev1_bibj,
                0502 CADJ &     key=ikey, byte=isbyte
                0503 CADJ STORE shiTransCoeffT(:,:,bi,bj) = comlev1_bibj,
                0504 CADJ &     key=ikey, byte=isbyte
                0505 # endif /* SHI_ALLOW_GAMMAFRICT */
                0506 #endif /* ALLOW_AUTODIFF_TAMC */
04ba99fb18 Mart*0507 #ifdef ALLOW_ISOMIP_TD
                0508         IF ( useISOMIPTD ) THEN
f39e9f371c Mart*0509          DO j = 1, sNy
                0510           DO i = 1, sNx
                0511            k = kTopC(i,j,bi,bj)
                0512            IF ( k .NE. 0 .AND. pLoc(i,j) .GT. 0. _d 0 ) THEN
ffe464dc7d Mart*0513 C--   Calculate freezing temperature as a function of salinity and pressure
ead4e7560e Jean*0514             thetaFreeze =
f39e9f371c Mart*0515      &           sLoc(i,j) * ( a0 + a1*sqrt(sLoc(i,j)) + a2*sLoc(i,j) )
                0516      &           + b0*pLoc(i,j) + c0
ffe464dc7d Mart*0517 C--   Calculate the upward heat and  fresh water fluxes
f39e9f371c Mart*0518             shelfIceHeatFlux(i,j,bi,bj) = maskC(i,j,k,bi,bj)
e4305b0f18 Patr*0519      &           * shiTransCoeffT(i,j,bi,bj)
f39e9f371c Mart*0520      &           * ( tLoc(i,j) - thetaFreeze )
e4305b0f18 Patr*0521      &           * HeatCapacity_Cp*rUnit2mass
3bafcf6020 Timo*0522 #ifdef ALLOW_GENTIM2D_CONTROL
f39e9f371c Mart*0523      &           - xx_shifwflx_loc(i,j,bi,bj)*SHELFICElatentHeat
3bafcf6020 Timo*0524 #endif /*  ALLOW_GENTIM2D_CONTROL */
ffe464dc7d Mart*0525 C     upward heat flux into the shelf-ice implies basal melting,
2bf81d7d9d Mart*0526 C     thus a downward (negative upward) fresh water flux (as a mass flux),
                0527 C     and vice versa
f39e9f371c Mart*0528             shelfIceFreshWaterFlux(i,j,bi,bj) =
                0529      &           - shelfIceHeatFlux(i,j,bi,bj)
6c8ffc8ddb Mart*0530      &           *recip_latentHeat
ffe464dc7d Mart*0531 C--   compute surface tendencies
                0532             shelficeForcingT(i,j,bi,bj) =
f39e9f371c Mart*0533      &           - shelfIceHeatFlux(i,j,bi,bj)
0b1017b546 Jean*0534      &           *recip_Cp*mass2rUnit
f39e9f371c Mart*0535      &           - cFac * shelfIceFreshWaterFlux(i,j,bi,bj)*mass2rUnit
                0536      &           * ( thetaFreeze - tLoc(i,j) )
ead4e7560e Jean*0537             shelficeForcingS(i,j,bi,bj) =
f39e9f371c Mart*0538      &           shelfIceFreshWaterFlux(i,j,bi,bj) * mass2rUnit
                0539      &           * ( cFac*sLoc(i,j) + (1. _d 0-cFac)*convertFW2SaltLoc )
ffe464dc7d Mart*0540 C--   stress at the ice/water interface is computed in separate
                0541 C     routines that are called from mom_fluxform/mom_vecinv
71125c711d Mart*0542            ELSE
f39e9f371c Mart*0543             shelfIceHeatFlux      (i,j,bi,bj) = 0. _d 0
                0544             shelfIceFreshWaterFlux(i,j,bi,bj) = 0. _d 0
                0545             shelficeForcingT      (i,j,bi,bj) = 0. _d 0
                0546             shelficeForcingS      (i,j,bi,bj) = 0. _d 0
ffe464dc7d Mart*0547            ENDIF
                0548           ENDDO
                0549          ENDDO
                0550         ELSE
ead4e7560e Jean*0551 #else
ffe464dc7d Mart*0552         IF ( .TRUE. ) THEN
                0553 #endif /* ALLOW_ISOMIP_TD */
17292dde13 Mart*0554 C     use BRIOS thermodynamics, following Hellmers PhD thesis:
ead4e7560e Jean*0555 C     Hellmer, H., 1989, A two-dimensional model for the thermohaline
17292dde13 Mart*0556 C     circulation under an ice shelf, Reports on Polar Research, No. 60
                0557 C     (in German).
a8eec31802 Mart*0558
f39e9f371c Mart*0559          DO j = 1, sNy
                0560           DO i = 1, sNx
                0561            k    = kTopC(i,j,bi,bj)
                0562            IF ( k .NE. 0 .AND. pLoc(i,j) .GT. 0. _d 0 ) THEN
823dd639fd Jean*0563 C     heat flux into the ice shelf, default is diffusive flux
37e32ac0e0 Mart*0564 C     (Holland and Jenkins, 1999, eq.21)
f39e9f371c Mart*0565             thetaFreeze = a0*sLoc(i,j)+c0+b0*pLoc(i,j)
37e32ac0e0 Mart*0566             fwflxFac    = 0. _d 0
f39e9f371c Mart*0567             IF ( tLoc(i,j) .GT. thetaFreeze ) fwflxFac = dFac
e4305b0f18 Patr*0568 C     a few abbreviations
                0569             eps1 = rUnit2mass*HeatCapacity_Cp
                0570      &           *shiTransCoeffT(i,j,bi,bj)
                0571             eps2 = rUnit2mass*SHELFICElatentHeat
                0572      &           *shiTransCoeffS(i,j,bi,bj)
                0573
                0574 C     solve quadratic equation for salinity at shelfice-ocean interface
71125c711d Mart*0575 C     note: this part of the code is not very intuitive as it involves
                0576 C     many arbitrary abbreviations that were introduced to derive the
                0577 C     correct form of the quadratic equation for salinity. The abbreviations
00f81e6785 Ou W*0578 C     are explained in the documentation.
823dd639fd Jean*0579 C
296c868ef5 Mart*0580 C     eps3a was introduced as a constant variant of eps3 to avoid AD of
00f81e6785 Ou W*0581 C     code of type (pLoc-const)/pLoc
397d7a6973 Patr*0582             eps3a = rhoShelfIce*SHELFICEheatCapacity_Cp
37e32ac0e0 Mart*0583      &           * SHELFICEkappa *  ( 1. _d 0 - dFac )
f39e9f371c Mart*0584             eps3 = eps3a/pLoc(i,j)
                0585             eps4 = b0*pLoc(i,j) + c0
                0586             eps6 = eps4 - tLoc(i,j)
71125c711d Mart*0587             eps7 = eps4 - SHELFICEthetaSurface
37e32ac0e0 Mart*0588             eps8 = rUnit2mass*SHELFICEheatCapacity_Cp
                0589      &           *shiTransCoeffS(i,j,bi,bj) * fwflxFac
f39e9f371c Mart*0590             epsq = eps1*eps6 + eps3*eps7
                0591 C     coefficients of the quadratic equations for saltFreeze
                0592             aqe = a0*(eps1+eps3-eps8)
17292dde13 Mart*0593             recip_aqe = 0. _d 0
71125c711d Mart*0594             IF ( aqe .NE. 0. _d 0 ) recip_aqe = 0.5 _d 0/aqe
f39e9f371c Mart*0595             bqe = epsq - eps2 + eps8*( a0*sLoc(i,j) - eps7 )
                0596      &           - a0*(eps1+eps3)*SHELFICEsalinity
                0597             cqe = ( eps2 + eps8*eps7 )*sLoc(i,j) - epsq*SHELFICEsalinity
17292dde13 Mart*0598             discrim = bqe*bqe - 4. _d 0*aqe*cqe
71125c711d Mart*0599 #undef ALLOW_SHELFICE_DEBUG
17292dde13 Mart*0600 #ifdef ALLOW_SHELFICE_DEBUG
                0601             IF ( discrim .LT. 0. _d 0 ) THEN
                0602              print *, 'ml-shelfice: discrim = ', discrim,aqe,bqe,cqe
f39e9f371c Mart*0603              print *, 'ml-shelfice: pLoc    = ', pLoc(i,j)
                0604              print *, 'ml-shelfice: tLoc    = ', tLoc(i,j)
                0605              print *, 'ml-shelfice: sLoc    = ', sLoc(i,j)
ead4e7560e Jean*0606              print *, 'ml-shelfice: tsurface= ',
17292dde13 Mart*0607      &            SHELFICEthetaSurface
                0608              print *, 'ml-shelfice: eps1    = ', eps1
                0609              print *, 'ml-shelfice: eps2    = ', eps2
                0610              print *, 'ml-shelfice: eps3    = ', eps3
                0611              print *, 'ml-shelfice: eps4    = ', eps4
71125c711d Mart*0612              print *, 'ml-shelfice: eps6    = ', eps6
                0613              print *, 'ml-shelfice: eps7    = ', eps7
296c868ef5 Mart*0614              print *, 'ml-shelfice: eps8    = ', eps8
2bf81d7d9d Mart*0615              print *, 'ml-shelfice: rU2mass = ', rUnit2mass
17292dde13 Mart*0616              print *, 'ml-shelfice: rhoIce  = ', rhoShelfIce
a8eec31802 Mart*0617              print *, 'ml-shelfice: cFac    = ', cFac
17292dde13 Mart*0618              print *, 'ml-shelfice: Cp_W    = ', HeatCapacity_Cp
                0619              print *, 'ml-shelfice: Cp_I    = ',
                0620      &            SHELFICEHeatCapacity_Cp
ead4e7560e Jean*0621              print *, 'ml-shelfice: gammaT  = ',
17292dde13 Mart*0622      &            SHELFICEheatTransCoeff
ead4e7560e Jean*0623              print *, 'ml-shelfice: gammaS  = ',
17292dde13 Mart*0624      &            SHELFICEsaltTransCoeff
ead4e7560e Jean*0625              print *, 'ml-shelfice: lat.heat= ',
17292dde13 Mart*0626      &            SHELFICElatentHeat
                0627              STOP 'ABNORMAL END in S/R SHELFICE_THERMODYNAMICS'
                0628             ENDIF
                0629 #endif /* ALLOW_SHELFICE_DEBUG */
                0630             saltFreeze = (- bqe - SQRT(discrim))*recip_aqe
f4aa450f28 Mart*0631             IF ( saltFreeze .LT. 0. _d 0 )
17292dde13 Mart*0632      &           saltFreeze = (- bqe + SQRT(discrim))*recip_aqe
7108c05fdf Mart*0633             thetaFreeze = a0*saltFreeze + eps4
f39e9f371c Mart*0634 #ifdef SHI_SALTBAL_FWFLX
2bf81d7d9d Mart*0635 C--   upward fresh water flux due to melting (in kg/m^2/s)
397d7a6973 Patr*0636 cph change to identical form
e4305b0f18 Patr*0637 cph            freshWaterFlux = rUnit2mass
                0638 cph     &           * shiTransCoeffS(i,j,bi,bj)
f39e9f371c Mart*0639 cph     &           * ( saltFreeze - sLoc(i,j) ) / saltFreeze
                0640             freshWaterFlux = rUnit2mass * shiTransCoeffS(i,j,bi,bj)
                0641 C     Catch unlikely case of saltFreeze = 0 (often when sLoc=0)
                0642             IF ( saltFreeze .NE. 0. _d 0 )
                0643      &           freshWaterFlux = freshWaterFlux
                0644      &           * ( 1. _d 0 - sLoc(i,j) / saltFreeze )
                0645 #else /* ndef SHI_SALTBAL_FWFLX */
                0646 C     This formulation of fwflux, derived from the heat balance equation
                0647 C     instead of the salt balance equation, can handle the case when the
                0648 C     salinity of the ocean, boundary layer, and ice are identical.
                0649 C     Default thermodynamics specify a linear temperature gradient
                0650 C     through the ice (Holland and Jenkins, 1999, Section 2).
                0651 C     Here, via eps3=0 and fwflxFac, we also cover the case of a nonlinear
                0652 C     temperature gradient through the ice (Holland and Jenkins, 1999,
                0653 C     Section 3), but only for melting (Eq. 31 of Holland and Jenkins,
                0654 C     1999). fwflxFac = 0 for freezing and 1 for melting conditions.
                0655             freshWaterFlux = (
                0656      &             eps1*( thetaFreeze - tLoc(i,j) )
                0657      &           + eps3*( thetaFreeze - SHELFICEthetaSurface )
                0658      &           ) / ( SHELFICElatentHeat
                0659      &                + fwflxFac * SHELFICEheatCapacity_Cp*
                0660      &                  (thetaFreeze - SHELFICEthetaSurface) )
                0661 #endif /* ndef SHI_SALTBAL_FWFLX */
                0662
3bafcf6020 Timo*0663 #ifdef ALLOW_GENTIM2D_CONTROL
f39e9f371c Mart*0664             freshWaterFlux = freshWaterFlux + xx_shifwflx_loc(i,j,bi,bj)
                0665 #endif /* ALLOW_GENTIM2D_CONTROL */
71125c711d Mart*0666 C--   Calculate the upward heat and fresh water fluxes;
ead4e7560e Jean*0667 C--   MITgcm sign conventions: downward (negative) fresh water flux
71125c711d Mart*0668 C--   implies melting and due to upward (positive) heat flux
f39e9f371c Mart*0669             shelfIceHeatFlux(i,j,bi,bj) =
823dd639fd Jean*0670      &           ( eps3
37e32ac0e0 Mart*0671      &           - freshWaterFlux*SHELFICEheatCapacity_Cp*fwflxFac )
                0672      &           * ( thetaFreeze - SHELFICEthetaSurface )
ead4e7560e Jean*0673      &           -  cFac*freshWaterFlux*( SHELFICElatentHeat
f39e9f371c Mart*0674      &             - HeatCapacity_Cp*( thetaFreeze - rFac*tLoc(i,j) ) )
                0675             shelfIceFreshWaterFlux(i,j,bi,bj) = freshWaterFlux
17292dde13 Mart*0676 C--   compute surface tendencies
9952f046d7 dngo*0677 C--   DNG: correction to use cell value for flux rather than BL values
                0678 C--        in order to conserve salt and temp even with real FW Flux
17292dde13 Mart*0679             shelficeForcingT(i,j,bi,bj) =
e4305b0f18 Patr*0680      &           ( shiTransCoeffT(i,j,bi,bj)
9952f046d7 dngo*0681      &           - cFac*freshWaterFlux*mass2rUnit
f39e9f371c Mart*0682      &           )*( thetaFreeze - tLoc(i,j) )
9952f046d7 dngo*0683      &           - rFWinBL*freshWaterFlux*mass2rUnit
f39e9f371c Mart*0684      &            *( tLoc(i,j) - theta(i,j,k,bi,bj) )
ead4e7560e Jean*0685             shelficeForcingS(i,j,bi,bj) =
e4305b0f18 Patr*0686      &           ( shiTransCoeffS(i,j,bi,bj)
9952f046d7 dngo*0687      &           - cFac*freshWaterFlux*mass2rUnit
f39e9f371c Mart*0688      &           )*( saltFreeze - sLoc(i,j) )
9952f046d7 dngo*0689      &           - rFWinBL*freshWaterFlux*mass2rUnit
f39e9f371c Mart*0690      &            *( sLoc(i,j) - salt(i,j,k,bi,bj) )
71125c711d Mart*0691            ELSE
f39e9f371c Mart*0692             shelfIceHeatFlux      (i,j,bi,bj) = 0. _d 0
                0693             shelfIceFreshWaterFlux(i,j,bi,bj) = 0. _d 0
                0694             shelficeForcingT      (i,j,bi,bj) = 0. _d 0
                0695             shelficeForcingS      (i,j,bi,bj) = 0. _d 0
17292dde13 Mart*0696            ENDIF
ead4e7560e Jean*0697           ENDDO
17292dde13 Mart*0698          ENDDO
                0699 C     endif (not) useISOMIPTD
9952f046d7 dngo*0700         ENDIF
                0701 C     end bi,bj loops
ffe464dc7d Mart*0702        ENDDO
                0703       ENDDO
00f81e6785 Ou W*0704 #ifdef ALLOW_AUTODIFF_TAMC
                0705 C     Avoid recomputing most of the routine body in the AD-code
                0706 CADJ STORE shelfIceFreshWaterFlux
                0707 CADJ &     = comlev1, key=ikey_dynamics, kind=isbyte
                0708 #endif
ead4e7560e Jean*0709
a0178c5a01 Jean*0710       IF (SHELFICEMassStepping) THEN
3d2f509a67 Dani*0711        CALL SHELFICE_STEP_ICEMASS( myTime, myIter, myThid )
                0712       ENDIF
                0713
4c0efb30d8 Mart*0714 #ifdef ALLOW_DIAGNOSTICS
                0715       IF ( useDiagnostics ) THEN
ead4e7560e Jean*0716        CALL DIAGNOSTICS_FILL_RS(shelfIceFreshWaterFlux,'SHIfwFlx',
4c0efb30d8 Mart*0717      &      0,1,0,1,1,myThid)
ead4e7560e Jean*0718        CALL DIAGNOSTICS_FILL_RS(shelfIceHeatFlux,      'SHIhtFlx',
4c0efb30d8 Mart*0719      &      0,1,0,1,1,myThid)
a86b6152f4 Dimi*0720 C     SHIForcT (Ice shelf forcing for theta [W/m2], >0 increases theta)
                0721        tmpFac = HeatCapacity_Cp*rUnit2mass
                0722        CALL DIAGNOSTICS_SCALE_FILL(shelficeForcingT,tmpFac,1,
7462d20066 Mart*0723      &      'SHIForcT',0,1,0,1,1,myThid)
a86b6152f4 Dimi*0724 C     SHIForcS (Ice shelf forcing for salt [g/m2/s], >0 increases salt)
                0725        tmpFac = rUnit2mass
                0726        CALL DIAGNOSTICS_SCALE_FILL(shelficeForcingS,tmpFac,1,
7462d20066 Mart*0727      &      'SHIForcS',0,1,0,1,1,myThid)
                0728 C     Transfer coefficients
                0729        CALL DIAGNOSTICS_FILL(shiTransCoeffT,'SHIgammT',
                0730      &      0,1,0,1,1,myThid)
                0731        CALL DIAGNOSTICS_FILL(shiTransCoeffS,'SHIgammS',
                0732      &      0,1,0,1,1,myThid)
7b8b86ab99 Timo*0733
5744811f23 Mart*0734 C     Friction velocity
fbf1101a4a Jean*0735 #ifdef SHI_ALLOW_GAMMAFRICT
7b8b86ab99 Timo*0736        IF ( SHELFICEuseGammaFrict ) THEN
                0737         CALL DIAGNOSTICS_FILL(uStarDiag,'SHIuStar',0,1,0,1,1,myThid)
                0738         CALL DIAGNOSTICS_FILL(shiCDragFld,'SHICDrag',
                0739      &      0,1,0,1,1,myThid)
                0740        ENDIF
fbf1101a4a Jean*0741 #endif /* SHI_ALLOW_GAMMAFRICT */
9952f046d7 dngo*0742 #ifdef ALLOW_SHELFICE_REMESHING
3b86795949 Jean*0743        CALL DIAGNOSTICS_FILL_RS( R_shelfIce, 'SHIRshel',
9952f046d7 dngo*0744      &                           0, 1, 0, 1, 1, myThid )
                0745 #endif
4c0efb30d8 Mart*0746       ENDIF
                0747 #endif /* ALLOW_DIAGNOSTICS */
a86b6152f4 Dimi*0748
ffe464dc7d Mart*0749 #endif /* ALLOW_SHELFICE */
                0750       RETURN
                0751       END