pkg/thsice/thsice_extend.F

87ea84cac6 Jean*0001 #include "THSICE_OPTIONS.h"
6b47d550f4 Mart*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
87ea84cac6 Jean*0005
                0006 CBOP
                0007 C     !ROUTINE: THSICE_EXTEND
                0008 C     !INTERFACE:
                0009       SUBROUTINE THSICE_EXTEND(
6dc8890c80 Patr*0010      I                  bi, bj,
7269783f6f Jean*0011      I                  iMin,iMax, jMin,jMax, dBugFlag,
                0012      I                  fzMlOc, tFrz, tOce,
                0013      U                  icFrac, hIce, hSnow,
                0014      U                  tSrf, tIc1, tIc2, qIc1, qIc2,
                0015      O                  flx2oc, frw2oc, fsalt,
                0016      I                  myTime, myIter, myThid )
87ea84cac6 Jean*0017 C     !DESCRIPTION: \bv
                0018 C     *==========================================================*
7269783f6f Jean*0019 C     | S/R  THSICE_EXTEND
87ea84cac6 Jean*0020 C     | o Extend sea-ice area incresing ice fraction
                0021 C     *==========================================================*
7269783f6f Jean*0022 C     | o incorporate surplus of energy to
                0023 C     |   make new ice or make ice grow laterally
87ea84cac6 Jean*0024 C     *==========================================================*
                0025 C     \ev
                0026
                0027 C     !USES:
                0028       IMPLICIT NONE
                0029
                0030 C     == Global variables ==
dbce8fc2d4 Jean*0031 #include "EEPARAMS.h"
f8d7459e30 Patr*0032 #include "SIZE.h"
87ea84cac6 Jean*0033 #include "THSICE_SIZE.h"
                0034 #include "THSICE_PARAMS.h"
d6f06800ae Patr*0035 #ifdef ALLOW_AUTODIFF_TAMC
                0036 # include "tamc.h"
                0037 #endif
87ea84cac6 Jean*0038
                0039 C     !INPUT/OUTPUT PARAMETERS:
                0040 C     == Routine Arguments ==
7269783f6f Jean*0041 C     bi,bj       :: tile indices
                0042 C     iMin,iMax   :: computation domain: 1rst index range
                0043 C     jMin,jMax   :: computation domain: 2nd  index range
                0044 C     dBugFlag    :: allow to print debugging stuff (e.g. on 1 grid point).
                0045 C---  Input:
                0046 C         iceMask :: sea-ice fractional mask [0-1]
                0047 C  fzMlOc (esurp) :: ocean mixed-layer freezing/melting potential [W/m2]
                0048 C  tFrz    (Tf)   :: sea-water freezing temperature [oC] (function of S)
790347483b Jean*0049 C  tOce           :: surface level oceanic temperature [oC]
7269783f6f Jean*0050 C---  Modified (input&output):
                0051 C  icFrac(iceFrac):: fraction of grid area covered in ice
                0052 C  hIce (iceThick):: ice height [m]
790347483b Jean*0053 C  hSnow          :: snow height [m]
7269783f6f Jean*0054 C  tSrf           :: surface (ice or snow) temperature [oC]
                0055 C  tIc1           :: temperature of ice layer 1 [oC]
                0056 C  tIc2           :: temperature of ice layer 2 [oC]
                0057 C  qIc1   (qicen) :: ice enthalpy (J/kg), 1rst level
                0058 C  qIc2   (qicen) :: ice enthalpy (J/kg), 2nd level
                0059 C---  Output
                0060 C  flx2oc   (=)   :: (additional) heat flux to ocean    [W/m2]        (+=dwn)
                0061 C  frw2oc   (=)   :: (additional) fresh water flux to ocean [kg/m2/s] (+=dwn)
                0062 C  fsalt    (=)   :: (additional) salt flux to ocean        [g/m2/s]  (+=dwn)
                0063 C---  Input:
                0064 C     myTime      :: current Time of simulation [s]
                0065 C     myIter      :: current Iteration number in simulation
                0066 C     myThid      :: my Thread Id number
                0067       INTEGER bi,bj
                0068       INTEGER iMin, iMax
                0069       INTEGER jMin, jMax
                0070       LOGICAL dBugFlag
