pkg/steep_icecavity/stic_thermodynamics.F

00f81e6785 Ou W*0001 #include "STIC_OPTIONS.h"
                0002 #ifdef ALLOW_SHELFICE
                0003 # include "SHELFICE_OPTIONS.h"
                0004 #endif
                0005 #ifdef ALLOW_AUTODIFF
                0006 # include "AUTODIFF_OPTIONS.h"
                0007 #endif
                0008 #ifdef ALLOW_CTRL
                0009 # include "CTRL_OPTIONS.h"
                0010 #endif
                0011
                0012 CBOP
                0013 C     !ROUTINE: STIC_THERMODYNAMICS
                0014 C     !INTERFACE:
                0015       SUBROUTINE STIC_THERMODYNAMICS(
                0016      I                        myTime, myIter, myThid )
                0017 C     !DESCRIPTION: \bv
                0018 C     *=============================================================*
                0019 C     | S/R  STIC_THERMODYNAMICS
                0020 C     | o shelf-ice main routine.
                0021 C     |   compute temperature and (virtual) salt flux at the
                0022 C     |   shelf-ice ocean interface
                0023 C     |
                0024 C     | stresses at the ice/water interface are computed in separate
                0025 C     | routines that are called from mom_fluxform/mom_vecinv
54f5c8e380 Jean*0026 C     |
00f81e6785 Ou W*0027 CIGF  | ASSUMES
                0028 C---  |   * SHELFICEconserve = true
                0029 C     *=============================================================*
                0030 C     \ev
                0031
                0032 C     !USES:
                0033       IMPLICIT NONE
                0034
                0035 C     === Global variables ===
                0036 #include "SIZE.h"
                0037 #include "EEPARAMS.h"
                0038 #include "PARAMS.h"
                0039 #include "EOS.h"
                0040 #include "GRID.h"
                0041 #include "DYNVARS.h"
                0042 #include "FFIELDS.h"
                0043 #ifdef ALLOW_SHELFICE
                0044 # include "SHELFICE.h"
                0045 # include "SHELFICE_COST.h"
                0046 #endif
                0047 #include "STIC.h"
                0048 #ifdef ALLOW_CTRL
                0049 # include "CTRL_SIZE.h"
                0050 # include "CTRL.h"
                0051 # include "CTRL_DUMMY.h"
                0052 # ifdef ALLOW_GENTIM2D_CONTROL
                0053 #  include "CTRL_GENARR.h"
                0054 # endif
                0055 #endif /* ALLOW_CTRL */
                0056 #ifdef ALLOW_AUTODIFF_TAMC
                0057 # include "tamc.h"
                0058 #endif /* ALLOW_AUTODIFF_TAMC */
                0059
                0060 C     !INPUT/OUTPUT PARAMETERS:
                0061 C     === Routine arguments ===
                0062 C     myTime :: time counter for this thread
54f5c8e380 Jean*0063 C     myIter :: iteration counter for this thread
00f81e6785 Ou W*0064 C     myThid :: thread number for this instance of the routine.
                0065       _RL  myTime
                0066       INTEGER myIter
                0067       INTEGER myThid
                0068
                0069 C     !LOCAL VARIABLES :
                0070 C     === Local variables ===
                0071 C     i,j,k,bi,bj      :: loop counters
                0072 C     tLoc, sLoc, pLoc :: local potential temperature, salinity, pressure
                0073 C     gammaT/SLoc      :: local heat and salt transfer coefficients
                0074 C     theta/saltFreeze :: temperature and salinity of water at the
                0075 C                         ice-ocean interface (at the freezing point)
                0076 C     iceFrontCellThickness   :: the ratio of the grid cell area to
                0077 C                         the horizontal length of the ice front.
                0078 C                         unit meters.  Approximately the length of the
                0079 C                         column perpendicular to the ice front extended
                0080 C                         to the far side of the tracer cell.
                0081 C     iceFrontWidth    :: the width of the ice front.  unit meters.
                0082
54f5c8e380 Jean*0083       INTEGER i, j, k
                0084       INTEGER bi, bj
