pkg/thsice/thsice_solve4temp.F

87ea84cac6 Jean*0001 #include "THSICE_OPTIONS.h"
6b47d550f4 Mart*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
e3123da85c Patr*0004 # ifdef ALLOW_EXF
                0005 #  include "EXF_OPTIONS.h"
                0006 # endif
a85293d087 Mart*0007 # define ALLOW_AUTODIFF_TAMC_MORE
e3123da85c Patr*0008 #endif
87ea84cac6 Jean*0009
                0010 CBOP
                0011 C     !ROUTINE: THSICE_SOLVE4TEMP
                0012 C     !INTERFACE:
                0013       SUBROUTINE THSICE_SOLVE4TEMP(
6dc8890c80 Patr*0014      I                  bi, bj,
72e380ddb5 Jean*0015      I                  iMin,iMax, jMin,jMax, dBugFlag,
2a9474d935 Mart*0016      I                  useBulkForce, useEXF,
a85293d087 Mart*0017      I                  icMask, hIce, hSnow, tFrz, flxExSW,
                0018      U                  flxSW, tSrf, qIc1, qIc2,
                0019      O                  tIc1, tIc2, dTsrf,
7269783f6f Jean*0020      O                  sHeat, flxCnB, flxAtm, evpAtm,
                0021      I                  myTime, myIter, myThid )
87ea84cac6 Jean*0022 C     !DESCRIPTION: \bv
                0023 C     *==========================================================*
                0024 C     | S/R  THSICE_SOLVE4TEMP
                0025 C     *==========================================================*
                0026 C     | Solve (implicitly) for sea-ice and surface temperature
                0027 C     *==========================================================*
                0028 C     \ev
                0029
7269783f6f Jean*0030 C ADAPTED FROM:
                0031 C LANL CICE.v2.0.2
                0032 C-----------------------------------------------------------------------
                0033 C.. thermodynamics (vertical physics) based on M. Winton 3-layer model
88f72205aa Jean*0034 C.. See Bitz, C. M. and W. H. Lipscomb, 1999:  An energy-conserving
                0035 C..       thermodynamic sea ice model for climate study.
                0036 C..       J. Geophys. Res., 104, 15669 - 15677.
7269783f6f Jean*0037 C..     Winton, M., 1999:  "A reformulated three-layer sea ice model."
                0038 C..       Submitted to J. Atmos. Ocean. Technol.
                0039 C.. authors Elizabeth C. Hunke and William Lipscomb
                0040 C..         Fluid Dynamics Group, Los Alamos National Laboratory
                0041 C-----------------------------------------------------------------------
                0042 Cc****subroutine thermo_winton(n,fice,fsnow,dqice,dTsfc)
                0043 C.. Compute temperature change using Winton model with 2 ice layers, of
                0044 C.. which only the top layer has a variable heat capacity.
                0045
87ea84cac6 Jean*0046 C     !USES:
                0047       IMPLICIT NONE
                0048
                0049 C     == Global variables ===
dbce8fc2d4 Jean*0050 #include "EEPARAMS.h"
9dcf02c6ac Jean*0051 #include "SIZE.h"
87ea84cac6 Jean*0052 #include "THSICE_PARAMS.h"
6b47d550f4 Mart*0053 #ifdef ALLOW_AUTODIFF
                0054 # include "THSICE_SIZE.h"
1d053ba71b Jean*0055 c#include "THSICE_VARS.h"
e3123da85c Patr*0056 # ifdef ALLOW_EXF
c1c3d0f9d7 Patr*0057 #  include "EXF_PARAM.h"
9623863f82 Jean*0058 #  include "EXF_CONSTANTS.h"
6b47d550f4 Mart*0059 #  include "EXF_FIELDS.h"
e3123da85c Patr*0060 # endif
d6f06800ae Patr*0061 #endif
6b47d550f4 Mart*0062 #ifdef ALLOW_AUTODIFF_TAMC
                0063 # include "tamc.h"
                0064 #endif
87ea84cac6 Jean*0065
                0066 C     !INPUT/OUTPUT PARAMETERS:
                0067 C     == Routine Arguments ==
7269783f6f Jean*0068 C     bi,bj       :: tile indices
                0069 C     iMin,iMax   :: computation domain: 1rst index range
                0070 C     jMin,jMax   :: computation domain: 2nd  index range
                0071 C     dBugFlag    :: allow to print debugging stuff (e.g. on 1 grid point).
2a9474d935 Mart*0072 C     useBulkForce:: use surf. fluxes from bulk-forcing external S/R
                0073 C     useEXF      :: use surf. fluxes from exf          external S/R
7269783f6f Jean*0074 C---  Input:
9623863f82 Jean*0075 C     icMask      :: sea-ice fractional mask [0-1]
9dcf02c6ac Jean*0076 C     hIce        :: ice height [m]
a85293d087 Mart*0077 C     hSnow       :: snow height [m]
9dcf02c6ac Jean*0078 C     tFrz        :: sea-water freezing temperature [oC] (function of S)
a85293d087 Mart*0079 C     flxExSW     :: Net (minus SW) surf. heat flux (+=down) fct of surf temp Ts
7269783f6f Jean*0080 C                    0: Flx(Ts=0) ; 1: Flx(Ts=Ts^n) ; 2: d.Flx/dTs(Ts=Ts^n)
                0081 C---  Modified (input&output):
9dcf02c6ac Jean*0082 C     flxSW       :: net Short-Wave flux (+=down) [W/m2]: input= at surface
                0083 C                 ::                 output= below sea-ice, into the ocean
a85293d087 Mart*0084 C     tSrf        :: surface (ice or snow) temperature
                0085 C     qIc1        :: ice enthalpy [J/kg], 1rst level
                0086 C     qIc2        :: ice enthalpy [J/kg], 2nd level
