Back to home page

MITgcm

 
 

    


File indexing completed on 2025-11-07 06:08:50 UTC

view on githubraw file Latest commit b7411f1a on 2025-11-06 19:05:26 UTC
ffe464dc7d Mart*0001 #include "SHELFICE_OPTIONS.h"
ceca9ed7c5 Jean*0002 #ifdef ALLOW_COST
                0003 # include "COST_OPTIONS.h"
                0004 #endif
                0005 #ifdef ALLOW_CTRL
                0006 # include "CTRL_OPTIONS.h"
                0007 #endif
a340904e5a Ou W*0008 #ifdef ALLOW_MOM_COMMON
                0009 # include "MOM_COMMON_OPTIONS.h"
                0010 #endif
ffe464dc7d Mart*0011 
                0012       SUBROUTINE SHELFICE_INIT_FIXED( myThid )
afaeb4fe62 Jean*0013 C     *============================================================*
                0014 C     | SUBROUTINE SHELFICE_INIT_FIXED
                0015 C     | o Routine to initialize SHELFICE parameters and variables.
                0016 C     *============================================================*
                0017 C     | Initialize SHELFICE parameters and variables.
                0018 C     *============================================================*
ffe464dc7d Mart*0019       IMPLICIT NONE
                0020 
                0021 C     === Global variables ===
                0022 #include "SIZE.h"
                0023 #include "EEPARAMS.h"
                0024 #include "PARAMS.h"
                0025 #include "GRID.h"
                0026 #include "SHELFICE.h"
abed653eb5 Mart*0027 #ifdef ALLOW_COST
                0028 # include "cost.h"
                0029 # include "SHELFICE_COST.h"
                0030 #endif /* ALLOW_COST */
ffe464dc7d Mart*0031 
                0032 C     === Routine arguments ===
4354afafa2 Jean*0033 C     myThid ::  Number of this instance of SHELFICE_INIT_FIXED
ffe464dc7d Mart*0034       INTEGER myThid
                0035 
                0036 #ifdef ALLOW_SHELFICE
                0037 C     === Local variables ===
4354afafa2 Jean*0038 C     i, j, bi, bj :: Loop counters
005af54e38 Jean*0039 C     SHI_minDepth :: minimum Ice-Shelf depth
                0040 C     msgBuf       :: Informational/error message buffer
4354afafa2 Jean*0041       INTEGER i, j, bi, bj
005af54e38 Jean*0042       _RS SHI_minDepth
                0043       CHARACTER*(MAX_LEN_MBUF) msgBuf
4c0efb30d8 Mart*0044 #ifdef ALLOW_DIAGNOSTICS
                0045       INTEGER       diagNum
8e7e785cad Jean*0046       INTEGER       diagMate
4c0efb30d8 Mart*0047       CHARACTER*8   diagName
                0048       CHARACTER*16  diagCode
                0049       CHARACTER*16  diagUnits
                0050       CHARACTER*(80) diagTitle
                0051 #endif /* ALLOW_DIAGNOSTICS */
6b47d550f4 Mart*0052 #ifdef ALLOW_CTRL
6861e183cf Jean*0053       INTEGER k
abed653eb5 Mart*0054 #endif
ffe464dc7d Mart*0055 
afaeb4fe62 Jean*0056 #ifdef ALLOW_MNC
                0057 C     Initialize MNC variable information for SHELFICE
b7411f1a84 Jean*0058       IF ( useMNC .AND. shelfice_dump_mnc ) THEN
afaeb4fe62 Jean*0059         CALL SHELFICE_MNC_INIT( myThid )
                0060       ENDIF
                0061 #endif /* ALLOW_MNC */
                0062 
