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