7269783f6f Jean*0087 C---  Output
9dcf02c6ac Jean*0088 C     tIc1        :: temperature of ice layer 1 [oC]
                0089 C     tIc2        :: temperature of ice layer 2 [oC]
a85293d087 Mart*0090 C     dTsrf       :: surface temp adjusment: Ts^n+1 - Ts^n
                0091 C     sHeat       :: surf heat flux left to melt snow or ice (=Atmos-conduction)
9dcf02c6ac Jean*0092 C     flxCnB      :: heat flux conducted through the ice to bottom surface
                0093 C     flxAtm      :: net flux of energy from the atmosphere [W/m2] (+=down)
7269783f6f Jean*0094 C                    without snow precip. (energy=0 for liquid water at 0.oC)
9dcf02c6ac Jean*0095 C     evpAtm      :: evaporation to the atmosphere [kg/m2/s] (>0 if evaporate)
7269783f6f Jean*0096 C---  Input:
                0097 C     myTime      :: current Time of simulation [s]
                0098 C     myIter      :: current Iteration number in simulation
                0099 C     myThid      :: my Thread Id number
                0100       INTEGER bi,bj
                0101       INTEGER iMin, iMax
                0102       INTEGER jMin, jMax
                0103       LOGICAL dBugFlag
2a9474d935 Mart*0104       LOGICAL useBulkForce
                0105       LOGICAL useEXF