005af54e38 Jean*0063       DO bj = myByLo(myThid), myByHi(myThid)
                0064        DO bi = myBxLo(myThid), myBxHi(myThid)
                0065         DO j = 1-OLy, sNy+OLy
                0066          DO i = 1-OLx, sNx+OLx
                0067           shelficeMassInit   (i,j,bi,bj) = 0. _d 0
                0068           shelficeLoadAnomaly(i,j,bi,bj) = 0. _d 0
                0069           shelfIceMassDynTendency(i,j,bi,bj) = 0. _d 0
                0070          ENDDO
                0071         ENDDO
                0072        ENDDO
                0073       ENDDO
                0074 
                0075       IF ( SHELFICEloadAnomalyFile .NE. ' ' ) THEN
                0076        CALL READ_FLD_XY_RL( SHELFICEloadAnomalyFile, ' ',
                0077      &                      shelficeLoadAnomaly, 0, myThid )
                0078       ENDIF
                0079       IF ( SHELFICEmassFile.NE.' ' ) THEN
                0080        CALL READ_FLD_XY_RL( SHELFICEmassFile, ' ',
                0081      &                      shelficeMassInit, 0, myThid )
                0082       ELSE
                0083        DO bj = myByLo(myThid), myByHi(myThid)
                0084         DO bi = myBxLo(myThid), myBxHi(myThid)
                0085          DO j = 1, sNy
                0086           DO i = 1, sNx
                0087            shelficeMassInit(i,j,bi,bj) =
                0088      &         shelficeLoadAnomaly(i,j,bi,bj)*recip_gravity
                0089      &       - rhoConst*Ro_surf(i,j,bi,bj)
                0090           ENDDO
                0091          ENDDO
                0092         ENDDO
                0093        ENDDO
                0094       ENDIF
                0095       _EXCH_XY_RL( shelficeMassInit, myThid )
                0096 
ffe464dc7d Mart*0097 C-----------------------------------------------------------------------
8f81ce1bc7 Mart*0098 C--   Initialize SHELFICE variables kTopC
9f5091f233 Jean*0099 C--   kTopC is the same as kSurfC, except outside ice-shelf area:
                0100 C--   kTop = 0 where there is no ice-shelf (where kSurfC=1)
                0101 C--   and over land (completely dry column) where kSurfC = Nr+1
ffe464dc7d Mart*0102 C-----------------------------------------------------------------------
                0103 
005af54e38 Jean*0104 C--   Currently, there is no specific pkg/shelfice criteria to decide when very
                0105 C     thin ice-shelf vanishes: some limits come from ocean hFacMin/hFacMinDr;
                0106 C     but should depend on shelficeMass if SHELFICEMassStepping is used.
                0107 C     Here the SHI_minDepth threshold is just to get around machine round-off
                0108 C     isssues so does not need to be a run-time params.
                0109       SHI_minDepth = 1. _d -6
                0110       SHI_minDepth = rF(1) - drF(1)*SHI_minDepth
ffe464dc7d Mart*0111       DO bj = myByLo(myThid), myByHi(myThid)
                0112        DO bi = myBxLo(myThid), myBxHi(myThid)
005af54e38 Jean*0113         IF ( SHELFICEMassStepping ) THEN
                0114 C-      In this case, it is recommended to initialise ice-shelf Mass (from file)
                0115 C       And if updating "kTopC" is permitted (SHI_update_kTopC=T), kTopC will
                0116 C       be reset in SHELFICE_INIT_VARIA after loading shelficeMass from pickup
                0117          DO j = 1-OLy, sNy+OLy
                0118           DO i = 1-OLx, sNx+OLx
                0119            IF ( kSurfC(i,j,bi,bj).LE.Nr .AND.
                0120      &          shelficeMassInit(i,j,bi,bj).GT.zeroRL ) THEN
4354afafa2 Jean*0121             kTopC(i,j,bi,bj) = kSurfC(i,j,bi,bj)
005af54e38 Jean*0122            ELSE
4354afafa2 Jean*0123             kTopC(i,j,bi,bj) = 0
005af54e38 Jean*0124            ENDIF
                0125           ENDDO
