| 1 | dgoldberg | 1.1 | C $Header: /u/gcmpack/MITgcm_contrib/shelfice_remeshing/code/shelfice_thermodynamics.F,v 1.2 2015/08/07 10:35:32 dgoldberg Exp $ | 
| 2 |  |  | C $Name:  $ | 
| 3 |  |  |  | 
| 4 |  |  | #include "SHELFICE_OPTIONS.h" | 
| 5 |  |  | #ifdef ALLOW_AUTODIFF | 
| 6 |  |  | # include "AUTODIFF_OPTIONS.h" | 
| 7 |  |  | #endif | 
| 8 |  |  | #ifdef ALLOW_CTRL | 
| 9 |  |  | # include "CTRL_OPTIONS.h" | 
| 10 |  |  | #endif | 
| 11 |  |  |  | 
| 12 |  |  | CBOP | 
| 13 |  |  | C     !ROUTINE: SHELFICE_THERMODYNAMICS | 
| 14 |  |  | C     !INTERFACE: | 
| 15 |  |  | SUBROUTINE SHELFICE_THERMODYNAMICS( | 
| 16 |  |  | I                        myTime, myIter, myThid ) | 
| 17 |  |  | C     !DESCRIPTION: \bv | 
| 18 |  |  | C     *=============================================================* | 
| 19 |  |  | C     | S/R  SHELFICE_THERMODYNAMICS | 
| 20 |  |  | C     | o shelf-ice main routine. | 
| 21 |  |  | C     |   compute temperature and (virtual) salt flux at the | 
| 22 |  |  | C     |   shelf-ice ocean interface | 
| 23 |  |  | C     | | 
| 24 |  |  | C     | stresses at the ice/water interface are computed in separate | 
| 25 |  |  | C     | routines that are called from mom_fluxform/mom_vecinv | 
| 26 |  |  | C     *=============================================================* | 
| 27 |  |  | C     \ev | 
| 28 |  |  |  | 
| 29 |  |  | C     !USES: | 
| 30 |  |  | IMPLICIT NONE | 
| 31 |  |  |  | 
| 32 |  |  | C     === Global variables === | 
| 33 |  |  | #include "SIZE.h" | 
| 34 |  |  | #include "EEPARAMS.h" | 
| 35 |  |  | #include "PARAMS.h" | 
| 36 |  |  | #include "GRID.h" | 
| 37 |  |  | #include "DYNVARS.h" | 
| 38 |  |  | #include "FFIELDS.h" | 
| 39 |  |  | #include "SHELFICE.h" | 
| 40 |  |  | #include "SHELFICE_COST.h" | 
| 41 |  |  | #ifdef ALLOW_AUTODIFF | 
| 42 |  |  | # include "CTRL_SIZE.h" | 
| 43 |  |  | # include "ctrl.h" | 
| 44 |  |  | # include "ctrl_dummy.h" | 
| 45 |  |  | #endif /* ALLOW_AUTODIFF */ | 
| 46 |  |  | #ifdef ALLOW_AUTODIFF_TAMC | 
| 47 |  |  | # ifdef SHI_ALLOW_GAMMAFRICT | 
| 48 |  |  | #  include "tamc.h" | 
| 49 |  |  | #  include "tamc_keys.h" | 
| 50 |  |  | # endif /* SHI_ALLOW_GAMMAFRICT */ | 
| 51 |  |  | #endif /* ALLOW_AUTODIFF_TAMC */ | 
| 52 |  |  |  | 
| 53 |  |  | C     !INPUT/OUTPUT PARAMETERS: | 
| 54 |  |  | C     === Routine arguments === | 
| 55 |  |  | C     myIter :: iteration counter for this thread | 
| 56 |  |  | C     myTime :: time counter for this thread | 
| 57 |  |  | C     myThid :: thread number for this instance of the routine. | 
| 58 |  |  | _RL  myTime | 
| 59 |  |  | INTEGER myIter | 
| 60 |  |  | INTEGER myThid | 
| 61 |  |  |  | 
| 62 |  |  | #ifdef ALLOW_SHELFICE | 
| 63 |  |  | C     !LOCAL VARIABLES : | 
| 64 |  |  | C     === Local variables === | 
| 65 |  |  | C     I,J,K,Kp1,bi,bj  :: loop counters | 
| 66 |  |  | C     tLoc, sLoc, pLoc :: local in-situ temperature, salinity, pressure | 
| 67 |  |  | C     theta/saltFreeze :: temperature and salinity of water at the | 
| 68 |  |  | C                         ice-ocean interface (at the freezing point) | 
| 69 |  |  | C     freshWaterFlux   :: local variable for fresh water melt flux due | 
| 70 |  |  | C                         to melting in kg/m^2/s | 
| 71 |  |  | C                         (negative density x melt rate) | 
| 72 |  |  | C     convertFW2SaltLoc:: local copy of convertFW2Salt | 
| 73 |  |  | C     cFac             :: 1 for conservative form, 0, otherwise | 
| 74 |  |  | C     rFac             :: realFreshWaterFlux factor | 
| 75 |  |  | C     dFac             :: 0 for diffusive heat flux (Holland and Jenkins, 1999, | 
| 76 |  |  | C                           eq21) | 
| 77 |  |  | C                         1 for advective and diffusive heat flux (eq22, 26, 31) | 
| 78 |  |  | C     fwflxFac         :: only effective for dFac=1, 1 if we expect a melting | 
| 79 |  |  | C                         fresh water flux, 0 otherwise | 
| 80 |  |  | C     auxiliary variables and abbreviations: | 
| 81 |  |  | C     a0, a1, a2, b, c0 | 
| 82 |  |  | C     eps1, eps2, eps3, eps3a, eps4, eps5, eps6, eps7, eps8 | 
| 83 |  |  | C     aqe, bqe, cqe, discrim, recip_aqe | 
| 84 |  |  | C     drKp1, recip_drLoc | 
| 85 |  |  | INTEGER I,J,K,Kp1 | 
| 86 |  |  | INTEGER bi,bj | 
| 87 |  |  | _RL tLoc(1:sNx,1:sNy) | 
| 88 |  |  | _RL sLoc(1:sNx,1:sNy) | 
| 89 |  |  | _RL pLoc(1:sNx,1:sNy) | 
| 90 |  |  | _RL uLoc(1:sNx,1:sNy) | 
| 91 |  |  | _RL vLoc(1:sNx,1:sNy) | 
| 92 |  |  | _RL u_topdr(1:sNx+1,1:sNy+1,nSx,nSy) | 
| 93 |  |  | _RL v_topdr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) | 
| 94 |  |  | _RL thetaFreeze, saltFreeze, recip_Cp | 
| 95 |  |  | _RL freshWaterFlux, convertFW2SaltLoc | 
| 96 |  |  | _RL a0, a1, a2, b, c0 | 
| 97 |  |  | _RL eps1, eps2, eps3, eps3a, eps4, eps5, eps6, eps7, eps8 | 
| 98 |  |  | _RL cFac, rFac, dFac, fwflxFac | 
| 99 |  |  | _RL aqe, bqe, cqe, discrim, recip_aqe | 
| 100 |  |  | _RL drKp1, recip_drLoc | 
| 101 |  |  | _RL recip_latentHeat | 
| 102 |  |  | _RL tmpFac | 
| 103 |  |  |  | 
| 104 |  |  | #ifdef SHI_ALLOW_GAMMAFRICT | 
| 105 |  |  | _RL shiPr, shiSc, shiLo, recip_shiKarman, shiTwoThirds | 
| 106 |  |  | _RL gammaTmoleT, gammaTmoleS, gammaTurb, gammaTurbConst | 
| 107 |  |  | _RL ustar, ustarSq, etastar | 
| 108 |  |  | PARAMETER ( shiTwoThirds = 0.66666666666666666666666666667D0 ) | 
| 109 |  |  | #ifdef ALLOW_DIAGNOSTICS | 
| 110 |  |  | _RL uStarDiag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) | 
| 111 |  |  | #endif /* ALLOW_DIAGNOSTICS */ | 
| 112 |  |  | #endif | 
| 113 |  |  |  | 
| 114 |  |  | #ifndef ALLOW_OPENAD | 
| 115 |  |  | _RL SW_TEMP | 
| 116 |  |  | EXTERNAL SW_TEMP | 
| 117 |  |  | #endif | 
| 118 |  |  |  | 
| 119 |  |  | #ifdef ALLOW_SHIFWFLX_CONTROL | 
| 120 |  |  | _RL xx_shifwflx_loc(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) | 
| 121 |  |  | #endif | 
| 122 |  |  | CEOP | 
| 123 |  |  | C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| | 
| 124 |  |  |  | 
| 125 |  |  | #ifdef SHI_ALLOW_GAMMAFRICT | 
| 126 |  |  | #ifdef ALLOW_AUTODIFF | 
| 127 |  |  | C     re-initialize here again, curtesy to TAF | 
| 128 |  |  | DO bj = myByLo(myThid), myByHi(myThid) | 
| 129 |  |  | DO bi = myBxLo(myThid), myBxHi(myThid) | 
| 130 |  |  | DO J = 1-OLy,sNy+OLy | 
| 131 |  |  | DO I = 1-OLx,sNx+OLx | 
| 132 |  |  | shiTransCoeffT(i,j,bi,bj) = SHELFICEheatTransCoeff | 
| 133 |  |  | shiTransCoeffS(i,j,bi,bj) = SHELFICEsaltTransCoeff | 
| 134 |  |  | ENDDO | 
| 135 |  |  | ENDDO | 
| 136 |  |  | ENDDO | 
| 137 |  |  | ENDDO | 
| 138 |  |  | #endif /* ALLOW_AUTODIFF */ | 
| 139 |  |  | IF ( SHELFICEuseGammaFrict ) THEN | 
| 140 |  |  | C     Implement friction velocity-dependent transfer coefficient | 
| 141 |  |  | C     of Holland and Jenkins, JPO, 1999 | 
| 142 |  |  | recip_shiKarman= 1. _d 0 / 0.4 _d 0 | 
| 143 |  |  | shiLo = 0. _d 0 | 
| 144 |  |  | shiPr = shiPrandtl**shiTwoThirds | 
| 145 |  |  | shiSc = shiSchmidt**shiTwoThirds | 
| 146 |  |  | cph      shiPr = (viscArNr(1)/diffKrNrT(1))**shiTwoThirds | 
| 147 |  |  | cph      shiSc = (viscArNr(1)/diffKrNrS(1))**shiTwoThirds | 
| 148 |  |  | gammaTmoleT = 12.5 _d 0 * shiPr - 6. _d 0 | 
| 149 |  |  | gammaTmoleS = 12.5 _d 0 * shiSc - 6. _d 0 | 
| 150 |  |  | C     instead of etastar = sqrt(1+zetaN*ustar./(f*Lo*Rc)) | 
| 151 |  |  | etastar = 1. _d 0 | 
| 152 |  |  | gammaTurbConst  = 1. _d 0 / (2. _d 0 * shiZetaN*etastar) | 
| 153 |  |  | &      - recip_shiKarman | 
| 154 |  |  | #ifdef ALLOW_AUTODIFF | 
| 155 |  |  | DO bj = myByLo(myThid), myByHi(myThid) | 
| 156 |  |  | DO bi = myBxLo(myThid), myBxHi(myThid) | 
| 157 |  |  | DO J = 1-OLy,sNy+OLy | 
| 158 |  |  | DO I = 1-OLx,sNx+OLx | 
| 159 |  |  | shiTransCoeffT(i,j,bi,bj) = 0. _d 0 | 
| 160 |  |  | shiTransCoeffS(i,j,bi,bj) = 0. _d 0 | 
| 161 |  |  | ENDDO | 
| 162 |  |  | ENDDO | 
| 163 |  |  | ENDDO | 
| 164 |  |  | ENDDO | 
| 165 |  |  | #endif /* ALLOW_AUTODIFF */ | 
| 166 |  |  | ENDIF | 
| 167 |  |  | #endif /* SHI_ALLOW_GAMMAFRICT */ | 
| 168 |  |  |  | 
| 169 |  |  | recip_latentHeat = 0. _d 0 | 
| 170 |  |  | IF ( SHELFICElatentHeat .NE. 0. _d 0 ) | 
| 171 |  |  | &     recip_latentHeat = 1. _d 0/SHELFICElatentHeat | 
| 172 |  |  | C     are we doing the conservative form of Jenkins et al. (2001)? | 
| 173 |  |  | recip_Cp = 1. _d 0 / HeatCapacity_Cp | 
| 174 |  |  | cFac = 0. _d 0 | 
| 175 |  |  | IF ( SHELFICEconserve ) cFac = 1. _d 0 | 
| 176 |  |  | C     with "real fresh water flux" (affecting ETAN), | 
| 177 |  |  | C     there is more to modify | 
| 178 |  |  | rFac = 1. _d 0 | 
| 179 |  |  | IF ( SHELFICEconserve .AND. useRealFreshWaterFlux ) rFac = 0. _d 0 | 
| 180 |  |  | C     heat flux into the ice shelf, default is diffusive flux | 
| 181 |  |  | C     (Holland and Jenkins, 1999, eq.21) | 
| 182 |  |  | dFac = 0. _d 0 | 
| 183 |  |  | IF ( SHELFICEadvDiffHeatFlux ) dFac = 1. _d 0 | 
| 184 |  |  | fwflxFac = 0. _d 0 | 
| 185 |  |  | C     linear dependence of freezing point on salinity | 
| 186 |  |  | a0 = -0.0575   _d  0 | 
| 187 |  |  | a1 =  0.0      _d -0 | 
| 188 |  |  | a2 =  0.0      _d -0 | 
| 189 |  |  | c0 =  0.0901   _d  0 | 
| 190 |  |  | b  =  -7.61    _d -4 | 
| 191 |  |  | #ifdef ALLOW_ISOMIP_TD | 
| 192 |  |  | IF ( useISOMIPTD ) THEN | 
| 193 |  |  | C     non-linear dependence of freezing point on salinity | 
| 194 |  |  | a0 = -0.0575   _d  0 | 
| 195 |  |  | a1 = 1.710523  _d -3 | 
| 196 |  |  | a2 = -2.154996 _d -4 | 
| 197 |  |  | b  = -7.53     _d -4 | 
| 198 |  |  | c0 = 0. _d 0 | 
| 199 |  |  | ENDIF | 
| 200 |  |  | convertFW2SaltLoc = convertFW2Salt | 
| 201 |  |  | C     hardcoding this value here is OK because it only applies to ISOMIP | 
| 202 |  |  | C     where this value is part of the protocol | 
| 203 |  |  | IF ( convertFW2SaltLoc .EQ. -1. ) convertFW2SaltLoc = 33.4 _d 0 | 
| 204 |  |  | #endif /* ALLOW_ISOMIP_TD */ | 
| 205 |  |  |  | 
| 206 |  |  | DO bj = myByLo(myThid), myByHi(myThid) | 
| 207 |  |  | DO bi = myBxLo(myThid), myBxHi(myThid) | 
| 208 |  |  | DO J = 1-OLy,sNy+OLy | 
| 209 |  |  | DO I = 1-OLx,sNx+OLx | 
| 210 |  |  | shelfIceHeatFlux      (I,J,bi,bj) = 0. _d 0 | 
| 211 |  |  | shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0 | 
| 212 |  |  | shelficeForcingT      (I,J,bi,bj) = 0. _d 0 | 
| 213 |  |  | shelficeForcingS      (I,J,bi,bj) = 0. _d 0 | 
| 214 |  |  | #if (defined SHI_ALLOW_GAMMAFRICT && defined ALLOW_DIAGNOSTICS) | 
| 215 |  |  | uStarDiag             (I,J,bi,bj) = 0. _d 0 | 
| 216 |  |  | #endif /* SHI_ALLOW_GAMMAFRICT and ALLOW_DIAGNOSTICS */ | 
| 217 |  |  | ENDDO | 
| 218 |  |  | ENDDO | 
| 219 |  |  | ENDDO | 
| 220 |  |  | ENDDO | 
| 221 |  |  | #ifdef ALLOW_SHIFWFLX_CONTROL | 
| 222 |  |  | DO bj = myByLo(myThid), myByHi(myThid) | 
| 223 |  |  | DO bi = myBxLo(myThid), myBxHi(myThid) | 
| 224 |  |  | DO J = 1-OLy,sNy+OLy | 
| 225 |  |  | DO I = 1-OLx,sNx+OLx | 
| 226 |  |  | xx_shifwflx_loc(I,J,bi,bj) = 0. _d 0 | 
| 227 |  |  | ENDDO | 
| 228 |  |  | ENDDO | 
| 229 |  |  | ENDDO | 
| 230 |  |  | ENDDO | 
| 231 |  |  | #ifdef ALLOW_CTRL | 
| 232 |  |  | if (useCTRL) CALL CTRL_GET_GEN ( | 
| 233 |  |  | &     xx_shifwflx_file, xx_shifwflxstartdate, xx_shifwflxperiod, | 
| 234 |  |  | &     maskSHI, xx_shifwflx_loc, xx_shifwflx0, xx_shifwflx1, | 
| 235 |  |  | &     xx_shifwflx_dummy, | 
| 236 |  |  | &     xx_shifwflx_remo_intercept, xx_shifwflx_remo_slope, | 
| 237 |  |  | &     wshifwflx, | 
| 238 |  |  | &     myTime, myIter, myThid ) | 
| 239 |  |  | #endif | 
| 240 |  |  | #endif /* ALLOW_SHIFWFLX_CONTROL */ | 
| 241 |  |  | DO bj = myByLo(myThid), myByHi(myThid) | 
| 242 |  |  | DO bi = myBxLo(myThid), myBxHi(myThid) | 
| 243 |  |  |  | 
| 244 |  |  | IF ( SHELFICEBoundaryLayer ) THEN | 
| 245 |  |  | C--   average over boundary layer width | 
| 246 |  |  | DO J = 1, sNy+1 | 
| 247 |  |  | DO I = 1, sNx+1 | 
| 248 |  |  | u_topdr(I,J,bi,bj) = 0.0 | 
| 249 |  |  | v_topdr(I,J,bi,bj) = 0.0 | 
| 250 |  |  | ENDDO | 
| 251 |  |  | ENDDO | 
| 252 |  |  | ENDIF | 
| 253 |  |  |  | 
| 254 |  |  | #ifdef ALLOW_AUTODIFF_TAMC | 
| 255 |  |  | # ifdef SHI_ALLOW_GAMMAFRICT | 
| 256 |  |  | act1 = bi - myBxLo(myThid) | 
| 257 |  |  | max1 = myBxHi(myThid) - myBxLo(myThid) + 1 | 
| 258 |  |  | act2 = bj - myByLo(myThid) | 
| 259 |  |  | max2 = myByHi(myThid) - myByLo(myThid) + 1 | 
| 260 |  |  | act3 = myThid - 1 | 
| 261 |  |  | max3 = nTx*nTy | 
| 262 |  |  | act4 = ikey_dynamics - 1 | 
| 263 |  |  | ikey = (act1 + 1) + act2*max1 | 
| 264 |  |  | &                    + act3*max1*max2 | 
| 265 |  |  | &                    + act4*max1*max2*max3 | 
| 266 |  |  | # endif /* SHI_ALLOW_GAMMAFRICT */ | 
| 267 |  |  | #endif /* ALLOW_AUTODIFF_TAMC */ | 
| 268 |  |  | DO J = 1, sNy | 
| 269 |  |  | DO I = 1, sNx | 
| 270 |  |  | C--   make local copies of temperature, salinity and depth (pressure in deci-bar) | 
| 271 |  |  | C--   underneath the ice | 
| 272 |  |  | K         = MAX(1,kTopC(I,J,bi,bj)) | 
| 273 |  |  | pLoc(I,J) = ABS(R_shelfIce(I,J,bi,bj)) | 
| 274 |  |  | c         pLoc(I,J) = shelficeMass(I,J,bi,bj)*gravity*1. _d -4 | 
| 275 |  |  | tLoc(I,J) = theta(I,J,K,bi,bj) | 
| 276 |  |  | sLoc(I,J) = MAX(salt(I,J,K,bi,bj), zeroRL) | 
| 277 |  |  | c              print *, 'JJHA' , hFacC(I,J,K,bi,bj) | 
| 278 |  |  | IF ( .not.SHELFICEBoundaryLayer ) THEN | 
| 279 |  |  | uLoc(I,J) = recip_hFacC(I,J,K,bi,bj) * | 
| 280 |  |  | &         ( uVel(I,  J,K,bi,bj) * _hFacW(I,  J,K,bi,bj) | 
| 281 |  |  | &         + uVel(I+1,J,K,bi,bj) * _hFacW(I+1,J,K,bi,bj) ) | 
| 282 |  |  | vLoc(I,J) = recip_hFacC(I,J,K,bi,bj) * | 
| 283 |  |  | &         ( vVel(I,  J,K,bi,bj) * _hFacS(I,  J,K,bi,bj) | 
| 284 |  |  | &         + vVel(I,J+1,K,bi,bj) * _hFacS(I,J+1,K,bi,bj) ) | 
| 285 |  |  | ENDIF | 
| 286 |  |  | !           tmpfac = | 
| 287 |  |  | !      &         ( _hFacS(I,J,  K,bi,bj) + | 
| 288 |  |  | !      &           _hFacS(I,J+1,K,bi,bj) ) / 2.0 | 
| 289 |  |  | !           vLoc(I,J) = | 
| 290 |  |  | !      &         ( vVel(I,J,  K,bi,bj) * _hFacS(I,J,  K,bi,bj) | 
| 291 |  |  | !      &         + vVel(I,J+1,K,bi,bj) * _hFacS(I,J+1,K,bi,bj) ) | 
| 292 |  |  | !           if (tmpfac .gt. 0.0) then | 
| 293 |  |  | !                if (j.lt.100) then | 
| 294 |  |  | !                print *, 1./tmpfac, recip_hFacC(I,J,K,bi,bj) | 
| 295 |  |  | !                endif | 
| 296 |  |  | !                vLoc(I,J) = vLoc(I,J) / tmpfac | 
| 297 |  |  | !           else | 
| 298 |  |  | !                vLoc(I,J) = 0.0 | 
| 299 |  |  | !           endif | 
| 300 |  |  | ENDDO | 
| 301 |  |  | ENDDO | 
| 302 |  |  |  | 
| 303 |  |  | IF ( SHELFICEBoundaryLayer ) THEN | 
| 304 |  |  | DO J = 1, sNy+1 | 
| 305 |  |  | DO I = 1, sNx+1 | 
| 306 |  |  | K = ksurfW(I,J,bi,bj) | 
| 307 |  |  | Kp1 = K+1 | 
| 308 |  |  | IF (K.lt.Nr) then | 
| 309 |  |  | drKp1 = drF(K)*(1. _d 0-_hFacW(I,J,K,bi,bj)) | 
| 310 |  |  | drKp1 = max (drKp1, 0. _d 0) | 
| 311 |  |  | recip_drLoc = 1.0 / | 
| 312 |  |  | &       (drF(K)*_hFacW(I,J,K,bi,bj)+drKp1) | 
| 313 |  |  | u_topdr(I,J,bi,bj) = | 
| 314 |  |  | &       (drF(K)*_hFacW(I,J,K,bi,bj)*uVel(I,J,K,bi,bj) + | 
| 315 |  |  | &        drKp1*uVel(I,J,Kp1,bi,bj)) | 
| 316 |  |  | &      * recip_drLoc | 
| 317 |  |  | ELSE | 
| 318 |  |  | u_topdr(I,J,bi,bj) = 0. _d 0 | 
| 319 |  |  | ENDIF | 
| 320 |  |  |  | 
| 321 |  |  | K = ksurfS(I,J,bi,bj) | 
| 322 |  |  | Kp1 = K+1 | 
| 323 |  |  | IF (K.lt.Nr) then | 
| 324 |  |  | drKp1 = drF(K)*(1. _d 0-_hFacS(I,J,K,bi,bj)) | 
| 325 |  |  | drKp1 = max (drKp1, 0. _d 0) | 
| 326 |  |  | recip_drLoc = 1.0 / | 
| 327 |  |  | &       (drF(K)*_hFacS(I,J,K,bi,bj)+drKp1) | 
| 328 |  |  | v_topdr(I,J,bi,bj) = | 
| 329 |  |  | &       (drF(K)*_hFacS(I,J,K,bi,bj)*vVel(I,J,K,bi,bj) + | 
| 330 |  |  | &        drKp1*vVel(I,J,Kp1,bi,bj)) | 
| 331 |  |  | &      * recip_drLoc | 
| 332 |  |  | ELSE | 
| 333 |  |  | v_topdr(I,J,bi,bj) = 0. _d 0 | 
| 334 |  |  | ENDIF | 
| 335 |  |  |  | 
| 336 |  |  | ENDDO | 
| 337 |  |  | ENDDO | 
| 338 |  |  | ENDIF | 
| 339 |  |  |  | 
| 340 |  |  | IF ( SHELFICEBoundaryLayer ) THEN | 
| 341 |  |  | C--   average over boundary layer width | 
| 342 |  |  | DO J = 1, sNy | 
| 343 |  |  | DO I = 1, sNx | 
| 344 |  |  | K   = kTopC(I,J,bi,bj) | 
| 345 |  |  | IF ( K .NE. 0 .AND. K .LT. Nr ) THEN | 
| 346 |  |  | Kp1 = MIN(Nr,K+1) | 
| 347 |  |  | C--   overlap into lower cell | 
| 348 |  |  | drKp1 = drF(K)*( 1. _d 0 - _hFacC(I,J,K,bi,bj) ) | 
| 349 |  |  | C-- Dans fix | 
| 350 |  |  | drKp1 = MAX(drKp1, 0.) | 
| 351 |  |  | C--   lower cell may not be as thick as required | 
| 352 |  |  | drKp1 = MIN( drKp1, drF(Kp1) * _hFacC(I,J,Kp1,bi,bj) ) | 
| 353 |  |  | recip_drLoc = 1. _d 0 / | 
| 354 |  |  | &           ( drF(K)*_hFacC(I,J,K,bi,bj) + drKp1 ) | 
| 355 |  |  | tLoc(I,J) = ( tLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj) | 
| 356 |  |  | &           + theta(I,J,Kp1,bi,bj) *drKp1 ) | 
| 357 |  |  | &           * recip_drLoc | 
| 358 |  |  | sLoc(I,J) = ( sLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj) | 
| 359 |  |  | &           + MAX(salt(I,J,Kp1,bi,bj), zeroRL) * drKp1 ) | 
| 360 |  |  | &           * recip_drLoc | 
| 361 |  |  |  | 
| 362 |  |  | !              uLoc(I,J) = ( uLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj) | 
| 363 |  |  | !      &           + drKp1 * recip_hFacC(I,J,Kp1,bi,bj) * | 
| 364 |  |  | !      &           ( uVel(I,  J,Kp1,bi,bj) * _hFacW(I,  J,Kp1,bi,bj) | 
| 365 |  |  | !      &           + uVel(I+1,J,Kp1,bi,bj) * _hFacW(I+1,J,Kp1,bi,bj) ) | 
| 366 |  |  | !      &           ) * recip_drLoc | 
| 367 |  |  | !             vLoc(I,J) = ( vLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj) | 
| 368 |  |  | !      &           + drKp1 * recip_hFacC(I,J,Kp1,bi,bj) * | 
| 369 |  |  | !      &           ( vVel(I,J,  Kp1,bi,bj) * _hFacS(I,J,  Kp1,bi,bj) | 
| 370 |  |  | !      &           + vVel(I,J+1,Kp1,bi,bj) * _hFacS(I,J+1,Kp1,bi,bj) ) | 
| 371 |  |  | !      &           ) * recip_drLoc | 
| 372 |  |  | ENDIF | 
| 373 |  |  | ENDDO | 
| 374 |  |  | ENDDO | 
| 375 |  |  | ENDIF | 
| 376 |  |  |  | 
| 377 |  |  | IF ( SHELFICEBoundaryLayer ) THEN | 
| 378 |  |  | DO J = 1, sNy | 
| 379 |  |  | DO I = 1, sNx | 
| 380 |  |  | uLoc(I,J) = | 
| 381 |  |  | &      u_topdr(I,J,bi,bj) + u_topdr(I+1,J,bi,bj) | 
| 382 |  |  | vLoc(I,J) = | 
| 383 |  |  | &      v_topdr(I,J,bi,bj) + v_topdr(I,J+1,bi,bj) | 
| 384 |  |  | ENDDO | 
| 385 |  |  | ENDDO | 
| 386 |  |  | ENDIF | 
| 387 |  |  |  | 
| 388 |  |  | C--   turn potential temperature into in-situ temperature relative | 
| 389 |  |  | C--   to the surface | 
| 390 |  |  | DO J = 1, sNy | 
| 391 |  |  | DO I = 1, sNx | 
| 392 |  |  | #ifndef ALLOW_OPENAD | 
| 393 |  |  | tLoc(I,J) = SW_TEMP(sLoc(I,J),tLoc(I,J),pLoc(I,J),zeroRL) | 
| 394 |  |  | #else | 
| 395 |  |  | CALL SW_TEMP(sLoc(I,J),tLoc(I,J),pLoc(I,J),zeroRL,tLoc(I,J)) | 
| 396 |  |  | #endif | 
| 397 |  |  | ENDDO | 
| 398 |  |  | ENDDO | 
| 399 |  |  |  | 
| 400 |  |  | #ifdef SHI_ALLOW_GAMMAFRICT | 
| 401 |  |  | IF ( SHELFICEuseGammaFrict ) THEN | 
| 402 |  |  | DO J = 1, sNy | 
| 403 |  |  | DO I = 1, sNx | 
| 404 |  |  | K = kTopC(I,J,bi,bj) | 
| 405 |  |  | IF ( K .NE. 0 .AND. pLoc(I,J) .GT. 0. _d 0 ) THEN | 
| 406 |  |  | ustarSq = shiCdrag * MAX( 1.D-6, | 
| 407 |  |  | &           0.25 _d 0 *(uLoc(I,J)*uLoc(I,J)+vLoc(I,J)*vLoc(I,J)) ) | 
| 408 |  |  | ustar   = SQRT(ustarSq) | 
| 409 |  |  | #ifdef ALLOW_DIAGNOSTICS | 
| 410 |  |  | uStarDiag(I,J,bi,bj) = ustar | 
| 411 |  |  | #endif /* ALLOW_DIAGNOSTICS */ | 
| 412 |  |  | C     instead of etastar = sqrt(1+zetaN*ustar./(f*Lo*Rc)) | 
| 413 |  |  | C           etastar = 1. _d 0 | 
| 414 |  |  | C           gammaTurbConst  = 1. _d 0 / (2. _d 0 * shiZetaN*etastar) | 
| 415 |  |  | C    &           - recip_shiKarman | 
| 416 |  |  | IF ( fCori(I,J,bi,bj) .NE. 0. _d 0 ) THEN | 
| 417 |  |  | gammaTurb = LOG( ustarSq * shiZetaN * etastar**2 | 
| 418 |  |  | &            / ABS(fCori(I,J,bi,bj) * 5.0 _d 0 * shiKinVisc)) | 
| 419 |  |  | &            * recip_shiKarman | 
| 420 |  |  | &            + gammaTurbConst | 
| 421 |  |  | C     Do we need to catch the unlikely case of very small ustar | 
| 422 |  |  | C     that can lead to negative gammaTurb? | 
| 423 |  |  | C            gammaTurb = MAX(0.D0, gammaTurb) | 
| 424 |  |  | ELSE | 
| 425 |  |  | gammaTurb = gammaTurbConst | 
| 426 |  |  | ENDIF | 
| 427 |  |  | shiTransCoeffT(i,j,bi,bj) = MAX( zeroRL, | 
| 428 |  |  | &           ustar/(gammaTurb + gammaTmoleT) ) | 
| 429 |  |  | shiTransCoeffS(i,j,bi,bj) = MAX( zeroRL, | 
| 430 |  |  | &           ustar/(gammaTurb + gammaTmoleS) ) | 
| 431 |  |  | ENDIF | 
| 432 |  |  | ENDDO | 
| 433 |  |  | ENDDO | 
| 434 |  |  | ENDIF | 
| 435 |  |  | #endif /* SHI_ALLOW_GAMMAFRICT */ | 
| 436 |  |  |  | 
| 437 |  |  | #ifdef ALLOW_AUTODIFF_TAMC | 
| 438 |  |  | # ifdef SHI_ALLOW_GAMMAFRICT | 
| 439 |  |  | CADJ STORE shiTransCoeffS(:,:,bi,bj) = comlev1_bibj, | 
| 440 |  |  | CADJ &     key=ikey, byte=isbyte | 
| 441 |  |  | CADJ STORE shiTransCoeffT(:,:,bi,bj) = comlev1_bibj, | 
| 442 |  |  | CADJ &     key=ikey, byte=isbyte | 
| 443 |  |  | # endif /* SHI_ALLOW_GAMMAFRICT */ | 
| 444 |  |  | #endif /* ALLOW_AUTODIFF_TAMC */ | 
| 445 |  |  | #ifdef ALLOW_ISOMIP_TD | 
| 446 |  |  | IF ( useISOMIPTD ) THEN | 
| 447 |  |  | DO J = 1, sNy | 
| 448 |  |  | DO I = 1, sNx | 
| 449 |  |  | K = kTopC(I,J,bi,bj) | 
| 450 |  |  | IF ( K .NE. 0 .AND. pLoc(I,J) .GT. 0. _d 0 ) THEN | 
| 451 |  |  | C--   Calculate freezing temperature as a function of salinity and pressure | 
| 452 |  |  | thetaFreeze = | 
| 453 |  |  | &           sLoc(I,J) * ( a0 + a1*sqrt(sLoc(I,J)) + a2*sLoc(I,J) ) | 
| 454 |  |  | &           + b*pLoc(I,J) + c0 | 
| 455 |  |  | C--   Calculate the upward heat and  fresh water fluxes | 
| 456 |  |  | shelfIceHeatFlux(I,J,bi,bj) = maskC(I,J,K,bi,bj) | 
| 457 |  |  | &           * shiTransCoeffT(i,j,bi,bj) | 
| 458 |  |  | &           * ( tLoc(I,J) - thetaFreeze ) | 
| 459 |  |  | &           * HeatCapacity_Cp*rUnit2mass | 
| 460 |  |  | #ifdef ALLOW_SHIFWFLX_CONTROL | 
| 461 |  |  | &           - xx_shifwflx_loc(I,J,bi,bj)*SHELFICElatentHeat | 
| 462 |  |  | #endif /*  ALLOW_SHIFWFLX_CONTROL */ | 
| 463 |  |  | C     upward heat flux into the shelf-ice implies basal melting, | 
| 464 |  |  | C     thus a downward (negative upward) fresh water flux (as a mass flux), | 
| 465 |  |  | C     and vice versa | 
| 466 |  |  | shelfIceFreshWaterFlux(I,J,bi,bj) = | 
| 467 |  |  | &           - shelfIceHeatFlux(I,J,bi,bj) | 
| 468 |  |  | &           *recip_latentHeat | 
| 469 |  |  | C--   compute surface tendencies | 
| 470 |  |  | shelficeForcingT(i,j,bi,bj) = | 
| 471 |  |  | &           - shelfIceHeatFlux(I,J,bi,bj) | 
| 472 |  |  | &           *recip_Cp*mass2rUnit | 
| 473 |  |  | &           - cFac * shelfIceFreshWaterFlux(I,J,bi,bj)*mass2rUnit | 
| 474 |  |  | &           * ( thetaFreeze - tLoc(I,J) ) | 
| 475 |  |  | shelficeForcingS(i,j,bi,bj) = | 
| 476 |  |  | &           shelfIceFreshWaterFlux(I,J,bi,bj) * mass2rUnit | 
| 477 |  |  | &           * ( cFac*sLoc(I,J) + (1. _d 0-cFac)*convertFW2SaltLoc ) | 
| 478 |  |  | C--   stress at the ice/water interface is computed in separate | 
| 479 |  |  | C     routines that are called from mom_fluxform/mom_vecinv | 
| 480 |  |  | ELSE | 
| 481 |  |  | shelfIceHeatFlux      (I,J,bi,bj) = 0. _d 0 | 
| 482 |  |  | shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0 | 
| 483 |  |  | shelficeForcingT      (I,J,bi,bj) = 0. _d 0 | 
| 484 |  |  | shelficeForcingS      (I,J,bi,bj) = 0. _d 0 | 
| 485 |  |  | ENDIF | 
| 486 |  |  | ENDDO | 
| 487 |  |  | ENDDO | 
| 488 |  |  | ELSE | 
| 489 |  |  | #else | 
| 490 |  |  | IF ( .TRUE. ) THEN | 
| 491 |  |  | #endif /* ALLOW_ISOMIP_TD */ | 
| 492 |  |  | C     use BRIOS thermodynamics, following Hellmers PhD thesis: | 
| 493 |  |  | C     Hellmer, H., 1989, A two-dimensional model for the thermohaline | 
| 494 |  |  | C     circulation under an ice shelf, Reports on Polar Research, No. 60 | 
| 495 |  |  | C     (in German). | 
| 496 |  |  |  | 
| 497 |  |  | DO J = 1, sNy | 
| 498 |  |  | DO I = 1, sNx | 
| 499 |  |  | K    = kTopC(I,J,bi,bj) | 
| 500 |  |  | IF ( K .NE. 0 .AND. pLoc(I,J) .GT. 0. _d 0 ) THEN | 
| 501 |  |  | C     heat flux into the ice shelf, default is diffusive flux | 
| 502 |  |  | C     (Holland and Jenkins, 1999, eq.21) | 
| 503 |  |  | thetaFreeze = a0*sLoc(I,J)+c0+b*pLoc(I,J) | 
| 504 |  |  | fwflxFac    = 0. _d 0 | 
| 505 |  |  | IF ( tLoc(I,J) .GT. thetaFreeze ) fwflxFac = dFac | 
| 506 |  |  | C     a few abbreviations | 
| 507 |  |  | eps1 = rUnit2mass*HeatCapacity_Cp | 
| 508 |  |  | &           *shiTransCoeffT(i,j,bi,bj) | 
| 509 |  |  | eps2 = rUnit2mass*SHELFICElatentHeat | 
| 510 |  |  | &           *shiTransCoeffS(i,j,bi,bj) | 
| 511 |  |  | eps5 = rUnit2mass*HeatCapacity_Cp | 
| 512 |  |  | &           *shiTransCoeffS(i,j,bi,bj) | 
| 513 |  |  |  | 
| 514 |  |  | C     solve quadratic equation for salinity at shelfice-ocean interface | 
| 515 |  |  | C     note: this part of the code is not very intuitive as it involves | 
| 516 |  |  | C     many arbitrary abbreviations that were introduced to derive the | 
| 517 |  |  | C     correct form of the quadratic equation for salinity. The abbreviations | 
| 518 |  |  | C     only make sense in connection with my notes on this (M.Losch) | 
| 519 |  |  | C | 
| 520 |  |  | C     eps3a was introduced as a constant variant of eps3 to avoid AD of | 
| 521 |  |  | C     code of typ (pLoc-const)/pLoc | 
| 522 |  |  | eps3a = rhoShelfIce*SHELFICEheatCapacity_Cp | 
| 523 |  |  | &           * SHELFICEkappa *  ( 1. _d 0 - dFac ) | 
| 524 |  |  | eps3 = eps3a/pLoc(I,J) | 
| 525 |  |  | eps4 = b*pLoc(I,J) + c0 | 
| 526 |  |  | eps6 = eps4 - tLoc(I,J) | 
| 527 |  |  | eps7 = eps4 - SHELFICEthetaSurface | 
| 528 |  |  | eps8 = rUnit2mass*SHELFICEheatCapacity_Cp | 
| 529 |  |  | &           *shiTransCoeffS(i,j,bi,bj) * fwflxFac | 
| 530 |  |  | aqe = a0  *(eps1+eps3-eps8) | 
| 531 |  |  | recip_aqe = 0. _d 0 | 
| 532 |  |  | IF ( aqe .NE. 0. _d 0 ) recip_aqe = 0.5 _d 0/aqe | 
| 533 |  |  | c           bqe = eps1*eps6 + eps3*eps7 - eps2 | 
| 534 |  |  | bqe = eps1*eps6 | 
| 535 |  |  | &           + eps3a*( b | 
| 536 |  |  | &                   + ( c0 - SHELFICEthetaSurface )/pLoc(I,J) ) | 
| 537 |  |  | &           - eps2 | 
| 538 |  |  | &           + eps8*( a0*sLoc(I,J) - eps7 ) | 
| 539 |  |  | cqe = ( eps2 + eps8*eps7 )*sLoc(I,J) | 
| 540 |  |  | discrim = bqe*bqe - 4. _d 0*aqe*cqe | 
| 541 |  |  | #undef ALLOW_SHELFICE_DEBUG | 
| 542 |  |  | #ifdef ALLOW_SHELFICE_DEBUG | 
| 543 |  |  | IF ( discrim .LT. 0. _d 0 ) THEN | 
| 544 |  |  | print *, 'ml-shelfice: discrim = ', discrim,aqe,bqe,cqe | 
| 545 |  |  | print *, 'ml-shelfice: pLoc    = ', pLoc(I,J) | 
| 546 |  |  | print *, 'ml-shelfice: tLoc    = ', tLoc(I,J) | 
| 547 |  |  | print *, 'ml-shelfice: sLoc    = ', sLoc(I,J) | 
| 548 |  |  | print *, 'ml-shelfice: tsurface= ', | 
| 549 |  |  | &            SHELFICEthetaSurface | 
| 550 |  |  | print *, 'ml-shelfice: eps1    = ', eps1 | 
| 551 |  |  | print *, 'ml-shelfice: eps2    = ', eps2 | 
| 552 |  |  | print *, 'ml-shelfice: eps3    = ', eps3 | 
| 553 |  |  | print *, 'ml-shelfice: eps4    = ', eps4 | 
| 554 |  |  | print *, 'ml-shelfice: eps5    = ', eps5 | 
| 555 |  |  | print *, 'ml-shelfice: eps6    = ', eps6 | 
| 556 |  |  | print *, 'ml-shelfice: eps7    = ', eps7 | 
| 557 |  |  | print *, 'ml-shelfice: eps8    = ', eps8 | 
| 558 |  |  | print *, 'ml-shelfice: rU2mass = ', rUnit2mass | 
| 559 |  |  | print *, 'ml-shelfice: rhoIce  = ', rhoShelfIce | 
| 560 |  |  | print *, 'ml-shelfice: cFac    = ', cFac | 
| 561 |  |  | print *, 'ml-shelfice: Cp_W    = ', HeatCapacity_Cp | 
| 562 |  |  | print *, 'ml-shelfice: Cp_I    = ', | 
| 563 |  |  | &            SHELFICEHeatCapacity_Cp | 
| 564 |  |  | print *, 'ml-shelfice: gammaT  = ', | 
| 565 |  |  | &            SHELFICEheatTransCoeff | 
| 566 |  |  | print *, 'ml-shelfice: gammaS  = ', | 
| 567 |  |  | &            SHELFICEsaltTransCoeff | 
| 568 |  |  | print *, 'ml-shelfice: lat.heat= ', | 
| 569 |  |  | &            SHELFICElatentHeat | 
| 570 |  |  | STOP 'ABNORMAL END in S/R SHELFICE_THERMODYNAMICS' | 
| 571 |  |  | ENDIF | 
| 572 |  |  | #endif /* ALLOW_SHELFICE_DEBUG */ | 
| 573 |  |  | saltFreeze = (- bqe - SQRT(discrim))*recip_aqe | 
| 574 |  |  | IF ( saltFreeze .LT. 0. _d 0 ) | 
| 575 |  |  | &           saltFreeze = (- bqe + SQRT(discrim))*recip_aqe | 
| 576 |  |  | thetaFreeze = a0*saltFreeze + eps4 | 
| 577 |  |  | C--   upward fresh water flux due to melting (in kg/m^2/s) | 
| 578 |  |  | cph change to identical form | 
| 579 |  |  | cph            freshWaterFlux = rUnit2mass | 
| 580 |  |  | cph     &           * shiTransCoeffS(i,j,bi,bj) | 
| 581 |  |  | cph     &           * ( saltFreeze - sLoc(I,J) ) / saltFreeze | 
| 582 |  |  | freshWaterFlux = rUnit2mass | 
| 583 |  |  | &           * shiTransCoeffS(i,j,bi,bj) | 
| 584 |  |  | &           * ( 1. _d 0 - sLoc(I,J) / saltFreeze ) | 
| 585 |  |  | #ifdef ALLOW_SHIFWFLX_CONTROL | 
| 586 |  |  | &           + xx_shifwflx_loc(I,J,bi,bj) | 
| 587 |  |  | #endif /*  ALLOW_SHIFWFLX_CONTROL */ | 
| 588 |  |  | C--   Calculate the upward heat and fresh water fluxes; | 
| 589 |  |  | C--   MITgcm sign conventions: downward (negative) fresh water flux | 
| 590 |  |  | C--   implies melting and due to upward (positive) heat flux | 
| 591 |  |  | shelfIceHeatFlux(I,J,bi,bj) = | 
| 592 |  |  | &           ( eps3 | 
| 593 |  |  | &           - freshWaterFlux*SHELFICEheatCapacity_Cp*fwflxFac ) | 
| 594 |  |  | &           * ( thetaFreeze - SHELFICEthetaSurface ) | 
| 595 |  |  | &           -  cFac*freshWaterFlux*( SHELFICElatentHeat | 
| 596 |  |  | &             - HeatCapacity_Cp*( thetaFreeze - rFac*tLoc(I,J) ) ) | 
| 597 |  |  | shelfIceFreshWaterFlux(I,J,bi,bj) = freshWaterFlux | 
| 598 |  |  | C--   compute surface tendencies | 
| 599 |  |  | shelficeForcingT(i,j,bi,bj) = | 
| 600 |  |  | &           ( shiTransCoeffT(i,j,bi,bj) | 
| 601 |  |  | &           - cFac*shelfIceFreshWaterFlux(I,J,bi,bj)*mass2rUnit ) | 
| 602 |  |  | &           * ( thetaFreeze - tLoc(I,J) ) | 
| 603 |  |  | shelficeForcingS(i,j,bi,bj) = | 
| 604 |  |  | &           ( shiTransCoeffS(i,j,bi,bj) | 
| 605 |  |  | &           - cFac*shelfIceFreshWaterFlux(I,J,bi,bj)*mass2rUnit ) | 
| 606 |  |  | &           * ( saltFreeze - sLoc(I,J) ) | 
| 607 |  |  | ELSE | 
| 608 |  |  | shelfIceHeatFlux      (I,J,bi,bj) = 0. _d 0 | 
| 609 |  |  | shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0 | 
| 610 |  |  | shelficeForcingT      (I,J,bi,bj) = 0. _d 0 | 
| 611 |  |  | shelficeForcingS      (I,J,bi,bj) = 0. _d 0 | 
| 612 |  |  | ENDIF | 
| 613 |  |  | ENDDO | 
| 614 |  |  | ENDDO | 
| 615 |  |  | ENDIF | 
| 616 |  |  | C     endif (not) useISOMIPTD | 
| 617 |  |  | ENDDO | 
| 618 |  |  | ENDDO | 
| 619 |  |  |  | 
| 620 |  |  | IF (SHELFICEMassStepping) THEN | 
| 621 |  |  | CALL SHELFICE_STEP_ICEMASS( myTime, myIter, myThid ) | 
| 622 |  |  | ENDIF | 
| 623 |  |  |  | 
| 624 |  |  | C--  Calculate new loading anomaly (in case the ice-shelf mass was updated) | 
| 625 |  |  | #ifndef ALLOW_AUTODIFF | 
| 626 |  |  | c     IF ( SHELFICEloadAnomalyFile .EQ. ' ' ) THEN | 
| 627 |  |  | DO bj = myByLo(myThid), myByHi(myThid) | 
| 628 |  |  | DO bi = myBxLo(myThid), myBxHi(myThid) | 
| 629 |  |  | DO j = 1-OLy, sNy+OLy | 
| 630 |  |  | DO i = 1-OLx, sNx+OLx | 
| 631 |  |  | shelficeLoadAnomaly(i,j,bi,bj) = gravity | 
| 632 |  |  | &      *( shelficeMass(i,j,bi,bj) + rhoConst*Ro_surf(i,j,bi,bj) ) | 
| 633 |  |  | ENDDO | 
| 634 |  |  | ENDDO | 
| 635 |  |  | ENDDO | 
| 636 |  |  | ENDDO | 
| 637 |  |  | c     ENDIF | 
| 638 |  |  | #endif /* ndef ALLOW_AUTODIFF */ | 
| 639 |  |  |  | 
| 640 |  |  | #ifdef ALLOW_DIAGNOSTICS | 
| 641 |  |  | IF ( useDiagnostics ) THEN | 
| 642 |  |  | CALL DIAGNOSTICS_FILL_RS(shelfIceFreshWaterFlux,'SHIfwFlx', | 
| 643 |  |  | &      0,1,0,1,1,myThid) | 
| 644 |  |  | CALL DIAGNOSTICS_FILL_RS(shelfIceHeatFlux,      'SHIhtFlx', | 
| 645 |  |  | &      0,1,0,1,1,myThid) | 
| 646 |  |  | C     SHIForcT (Ice shelf forcing for theta [W/m2], >0 increases theta) | 
| 647 |  |  | tmpFac = HeatCapacity_Cp*rUnit2mass | 
| 648 |  |  | CALL DIAGNOSTICS_SCALE_FILL(shelficeForcingT,tmpFac,1, | 
| 649 |  |  | &      'SHIForcT',0,1,0,1,1,myThid) | 
| 650 |  |  | C     SHIForcS (Ice shelf forcing for salt [g/m2/s], >0 increases salt) | 
| 651 |  |  | tmpFac = rUnit2mass | 
| 652 |  |  | CALL DIAGNOSTICS_SCALE_FILL(shelficeForcingS,tmpFac,1, | 
| 653 |  |  | &      'SHIForcS',0,1,0,1,1,myThid) | 
| 654 |  |  | C     Transfer coefficients | 
| 655 |  |  | CALL DIAGNOSTICS_FILL(shiTransCoeffT,'SHIgammT', | 
| 656 |  |  | &      0,1,0,1,1,myThid) | 
| 657 |  |  | CALL DIAGNOSTICS_FILL(shiTransCoeffS,'SHIgammS', | 
| 658 |  |  | &      0,1,0,1,1,myThid) | 
| 659 |  |  | C     Friction velocity | 
| 660 |  |  | #ifdef SHI_ALLOW_GAMMAFRICT | 
| 661 |  |  | IF ( SHELFICEuseGammaFrict ) | 
| 662 |  |  | &  CALL DIAGNOSTICS_FILL(uStarDiag,'SHIuStar',0,1,0,1,1,myThid) | 
| 663 |  |  | #endif /* SHI_ALLOW_GAMMAFRICT */ | 
| 664 |  |  | ENDIF | 
| 665 |  |  | #endif /* ALLOW_DIAGNOSTICS */ | 
| 666 |  |  |  | 
| 667 |  |  | #endif /* ALLOW_SHELFICE */ | 
| 668 |  |  |  | 
| 669 |  |  | RETURN | 
| 670 |  |  |  | 
| 671 |  |  | END |