a85293d087 Mart*0106       _RL icMask (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
6dc8890c80 Patr*0107       _RL hIce   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
a85293d087 Mart*0108       _RL hSnow  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
6dc8890c80 Patr*0109       _RL tFrz   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0110       _RL flxSW  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
a85293d087 Mart*0111       _RL tSrf   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
6dc8890c80 Patr*0112       _RL qIc1   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0113       _RL qIc2   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0114       _RL tIc1   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0115       _RL tIc2   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0116       _RL sHeat  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0117       _RL flxCnB (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0118       _RL flxAtm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0119       _RL evpAtm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
a85293d087 Mart*0120 C
                0121       _RL flxExSW(1:sNx,1:sNy,0:2)
                0122       _RL dTsrf  (1:sNx,1:sNy)
d99870c051 Patr*0123       _RL myTime
7269783f6f Jean*0124       INTEGER myIter
                0125       INTEGER myThid
                0126 CEOP
                0127
                0128 #ifdef ALLOW_THSICE
                0129 C     !LOCAL VARIABLES:
87ea84cac6 Jean*0130 C     == Local Variables ==
9dcf02c6ac Jean*0131 C     useBlkFlx    :: use some bulk-formulae to compute fluxes over sea-ice
                0132 C                  :: (otherwise, fluxes are passed as argument from AIM)
                0133 C     iterate4Tsf  :: to stop to iterate when all icy grid pts Tsf did converged
                0134 C     iceFlag      :: True= do iterate for Surf.Temp ; False= do nothing
98b4e0ca2d Jean*0135 C     frsnow       :: fractional snow cover
a85293d087 Mart*0136 C     fswpen       :: SW penetrating beneath surface [W/m2]
                0137 C     fswdn        :: SW absorbed at surface [W/m2]
                0138 C     fswint       :: SW absorbed in ice [W/m2]
                0139 C     fswocn       :: SW passed through ice to ocean [W/m2]
                0140 C     Tsf_mlt      :: surface temperature of snow or ice at melting point
                0141 C     flx0exSW     :: net (minus SW) surface heat flux over melting snow/ice
                0142 C     flxTexSW     :: net (minus SW) surface heat flux (+=down) [W/m2]
                0143 C     dFlxdT       :: derivative of flxNet wrt tSrf [W/m2/K]
                0144 C     evap_0       :: evaporation over melting snow/ice (>0 evaporate) [kg/m2/s]
                0145 C                     renamed to avoid conflict with EXF evap0 (AUTODIFF)
                0146 C     evapT        :: evaporation over snow/ice (>0 if evaporate) [kg/m2/s]
                0147 C     dEvdT        :: derivative of evap. with respect to tSrf [kg/m2/s/K]
                0148 C     flxNet       :: net surf heat flux (+=down), from atmos to sea-ice [W/m2]
9dcf02c6ac Jean*0149 C     netSW        :: net Short-Wave flux at surface (+=down) [W/m2]
98b4e0ca2d Jean*0150 C     fct          :: heat conducted to top surface
                0151 C     k12, k32     :: thermal conductivity terms
9dcf02c6ac Jean*0152 C     a10,b10,c10  :: coefficients in quadratic eqn for T1
98b4e0ca2d Jean*0153 C     a1, b1, c1   :: coefficients in quadratic eqn for T1
                0154 C     dt           :: timestep
9dcf02c6ac Jean*0155       LOGICAL useBlkFlx
                0156       LOGICAL iterate4Tsf
389b16a470 Jean*0157       INTEGER i, j, k, iterMax
                0158       INTEGER ii, jj, icount, stdUnit
9f5240b52a Jean*0159 #if ( defined ALLOW_DBUG_THSICE || !defined ALLOW_AUTODIFF )
389b16a470 Jean*0160       CHARACTER*(MAX_LEN_MBUF) msgBuf
9f5240b52a Jean*0161 #endif
98b4e0ca2d Jean*0162       _RL  frsnow
                0163       _RL  fswpen
                0164       _RL  fswdn
                0165       _RL  fswint
                0166       _RL  fswocn
c1c3d0f9d7 Patr*0167       _RL  iceFlag (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
a85293d087 Mart*0168       _RL  Tsf_mlt (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
9dcf02c6ac Jean*0169       _RL  flx0exSW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0170       _RL  flxTexSW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0171       _RL  dFlxdT  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
a85293d087 Mart*0172       _RL  evap_0  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
9dcf02c6ac Jean*0173       _RL  evapT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0174       _RL  dEvdT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
170766e9fd Jean*0175       _RL  flxNet
98b4e0ca2d Jean*0176       _RL  fct
9dcf02c6ac Jean*0177       _RL  k12     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0178       _RL  k32
                0179       _RL  a10     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0180       _RL  b10     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0181       _RL  c10     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98b4e0ca2d Jean*0182       _RL  a1, b1, c1
                0183       _RL  dt
72e380ddb5 Jean*0184       _RL  recip_dhSnowLin
a85293d087 Mart*0185       _RL  TsfTmp, tIc2Tmp
9dcf02c6ac Jean*0186 #ifdef ALLOW_DBUG_THSICE
                0187       _RL  netSW
                0188 #endif
7c50f07931 Mart*0189 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0190 C     tkey  :: tape key (depends on tiles)
                0191 C     ikey  :: tape key (depends on iteration index and tiles)
                0192       INTEGER tkey, ikey
7c50f07931 Mart*0193 #endif
7269783f6f Jean*0194
a85293d087 Mart*0195 #ifdef ALLOW_DBUG_THSICE
9dcf02c6ac Jean*0196 C-    Define grid-point location where to print debugging values
7269783f6f Jean*0197 #include "THSICE_DEBUG.h"
9865bb1ac3 Jean*0198  1010 FORMAT(A,I3,3F11.6)
                0199  1020 FORMAT(A,1P4E14.6)
a85293d087 Mart*0200 #endif
                0201 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
87ea84cac6 Jean*0202
d6f06800ae Patr*0203 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0204       tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0205 CADJ STORE flxsw  = comlev1_bibj, key = tkey, kind = isbyte
                0206 CADJ STORE hSnow  = comlev1_bibj, key = tkey, kind = isbyte
6b47d550f4 Mart*0207 #endif /* ALLOW_AUTODIFF_TAMC */
                0208
                0209 #ifdef ALLOW_AUTODIFF
d4b7a0d2bc Patr*0210       DO j = 1-OLy, sNy+OLy
                0211        DO i = 1-OLx, sNx+OLx
389b16a470 Jean*0212         tIc1(i,j) = 0. _d 0
                0213         tIc2(i,j) = 0. _d 0
0707ba3b3d Jean*0214 C-   set these arrays everywhere: overlap are not set and not used,
                0215 C     but some arrays are stored and storage includes overlap.
                0216         flx0exSW(i,j) = 0. _d 0
                0217         flxTexSW(i,j) = 0. _d 0
                0218         dFlxdT  (i,j) = 0. _d 0
a85293d087 Mart*0219         evap_0  (i,j) = 0. _d 0
0707ba3b3d Jean*0220         evapT   (i,j) = 0. _d 0
                0221         dEvdT   (i,j) = 0. _d 0
                0222         iceFlag (i,j) = 0. _d 0
a85293d087 Mart*0223         Tsf_mlt (i,j) = 0. _d 0
389b16a470 Jean*0224        ENDDO
                0225       ENDDO
6dc8890c80 Patr*0226 #endif
                0227
389b16a470 Jean*0228       stdUnit = standardMessageUnit
2a9474d935 Mart*0229       useBlkFlx = useEXF .OR. useBulkForce
9dcf02c6ac Jean*0230       dt  = thSIce_dtTemp
72e380ddb5 Jean*0231       IF ( dhSnowLin.GT.0. _d 0 ) THEN
                0232         recip_dhSnowLin = 1. _d 0 / dhSnowLin
                0233       ELSE
                0234         recip_dhSnowLin = 0. _d 0
                0235       ENDIF
2a9474d935 Mart*0236
9dcf02c6ac Jean*0237       iterate4Tsf = .FALSE.
9c4b3393ef Mart*0238       icount = 0
389b16a470 Jean*0239
7269783f6f Jean*0240       DO j = jMin, jMax
                0241        DO i = iMin, iMax
c1c3d0f9d7 Patr*0242         IF ( icMask(i,j).GT.0. _d 0) THEN
9dcf02c6ac Jean*0243           iterate4Tsf  = .TRUE.
c1c3d0f9d7 Patr*0244           iceFlag(i,j) = 1. _d 0
7269783f6f Jean*0245 #ifdef ALLOW_DBUG_THSICE
389b16a470 Jean*0246           IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,'(A,2I4,2I2)')
7269783f6f Jean*0247      &         'ThSI_SOLVE4T: i,j=',i,j,bi,bj
                0248 #endif
9dcf02c6ac Jean*0249           IF ( hIce(i,j).LT.hIceMin ) THEN
9c4b3393ef Mart*0250 C     if hi < hIceMin, melt the ice.
                0251 C     keep the position of this problem but do the stop later
                0252            ii = i
                0253            jj = j
                0254            icount = icount + 1
77008a74ba Patr*0255           ENDIF
87ea84cac6 Jean*0256
9dcf02c6ac Jean*0257 C--   Fractional snow cover:
                0258 C     assume a linear distribution of snow thickness, with dhSnowLin slope,
                0259 C      from hs-dhSnowLin to hs+dhSnowLin if full ice & snow cover.
                0260 C     frsnow = fraction of snow over the ice-covered part of the grid cell
a85293d087 Mart*0261           IF ( hSnow(i,j) .GT. icMask(i,j)*dhSnowLin ) THEN
9dcf02c6ac Jean*0262            frsnow = 1. _d 0
                0263           ELSE
a85293d087 Mart*0264            frsnow = hSnow(i,j)*recip_dhSnowLin/icMask(i,j)
9dcf02c6ac Jean*0265            IF ( frsnow.GT.0. _d 0 ) frsnow = SQRT(frsnow)
                0266           ENDIF
87ea84cac6 Jean*0267
9dcf02c6ac Jean*0268 C--   Compute SW flux absorbed at surface and penetrating to layer 1.
                0269           fswpen = flxSW(i,j) * (1. _d 0 - frsnow) * i0swFrac
                0270           fswocn = fswpen * exp(-ksolar*hIce(i,j))
                0271           fswint = fswpen - fswocn
                0272           fswdn  = flxSW(i,j) - fswpen
                0273
                0274 C     Initialise Atmospheric surf. heat flux with net SW flux at the surface
                0275           flxAtm(i,j) = flxSW(i,j)
                0276 C     SW flux at sea-ice base left to the ocean
                0277           flxSW(i,j) = fswocn
                0278 C     Initialise surface Heating with SW contribution
                0279           sHeat(i,j) = fswdn
                0280
                0281 C--   Compute conductivity terms at layer interfaces.
                0282           k12(i,j) = 4. _d 0*kIce*kSnow
a85293d087 Mart*0283      &             / (kSnow*hIce(i,j) + 4. _d 0*kIce*hSnow(i,j))
9dcf02c6ac Jean*0284           k32      = 2. _d 0*kIce  / hIce(i,j)
                0285
                0286 C--   Compute ice temperatures
                0287           a1 = cpIce
                0288           b1 = qIc1(i,j) + (cpWater-cpIce )*Tmlt1 - Lfresh
                0289           c1 = Lfresh * Tmlt1
                0290           tIc1(i,j) = 0.5 _d 0 *(-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/a1
                0291           tIc2(i,j) = (Lfresh-qIc2(i,j)) / cpIce
                0292
9c4b3393ef Mart*0293 #ifdef ALLOW_DBUG_THSICE
9dcf02c6ac Jean*0294           IF (tIc1(i,j).GT.0. _d 0 ) THEN
389b16a470 Jean*0295            WRITE(stdUnit,'(A,I12,1PE14.6)')
9dcf02c6ac Jean*0296      &       ' BBerr: Tice(1) > 0 ; it=', myIter, qIc1(i,j)
389b16a470 Jean*0297            WRITE(stdUnit,'(A,4I5,2F11.4)')
9dcf02c6ac Jean*0298      &      ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,tIc1(i,j),tIc2(i,j)
                0299           ENDIF
                0300           IF ( tIc2(i,j).GT.0. _d 0) THEN
389b16a470 Jean*0301            WRITE(stdUnit,'(A,I12,1PE14.6)')
9dcf02c6ac Jean*0302      &       ' BBerr: Tice(2) > 0 ; it=', myIter, qIc2(i,j)
389b16a470 Jean*0303            WRITE(stdUnit,'(A,4I5,2F11.4)')
9dcf02c6ac Jean*0304      &      ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,tIc1(i,j),tIc2(i,j)
                0305           ENDIF
389b16a470 Jean*0306           IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1010)
a85293d087 Mart*0307      &     'ThSI_SOLVE4T: k, Ts, Tice=',0,tSrf(i,j),tIc1(i,j),tIc2(i,j)
7269783f6f Jean*0308 #endif
87ea84cac6 Jean*0309
9dcf02c6ac Jean*0310 C--   Compute coefficients used in quadratic formula.
87ea84cac6 Jean*0311
9dcf02c6ac Jean*0312           a10(i,j) = rhoi*cpIce *hIce(i,j)/(2. _d 0*dt) +
                0313      &          k32*( 4. _d 0*dt*k32 + rhoi*cpIce *hIce(i,j) )
                0314      &           / ( 6. _d 0*dt*k32 + rhoi*cpIce *hIce(i,j) )
                0315           b10(i,j) = -hIce(i,j)*
                0316      &          ( rhoi*cpIce*tIc1(i,j) + rhoi*Lfresh*Tmlt1/tIc1(i,j) )
                0317      &           /(2. _d 0*dt)
                0318      &        - k32*( 4. _d 0*dt*k32*tFrz(i,j)
                0319      &               +rhoi*cpIce*hIce(i,j)*tIc2(i,j) )
                0320      &           / ( 6. _d 0*dt*k32 + rhoi*cpIce *hIce(i,j) )
                0321      &        - fswint
                0322           c10(i,j) = rhoi*Lfresh*hIce(i,j)*Tmlt1 / (2. _d 0*dt)
                0323
                0324         ELSE
c1c3d0f9d7 Patr*0325           iceFlag(i,j) = 0. _d 0
9dcf02c6ac Jean*0326         ENDIF
                0327        ENDDO
                0328       ENDDO
c1c3d0f9d7 Patr*0329 #ifndef ALLOW_AUTODIFF
9c4b3393ef Mart*0330       IF ( icount .gt. 0 ) THEN
389b16a470 Jean*0331        WRITE(stdUnit,'(A,I5,A)')
9c4b3393ef Mart*0332      &      'THSICE_SOLVE4TEMP: there are ',icount,
                0333      &      ' case(s) where hIce<hIceMin;'
389b16a470 Jean*0334        WRITE(stdUnit,'(A,I3,A1,I3,A)')
9c4b3393ef Mart*0335      &      'THSICE_SOLVE4TEMP: the last one was at (',ii,',',jj,
486e9de820 Jean*0336      &      ') with hIce = ', hIce(ii,jj)
389b16a470 Jean*0337        WRITE( msgBuf, '(A)')
                0338      &      'THSICE_SOLVE4TEMP: should not enter if hIce<hIceMin'
                0339        CALL PRINT_ERROR( msgBuf , myThid )
                0340        STOP 'ABNORMAL END: S/R THSICE_SOLVE4TEMP'
9c4b3393ef Mart*0341       ENDIF
c1c3d0f9d7 Patr*0342 #endif /* ALLOW_AUTODIFF */
87ea84cac6 Jean*0343
a85293d087 Mart*0344 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
6dc8890c80 Patr*0345
a85293d087 Mart*0346 C--   Get surface fluxes over melting surface (assumed to be @ 0 oC)
c1c3d0f9d7 Patr*0347 #ifndef ALLOW_AUTODIFF
9dcf02c6ac Jean*0348       IF ( useBlkFlx .AND. iterate4Tsf  ) THEN
                0349 #endif
                0350         DO j = jMin, jMax
                0351          DO i = iMin, iMax
a85293d087 Mart*0352            Tsf_mlt(i,j) = 0.
9dcf02c6ac Jean*0353          ENDDO
                0354         ENDDO
c1c3d0f9d7 Patr*0355 #ifndef ALLOW_AUTODIFF
9dcf02c6ac Jean*0356         IF ( useEXF ) THEN
c1c3d0f9d7 Patr*0357 #endif
a163d11794 Patr*0358            k = 1
                0359            CALL THSICE_GET_EXF(   bi, bj, k,
9dcf02c6ac Jean*0360      I                            iMin, iMax, jMin, jMax,
a85293d087 Mart*0361      I                            iceFlag, hSnow, Tsf_mlt,
                0362      O                            flx0exSW, dFlxdT, evap_0, dEvdT,
9dcf02c6ac Jean*0363      I                            myTime, myIter, myThid )
9623863f82 Jean*0364 #ifndef ALLOW_AUTODIFF
9dcf02c6ac Jean*0365 C-    could add this "ifdef" to hide THSICE_GET_BULKF from TAF
                0366         ELSEIF ( useBulkForce ) THEN
6dc8890c80 Patr*0367            CALL THSICE_GET_BULKF( bi, bj,
9dcf02c6ac Jean*0368      I                            iMin, iMax, jMin, jMax,
a85293d087 Mart*0369      I                            iceFlag, hSnow, Tsf_mlt,
                0370      O                            flx0exSW, dFlxdT, evap_0, dEvdT,
9dcf02c6ac Jean*0371      I                            myTime, myIter, myThid )
c1c3d0f9d7 Patr*0372 C--- end if: IF ( useEXF ) THEN
9dcf02c6ac Jean*0373         ENDIF
c1c3d0f9d7 Patr*0374 C--- end if: IF ( useBlkFlx .AND. iterate4Tsf  ) THEN
9dcf02c6ac Jean*0375       ENDIF
c1c3d0f9d7 Patr*0376 #endif
170766e9fd Jean*0377
a85293d087 Mart*0378 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0379 C---  Compute new surface and ice temperatures; iterate until tSrf converges.
9dcf02c6ac Jean*0380       DO j = jMin, jMax
                0381         DO i = iMin, iMax
a85293d087 Mart*0382           dTsrf(i,j) = Terrmax
9dcf02c6ac Jean*0383         ENDDO
                0384       ENDDO
9a68c0a761 Jean*0385       IF ( useBlkFlx ) THEN
                0386         iterMax = nitMaxTsf
                0387       ELSE
                0388         iterMax = 1
                0389       ENDIF
                0390
87ea84cac6 Jean*0391 C ----- begin iteration  -----
9a68c0a761 Jean*0392       DO k = 1,iterMax
c1c3d0f9d7 Patr*0393 #ifndef ALLOW_AUTODIFF
9dcf02c6ac Jean*0394        IF ( iterate4Tsf ) THEN
6b47d550f4 Mart*0395         iterate4Tsf = .FALSE.
c1c3d0f9d7 Patr*0396 C
6b47d550f4 Mart*0397         IF ( useBlkFlx ) THEN
                0398 #endif
                0399 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0400          ikey = (tkey-1)*MaxTsf + k
                0401 CADJ STORE iceflag  = comlev1_thsice_s4t, key = ikey, kind = isbyte
                0402 CADJ STORE tSrf     = comlev1_thsice_s4t, key = ikey, kind = isbyte
                0403 CADJ STORE flxTexSW = comlev1_thsice_s4t, key = ikey, kind = isbyte
                0404 CADJ STORE dFlxdT   = comlev1_thsice_s4t, key = ikey, kind = isbyte
                0405 CADJ STORE dEvdT    = comlev1_thsice_s4t, key = ikey, kind = isbyte
6b47d550f4 Mart*0406 #endif
c1c3d0f9d7 Patr*0407
9dcf02c6ac Jean*0408 C--   Compute top surface flux.
c1c3d0f9d7 Patr*0409 #ifndef ALLOW_AUTODIFF
170766e9fd Jean*0410          IF ( useEXF ) THEN
c1c3d0f9d7 Patr*0411 #endif
a163d11794 Patr*0412            CALL THSICE_GET_EXF(   bi, bj, k+1,
9dcf02c6ac Jean*0413      I                            iMin, iMax, jMin, jMax,
a85293d087 Mart*0414      I                            iceFlag, hSnow, tSrf,
9dcf02c6ac Jean*0415      O                            flxTexSW, dFlxdT, evapT, dEvdT,
                0416      I                            myTime, myIter, myThid )
                0417 C-    could add this "ifdef" to hide THSICE_GET_BULKF from TAF
9623863f82 Jean*0418 #ifndef ALLOW_AUTODIFF
170766e9fd Jean*0419          ELSEIF ( useBulkForce ) THEN
6dc8890c80 Patr*0420            CALL THSICE_GET_BULKF( bi, bj,
9dcf02c6ac Jean*0421      I                            iMin, iMax, jMin, jMax,
a85293d087 Mart*0422      I                            iceFlag, hSnow, tSrf,
9dcf02c6ac Jean*0423      O                            flxTexSW, dFlxdT, evapT, dEvdT,
                0424      I                            myTime, myIter, myThid )
c1c3d0f9d7 Patr*0425 C--- end if: IF ( useEXF ) THEN
2a9474d935 Mart*0426          ENDIF
6b47d550f4 Mart*0427         ELSE
9dcf02c6ac Jean*0428          DO j = jMin, jMax
                0429           DO i = iMin, iMax
c1c3d0f9d7 Patr*0430            IF ( iceFlag(i,j).GT.0. _d 0 ) THEN
9dcf02c6ac Jean*0431              flxTexSW(i,j) = flxExSW(i,j,1)
                0432              dFlxdT(i,j) = flxExSW(i,j,2)
                0433            ENDIF
                0434           ENDDO
                0435          ENDDO
c1c3d0f9d7 Patr*0436 C--- end if: IF ( useBlkFlx ) THEN
6b47d550f4 Mart*0437         ENDIF
c1c3d0f9d7 Patr*0438 #endif /* ALLOW_AUTODIFF */
9dcf02c6ac Jean*0439
a85293d087 Mart*0440 #ifdef ALLOW_AUTODIFF_TAMC_MORE
                0441 C     avoid calling thsice_get_exf in ad-code
edb6656069 Mart*0442 CADJ STORE flxTexSW = comlev1_thsice_s4t, key = ikey, kind = isbyte
                0443 CADJ STORE dFlxdT   = comlev1_thsice_s4t, key = ikey, kind = isbyte
1818702d6f Patr*0444 #endif
9dcf02c6ac Jean*0445 C--   Compute new top layer and surface temperatures.
a85293d087 Mart*0446 C     If tSrf is computed to be > 0 C, fix tSrf = 0 and recompute tIc1
9dcf02c6ac Jean*0447 C     with different coefficients.
6b47d550f4 Mart*0448         DO j = jMin, jMax
                0449          DO i = iMin, iMax
                0450           IF ( iceFlag(i,j).GT.0. _d 0 ) THEN
                0451            flxNet = sHeat(i,j) + flxTexSW(i,j)
7269783f6f Jean*0452 #ifdef ALLOW_DBUG_THSICE
6b47d550f4 Mart*0453            IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1020)
9dcf02c6ac Jean*0454      &     'ThSI_SOLVE4T: flxNet,dFlxdT,k12,D=',
                0455      &      flxNet, dFlxdT(i,j), k12(i,j), k12(i,j)-dFlxdT(i,j)
7269783f6f Jean*0456 #endif
87ea84cac6 Jean*0457
6b47d550f4 Mart*0458            a1 = a10(i,j) - k12(i,j)*dFlxdT(i,j) / (k12(i,j)-dFlxdT(i,j))
a85293d087 Mart*0459            b1 = b10(i,j) - k12(i,j)*(flxNet-dFlxdT(i,j)*tSrf(i,j))
6b47d550f4 Mart*0460      &                   /(k12(i,j)-dFlxdT(i,j))
                0461            c1 = c10(i,j)
                0462            tIc1(i,j)  = -(b1 + SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)
a85293d087 Mart*0463            dTsrf(i,j) = (flxNet + k12(i,j)*(tIc1(i,j)-tSrf(i,j)))
6b47d550f4 Mart*0464      &                   /(k12(i,j)-dFlxdT(i,j))
a85293d087 Mart*0465 C     this line causes 3 recomputation warings for TAF because tSrf is
                0466 C     updated and the updated value is used to update tSrf (and other
                0467 C     variables) again; with this intermediate variable, the warnings go
                0468 C     away.
                0469 C          tSrf(i,j) = tSrf(i,j) + dTsrf(i,j)
                0470            TsfTmp = tSrf(i,j) + dTsrf(i,j)
c1c3d0f9d7 Patr*0471 C
a85293d087 Mart*0472            IF ( TsfTmp .GT. 0. _d 0 ) THEN
7269783f6f Jean*0473 #ifdef ALLOW_DBUG_THSICE
6b47d550f4 Mart*0474             IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1010)
a85293d087 Mart*0475      &     'ThSI_SOLVE4T: k,ts,t1,dTs=',k,TsfTmp,tIc1(i,j),dTsrf(i,j)
7269783f6f Jean*0476 #endif
6b47d550f4 Mart*0477             a1 = a10(i,j) + k12(i,j)
98b4e0ca2d Jean*0478 C      note: b1 = b10 - k12*Tf0
6b47d550f4 Mart*0479             b1 = b10(i,j)
                0480             tIc1(i,j) = (-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)
a85293d087 Mart*0481             tSrf(i,j) = 0. _d 0
c1c3d0f9d7 Patr*0482 #ifndef ALLOW_AUTODIFF
6b47d550f4 Mart*0483             IF ( useBlkFlx ) THEN
c1c3d0f9d7 Patr*0484 #endif /* ALLOW_AUTODIFF */
6b47d550f4 Mart*0485              flxTexSW(i,j) = flx0exSW(i,j)
a85293d087 Mart*0486              evapT(i,j) = evap_0(i,j)
                0487              dTsrf(i,j) = 0. _d 0
c1c3d0f9d7 Patr*0488 #ifndef ALLOW_AUTODIFF
6b47d550f4 Mart*0489             ELSE
                0490              flxTexSW(i,j) = flxExSW(i,j,0)
a85293d087 Mart*0491              dTsrf(i,j) = 1000.
6b47d550f4 Mart*0492              dFlxdT(i,j) = 0.
                0493             ENDIF
c1c3d0f9d7 Patr*0494 #endif /* ALLOW_AUTODIFF */
a85293d087 Mart*0495            ELSE
                0496             tSrf(i,j) = TsfTmp
6b47d550f4 Mart*0497            ENDIF
9dcf02c6ac Jean*0498
                0499 C--   Check for convergence.  If no convergence, then repeat.
a85293d087 Mart*0500            IF (ABS(dTsrf(i,j)).GE.Terrmax) THEN
6b47d550f4 Mart*0501             iceFlag(i,j) = 1. _d 0
                0502            ELSE
                0503             iceFlag(i,j) = 0. _d 0
                0504            ENDIF
a85293d087 Mart*0505 #ifndef ALLOW_AUTODIFF
6b47d550f4 Mart*0506            iterate4Tsf = iterate4Tsf .OR. (iceFlag(i,j).GT.0. _d 0)
a85293d087 Mart*0507 #endif
87ea84cac6 Jean*0508
a85293d087 Mart*0509 C     Convergence test: Make sure tSrf has converged, within prescribed error.
9dcf02c6ac Jean*0510 C     (Energy conservation is guaranteed within machine roundoff, even
a85293d087 Mart*0511 C      if tSrf has not converged.)
9dcf02c6ac Jean*0512 C     If no convergence, then repeat.
87ea84cac6 Jean*0513
7269783f6f Jean*0514 #ifdef ALLOW_DBUG_THSICE
6b47d550f4 Mart*0515            IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1010)
a85293d087 Mart*0516      &  'ThSI_SOLVE4T: k,ts,t1,dTs=', k,tSrf(i,j),tIc1(i,j),dTsrf(i,j)
6b47d550f4 Mart*0517            IF ( useBlkFlx .AND. k.EQ.nitMaxTsf .AND.
                0518      &                     (iceFlag(i,j).GT.0. _d 0) ) THEN
                0519             WRITE(stdUnit,'(A,4I4,I12,F15.9)')