ffe464dc7d Mart*0126          ENDDO
005af54e38 Jean*0127         ELSE
                0128 C-      no time-evolving Ice-Shelf: set kTopC according to Ro_surf
                0129 C       Note: use Ro_surf instead of R_shelfIce to stay consistent
                0130 C             with ocean mask & hFac, see issue #99
                0131          DO j = 1-OLy, sNy+OLy
                0132           DO i = 1-OLx, sNx+OLx
                0133            IF ( kSurfC(i,j,bi,bj).LE.Nr .AND.
                0134      &          Ro_surf(i,j,bi,bj).LT.SHI_minDepth ) THEN
                0135 c    &          R_shelfIce(i,j,bi,bj).LT.rF(1) ) THEN
                0136             kTopC(i,j,bi,bj) = kSurfC(i,j,bi,bj)
                0137            ELSE
                0138             kTopC(i,j,bi,bj) = 0
                0139            ENDIF
                0140           ENDDO
                0141          ENDDO
                0142         ENDIF
ffe464dc7d Mart*0143        ENDDO
                0144       ENDDO
6861e183cf Jean*0145 
6b47d550f4 Mart*0146 #ifdef ALLOW_CTRL
6861e183cf Jean*0147 C     maskSHI is a hack to play along with the general ctrl-package
                0148 C     infrastructure, where only the k=1 layer of a 3D mask is used
                0149 C     for 2D fields. We cannot use maskInC instead, because routines
                0150 C     like ctrl_get_gen and ctrl_set_unpack_xy require 3D masks.
                0151       DO bj = myByLo(myThid), myByHi(myThid)
                0152        DO bi = myBxLo(myThid), myBxHi(myThid)
                0153         DO k=1,Nr
                0154          DO j=1-OLy,sNy+OLy
                0155           DO i=1-OLx,sNx+OLx
                0156            maskSHI(i,j,k,bi,bj) = 0. _d 0
                0157           ENDDO
                0158          ENDDO
                0159         ENDDO
                0160         DO k=1,Nr
                0161          DO j=1-OLy,sNy+OLy
                0162           DO i=1-OLx,sNx+OLx
005af54e38 Jean*0163            IF ( kTopC(i,j,bi,bj).GE.1
f8ef4d8fd6 Oliv*0164      &          .AND. maskC(i,j,k,bi,bj).NE.zeroRS ) THEN
6861e183cf Jean*0165             maskSHI(i,j,k,bi,bj) = 1. _d 0
                0166             maskSHI(i,j,1,bi,bj) = 1. _d 0
                0167            ENDIF
                0168           ENDDO
                0169          ENDDO
                0170         ENDDO
                0171        ENDDO
                0172       ENDDO
6b47d550f4 Mart*0173 #endif /* ALLOW_CTRL */
6861e183cf Jean*0174 
005af54e38 Jean*0175       IF ( debugLevel.GE.debLevC ) THEN
                0176        WRITE(msgBuf,'(A)')
                0177      &  'SHELFICE_INIT_FIXED: checking Ice-Shelf extension & kTopC'
                0178        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0179      &                     SQUEEZE_RIGHT, myThid )
                0180        WRITE(msgBuf,'(A,1PE18.10)')
                0181      &  'SHELFICE_INIT_FIXED:  SHI_minDepth =', SHI_minDepth
                0182        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0183      &                     SQUEEZE_RIGHT, myThid )
4354afafa2 Jean*0184        DO bj = myByLo(myThid), myByHi(myThid)
                0185         DO bi = myBxLo(myThid), myBxHi(myThid)
005af54e38 Jean*0186 c        DO j = 1-OLy, sNy+OLy
                0187 c         DO i = 1-OLx, sNx+OLx
4354afafa2 Jean*0188          DO j = 1, sNy
                0189           DO i = 1, sNx