00f81e6785 Ou W*0085       INTEGER CURI, CURJ, FRONT_K
                0086
                0087       _RL tLoc
                0088       _RL sLoc
                0089       _RL pLoc
                0090       _RL gammaTLoc
                0091       _RL gammaSLoc
                0092
                0093       _RL ice_bottom_Z_C
                0094       _RL wet_top_Z_N, wet_bottom_Z_N
                0095       _RL iceFrontWetContact_Z_max
                0096       _RL iceFrontVertContactFrac, iceFrontCellThickness
                0097       _RL iceFrontWidth, iceFrontFaceArea
                0098       _RL thermalConductionDistance, thermalConductionTemp
                0099       _RL tmpHeatFlux, tmpFWFLX
                0100       _RL tmpForcingT, tmpForcingS
                0101       _RL tmpFac, icfgridareaFrac
                0102       _RL tmpHeatFluxscaled, tmpFWFLXscaled
                0103       INTEGER SI
                0104       _RL insituT
                0105
                0106 #ifdef ALLOW_DIAGNOSTICS
                0107 # ifdef SHI_ALLOW_GAMMAFRICT
                0108       _RL uStarDiag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0109 # endif
                0110 #ifdef NONLIN_FRSURF
                0111       _RL tmpDiag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0112       _RL tmpDiagT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0113 #endif
                0114       _RL tmpDiagIcfForcingT(1-OLx:sNx+OLx,
                0115      & 1-OLy:sNy+OLy,Nr,nSx,nSy)
                0116       _RL tmpDiagIcfForcingS(1-OLx:sNx+OLx,
                0117      & 1-OLy:sNy+OLy,Nr,nSx,nSy)
                0118       _RL tmpDiagShelficeForcingT(1-OLx:sNx+OLx,
                0119      & 1-OLy:sNy+OLy,nSx,nSy)
                0120       _RL tmpDiagShelficeForcingS(1-OLx:sNx+OLx,
                0121      & 1-OLy:sNy+OLy,nSx,nSy)
                0122 #endif /* ALLOW_DIAGNOSTICS */
                0123
                0124       _RL epsilon_H
                0125       _RL stic_addMass(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0126
                0127 #ifdef ALLOW_GENTIM2D_CONTROL
                0128       INTEGER iarr
                0129       _RL xx_shifwflx_loc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0130 #endif
                0131
                0132 #ifdef ALLOW_AUTODIFF_TAMC
                0133 C     tkey :: tape key (depends on tiles)
                0134 C     ikey :: tape key (depends on ij-indices and tiles)
                0135 C     kkey :: tape key (depends on levels, ij-indices, and tiles)
                0136       INTEGER tkey, ikey, kkey
                0137 #endif
                0138
                0139 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0140
                0141 #ifdef ALLOW_AUTODIFF_TAMC
                0142 CADJ STORE kTopC          = comlev1, key = ikey_dynamics
                0143 #endif
                0144
                0145 C--   minimum fraction of a cell adjacent to an ice front that must be
                0146 C--   wet for exchange to happen
                0147       epsilon_H = 1. _d -03
                0148
                0149 C--   hard coded for now.
                0150       thermalConductionDistance = 100.0 _d 0
                0151       thermalConductionTemp     = -20.0 _d 0
                0152       icfgridareaFrac = 1.0 _d 0
                0153       insituT = 0. _d 0
                0154
                0155 C     heat flux into the ice shelf, default is diffusive flux
                0156 C     (Holland and Jenkins, 1999, eq.21)
                0157
                0158       DO bj = myByLo(myThid), myByHi(myThid)
                0159        DO bi = myBxLo(myThid), myBxHi(myThid)
                0160         DO j = 1-OLy,sNy+OLy
                0161          DO i = 1-OLx,sNx+OLx
                0162           shelfIceHeatFlux        (i,j,bi,bj) = 0. _d 0
                0163           shelfIceFreshWaterFlux  (i,j,bi,bj) = 0. _d 0
                0164           sticfHeatFlux           (i,j,bi,bj) = 0. _d 0
                0165           sticfFreshWaterFlux     (i,j,bi,bj) = 0. _d 0
                0166 #ifdef ALLOW_GENTIM2D_CONTROL
                0167           xx_shifwflx_loc         (i,j,bi,bj) = 0. _d 0
                0168 #endif /* ALLOW_GENTIM2D_CONTROL */
                0169 #ifndef ALLOW_SHITRANSCOEFF_3D
                0170           shiTransCoeffS(i,j,bi,bj) = SHELFICEsaltToHeatRatio *
                0171      &         shiTransCoeffT(i,j,bi,bj)
                0172 #endif
                0173          ENDDO
                0174         ENDDO
                0175         DO k = 1, Nr
                0176          DO j = 1-OLy,sNy+OLy
                0177           DO i = 1-OLx,sNx+OLx
                0178 #ifdef ALLOW_SHITRANSCOEFF_3D
                0179            shiTransCoeffS3d(i,j,k,bi,bj) = SHELFICEsaltToHeatRatio *
                0180      &          shiTransCoeffT3d(i,j,k,bi,bj)
                0181 #endif
                0182            icfHeatFlux(i,j,k,bi,bj)       = 0. _d 0
                0183            icfFreshWaterFlux(i,j,k,bi,bj) = 0. _d 0
                0184            stic_gT(i,j,k,bi,bj)           = 0. _d 0
                0185            stic_gS(i,j,k,bi,bj)           = 0. _d 0
                0186            stic_addMass(i,j,k,bi,bj)      = 0. _d 0
                0187           ENDDO
                0188          ENDDO
                0189         ENDDO
                0190 #ifdef ALLOW_DIAGNOSTICS
                0191         IF ( useDiagnostics ) THEN
                0192          DO j = 1-OLy,sNy+OLy
                0193           DO i = 1-OLx,sNx+OLx
                0194            tmpDiagShelficeForcingT(i,j,bi,bj) = 0. _d 0
                0195            tmpDiagShelficeForcingS(i,j,bi,bj) = 0. _d 0
                0196           ENDDO
                0197          ENDDO
                0198          DO k = 1, Nr
                0199           DO j = 1-OLy,sNy+OLy
                0200            DO i = 1-OLx,sNx+OLx
                0201             tmpDiagIcfForcingT(i,j,k,bi,bj) = 0. _d 0
                0202             tmpDiagIcfForcingS(i,j,k,bi,bj) = 0. _d 0
                0203 #ifdef NONLIN_FRSURF
                0204             tmpDiag           (i,j,k,bi,bj) = 0. _d 0
                0205             tmpDiagT          (i,j,k,bi,bj) = 0. _d 0
                0206 #endif
                0207            ENDDO
                0208           ENDDO
                0209          ENDDO
                0210         ENDIF
                0211 #endif /* ALLOW_DIAGNOSTICS */
                0212        ENDDO
                0213       ENDDO
                0214
                0215 #if (defined ALLOW_CTRL && defined ALLOW_GENTIM2D_CONTROL)
                0216       IF ( useCTRL ) THEN
                0217        DO iarr = 1, maxCtrlTim2D
                0218         IF (xx_gentim2d_file(iarr)(1:11).EQ.'xx_shifwflx') THEN
                0219          DO bj = myByLo(myThid),myByHi(myThid)
                0220           DO bi = myBxLo(myThid),myBxHi(myThid)
                0221            DO j = 1,sNy
                0222             DO i = 1,sNx
                0223              xx_shifwflx_loc(i,j,bi,bj)=xx_gentim2d(i,j,bi,bj,iarr)
                0224             ENDDO
                0225            ENDDO
                0226           ENDDO
                0227          ENDDO
                0228         ENDIF
                0229        ENDDO
                0230       ENDIF
                0231 #endif
                0232
                0233       DO bj = myByLo(myThid), myByHi(myThid)
                0234         DO bi = myBxLo(myThid), myBxHi(myThid)
                0235 #ifdef ALLOW_AUTODIFF_TAMC
                0236           tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0237 #endif
                0238 C--   First ice front then ice shelf.  Loop through each i,j point
                0239 C--   process ice fronts in k, then process ice shelf.
                0240           DO j = 1-OLy+1,sNy+OLy-1
                0241             DO i = 1-OLx+1,sNx+OLx-1
                0242
                0243 #ifdef ALLOW_AUTODIFF_TAMC
                0244               ikey = i + OLx + (j+OLy-1)*(sNx+2*OLx)
                0245      &             + (tkey-1)*(sNx+2*OLx)*(sNy+2*OLy)
                0246 #endif
                0247
                0248 C--   The K index where the ice front ends (0 if no ice front)
                0249               FRONT_K = kIcf(i,j,bi,bj)
                0250
                0251 C--   If there is an ice front at this (i,j) continue
                0252               IF (FRONT_K .GT. 0) THEN
                0253
                0254 C--   Loop through all depths where the ice front is fround
                0255                 DO k = 1, FRONT_K
                0256 C--   Loop around the four laterally neighboring cells of the ice front.
                0257 C--   If any neighboring points has wet volume in contact with the ice
                0258 C--   front at (i,j) then calculate ice-ocean exchanges.
                0259 C--   The four laterally neighboring point are at (CURI,CURJ)
                0260
                0261                   DO SI = 1,4
                0262                     CURI=CURI_ARR(i,j,bi,bj,SI)
                0263                     CURJ=CURJ_ARR(i,j,bi,bj,SI)
                0264                     iceFrontWidth=sticfWidth_arr(i,j,bi,bj,SI)
                0265
                0266 C--                 cell depth describes the average distance
                0267 C--                 perpendicular to the ice front fact
                0268
                0269                     iceFrontCellThickness = 0. _d 0
                0270                     IF(iceFrontWidth.NE.0. _d 0)
                0271      &                iceFrontCellThickness = RA(CURI,CURJ,bi,bj)
                0272      &                                  /iceFrontWidth
                0273                     iceFrontFaceArea  = drF(k)*iceFrontWidth
                0274
                0275 C--   First, make sure the adjacent point has at least some water in it.
                0276                     IF (_hFacC(CURI,CURJ,k,bi,bj) .GT. zeroRL) THEN
                0277
                0278 C--   we need to determine how much of the ice front is in contact with
                0279 C--   water in the neighboring grid cell at this depth level.
                0280
                0281 C--   1. Determine the top depth with water in the current cell
                0282 C--   2. Determine the top depth with water in the neighbor cell
                0283 C--   3. Determine the depth where water  gap between (1) and (2).
                0284 C--   4. If there is a gap then ice front is in contact with water in
                0285 C--      the neighboring cell
                0286
                0287 C--   ice_bottom_Z_C: the depth (m) of the bottom of the ice in the
                0288 C--               current cell.  Bounded between rF(k) and rF(k+1).
                0289 C--               * If the ice extends past the bottom of the cell then
                0290 C--                 ice_bottom_Z_C = rF(k+1)
                0291 C--               [rF(k) >= ice_bottom_Z_C >= rF(k+1)]  (rF is negative)
                0292                       ice_bottom_Z_C = max(rF(k+1),
                0293      &                  min(Ro_surf(i,j, bi,bj), rF(k)))
                0294
                0295 C--   wet_top_Z_N: the depth (m) of the bottom of the ice in the
                0296 C--              neighboring grid.  If the neighboring cell has ice in
                0297 C--              (in the form of a shelf or front) then wet_top_Z_N is
                0298 C--              the depth of this neighboring ice.
                0299 C--
                0300 C--              * If neighbor cell has no ice, then Ro_surf = 0 and
                0301 C--                wet_top_Z_N = rF(k)
                0302 C--              [rF(k) >= wet_top_Z_N >= rF(k+1)]     (rF is negative)
                0303
                0304                       wet_top_Z_N = max(rF(k+1),
                0305      &                 min(Ro_surf(CURI,CURJ, bi,bj), rF(k)))
                0306
                0307 C--   wet_bottom_Z_N: the depth (m) of the bottom of the wet part of the
                0308 C--              neighboring cell.  If the seafloor reaches into
                0309 C--              the grid cell then the bottom of the wet part of the
                0310 C--              grid cell is at the seafloor.
                0311 C--
                0312 C--              * If the seafloor is deeper than this grid cell then
                0313 C--                wet_bottom_Z = rF(k+1)
                0314 C--              * If the seafloor is shallower than this grid cell then
                0315 C--                wet_bottom_Z = rF(k)
                0316 C--              * If the seafloor reaches partly into this grid cell
                0317 C--                then wet_bottom_Z = R_low
                0318
                0319 C--              [rF(k) >= wet_bottom_Z >= rF(k+1)]     (rF is negative)
                0320
                0321                       wet_bottom_Z_N = min(rF(k),
                0322      &                  max(R_low(CURI,CURJ, bi,bj), rF(k+1)))
                0323
                0324 C--   iceFrontWetContact_Z_max:  The deepest point where the
                0325 C--              the ice front at (i,j) is in contact with water
                0326 C--              in the neighboring cell.  The shallower of
                0327 C--              wet_bottom_Z_N (seafloor depth of neighboring point) and
                0328 C--              ice_bottom_Z_C (bottom of ice front in this center cell).
                0329
                0330 C--              * wet_bottom_Z_N if the seafloor of the neighboring
                0331 C--                cell is shallower than the ice draft at (i,j).
                0332 C--              * ice_bottom_Z_C if the ice draft at (i,j) is shallower
                0333 C--                than the seafloor of the neighboring cell.
                0334
                0335                       IF (ice_bottom_Z_C .GT. wet_bottom_Z_N) THEN
                0336                         iceFrontWetContact_Z_max = ice_bottom_Z_C
                0337                       ELSE
                0338                         iceFrontWetContact_Z_max = wet_bottom_Z_N
                0339                       ENDIF
                0340
                0341 C--   The shallowest depth where the ice front at (i,j) is in contact
                0342 C--   with water in the neighboring cell.  If the neighboring cell has
                0343 C--   no ice draft then wet_top_Z_N = rF(k), the top of the cell.
                0344 C--   Otherwise, the shallowest depth where the ice front at (i,j) can
                0345 C--   be in in contact with water (not ice) in (CURI, CURJ)
                0346 C--   is wet_top_Z_N.
                0347
                0348 C--   the fraction of the grid cell height that has ice draft in contact
                0349 C--   with water in the neighboring cell.
                0350                       iceFrontVertContactFrac =
                0351      &                  (wet_top_Z_N - iceFrontWetContact_Z_max)/ drF(k)
                0352
                0353 C--   Only proceed if iceFrontVertContactFrac is > 0, the
                0354 C--   ice draft at (i,j)
                0355 C--   is in contact with some water in the neighboring grid cell.
                0356                       IF (iceFrontVertContactFrac .GT. epsilon_H) THEN
                0357                         tLoc = theta(CURI,CURJ,k,bi,bj)
                0358                         sLoc = MAX(salt(CURI,CURJ,k,bi,bj), zeroRL)
                0359
                0360 C--   use pressure at the halfway point between the top and bottom of
                0361 C--   points of the ice front where the ice front is in contact with
                0362 C--   open water.
                0363                         pLoc = 0.5 _d 0 * ABS(wet_top_Z_N +
                0364      &                    iceFrontWetContact_Z_max)
                0365 #ifdef ALLOW_SHITRANSCOEFF_3D
                0366                         gammaTLoc = shiTransCoeffT3d(CURI,CURJ,k,bi,bj)
                0367                         gammaSLoc = shiTransCoeffS3d(CURI,CURJ,k,bi,bj)
                0368 #else
                0369                         gammaTLoc = shiTransCoeffT(CURI,CURJ,bi,bj)
                0370                         gammaSLoc = shiTransCoeffS(CURI,CURJ,bi,bj)
                0371 #endif
                0372
                0373                         CALL STIC_SOLVE4FLUXES(
                0374      I                    tLoc, sLoc, pLoc,
                0375      I                    gammaTLoc, gammaSLoc,
                0376      I                    thermalConductionDistance,
                0377      I                    thermalConductionTemp,
                0378      O                    tmpHeatFlux, tmpFWFLX,
                0379      O                    tmpForcingT, tmpForcingS,
                0380      O                    insituT,
                0381      I                    bi, bj, myTime, myIter, myThid )
                0382 #ifdef ALLOW_AUTODIFF_TAMC
                0383                         kkey = k + (ikey-1)*Nr
                0384 CADJ STORE tmpForcingT, tmpForcingS = comlev1_stic_bibj_ijk, key = kkey
                0385 #endif
                0386 C--   fluxes and forcing must be scaled by iceFrontVertContactFract and
                0387 C--   iceFrontContactFrac some fraction of the heigth and width of the
                0388 C--   grid cell face may not ice in contact with water.
                0389
                0390 C     tmpHeatFlux and tmpFWFLX come as W/m^2 and kg/m^2/s respectively
                0391 C--   but these rates only apply to the
                0392 C--   fraction of the grid cell that has ice in contact with seawater.
                0393 C--   we must scale by iceFrontVertContactFrac to get to the average
                0394 C--   fluxes in this grid cell.
                0395 C--   We also further scale by ratio of vertical to horizontal grid
                0396 C--   cell area so when comparing ice-front flux to ice-shelf flux we
                0397 C--   can just times them by the same area, i.e. horizontal grid cell area.
                0398
                0399 C--   ratio of vertical area to horizontal grid cell area
                0400                         icfgridareaFrac =
                0401      &                   iceFrontFaceArea/RA(CURI,CURJ,bi,bj)
                0402
                0403 C--   In units W/m^2
                0404                         tmpHeatFluxscaled =
                0405      &                    tmpHeatFlux*iceFrontVertContactFrac
                0406      &                    *icfgridareaFrac
                0407                         icfHeatFlux(CURI,CURJ,k,bi,bj) =
                0408      &                    icfHeatFlux(CURI,CURJ,k,bi,bj) +
                0409      &                    tmpHeatFluxscaled
                0410
                0411 C     In units of kg/s/m^2
                0412                         tmpFWFLXscaled =
                0413      &                    tmpFWFLX*iceFrontVertContactFrac
                0414      &                    *icfgridareaFrac
                0415                         icfFreshWaterFlux(CURI,CURJ,k,bi,bj) =
                0416      &                    icfFreshWaterFlux(CURI,CURJ,k,bi,bj) +
                0417      &                    tmpFWFLXscaled
                0418
                0419 C ow - 06/29/2018
                0420 C ow - Verticallly sum up the 3D icefront heat and freshwater fluxes to
                0421 C ow -  compute the total flux for the water column. The shelfice fluxes,
                0422 C ow -  which are 2D, will be added later. NOTE that only
                0423 C ow -  ice-front melts below shelf-ice are included to be consistent
                0424 C ow -  with Rignot's data
                0425                   if(k.GE.kTopC(i,j,bi,bj))then
                0426                    if(RA(CURI,CURJ,bi,bj).NE.0. _d 0)then
                0427                         sticfHeatFlux(CURI,CURJ,bi,bj) =
                0428      &                   sticfHeatFlux(CURI,CURJ,bi,bj) +
                0429      &                   tmpHeatFluxscaled
                0430                         sticfFreshWaterFlux(CURI,CURJ,bi,bj) =
                0431      &                   sticfFreshWaterFlux(CURI,CURJ,bi,bj) +
                0432      &                   tmpFWFLXscaled
                0433                    endif
                0434                   endif
                0435 C     stic_g[T,S] are tendency contributions due to ice front melt,
                0436 C--   but these constributions only apply to the
                0437 C--   fraction of the grid cell that has ice in contact with seawater.
                0438 C--   we must scale by iceFrontVertContactFrac to get to the average
                0439 C--   fluxes in this grid cell.  We must also divide the by the length
                0440 C--   of the grid cell perpendicular to the face.
                0441
                0442                        IF (iceFrontCellThickness .NE. 0. _d 0) THEN
                0443 C     In units of K / s
                0444                         stic_gT(CURI,CURJ,k,bi,bj) =
                0445      &                    stic_gT(CURI,CURJ,k,bi,bj) +
                0446      &                    tmpForcingT/iceFrontCellThickness*
                0447      &                    iceFrontVertContactFrac*
                0448      &                    _recip_hFacC(CURI,CURJ,k,bi,bj)
                0449 C     In units of psu /s
                0450                         stic_gS(CURI,CURJ,k,bi,bj) =
                0451      &                    stic_gS(CURI,CURJ,k,bi,bj) +
                0452      &                    tmpForcingS/iceFrontCellThickness*
                0453      &                    iceFrontVertContactFrac*
                0454      &                    _recip_hFacC(CURI,CURJ,k,bi,bj)
54f5c8e380 Jean*0455 #ifdef ALLOW_DIAGNOSTICS
                0456                         IF ( useDiagnostics ) THEN
                0457                          tmpDiagIcfForcingT(CURI,CURJ,k,bi,bj) =
                0458      &                     tmpDiagIcfForcingT(CURI,CURJ,k,bi,bj) +
                0459      &                     tmpForcingT/iceFrontCellThickness*
                0460      &                     iceFrontVertContactFrac*
                0461      &                     drF(k)
                0462                          tmpDiagIcfForcingS(CURI,CURJ,k,bi,bj) =
                0463      &                     tmpDiagIcfForcingS(CURI,CURJ,k,bi,bj) +
                0464      &                     tmpForcingS/iceFrontCellThickness*
                0465      &                     iceFrontVertContactFrac*
                0466      &                     drF(k)
                0467                         ENDIF
                0468 #endif /* ALLOW_DIAGNOSTICS */
00f81e6785 Ou W*0469                        ENDIF /* iceFrontCellThickness */
                0470 C     In units of kg /s
                0471                          stic_addMass(CURI,CURJ,k,bi,bj) =
                0472      &                     stic_addMass(CURI,CURJ,k,bi,bj) -
                0473      &                     tmpFWFLX*iceFrontFaceArea*
                0474      &                     iceFrontVertContactFrac
                0475 #ifdef ALLOW_DIAGNOSTICS
                0476                        IF ( useDiagnostics ) THEN
                0477 #ifdef NONLIN_FRSURF
                0478                         IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
                0479      &                   .AND. useRealFreshWaterFlux ) THEN
                0480                          tmpDiagT(CURI,CURJ,k,bi,bj) =
                0481      &                     stic_addMass(CURI,CURJ,k, bi,bj)/
                0482      &                     RA(CURI,CURJ,bi,bj)*
                0483      &                     insituT*HeatCapacity_Cp
                0484                         ENDIF
                0485 #endif /* NONLIN_FRSURF */
                0486                        ENDIF
                0487 #endif /* ALLOW_DIAGNOSTICS */
                0488                       ENDIF /* iceFrontVertContactFrac */
                0489                     ENDIF /* hFacC(CURI,CURJ,k,bi,bj) */
                0490                   ENDDO /* SI loop for adjacent cells */
                0491                 ENDDO /* k Loop */
                0492               ENDIF /* FRONT_K */
                0493
                0494 C--   ice shelf
                0495               k = kTopC(i,j,bi,bj)
                0496
                0497 C--   If there is an ice front at this (i,j) continue
                0498 C--   I am assuming K is only .GT. when there is at least some
                0499 C--   nonzero wet point below the shelf in the grid cell.
                0500               IF (k .GT. 0) THEN
                0501 C--   Initialize these values to zero
                0502                 pLoc = 0. _d 0
                0503                 tLoc = 0. _d 0
                0504                 sLoc = 0. _d 0
                0505                 gammaTLoc = 0. _d 0
                0506                 gammaSLoc = 0. _d 0
                0507
                0508 C--   make local copies of temperature, salinity and depth
                0509 C--   (pressure in deci-bar) underneath the ice
                0510 C--   for the ice shelf case we use hydrostatic pressure at the ice
                0511 C--   base of the ice shelf, top of the cavity.
                0512
                0513                 pLoc = ABS(R_shelfIce(i,j,bi,bj))
                0514                 tLoc = theta(i,j,k,bi,bj)
                0515                 sLoc = MAX(salt(i,j,k,bi,bj), zeroRL)
                0516 #ifdef ALLOW_SHITRANSCOEFF_3D
                0517                 gammaTLoc = shiTransCoeffT3d(i,j,k,bi,bj)
                0518                 gammaSLoc = shiTransCoeffS3d(i,j,k,bi,bj)
                0519 #else
                0520                 gammaTLoc = shiTransCoeffT(i,j,bi,bj)
                0521                 gammaSLoc = shiTransCoeffS(i,j,bi,bj)
                0522 #endif
                0523                 CALL STIC_SOLVE4FLUXES(
                0524      I            tLoc, sLoc, pLoc,
                0525      I            gammaTLoc, gammaSLoc,
                0526      I            pLoc, thermalConductionTemp,
                0527      O            tmpHeatFlux, tmpFWFLX,
                0528      O            tmpForcingT, tmpForcingS,
                0529      O            insituT,
                0530      I            bi, bj, myTime, myIter, myThid )
                0531 C     In units of W/m^2
                0532                 shelficeHeatFlux(i,j,bi,bj) = tmpHeatFlux
                0533 #ifdef ALLOW_GENTIM2D_CONTROL
                0534      &           - xx_shifwflx_loc(i,j,bi,bj)*SHELFICElatentHeat
                0535 #endif /*  ALLOW_GENTIM2D_CONTROL */
                0536 C     In units of kg/m^2/s
                0537                 shelfIceFreshWaterFlux(i,j,bi,bj) = tmpFWFLX
                0538 #ifdef ALLOW_GENTIM2D_CONTROL
                0539      &           + xx_shifwflx_loc(i,j,bi,bj)
                0540 #endif /*  ALLOW_GENTIM2D_CONTROL */
                0541 C ow - 06/29/2018
                0542 C ow - Now add shelfice heat and freshwater fluxes
                0543                         sticfHeatFlux(i,j,bi,bj) =
                0544      &                   sticfHeatFlux(i,j,bi,bj) +
                0545      &                   shelficeHeatFlux(i,j,bi,bj)
                0546                         sticfFreshWaterFlux(i,j,bi,bj) =
                0547      &                   sticfFreshWaterFlux(i,j,bi,bj) +
                0548      &                   shelfIceFreshWaterFlux(i,j,bi,bj)
                0549 C     In units of K/s
                0550                 stic_gT(i,j,k,bi,bj) = stic_gT(i,j,k,bi,bj) +
                0551      &              tmpForcingT*recip_drF(k)* _recip_hFacC(i,j,k,bi,bj)
                0552 C     In units of psu/s
                0553                 stic_gS(i,j,k,bi,bj) = stic_gS(i,j,k,bi,bj) +
                0554      &              tmpForcingS*recip_drF(k)* _recip_hFacC(i,j,k,bi,bj)
                0555 C     In units of kg/s  -- multiplication of area required first
                0556                 stic_addMass(i,j,k,bi,bj) = stic_addMass(i,j,k,bi,bj) -
                0557      &              tmpFWFLX*RA(i,j,bi,bj)
                0558 #ifdef ALLOW_DIAGNOSTICS
                0559                 IF ( useDiagnostics ) THEN