6dc8890c80 Patr*0071 c     _RL iceMask(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0072       _RL fzMlOc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0073       _RL tFrz   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0074       _RL tOce   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0075       _RL icFrac (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0076       _RL hIce   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0077       _RL hSnow  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0078       _RL tSrf   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0079       _RL tIc1   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0080       _RL tIc2   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0081       _RL qIc1   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0082       _RL qIc2   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0083       _RL flx2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0084       _RL frw2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0085       _RL fsalt  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
7269783f6f Jean*0086       _RL  myTime
                0087       INTEGER myIter
                0088       INTEGER myThid
                0089 CEOP
                0090
                0091 #ifdef ALLOW_THSICE
                0092 C     !LOCAL VARIABLES:
                0093 C---  local copy of input/output argument list variables (see description above)
87ea84cac6 Jean*0094       _RL esurp
                0095       _RL Tf
7269783f6f Jean*0096       _RL iceFrac
87ea84cac6 Jean*0097       _RL iceThick
                0098       _RL qicen(nlyr)
                0099
                0100 C     == Local variables ==
790347483b Jean*0101 C     iceVol    :: previous ice volume
                0102 C     newIce    :: new ice volume to produce
                0103 C     hNewIce   :: thickness of new ice to form
                0104 C     iceFormed :: ice-volume formed (new ice volume = iceVol+iceFormed )
                0105 C     qicAv     :: mean enthalpy of ice (layer 1 & 2) [J/m^3]
87ea84cac6 Jean*0106       _RL deltaTice ! time-step for ice model
790347483b Jean*0107       _RL iceVol
87ea84cac6 Jean*0108       _RL newIce
790347483b Jean*0109       _RL hNewIce
                0110       _RL iceFormed
87ea84cac6 Jean*0111       _RL qicAv
7269783f6f Jean*0112       INTEGER  i,j     ! loop indices
                0113
7c50f07931 Mart*0114 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0115 C     tkey :: tape key (depends on tiles)
                0116       INTEGER tkey
7c50f07931 Mart*0117 #endif
7269783f6f Jean*0118 C-    define grid-point location where to print debugging values
                0119 #include "THSICE_DEBUG.h"
87ea84cac6 Jean*0120
a85293d087 Mart*0121 #ifdef ALLOW_DBUG_THSICE
87ea84cac6 Jean*0122  1020 FORMAT(A,1P4E11.3)
a85293d087 Mart*0123 #endif
87ea84cac6 Jean*0124
7269783f6f Jean*0125 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0126
86e6447a82 Patr*0127       deltaTice = thSIce_deltaT
                0128
6b47d550f4 Mart*0129 #ifdef ALLOW_AUTODIFF
86e6447a82 Patr*0130       DO j = 1-OLy, sNy+OLy
                0131        DO i = 1-OLx, sNx+OLx
                0132          flx2oc(i,j) = 0. _d 0
                0133          frw2oc(i,j) = 0. _d 0
                0134          fsalt (i,j) = 0. _d 0
                0135        ENDDO
                0136       ENDDO
6b47d550f4 Mart*0137 #endif /* ALLOW_AUTODIFF */
86e6447a82 Patr*0138
d6f06800ae Patr*0139 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0140       tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0141 CADJ STORE hIce   = comlev1_bibj, key = tkey, kind = isbyte
                0142 CADJ STORE hSnow  = comlev1_bibj, key = tkey, kind = isbyte
                0143 CADJ STORE icFrac = comlev1_bibj, key = tkey, kind = isbyte
                0144 CADJ STORE qIc1   = comlev1_bibj, key = tkey, kind = isbyte
                0145 CADJ STORE qIc2   = comlev1_bibj, key = tkey, kind = isbyte
d6f06800ae Patr*0146 #endif
486e9de820 Jean*0147       DO j = jMin, jMax
                0148        DO i = iMin, iMax
d6f06800ae Patr*0149
7269783f6f Jean*0150         IF (fzMlOc(i,j).GT.0. _d 0) THEN
a85293d087 Mart*0151          esurp     = fzMlOc(i,j)
                0152          Tf        = tFrz(i,j)
                0153          iceFrac   = icFrac(i,j)
                0154          iceThick  = hIce(i,j)
                0155          qicen(1)  = qIc1(i,j)
                0156          qicen(2)  = qIc2(i,j)
7269783f6f Jean*0157 C---
87ea84cac6 Jean*0158 C--   start ice
a85293d087 Mart*0159          iceFormed = 0. _d 0
                0160          iceVol    = iceFrac*iceThick
87ea84cac6 Jean*0161
                0162 C-    enthalpy of new ice to form :
a85293d087 Mart*0163          IF ( iceFrac.LE.0. _d 0 ) THEN
                0164           qicen(1) = -cpWater*Tmlt1
                0165      &              + cpIce *(Tmlt1-Tf) + Lfresh*(1. _d 0-Tmlt1/Tf)
                0166           qicen(2) = -cpIce *Tf + Lfresh
                0167          ENDIF
                0168          qicAv     = rhoi*(qicen(1)+qicen(2))*0.5 _d 0
                0169          newIce    = esurp*deltaTice/qicAv
                0170
                0171          IF ( icFrac(i,j).EQ.0. _d 0 ) THEN
                0172 C-  to keep identical results (as it used to be):
                0173 c-old_v: IF ( newIce.GE.hThinIce*iceMaskMin ) THEN
790347483b Jean*0174 C-  here we allow to form ice earlier (as soon as min-ice-vol is reached)
5471f8df19 Jean*0175 c-new_v:
                0176           IF ( newIce.GT.hIceMin*iceMaskMin ) THEN
790347483b Jean*0177 C-    if there is no ice in grid-cell and enough ice to form:
a85293d087 Mart*0178 C-    make ice over iceMaskMin fraction, up to hThinIce,
                0179 C     and if more ice to form, then increase fraction
                0180            iceThick  = MIN(hThinIce,newIce/iceMaskMin)
                0181            iceThick  = MAX(iceThick,newIce/iceMaskMax)
                0182            iceFrac   = newIce/iceThick
                0183            iceFormed = newIce
87ea84cac6 Jean*0184           ENDIF
a85293d087 Mart*0185          ELSEIF ( iceVol.LT.hiMax*iceMaskMax ) THEN
87ea84cac6 Jean*0186 C-    if there is already some ice
790347483b Jean*0187 C     create ice with same thickness or hNewIceMax (the smallest of the 2)
a85293d087 Mart*0188            hNewIce = MIN(hIce(i,j),hNewIceMax)
                0189            iceFrac = MIN(icFrac(i,j)+newIce/hNewIce,iceMaskMax)
790347483b Jean*0190 C-    update thickness: area weighted average
5471f8df19 Jean*0191 c-new_v:
a85293d087 Mart*0192            iceThick = MIN(hiMax,(iceVol+newIce)/iceFrac)
790347483b Jean*0193 C-  to keep identical results: comment the line above and uncomment line below:
5471f8df19 Jean*0194 c-old_v     iceFrac = MIN(icFrac(i,j)+newIce/iceThick,iceMaskMax)
a85293d087 Mart*0195            iceFormed = iceThick*iceFrac - iceVol
87ea84cac6 Jean*0196 C-    spread snow out over ice
a85293d087 Mart*0197            hSnow(i,j) = hSnow(i,j)*icFrac(i,j)/iceFrac
                0198          ENDIF
7269783f6f Jean*0199 C-    oceanic fluxes:
a85293d087 Mart*0200          flx2oc(i,j) = qicAv*iceFormed/deltaTice
                0201          frw2oc(i,j) = -rhoi*iceFormed/deltaTice
                0202          fsalt (i,j) = -(rhoi*saltIce)*iceFormed/deltaTice
87ea84cac6 Jean*0203
7269783f6f Jean*0204 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0205 #ifdef ALLOW_DBUG_THSICE
a85293d087 Mart*0206          IF ( dBug(i,j,bi,bj) ) THEN
7269783f6f Jean*0207           WRITE(6,1020) 'ThSI_EXT: iceH, newIce, newIceFrac=',
790347483b Jean*0208      &                       iceThick, newIce, iceFrac-icFrac(i,j)
7269783f6f Jean*0209           WRITE(6,1020) 'ThSI_EXT: iceFrac,flx2oc,fsalt,frw2oc=',
                0210      &                  iceFrac,flx2oc(i,j),fsalt(i,j),frw2oc(i,j)
a85293d087 Mart*0211          ENDIF
7269783f6f Jean*0212 #endif
                0213 #ifdef CHECK_ENERGY_CONSERV
a85293d087 Mart*0214          CALL THSICE_CHECK_CONSERV( dBugFlag, i, j, bi, bj, 1,
790347483b Jean*0215      I            icFrac(i,j), iceFrac, iceThick, hSnow(i,j), qicen,
7269783f6f Jean*0216      I            flx2oc(i,j), frw2oc(i,j), fsalt(i,j),
                0217      I            myTime, myIter, myThid )
                0218 #endif /* CHECK_ENERGY_CONSERV */
                0219 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0220 C-- Update Sea-Ice state output:
a85293d087 Mart*0221          IF ( iceFrac.GT.0. _d 0 .AND. icFrac(i,j).EQ.0. _d 0) THEN
                0222 c         hSnow(i,j) = 0. _d 0
                0223           tSrf(i,j)  = tFrz(i,j)
                0224           tIc1(i,j)  = tFrz(i,j)
                0225           tIc2(i,j)  = tFrz(i,j)
                0226           qIc1(i,j)  = qicen(1)
                0227           qIc2(i,j)  = qicen(2)
                0228          ENDIF
                0229          icFrac(i,j) = iceFrac
                0230          hIce  (i,j) = iceThick
7269783f6f Jean*0231         ENDIF
                0232        ENDDO
                0233       ENDDO
87ea84cac6 Jean*0234
                0235 #endif /* ALLOW_THSICE */
                0236
                0237       RETURN
                0238       END