005af54e38 Jean*0190            IF ( kTopC(i,j,bi,bj).NE.0 .AND.
                0191      &          ( shelficeMassInit(i,j,bi,bj).LE.zeroRL
                0192      &            .OR. R_shelfIce(i,j,bi,bj).GE.rF(1) ) ) THEN
                0193              WRITE(standardMessageUnit,'(A,2I4,2I3,A,1P3E17.9)')
                0194      &        ' Ice-Shelf (kTopC>0) in i,j,bi,bj=', i, j, bi, bj,
                0195      &        ' but R_shelf,Ro_surf,MassIni=', R_shelfIce(i,j,bi,bj),
                0196      &         Ro_surf(i,j,bi,bj), shelficeMassInit(i,j,bi,bj)
                0197            ELSEIF ( kTopC(i,j,bi,bj).EQ.0 .AND.
                0198      &          ( shelficeMassInit(i,j,bi,bj).GT.zeroRL
                0199      &            .OR. R_shelfIce(i,j,bi,bj).LT.rF(1) ) ) THEN
                0200              WRITE(standardMessageUnit,'(A,2I4,2I3,A,1P3E17.9)')
                0201      &        ' no Ice-Shelf (kTopC=0) i,j,bi,bj=', i, j, bi, bj,
                0202      &        ' but R_shelf,Ro_surf,MassIni=', R_shelfIce(i,j,bi,bj),
                0203      &         Ro_surf(i,j,bi,bj), shelficeMassInit(i,j,bi,bj)
                0204            ENDIF
4354afafa2 Jean*0205           ENDDO
                0206          ENDDO
                0207         ENDDO
                0208        ENDDO
005af54e38 Jean*0209        WRITE(msgBuf,'(A)')
                0210      &  'SHELFICE_INIT_FIXED: checking Ice-Shelf extension: done'
                0211        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0212      &                     SQUEEZE_RIGHT, myThid )
4354afafa2 Jean*0213       ENDIF
                0214 c     IF ( SHELFICEloadAnomalyFile .EQ. ' ' ) THEN
                0215 C-   In case we need shelficeLoadAnomaly in phi0surf for initial pressure
0236599f51 Jean*0216 C    calculation (if using selectP_inEOS_Zc=2 or 3)
4354afafa2 Jean*0217        DO bj = myByLo(myThid), myByHi(myThid)
                0218         DO bi = myBxLo(myThid), myBxHi(myThid)
                0219          DO j = 1-OLy, sNy+OLy
                0220           DO i = 1-OLx, sNx+OLx
005af54e38 Jean*0221            IF ( kTopC(i,j,bi,bj).EQ.0 ) THEN
                0222             R_shelfIce         (i,j,bi,bj) = rF(1)
                0223             shelficeMassInit   (i,j,bi,bj) = 0. _d 0
                0224             shelficeLoadAnomaly(i,j,bi,bj) = 0. _d 0
                0225            ELSE
                0226             shelficeLoadAnomaly(i,j,bi,bj) = gravity
3d2f509a67 Dani*0227      &      *(shelficeMassInit(i,j,bi,bj)+rhoConst*Ro_surf(i,j,bi,bj))
005af54e38 Jean*0228            ENDIF
4354afafa2 Jean*0229           ENDDO
                0230          ENDDO
                0231         ENDDO
                0232        ENDDO
                0233 c     ELSE
                0234 c      _EXCH_XY_RS( shelficeLoadAnomaly, myThid )
                0235 c     ENDIF
005af54e38 Jean*0236       CALL WRITE_FLD_XY_RL ( 'shelficemassinit', ' ',
                0237      &                       shelficeMassInit, 0, myThid )
4354afafa2 Jean*0238       IF ( debugLevel.GE.debLevC ) THEN
4d294a5719 Jean*0239        CALL WRITE_FLD_XY_RL( 'SHICE_pLoadAnom', ' ',
4354afafa2 Jean*0240      I                       shelficeLoadAnomaly, -1, myThid )
                0241       ENDIF
                0242 
a0178c5a01 Jean*0243       IF ( SHELFICEMassStepping .AND.
                0244      &     SHELFICEMassDynTendFile .NE. ' ' ) THEN
                0245        CALL READ_FLD_XY_RS( SHELFICEMassDynTendFile, ' ',
                0246      &                      shelfIceMassDynTendency, 0, myThid )
