Back to home page

MITgcm
 
 

    

File indexing completed on 2024-01-13 06:10:46 UTC

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

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!