54f5c8e380 Jean*0560                  tmpDiagShelficeForcingT(i,j,bi,bj) = tmpForcingT
                0561                  tmpDiagShelficeForcingS(i,j,bi,bj) = tmpForcingS
00f81e6785 Ou W*0562 #ifdef NONLIN_FRSURF
                0563                  IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
                0564      &            .AND. useRealFreshWaterFlux ) THEN
                0565                   tmpDiagT(i,j,k,bi,bj) =
                0566      &                 stic_addMass(i,j,k,bi,bj) / rA(i,j,bi,bj)*
                0567      &                 insituT*HeatCapacity_Cp
                0568                  ENDIF
                0569 #endif /* NONLIN_FRSURF */
                0570                 ENDIF
                0571 #endif /* ALLOW_DIAGNOSTICS */
                0572               ENDIF /* Shelf k > 0 */
                0573             ENDDO /* i */
                0574           ENDDO /* j */
                0575         ENDDO /* bi */
                0576       ENDDO /* bj */
                0577
                0578       DO bj = myByLo(myThid), myByHi(myThid)
                0579         DO bi = myBxLo(myThid), myBxHi(myThid)
                0580           DO j = 1-OLy+1,sNy+OLy-1
                0581             DO i = 1-OLx+1,sNx+OLx-1
                0582               DO k = 1,Nr
                0583                  addMass(i,j,k,bi,bj) = addMass(i,j,k,bi,bj) +
                0584      &               stic_addMass(i,j,k,bi,bj)
                0585               ENDDO /* k */
                0586             ENDDO /* i */
                0587           ENDDO /* j */
                0588         ENDDO /* bi */
                0589       ENDDO /* bj */
                0590
                0591 #ifdef ALLOW_ADDFLUID
                0592       IF ( selectAddFluid.NE.0 .AND. salt_addMass.NE.UNSET_RL ) THEN
                0593         IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
                0594      &       .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
                0595           DO bj = myByLo(myThid), myByHi(myThid)
                0596             DO bi = myBxLo(myThid), myBxHi(myThid)
                0597               DO j=0,sNy+1
                0598                 DO i=0,sNx+1
                0599                   DO k=1,Nr
                0600                     stic_gT(i,j,k,bi,bj) = stic_gT(i,j,k,bi,bj)
                0601      &                - stic_addMass(i,j,k,bi,bj)*mass2rUnit
                0602      &                  *( temp_addMass - theta(i,j,k,bi,bj) )
                0603      &                  *recip_rA(i,j,bi,bj)
                0604      &                  *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0605 C    &                  *recip_deepFac2C(k)*recip_rhoFacC(k)
                0606                     stic_gS(i,j,k,bi,bj) = stic_gS(i,j,k,bi,bj)
                0607      &                - stic_addMass(i,j,k,bi,bj)*mass2rUnit
                0608      &                  *( salt_addMass - salt(i,j,k,bi,bj) )
                0609      &                  *recip_rA(i,j,bi,bj)
                0610      &                  *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0611 C    &                  *recip_deepFac2C(k)*recip_rhoFacC(k)
                0612                   ENDDO /* k */
                0613                 ENDDO /* i */
                0614               ENDDO /* j */
                0615             ENDDO /* bi */
                0616           ENDDO /* bj */
                0617         ELSE
                0618           DO bj = myByLo(myThid), myByHi(myThid)
                0619             DO bi = myBxLo(myThid), myBxHi(myThid)
                0620               DO j=0,sNy+1
                0621                 DO i=0,sNx+1
                0622                   DO k=1,Nr
                0623                     stic_gT(i,j,k,bi,bj) = stic_gT(i,j,k,bi,bj)
                0624      &                - stic_addMass(i,j,k,bi,bj)*mass2rUnit
                0625      &                  *( temp_addMass - tRef(k) )
                0626      &                  *recip_rA(i,j,bi,bj)
                0627      &                  *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0628 C    &                  *recip_deepFac2C(k)*recip_rhoFacC(k)
                0629                     stic_gS(i,j,k,bi,bj) = stic_gS(i,j,k,bi,bj)
                0630      &                - stic_addMass(i,j,k,bi,bj)*mass2rUnit
                0631      &                  *( salt_addMass - sRef(k) )
                0632      &                  *recip_rA(i,j,bi,bj)
                0633      &                  *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0634 C    &                  *recip_deepFac2C(k)*recip_rhoFacC(k)
                0635                   ENDDO /* k */
                0636                 ENDDO /* i */
                0637               ENDDO /* j */
                0638             ENDDO /* bi */
                0639           ENDDO /* bj */
                0640         ENDIF
                0641       ENDIF
                0642 #endif /* ALLOW_ADDFLUID */
                0643
                0644 C--  Calculate new loading anomaly (in case the ice-shelf mass was updated)
                0645 #ifndef ALLOW_AUTODIFF
                0646 c     IF ( SHELFICEloadAnomalyFile .EQ. ' ' ) THEN
                0647        DO bj = myByLo(myThid), myByHi(myThid)
                0648         DO bi = myBxLo(myThid), myBxHi(myThid)
                0649          DO j = 1-OLy, sNy+OLy
                0650           DO i = 1-OLx, sNx+OLx
                0651            shelficeLoadAnomaly(i,j,bi,bj) = gravity
                0652      &      *( shelficeMass(i,j,bi,bj) + rhoConst*Ro_surf(i,j,bi,bj) )
                0653           ENDDO
                0654          ENDDO
                0655         ENDDO
                0656        ENDDO
                0657 c     ENDIF
                0658 #endif /* ndef ALLOW_AUTODIFF */
                0659
                0660 #ifdef ALLOW_DIAGNOSTICS
                0661       IF ( useDiagnostics ) THEN