3d2f509a67 Dani*0247       ENDIF
                0248 
4c0efb30d8 Mart*0249 #ifdef ALLOW_DIAGNOSTICS
                0250       IF ( useDiagnostics ) THEN
                0251        diagName  = 'SHIfwFlx'
e84c33192a Mart*0252        diagTitle = 'Ice shelf fresh water flux (positive upward)'
2bf81d7d9d Mart*0253        diagUnits = 'kg/m^2/s        '
8e7e785cad Jean*0254        diagCode  = 'SM      L1      '
                0255        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0256      I      diagName, diagCode, diagUnits, diagTitle, 0, myThid )
4c0efb30d8 Mart*0257 
                0258        diagName  = 'SHIhtFlx'
e84c33192a Mart*0259        diagTitle = 'Ice shelf heat flux  (positive upward)'
4c0efb30d8 Mart*0260        diagUnits = 'W/m^2           '
8e7e785cad Jean*0261        diagCode  = 'SM      L1      '
                0262        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0263      I      diagName, diagCode, diagUnits, diagTitle, 0, myThid )
4c0efb30d8 Mart*0264 
1fa1235c63 Jean*0265        diagName  = 'SHI_TauX'
                0266        diagTitle =
                0267      &     'Ice shelf bottom stress, zonal  comp., >0 increases uVel'
                0268        diagUnits = 'N/m^2           '
                0269        diagCode  = 'UU      L1      '
8e7e785cad Jean*0270        diagMate  = diagNum + 2
                0271        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0272      I      diagName, diagCode, diagUnits, diagTitle, diagMate, myThid )
4c0efb30d8 Mart*0273 
1fa1235c63 Jean*0274        diagName  = 'SHI_TauY'
                0275        diagTitle =
                0276      &     'Ice shelf bottom stress, merid. comp., >0 increases vVel'
                0277        diagUnits = 'N/m^2           '
                0278        diagCode  = 'VV      L1      '
8e7e785cad Jean*0279        diagMate  = diagNum
                0280        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0281      I      diagName, diagCode, diagUnits, diagTitle, diagMate, myThid )
a86b6152f4 Dimi*0282 
                0283        diagName  = 'SHIForcT'
                0284        diagTitle = 'Ice shelf forcing for theta, >0 increases theta'
                0285        diagUnits = 'W/m^2           '
8e7e785cad Jean*0286        diagCode  = 'SM      L1      '
                0287        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0288      I      diagName, diagCode, diagUnits, diagTitle, 0, myThid )
a86b6152f4 Dimi*0289 
                0290        diagName  = 'SHIForcS'
                0291        diagTitle = 'Ice shelf forcing for salt, >0 increases salt'
                0292        diagUnits = 'g/m^2/s         '
8e7e785cad Jean*0293        diagCode  = 'SM      L1      '
                0294        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0295      I      diagName, diagCode, diagUnits, diagTitle, 0, myThid )
84ccdf230a Mart*0296 
                0297        diagName  = 'SHIgammT'
                0298        diagTitle = 'Ice shelf exchange coefficient for theta'
                0299        diagUnits = 'm/s             '
                0300        diagCode  = 'SM      L1      '
                0301        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0302      I      diagName, diagCode, diagUnits, diagTitle, 0, myThid )
                0303 
                0304        diagName  = 'SHIgammS'
                0305        diagTitle = 'Ice shelf exchange coefficient for salt'
                0306        diagUnits = 'm/s             '
                0307        diagCode  = 'SM      L1      '
                0308        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0309      I      diagName, diagCode, diagUnits, diagTitle, 0, myThid )
3d2f509a67 Dani*0310 
7b8b86ab99 Timo*0311        diagName  = 'SHI_mass'
                0312        diagTitle = 'dynamic ice shelf mass for surface load anomaly'
                0313        diagUnits = 'kg/m^2          '
                0314        diagCode  = 'SM      L1      '
                0315        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0316      I      diagName, diagCode, diagUnits, diagTitle, 0, myThid )
                0317 