9dcf02c6ac Jean*0520      &       ' BB: not converge: i,j,it,hi=',i,j,bi,bj,myIter,hIce(i,j)
6b47d550f4 Mart*0521             WRITE(stdUnit,*)
a85293d087 Mart*0522      &        'BB: not converge: tSrf, dTsrf=', tSrf(i,j), dTsrf(i,j)
6b47d550f4 Mart*0523             WRITE(stdUnit,*)
389b16a470 Jean*0524      &        'BB: not converge: flxNet,dFlxT=', flxNet, dFlxdT(i,j)
a85293d087 Mart*0525             IF ( tSrf(i,j).LT.-70. _d 0 ) THEN
389b16a470 Jean*0526              WRITE( msgBuf, '(A,2I4,2I3,I10,F12.3)')
a85293d087 Mart*0527      &        'THSICE_SOLVE4TEMP: Too low tSrf in', i, j, bi, bj,
                0528      &                            myIter, tSrf(i,j)
389b16a470 Jean*0529              CALL PRINT_ERROR( msgBuf , myThid )
                0530              STOP 'ABNORMAL END: S/R THSICE_SOLVE4TEMP'
6b47d550f4 Mart*0531             ENDIF
389b16a470 Jean*0532            ENDIF
9623863f82 Jean*0533 #endif /* ALLOW_DBUG_THSICE */
9dcf02c6ac Jean*0534
6b47d550f4 Mart*0535           ENDIF
                0536          ENDDO
