pkg/shelfice/shelfice_init_fixed.F

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