pkg/steep_icecavity/stic_init_fixed.F

00f81e6785 Ou W*0001 #include "STIC_OPTIONS.h"
                0002 #ifdef ALLOW_SHELFICE
                0003 # include "SHELFICE_OPTIONS.h"
                0004 #endif
                0005
                0006 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0007 CBOP
                0008 C !ROUTINE: STIC_INIT_FIXED
                0009
                0010 C !INTERFACE:
                0011       SUBROUTINE STIC_INIT_FIXED( myThid )
                0012
                0013 C     !DESCRIPTION:
                0014 C     Initialize STIC variables that are kept fixed during the run.
                0015
                0016 C     !USES:
                0017       IMPLICIT NONE
                0018 C     == Global variables ===
                0019 #include "EEPARAMS.h"
                0020 #include "SIZE.h"
                0021 #include "PARAMS.h"
                0022 #include "GRID.h"
                0023 #include "STIC.h"
                0024 #ifdef ALLOW_SHELFICE
                0025 # include "SHELFICE.h"
                0026 #endif
                0027
                0028 C     !INPUT/OUTPUT PARAMETERS:
                0029 C     myThid ::  my Thread Id number
                0030       INTEGER myThid
                0031 CEOP
                0032
                0033 C     !LOCAL VARIABLES:
                0034 C     i,j,bi,bj - Loop counters
                0035       INTEGER i, j, k, bi, bj
                0036       INTEGER CURI, CURJ, FRONT_K
                0037 C local variabls used to determine shelf-ice and ice-front masks
                0038 C     iceFrontCellThickness   :: the ratio of the horizontal length
                0039 C                         of the ice front in each model grid cell
                0040 C                         divided by the grid cell area.  The "thickness"
                0041 C                         of the colum perpendicular to the front
                0042 C     iceFrontWidth    :: the width of the ice front.
                0043       _RL ice_bottom_Z_C
                0044       _RL wet_top_Z_N, wet_bottom_Z_N
                0045       _RL iceFrontWetContact_Z_max
                0046       _RL iceFrontVertContactFrac, iceFrontCellThickness
                0047       _RL iceFrontWidth, iceFrontFaceArea
                0048       _RS fK_icefront (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0049
                0050       INTEGER SI
                0051       _RL epsilon_H
                0052
                0053 #ifdef ALLOW_DIAGNOSTICS
                0054       INTEGER       diagNum
                0055       CHARACTER*8   diagName
                0056       CHARACTER*16  diagCode
                0057       CHARACTER*16  diagUnits
                0058       CHARACTER*(80) diagTitle
                0059 #endif /* ALLOW_DIAGNOSTICS */
                0060
                0061 C ow - 06/29/2018
                0062 C     pkg/shelfice maskSHI above is not consistent with the spirit of gencost.
                0063 C     Use mask2dSHI and mask3dSHI below instead.
                0064 C     mask2dSHI and mask3dSHI are the 2d and 3d mask for shelfice. They are
                0065 C     zero if there is no shelfice and one if otherwise. For any i,j,
                0066 C     if there is at least one non-zero mask3dSHI in the vertical, then
                0067 C     mask2dSHI at i,j is one.
                0068       DO bj = myByLo(myThid), myByHi(myThid)
                0069        DO bi = myBxLo(myThid), myBxHi(myThid)
                0070         DO k=1,Nr
                0071          DO j=1-OLy,sNy+OLy
                0072           DO i=1-OLx,sNx+OLx
                0073            mask3dSHIICF(i,j,k,bi,bj) = 0. _d 0
                0074            mask3dSHI(i,j,k,bi,bj) = 0. _d 0
                0075            mask3dICF(i,j,k,bi,bj) = 0. _d 0
                0076            IF (k.EQ.1) THEN
                0077             mask2dSHIICF(i,j,bi,bj) = 0. _d 0
                0078             mask2dSHI(i,j,bi,bj) = 0. _d 0
                0079             mask2dICF(i,j,bi,bj) = 0. _d 0
                0080            ENDIF
                0081           ENDDO
                0082          ENDDO
                0083         ENDDO
                0084
                0085         DO j = 1-OLy,sNy+OLy
                0086          DO i = 1-OLx,sNx+OLx
                0087           sticfLength(i,j,bi,bj) = 0. _d 0
                0088           DO SI = 1,4
                0089            CURI_ARR(i,j,bi,bj,SI) = -9999
                0090            CURJ_ARR(i,j,bi,bj,SI) = -9999
                0091            sticfWidth_arr(i,j,bi,bj,SI) = 0. _d 0
                0092           ENDDO
                0093          ENDDO
                0094         ENDDO
                0095
                0096        ENDDO
                0097       ENDDO
                0098
                0099 C--  STEEP_ICECAVITY parameters (BEGIN)
                0100       DO bj = myByLo(myThid), myByHi(myThid)
                0101        DO bi = myBxLo(myThid), myBxHi(myThid)
                0102         DO j = 1-OLy, sNy+OLy
                0103          DO i = 1-OLx, sNx+OLx
                0104           kIcf(i,j,bi,bj) = 0
                0105           DO k = 1 , Nr
                0106            IF ( R_stic(i,j,bi,bj) .GT. ABS(rF(k))) kIcf(i,j,bi,bj) = k
                0107           ENDDO
                0108           fK_icefront(i,j,bi,bj) = FLOAT(kIcf(i,j,bi,bj))
                0109          ENDDO
                0110         ENDDO
                0111        ENDDO
                0112       ENDDO
                0113 C--  STEEP_ICECAVITY parameters (END)
                0114
                0115 C--   Create masks for shelf-ice and ice-front by modifyig code from
                0116 C     shelfice_thermodynamics.F minimum fraction of a cell adjacent to an
                0117 C     ice front that must be wet for exchange to happen
                0118       epsilon_H = 1. _d -03
                0119
                0120       DO bj = myByLo(myThid), myByHi(myThid)
                0121        DO bi = myBxLo(myThid), myBxHi(myThid)
                0122
                0123 C--   First ice front then ice shelf.  Loop through each i,j point
                0124 C--   process ice fronts in k, then process ice shelf.
                0125         DO j = 1-OLy+1,sNy+OLy-1
                0126          DO i = 1-OLx+1,sNx+OLx-1
                0127
                0128 C--   The k index where the ice front ends (0 if no ice front)
                0129           FRONT_K = kIcf(i,j,bi,bj)
                0130
                0131 C--   If there is an ice front at this (i,j) continue
                0132           IF (FRONT_K .GT. 0) THEN
                0133
                0134 C--   Loop through all depths where the ice front is fround
                0135            DO k = 1, FRONT_K
                0136 C--   Loop around the four laterally neighboring cells of the ice front.
                0137 C--   If any neighboring points has wet volume in contact with the ice
                0138 C--   front at (i,j) then calculate ice-ocean exchanges.
                0139 C--   The four laterally neighboring point are at (CURI,CURJ)
                0140             DO SI = 1,4
                0141              IF     (SI .EQ. 1) THEN
                0142 C--   Looking to right
                0143               CURI = i+1
                0144               CURJ = j
                0145               iceFrontWidth = dyG(i+1,j,bi,bj)
                0146              ELSEIF (SI .EQ. 2) THEN
                0147 C--   Looking to LEFT
                0148               CURI = i-1
                0149               CURJ = j
                0150               iceFrontWidth = dyG(i,j,bi,bj)
                0151              ELSEIF (SI .EQ. 3) THEN
                0152 C--   Looking to NORTH
                0153               CURI = i
                0154               CURJ = j+1
                0155               iceFrontWidth = dxG(i,j+1,bi,bj)
                0156              ELSEIF (SI .EQ. 4) THEN
                0157 C--   Looking to south
                0158               CURI = i
                0159               CURJ = j-1
                0160               iceFrontWidth = dxG(i,j,bi,bj)
                0161              ENDIF
                0162
                0163              CURI_ARR(i,j,bi,bj,SI) = CURI
                0164              CURJ_ARR(i,j,bi,bj,SI) = CURJ
                0165              sticfWidth_arr(i,j,bi,bj,SI) = iceFrontWidth
                0166
                0167 C--   cell depth describes the average distance perpendicular to the ice
                0168 C--   front fact
                0169
                0170              iceFrontCellThickness = 0. _d 0
                0171              IF(iceFrontWidth.NE.0. _d 0)
                0172      &         iceFrontCellThickness = RA(CURI,CURJ,bi,bj)/iceFrontWidth
                0173              iceFrontFaceArea  = drF(k)*iceFrontWidth
                0174
                0175 C--   First, make sure the adjacent point has at least some water in it.
                0176              IF (_hFacC(CURI,CURJ,k,bi,bj) .GT. zeroRL) THEN
                0177
                0178 C--   we need to determine how much of the ice front is in contact with
                0179 C--   water in the neighboring grid cell at this depth level.
                0180
                0181 C--   1. Determine the top depth with water in the current cell
                0182 C--   2. Determine the top depth with water in the neighbor cell
                0183 C--   3. Determine the depth where water  gap between (1) and (2).
                0184 C--   4. If there is a gap then ice front is in contact with water in
                0185 C--      the neighboring cell
                0186
                0187 C--   ice_bottom_Z_C: the depth (m) of the bottom of the ice in the
                0188 C--               current cell.  Bounded between rF(k) and rF(k+1).
                0189 C--               * If the ice extends past the bottom of the cell then
                0190 C--                 ice_bottom_Z_C = rF(k+1)
                0191 C--               [rF(k) >= ice_bottom_Z_C >= rF(k+1)]  (rF is negative)
                0192               ice_bottom_Z_C = MAX(rF(k+1),
                0193      &             MIN(Ro_surf(i,j, bi,bj), rF(k)))
                0194
                0195 C--   wet_top_Z_N: the depth (m) of the bottom of the ice in the
                0196 C--              neighboring grid.  If the neighboring cell has ice in
                0197 C--              (in the form of a shelf or front) then wet_top_Z_N is
                0198 C--              the depth of this neighboring ice.
                0199 C--
                0200 C--              * If neighbor cell has no ice, then Ro_surf = 0 and
                0201 C--                wet_top_Z_N = rF(k)
                0202 C--              [rF(k) >= wet_top_Z_N >= rF(k+1)]     (rF is negative)
                0203
                0204               wet_top_Z_N = MAX(rF(k+1),
                0205      &             MIN(Ro_surf(CURI,CURJ, bi,bj), rF(k)))
                0206
                0207 C--   wet_bottom_Z_N: the depth (m) of the bottom of the wet part of the
                0208 C--              neighboring cell.  If the seafloor reaches into
                0209 C--              the grid cell then the bottom of the wet part of the
                0210 C--              grid cell is at the seafloor.
                0211 C--
                0212 C--              * If the seafloor is deeper than this grid cell then
                0213 C--                wet_bottom_Z = rF(k+1)
                0214 C--              * If the seafloor is shallower than this grid cell then
                0215 C--                wet_bottom_Z = rF(k)
                0216 C--              * If the seafloor reaches partly into this grid cell
                0217 C--                then wet_bottom_Z = R_low
                0218
                0219 C--              [rF(k) >= wet_bottom_Z >= rF(k+1)]     (rF is negative)
                0220
                0221               wet_bottom_Z_N = MIN(rF(k),
                0222      &             MAX(R_low(CURI,CURJ, bi,bj), rF(k+1)))
                0223
                0224 C--   iceFrontWetContact_Z_max:  The deepest point where the
                0225 C--              the ice front at (i,j) is in contact with water
                0226 C--              in the neighboring cell.  The shallower of
                0227 C--              wet_bottom_Z_N (seafloor depth of neighboring point) and
                0228 C--              ice_bottom_Z_C (bottom of ice front in this center cell).
                0229
                0230 C--              * wet_bottom_Z_N if the seafloor of the neighboring
                0231 C--                cell is shallower than the ice draft at (i,j).
                0232 C--              * ice_bottom_Z_C if the ice draft at (i,j) is shallower
                0233 C--                than the seafloor of the neighboring cell.
                0234
                0235               IF (ice_bottom_Z_C .GT. wet_bottom_Z_N) THEN
                0236                iceFrontWetContact_Z_max = ice_bottom_Z_C
                0237               ELSE
                0238                iceFrontWetContact_Z_max = wet_bottom_Z_N
                0239               ENDIF
                0240
                0241 C--   The shallowest depth where the ice front at (i,j) is in contact
                0242 C--   with water in the neighboring cell.  If the neighboring cell has
                0243 C--   no ice draft then wet_top_Z_N = rF(k), the top of the cell.
                0244 C--   Otherwise, the shallowest depth where the ice front at (i,j) can
                0245 C--   be in in contact with water (not ice) in (CURI, CURJ)
                0246 C--   is wet_top_Z_N.
                0247
                0248 C--   the fraction of the grid cell height that has ice draft in contact
                0249 C--   with water in the neighboring cell.
                0250               iceFrontVertContactFrac =
                0251      &             (wet_top_Z_N - iceFrontWetContact_Z_max)/ drF(k)
                0252 C--   Only proceed if iceFrontVertContactFrac is > 0, the
                0253 C--   ice draft at (i,j)
                0254 C--   is in contact with some water in the neighboring grid cell.
                0255               IF (iceFrontVertContactFrac .GT. epsilon_H) THEN
                0256                mask3dSHIICF(CURI,CURJ,k,bi,bj) = 1. _d 0
                0257                mask2dSHIICF(CURI,CURJ,  bi,bj) = 1. _d 0
                0258                mask3dICF   (CURI,CURJ,k,bi,bj) = 1. _d 0
                0259                mask2dICF   (CURI,CURJ,  bi,bj) = 1. _d 0
                0260               ENDIF /* iceFrontVertContactFrac */
                0261              ENDIF /* hFacC(CURI,CURJ,k,bi,bj) */
                0262             ENDDO /* SI loop for adjacent cells */
                0263            ENDDO /* k loop */
                0264           ENDIF /* FRONT_K */
                0265
                0266 C--   ice shelf
                0267           k = kTopC(i,j,bi,bj)
                0268
                0269 C--   If there is an ice front at this (i,j) continue
                0270 C--   I am assuming k is only .GT. when there is at least some
                0271 C--   nonzero wet point below the shelf in the grid cell.
                0272           IF (k .GT. 0) THEN
                0273            mask3dSHIICF(i,j,k,bi,bj) = 1. _d 0
                0274            mask2dSHIICF(i,j,  bi,bj) = 1. _d 0
                0275            mask3dSHI   (i,j,k,bi,bj) = 1. _d 0
                0276            mask2dSHI   (i,j,  bi,bj) = 1. _d 0
                0277           ENDIF /* SHELF k > 0 */
                0278          ENDDO /* i */
                0279         ENDDO /* j */
                0280        ENDDO /* bi */
                0281       ENDDO /* bj */
                0282
                0283 C     fill in the halos
                0284       _EXCH_XY_RS(  mask2dSHIICF, myThid )
                0285       _EXCH_XY_RS(  mask2dICF   , myThid )
                0286       _EXCH_XY_RS(  mask2dSHI   , myThid )
                0287       _EXCH_XYZ_RS( mask3dSHIICF, myThid )
                0288       _EXCH_XYZ_RS( mask3dICF   , myThid )
                0289       _EXCH_XYZ_RS( mask3dSHI   , myThid )
                0290
                0291 C     output the masks
                0292       CALL WRITE_FLD_XY_RS ( 'mask2dSHIICF',' ',mask2dSHIICF,-1,myThid)
                0293       CALL WRITE_FLD_XYZ_RS( 'mask3dSHIICF',' ',mask3dSHIICF, 0,myThid)
                0294       CALL WRITE_FLD_XY_RS ( 'mask2dSHI',' ',mask2dSHI,-1,myThid)
                0295       CALL WRITE_FLD_XYZ_RS( 'mask3dSHI',' ',mask3dSHI, 0,myThid)
                0296       CALL WRITE_FLD_XY_RS ( 'mask2dICF',' ',mask2dICF,-1,myThid)
                0297       CALL WRITE_FLD_XYZ_RS( 'mask3dICF',' ',mask3dICF, 0,myThid)
                0298       CALL WRITE_FLD_XY_RS ( 'R_stic',' ',R_stic,-1,myThid)
                0299       CALL WRITE_FLD_XY_RS ( 'kIcf',' ',fK_icefront,-1,myThid)
                0300
                0301 #ifdef ALLOW_DIAGNOSTICS
                0302       IF ( useDiagnostics ) THEN
                0303
                0304 # ifdef ALLOW_SHITRANSCOEFF_3D
                0305        diagName  = 'SHIgam3T'
                0306        diagTitle = '3D-Ice shelf exchange coefficient for theta'
                0307        diagUnits = 'm/s             '
                0308        diagCode  = 'SMR     MR      '
                0309        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0310      I      diagName, diagCode, diagUnits, diagTitle, 0, myThid )
                0311
                0312        diagName  = 'SHIgam3S'
                0313        diagTitle = '3D-Ice shelf exchange coefficient for salt'
                0314        diagUnits = 'm/s             '
                0315        diagCode  = 'SMR     MR      '
                0316        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0317      I      diagName, diagCode, diagUnits, diagTitle, 0, myThid )
                0318 # endif
                0319
                0320        diagName  = 'SHIICFfw'
                0321        diagTitle = 'total ice shelf and front FW flux (+ upward)'
                0322        diagUnits = 'kg/m^2/s        '
                0323        diagCode  = 'SM      L1      '
                0324        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0325      I      diagName, diagCode, diagUnits, diagTitle, 0, myThid )
                0326
                0327        diagName  = 'SHIICFht'
                0328        diagTitle = 'total ice shelf and ice front heat flux (+ upward)'
                0329        diagUnits = 'W/m^2           '
                0330        diagCode  = 'SM      L1      '
                0331        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0332      I      diagName, diagCode, diagUnits, diagTitle, 0, myThid )
                0333
                0334        diagName  = 'STCfwFlx'
                0335        diagTitle ='Ice front freshwater flux (+ve increases ocean salt)'
                0336        diagUnits = 'kg/m^2/s        '
                0337        diagCode  = 'SM      MR      '
                0338        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0339      I      diagName, diagCode, diagUnits, diagTitle, 0, myThid )
                0340
                0341        diagName  = 'STChtFlx'
                0342        diagTitle = 'Ice front heat flux  (+ve cools ocean)'
                0343        diagUnits = 'W/m^2           '
                0344        diagCode  = 'SM      MR      '
                0345        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0346      I      diagName, diagCode, diagUnits, diagTitle, 0, myThid )
                0347
                0348        diagName  = 'STCForcT'
                0349        diagTitle = 'Ice front forcing for theta, >0 increases theta'
                0350        diagUnits = 'W/m^2           '
                0351        diagCode  = 'SM      MR      '
                0352        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0353      I      diagName, diagCode, diagUnits, diagTitle, 0, myThid )
                0354
                0355        diagName  = 'STCForcS'
                0356        diagTitle = 'Ice front forcing for salt, >0 increases salt'
                0357        diagUnits = 'g/m^2/s         '
                0358        diagCode  = 'SM      MR      '
                0359        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0360      I      diagName, diagCode, diagUnits, diagTitle, 0, myThid )
                0361
                0362 C     these diagnostics are never filled
                0363 c      diagName  = 'SHIICFFT'
                0364 c      diagTitle = 'total SHI and ICF forcing for T, >0 increases theta'
                0365 c      diagUnits = 'W/m^2           '
                0366 c      diagCode  = 'SM      L1      '
                0367 c      CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0368 c    I      diagName, diagCode, diagUnits, diagTitle, 0, myThid )
                0369
                0370 c      diagName  = 'SHIICFFS'
                0371 c      diagTitle = 'total SHI and ICF forcing for S, >0 increases salt'
                0372 c      diagUnits = 'g/m^2/s         '
                0373 c      diagCode  = 'SM      L1      '
                0374 c      CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0375 c    I      diagName, diagCode, diagUnits, diagTitle, 0, myThid )
                0376
                0377       ENDIF
                0378 #endif
                0379
                0380       RETURN
                0381       END