54f5c8e380 Jean*0662
00f81e6785 Ou W*0663 #ifdef NONLIN_FRSURF
                0664        IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
                0665      &         .AND. useRealFreshWaterFlux ) THEN
                0666         DO bj = myByLo(myThid), myByHi(myThid)
                0667          DO bi = myBxLo(myThid), myBxHi(myThid)
                0668           DO k = 1, Nr
                0669            DO j=1,sNy
                0670             DO i=1,sNx
                0671              tmpDiag(i,j,k,bi,bj) =
                0672      &        stic_addMass(i,j,k,bi,bj)/RA(i,j,bi,bj)*
                0673      &        theta(i,j,k,bi,bj)*HeatCapacity_Cp
                0674             ENDDO
                0675            ENDDO
                0676           ENDDO
                0677          ENDDO /* bi */
                0678         ENDDO /* bj */
                0679         CALL DIAGNOSTICS_FILL( tmpDiag,'TFLUXMLT',0,Nr,0,1,1,myThid )
                0680
                0681         DO bj = myByLo(myThid), myByHi(myThid)
                0682          DO bi = myBxLo(myThid), myBxHi(myThid)
                0683           DO k = 1, Nr
                0684            DO j=1,sNy
                0685             DO i=1,sNx
                0686              tmpDiag(i,j,k,bi,bj) =
                0687      &        stic_addMass(i,j,k,bi,bj)/RA(i,j,bi,bj)*salt(i,j,k,bi,bj)
                0688             ENDDO
                0689            ENDDO
                0690           ENDDO
                0691          ENDDO /* bi */
                0692         ENDDO /* bj */
                0693         CALL DIAGNOSTICS_FILL( tmpDiag,'SFLUXMLT',0,Nr,0,1,1,myThid )
                0694
                0695         CALL DIAGNOSTICS_FILL( tmpDiagT,'TFLUXMTI',0,Nr,0,1,1,myThid )
                0696        ENDIF
                0697 #endif /* NONLIN_FRSURF */
                0698
                0699        CALL DIAGNOSTICS_FILL_RS(shelfIceFreshWaterFlux,'SHIfwFlx',
                0700      &      0,1,0,1,1,myThid)
                0701        CALL DIAGNOSTICS_FILL_RS(shelfIceHeatFlux,      'SHIhtFlx',
                0702      &      0,1,0,1,1,myThid)
                0703
                0704 #ifdef ALLOW_STEEP_ICECAVITY