9dcf02c6ac Jean*0537         ENDDO
c1c3d0f9d7 Patr*0538 #ifndef ALLOW_AUTODIFF
                0539 C--- end if: IF ( iterate4Tsf ) THEN
9a68c0a761 Jean*0540        ENDIF
c1c3d0f9d7 Patr*0541 #endif /* ALLOW_AUTODIFF */
                0542 C--- end loop DO k = 1,iterMax
9a68c0a761 Jean*0543       ENDDO
87ea84cac6 Jean*0544 C ------ end iteration ------------
a85293d087 Mart*0545 #ifdef ALLOW_AUTODIFF_TAMC_MORE
                0546 C     avoid calling thsice_get_exf in ad-code
edb6656069 Mart*0547 CADJ STORE tIc1   = comlev1_bibj, key = tkey, kind = isbyte
                0548 CADJ STORE dTsrf  = comlev1_bibj, key = tkey, kind = isbyte
                0549 CADJ STORE tSrf   = comlev1_bibj, key = tkey, kind = isbyte
                0550 CADJ STORE dFlxdT = comlev1_bibj, key = tkey, kind = isbyte
                0551 CADJ STORE dEvdT  = comlev1_bibj, key = tkey, kind = isbyte
a85293d087 Mart*0552 #endif
87ea84cac6 Jean*0553
a85293d087 Mart*0554 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
87ea84cac6 Jean*0555
9dcf02c6ac Jean*0556       DO j = jMin, jMax
                0557        DO i = iMin, iMax