a340904e5a Ou W*0318 #if ( defined ALLOW_MOM_COMMON && defined ALLOW_MOM_TEND_EXTRA_DIAGS )
                0319 cHP adding diagnostic for drag from ice shelf to allow exp. dissip decomp
                0320       diagName  = 'UShIDrag'
                0321       diagTitle = 'U momentum tendency from ice shelf Drag'
                0322       diagCode  = 'UUR     MR      '
                0323       diagMate  = diagNum + 2
                0324       CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0325      I   diagName, diagCode, diagUnits, diagTitle, diagMate, myThid )
                0326       diagName  = 'VShIDrag'
                0327       diagTitle = 'V momentum tendency from ice shelf Drag'
                0328       diagCode  = 'VVR     MR      '
                0329       diagMate  = diagNum
                0330       CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0331      I   diagName, diagCode, diagUnits, diagTitle, diagMate, myThid )
                0332 #endif
                0333 
7b8b86ab99 Timo*0334 #ifdef SHI_ALLOW_GAMMAFRICT
5744811f23 Mart*0335        diagName  = 'SHIuStar'
                0336        diagTitle = 'Friction velocity at bottom of ice shelf'
                0337        diagUnits = 'm/s             '
                0338        diagCode  = 'SM      L1      '
                0339        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0340      I      diagName, diagCode, diagUnits, diagTitle, 0, myThid )
                0341 
7b8b86ab99 Timo*0342        diagName  = 'SHICDrag'
                0343        diagTitle = 'Shelfice drag coefficient for u* parameterization'
                0344        diagUnits = '1               '
3d2f509a67 Dani*0345        diagCode  = 'SM      L1      '
                0346        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
7b8b86ab99 Timo*0347      I          diagName, diagCode, diagUnits, diagTitle, 0, myThid )
                0348 #endif
9952f046d7 dngo*0349 
                0350 #ifdef ALLOW_SHELFICE_REMESHING
                0351        diagName  = 'SHIRshel'
                0352        diagTitle = 'depth of shelfice'
                0353        diagUnits = 'm               '
                0354        diagCode  = 'SM      L1      '
                0355        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0356      I      diagName, diagCode, diagUnits, diagTitle, 0, myThid )
                0357 #endif
                0358 
7b8b86ab99 Timo*0359 # ifdef ALLOW_AUTODIFF
                0360 #  ifndef SHI_ALLOW_GAMMAFRICT
                0361        diagName  = 'ADJshict'
                0362        diagTitle = 'dJ/dgammaT: Sens. to shelfice heat transfer coeff'
                0363        diagUnits = 'dJ/(m/s)        '
                0364        diagCode  = 'SM A    M1      '
                0365        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0366      I          diagName, diagCode, diagUnits, diagTitle, 0, myThid )
                0367 
                0368        diagName  = 'ADJshics'
                0369        diagTitle = 'dJ/dgammaS: Sens. to shelfice salt transfer coeff'
                0370        diagUnits = 'dJ/(m/s)        '
                0371        diagCode  = 'SM A    M1      '
                0372        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0373      I          diagName, diagCode, diagUnits, diagTitle, 0, myThid )
                0374 
                0375 #  else
                0376        diagName  = 'ADJshicd'
                0377        diagTitle = 'dJ/dcDrag: Sensitivity to shelfice u* drag coeff'
                0378        diagUnits = 'dJ/1            '
                0379        diagCode  = 'SM A    M1      '
                0380        CALL DIAGNOSTICS_ADDTOLIST( diagNum,
                0381      I          diagName, diagCode, diagUnits, diagTitle, 0, myThid )
                0382 #  endif
                0383 # endif
4c0efb30d8 Mart*0384       ENDIF
                0385 #endif /* ALLOW_DIAGNOSTICS */
ffe464dc7d Mart*0386 #endif /* ALLOW_SHELFICE */
                0387 
a93eea6699 Jean*0388       RETURN
                0389       END