Back to home page

MITgcm

 
 

    


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

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