a85293d087 Mart*0558 C     this intermediate variable avoids a recomputation warning
                0559         tIc2Tmp = tIc2(i,j)
c1c3d0f9d7 Patr*0560         IF ( icMask(i,j).GT.0. _d 0 ) THEN
9dcf02c6ac Jean*0561
                0562 C--   Compute new bottom layer temperature.
                0563           k32       = 2. _d 0*kIce  / hIce(i,j)
                0564           tIc2(i,j) = ( 2. _d 0*dt*k32*(tIc1(i,j)+2. _d 0*tFrz(i,j))
a85293d087 Mart*0565      &                 + rhoi*cpIce*hIce(i,j)*tIc2Tmp)
9dcf02c6ac Jean*0566      &               /(6. _d 0*dt*k32 + rhoi*cpIce*hIce(i,j))
7269783f6f Jean*0567 #ifdef ALLOW_DBUG_THSICE
389b16a470 Jean*0568           IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1010)
a85293d087 Mart*0569      &  'ThSI_SOLVE4T: k, Ts, Tice=',k,tSrf(i,j),tIc1(i,j),tIc2(i,j)
9dcf02c6ac Jean*0570           netSW   = flxAtm(i,j)
7269783f6f Jean*0571 #endif
87ea84cac6 Jean*0572
9dcf02c6ac Jean*0573 C--   Compute final flux values at surfaces.
a85293d087 Mart*0574           fct    = k12(i,j)*(tSrf(i,j)-tIc1(i,j))
9dcf02c6ac Jean*0575           flxCnB(i,j) = 4. _d 0*kIce *(tIc2(i,j)-tFrz(i,j))/hIce(i,j)
                0576           flxNet = sHeat(i,j) + flxTexSW(i,j)