54f5c8e380 Jean*0705        CALL DIAGNOSTICS_FILL( sticfFreshWaterFlux, 'SHIICFfw',
00f81e6785 Ou W*0706      &      0,1,0,1,1,myThid)
54f5c8e380 Jean*0707        CALL DIAGNOSTICS_FILL( sticfHeatFlux,       'SHIICFht',
00f81e6785 Ou W*0708      &      0,1,0,1,1,myThid)
                0709 #endif
                0710
                0711         CALL DIAGNOSTICS_FILL(icfFreshWaterFlux, 'STCfwFlx',
                0712      &      0,Nr,0,1,1,myThid)
                0713         CALL DIAGNOSTICS_FILL(icfHeatFlux, 'STChtFlx',
                0714      &      0,Nr,0,1,1,myThid)
                0715
                0716 C     SHIForcT (Ice shelf forcing for theta [W/m2], >0 increases theta)
                0717        tmpFac = HeatCapacity_Cp*rUnit2mass
                0718        CALL DIAGNOSTICS_SCALE_FILL(tmpDiagShelficeForcingT,tmpFac,1,
                0719      &      'SHIForcT',0,1,0,1,1,myThid)
                0720 C     SHIForcS (Ice shelf forcing for salt [g/m2/s], >0 increases salt)
                0721        tmpFac = rUnit2mass
                0722        CALL DIAGNOSTICS_SCALE_FILL(tmpDiagShelficeForcingS,tmpFac,1,
                0723      &      'SHIForcS',0,1,0,1,1,myThid)
                0724
                0725 C     STCForcT (Ice front forcing for theta [W/m2], >0 increases theta)
                0726        tmpFac = HeatCapacity_Cp*rUnit2mass
                0727        CALL DIAGNOSTICS_SCALE_FILL(tmpDiagIcfForcingT,tmpFac,1,
                0728      &      'STCForcT',0,Nr,0,1,1,myThid)
                0729 C     STCForcS (Ice front forcing for salt [g/m2/s], >0 increases salt)
                0730        tmpFac = rUnit2mass
                0731        CALL DIAGNOSTICS_SCALE_FILL(tmpDiagIcfForcingS,tmpFac,1,
                0732      &      'STCForcS',0,Nr,0,1,1,myThid)
                0733
                0734 C     Transfer coefficients
                0735 #ifdef ALLOW_SHITRANSCOEFF_3D
                0736        CALL DIAGNOSTICS_FILL(shiTransCoeffT3d,'SHIgam3T',
                0737      &      0,Nr,0,1,1,myThid)
                0738        CALL DIAGNOSTICS_FILL(shiTransCoeffS3d,'SHIgam3S',
                0739      &      0,Nr,0,1,1,myThid)
                0740 #else
                0741        CALL DIAGNOSTICS_FILL(shiTransCoeffT,'SHIgammT',
                0742      &      0,1,0,1,1,myThid)
                0743        CALL DIAGNOSTICS_FILL(shiTransCoeffS,'SHIgammS',
                0744      &      0,1,0,1,1,myThid)
                0745 #endif
                0746 C     Friction velocity
                0747 #ifdef SHI_ALLOW_GAMMAFRICT
                0748        IF ( SHELFICEuseGammaFrict )
                0749      &  CALL DIAGNOSTICS_FILL(uStarDiag,'SHIuStar',0,1,0,1,1,myThid)
                0750 #endif /* SHI_ALLOW_GAMMAFRICT */
                0751
                0752       ENDIF
54f5c8e380 Jean*0753 #endif /* ALLOW_DIAGNOSTICS */
00f81e6785 Ou W*0754
                0755       RETURN
                0756       END