a85293d087 Mart*0577           flxNet = flxNet + dFlxdT(i,j)*dTsrf(i,j)
c1c3d0f9d7 Patr*0578 #ifndef ALLOW_AUTODIFF
9dcf02c6ac Jean*0579           IF ( useBlkFlx ) THEN
c1c3d0f9d7 Patr*0580 #endif
a85293d087 Mart*0581 C--   need to update also Evap (tSrf changes) since Latent heat has
                0582 C     been updated
                0583             evpAtm(i,j) = evapT(i,j) + dEvdT(i,j)*dTsrf(i,j)
c1c3d0f9d7 Patr*0584 #ifndef ALLOW_AUTODIFF
9dcf02c6ac Jean*0585           ELSE
9865bb1ac3 Jean*0586 C- WARNING: Evap & +Evap*Lfresh are missing ! (but only affects Diagnostics)
9dcf02c6ac Jean*0587             evpAtm(i,j) = 0.
                0588           ENDIF
c1c3d0f9d7 Patr*0589 #endif
9dcf02c6ac Jean*0590 C-    Update energy flux to Atmos with other than SW contributions;
9865bb1ac3 Jean*0591 C     use latent heat = Lvap (energy=0 for liq. water at 0.oC)
9dcf02c6ac Jean*0592           flxAtm(i,j) = flxAtm(i,j) + flxTexSW(i,j)
a85293d087 Mart*0593      &                + dFlxdT(i,j)*dTsrf(i,j) + evpAtm(i,j)*Lfresh
9865bb1ac3 Jean*0594 C-    excess of energy @ surface (used for surface melting):
9dcf02c6ac Jean*0595           sHeat(i,j) = flxNet - fct
87ea84cac6 Jean*0596
7269783f6f Jean*0597 #ifdef ALLOW_DBUG_THSICE
389b16a470 Jean*0598           IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1020)
9dcf02c6ac Jean*0599      &     'ThSI_SOLVE4T: flxNet,fct,Dif,flxCnB=',
                0600      &                    flxNet,fct,flxNet-fct,flxCnB(i,j)
7269783f6f Jean*0601 #endif
87ea84cac6 Jean*0602
9dcf02c6ac Jean*0603 C--   Compute new enthalpy for each layer.
                0604           qIc1(i,j) = -cpWater*Tmlt1 + cpIce *(Tmlt1-tIc1(i,j))
                0605      &                + Lfresh*(1. _d 0-Tmlt1/tIc1(i,j))
                0606           qIc2(i,j) = -cpIce *tIc2(i,j) + Lfresh
87ea84cac6 Jean*0607
9c4b3393ef Mart*0608 #ifdef ALLOW_DBUG_THSICE
9dcf02c6ac Jean*0609 C--   Make sure internal ice temperatures do not exceed Tmlt.
                0610 C     (This should not happen for reasonable values of i0swFrac)
                0611           IF (tIc1(i,j) .GE. Tmlt1) THEN
389b16a470 Jean*0612            WRITE(stdUnit,'(A,2I4,2I3,1P2E14.6)')
9dcf02c6ac Jean*0613      &     ' BBerr - Bug: IceT(1) > Tmlt',i,j,bi,bj,tIc1(i,j),Tmlt1
                0614           ENDIF
                0615           IF (tIc2(i,j) .GE. 0. _d 0) THEN
389b16a470 Jean*0616            WRITE(stdUnit,'(A,2I4,2I3,1P2E14.6)')
9dcf02c6ac Jean*0617      &     ' BBerr - Bug: IceT(2) > 0',i,j,bi,bj,tIc2(i,j)
                0618           ENDIF
87ea84cac6 Jean*0619
7269783f6f Jean*0620           IF ( dBug(i,j,bi,bj) ) THEN
a85293d087 Mart*0621            WRITE(stdUnit,1020) 'ThSI_SOLV_4T: tSrf,Tice(1,2), dTsurf=',
                0622      &           tSrf(i,j), tIc1(i,j), tIc2(i,j), dTsrf(i,j)
389b16a470 Jean*0623            WRITE(stdUnit,1020) 'ThSI_SOLV_4T: sHeat, flxCndBt, Qice =',
9dcf02c6ac Jean*0624      &           sHeat(i,j), flxCnB(i,j), qIc1(i,j), qIc2(i,j)
389b16a470 Jean*0625            WRITE(stdUnit,1020) 'ThSI_SOLV_4T: flxA, evpA, fxSW_bf,af=',
7269783f6f Jean*0626      &           flxAtm(i,j), evpAtm(i,j), netSW, flxSW(i,j)
                0627           ENDIF
                0628 #endif
9dcf02c6ac Jean*0629
7269783f6f Jean*0630         ELSE
9dcf02c6ac Jean*0631 C--     ice-free grid point:
                0632 c         tIc1  (i,j) = 0. _d 0
                0633 c         tIc2  (i,j) = 0. _d 0
a85293d087 Mart*0634           dTsrf (i,j) = 0. _d 0
9dcf02c6ac Jean*0635 c         sHeat (i,j) = 0. _d 0
                0636 c         flxCnB(i,j) = 0. _d 0
                0637 c         flxAtm(i,j) = 0. _d 0
                0638 c         evpAtm(i,j) = 0. _d 0
                0639
7269783f6f Jean*0640         ENDIF
                0641        ENDDO
                0642       ENDDO
87ea84cac6 Jean*0643 #endif  /* ALLOW_THSICE */
                0644
a85293d087 Mart*0645 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
87ea84cac6 Jean*0646
                0647       RETURN
                0648       END