C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_solve4temp.F,v 1.37 2013/08/22 02:10:48 jmc Exp $
C $Name: $
#include "THSICE_OPTIONS.h"
#ifdef ALLOW_AUTODIFF_TAMC
# ifdef ALLOW_EXF
# include "EXF_OPTIONS.h"
# endif
#endif
CBOP
C !ROUTINE: THSICE_SOLVE4TEMP
C !INTERFACE:
SUBROUTINE THSICE_SOLVE4TEMP(
I bi, bj,
I iMin,iMax, jMin,jMax, dBugFlag,
I useBulkForce, useEXF,
I icMask, hIce, hSnow1, tFrz, flxExSW,
U flxSW, tSrf1, qIc1, qIc2,
O tIc1, tIc2, dTsrf1,
O sHeat, flxCnB, flxAtm, evpAtm,
I myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | S/R THSICE_SOLVE4TEMP
C *==========================================================*
C | Solve (implicitly) for sea-ice and surface temperature
C *==========================================================*
C \ev
C ADAPTED FROM:
C LANL CICE.v2.0.2
C-----------------------------------------------------------------------
C.. thermodynamics (vertical physics) based on M. Winton 3-layer model
C.. See Bitz, C. M. and W. H. Lipscomb, 1999: An energy-conserving
C.. thermodynamic sea ice model for climate study.
C.. J. Geophys. Res., 104, 15669 - 15677.
C.. Winton, M., 1999: "A reformulated three-layer sea ice model."
C.. Submitted to J. Atmos. Ocean. Technol.
C.. authors Elizabeth C. Hunke and William Lipscomb
C.. Fluid Dynamics Group, Los Alamos National Laboratory
C-----------------------------------------------------------------------
Cc****subroutine thermo_winton(n,fice,fsnow,dqice,dTsfc)
C.. Compute temperature change using Winton model with 2 ice layers, of
C.. which only the top layer has a variable heat capacity.
C !USES:
IMPLICIT NONE
C == Global variables ===
#include "EEPARAMS.h"
#include "SIZE.h"
#include "THSICE_PARAMS.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
# include "tamc_keys.h"
#include "THSICE_SIZE.h"
c#include "THSICE_VARS.h"
# ifdef ALLOW_EXF
# include "EXF_FIELDS.h"
# include "EXF_PARAM.h"
# include "EXF_CONSTANTS.h"
# endif
#endif
C !INPUT/OUTPUT PARAMETERS:
C == Routine Arguments ==
C bi,bj :: tile indices
C iMin,iMax :: computation domain: 1rst index range
C jMin,jMax :: computation domain: 2nd index range
C dBugFlag :: allow to print debugging stuff (e.g. on 1 grid point).
C useBulkForce:: use surf. fluxes from bulk-forcing external S/R
C useEXF :: use surf. fluxes from exf external S/R
C--- Input:
C icMask :: sea-ice fractional mask [0-1]
C hIce :: ice height [m]
C hSnow1 :: snow height [m]
C tFrz :: sea-water freezing temperature [oC] (function of S)
C flxExSW :: surf. heat flux (+=down) except SW, function of surf. temp Ts:
C 0: Flx(Ts=0) ; 1: Flx(Ts=Ts^n) ; 2: d.Flx/dTs(Ts=Ts^n)
C--- Modified (input&output):
C flxSW :: net Short-Wave flux (+=down) [W/m2]: input= at surface
C :: output= below sea-ice, into the ocean
C tSrf1 :: surface (ice or snow) temperature
C qIc1 :: ice enthalpy (J/kg), 1rst level
C qIc2 :: ice enthalpy (J/kg), 2nd level
C--- Output
C tIc1 :: temperature of ice layer 1 [oC]
C tIc2 :: temperature of ice layer 2 [oC]
C dTsrf1 :: surf. temp adjusment: Ts^n+1 - Ts^n
C sHeat :: surf heating flux left to melt snow or ice (= Atmos-conduction)
C flxCnB :: heat flux conducted through the ice to bottom surface
C flxAtm :: net flux of energy from the atmosphere [W/m2] (+=down)
C without snow precip. (energy=0 for liquid water at 0.oC)
C evpAtm :: evaporation to the atmosphere [kg/m2/s] (>0 if evaporate)
C--- Input:
C myTime :: current Time of simulation [s]
C myIter :: current Iteration number in simulation
C myThid :: my Thread Id number
INTEGER bi,bj
INTEGER iMin, iMax
INTEGER jMin, jMax
LOGICAL dBugFlag
LOGICAL useBulkForce
LOGICAL useEXF
_RL icMask(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL hIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL hSnow1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL tFrz (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL tSrf1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL qIc1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL qIc2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL tIc1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL tIc2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL sHeat (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL flxCnB (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL flxAtm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL evpAtm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL flxExSW(iMin:iMax,jMin:jMax,0:2)
_RL dTsrf1 (iMin:iMax,jMin:jMax)
_RL myTime
INTEGER myIter
INTEGER myThid
CEOP
#ifdef ALLOW_THSICE
C !LOCAL VARIABLES:
C == Local Variables ==
C useBlkFlx :: use some bulk-formulae to compute fluxes over sea-ice
C :: (otherwise, fluxes are passed as argument from AIM)
C iterate4Tsf :: to stop to iterate when all icy grid pts Tsf did converged
C iceFlag :: True= do iterate for Surf.Temp ; False= do nothing
C frsnow :: fractional snow cover
C fswpen :: SW penetrating beneath surface (W m-2)
C fswdn :: SW absorbed at surface (W m-2)
C fswint :: SW absorbed in ice (W m-2)
C fswocn :: SW passed through ice to ocean (W m-2)
C Tsf :: surface (ice or snow) temperature (local copy of tSrf1)
C flx0exSW :: net surface heat flux over melting snow/ice, except short-wave.
C flxTexSW :: net surface heat flux, except short-wave (W/m2)
C dFlxdT :: deriv of flxNet wrt Tsf (W m-2 deg-1)
C evap00 :: evaporation over melting snow/ice [kg/m2/s] (>0 if evaporate)
C renamed to evap00 because TAF confuses symbol with EXF evap0
C evapT :: evaporation over snow/ice [kg/m2/s] (>0 if evaporate)
C dEvdT :: derivative of evap. with respect to Tsf [kg/m2/s/K]
C flxNet :: net surf heat flux (+=down), from Atmos. to sea-ice (W m-2)
C netSW :: net Short-Wave flux at surface (+=down) [W/m2]
C fct :: heat conducted to top surface
C k12, k32 :: thermal conductivity terms
C a10,b10,c10 :: coefficients in quadratic eqn for T1
C a1, b1, c1 :: coefficients in quadratic eqn for T1
C dt :: timestep
LOGICAL useBlkFlx
LOGICAL iterate4Tsf
INTEGER i, j, k, iterMax
INTEGER ii, jj, icount, stdUnit
CHARACTER*(MAX_LEN_MBUF) msgBuf
_RL frsnow
_RL fswpen
_RL fswdn
_RL fswint
_RL fswocn
_RL iceFlag (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL Tsf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL flx0exSW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL flxTexSW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL dFlxdT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL evap00 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL evapT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL dEvdT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL flxNet
_RL fct
_RL k12 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL k32
_RL a10 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL b10 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL c10 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL a1, b1, c1
_RL dt
_RL recip_dhSnowLin
#ifdef ALLOW_DBUG_THSICE
_RL netSW
#endif
C- Define grid-point location where to print debugging values
#include "THSICE_DEBUG.h"
1010 FORMAT(A,I3,3F11.6)
1020 FORMAT(A,1P4E14.6)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
#ifdef ALLOW_AUTODIFF_TAMC
act1 = bi - myBxLo(myThid)
max1 = myBxHi(myThid) - myBxLo(myThid) + 1
act2 = bj - myByLo(myThid)
max2 = myByHi(myThid) - myByLo(myThid) + 1
act3 = myThid - 1
max3 = nTx*nTy
act4 = ikey_dynamics - 1
ticekey = (act1 + 1) + act2*max1
& + act3*max1*max2
& + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE flxsw(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
DO j = 1-OLy, sNy+OLy
DO i = 1-OLx, sNx+OLx
tIc1(i,j) = 0. _d 0
tIc2(i,j) = 0. _d 0
C- set these arrays everywhere: overlap are not set and not used,
C but some arrays are stored and storage includes overlap.
flx0exSW(i,j) = 0. _d 0
flxTexSW(i,j) = 0. _d 0
dFlxdT (i,j) = 0. _d 0
evap00 (i,j) = 0. _d 0
evapT (i,j) = 0. _d 0
dEvdT (i,j) = 0. _d 0
iceFlag (i,j) = 0. _d 0
Tsf (i,j) = 0. _d 0
ENDDO
ENDDO
#endif
stdUnit = standardMessageUnit
useBlkFlx = useEXF .OR. useBulkForce
dt = thSIce_dtTemp
IF ( dhSnowLin.GT.0. _d 0 ) THEN
recip_dhSnowLin = 1. _d 0 / dhSnowLin
ELSE
recip_dhSnowLin = 0. _d 0
ENDIF
iterate4Tsf = .FALSE.
icount = 0
DO j = jMin, jMax
DO i = iMin, iMax
#ifdef ALLOW_AUTODIFF_TAMC
ikey_1 = i
& + sNx*(j-1)
& + sNx*sNy*act1
& + sNx*sNy*max1*act2
& + sNx*sNy*max1*max2*act3
& + sNx*sNy*max1*max2*max3*act4
C--
CADJ STORE hSnow1(i,j) = comlev1_thsice_1, key=ikey_1
cCADJ STORE flxsw(i,j) = comlev1_thsice_1, key=ikey_1
cCADJ STORE qic1(i,j) = comlev1_thsice_1, key=ikey_1
cCADJ STORE qic2(i,j) = comlev1_thsice_1, key=ikey_1
#endif
IF ( icMask(i,j).GT.0. _d 0) THEN
iterate4Tsf = .TRUE.
iceFlag(i,j) = 1. _d 0
#ifdef ALLOW_DBUG_THSICE
IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,'(A,2I4,2I2)')
& 'ThSI_SOLVE4T: i,j=',i,j,bi,bj
#endif
IF ( hIce(i,j).LT.hIceMin ) THEN
C if hi < hIceMin, melt the ice.
C keep the position of this problem but do the stop later
ii = i
jj = j
icount = icount + 1
ENDIF
C-- Fractional snow cover:
C assume a linear distribution of snow thickness, with dhSnowLin slope,
C from hs-dhSnowLin to hs+dhSnowLin if full ice & snow cover.
C frsnow = fraction of snow over the ice-covered part of the grid cell
IF ( hSnow1(i,j) .GT. icMask(i,j)*dhSnowLin ) THEN
frsnow = 1. _d 0
ELSE
frsnow = hSnow1(i,j)*recip_dhSnowLin/icMask(i,j)
IF ( frsnow.GT.0. _d 0 ) frsnow = SQRT(frsnow)
ENDIF
C-- Compute SW flux absorbed at surface and penetrating to layer 1.
fswpen = flxSW(i,j) * (1. _d 0 - frsnow) * i0swFrac
fswocn = fswpen * exp(-ksolar*hIce(i,j))
fswint = fswpen - fswocn
fswdn = flxSW(i,j) - fswpen
C Initialise Atmospheric surf. heat flux with net SW flux at the surface
flxAtm(i,j) = flxSW(i,j)
C SW flux at sea-ice base left to the ocean
flxSW(i,j) = fswocn
C Initialise surface Heating with SW contribution
sHeat(i,j) = fswdn
C-- Compute conductivity terms at layer interfaces.
k12(i,j) = 4. _d 0*kIce*kSnow
& / (kSnow*hIce(i,j) + 4. _d 0*kIce*hSnow1(i,j))
k32 = 2. _d 0*kIce / hIce(i,j)
C-- Compute ice temperatures
a1 = cpIce
b1 = qIc1(i,j) + (cpWater-cpIce )*Tmlt1 - Lfresh
c1 = Lfresh * Tmlt1
tIc1(i,j) = 0.5 _d 0 *(-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/a1
tIc2(i,j) = (Lfresh-qIc2(i,j)) / cpIce
#ifdef ALLOW_DBUG_THSICE
IF (tIc1(i,j).GT.0. _d 0 ) THEN
WRITE(stdUnit,'(A,I12,1PE14.6)')
& ' BBerr: Tice(1) > 0 ; it=', myIter, qIc1(i,j)
WRITE(stdUnit,'(A,4I5,2F11.4)')
& ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,tIc1(i,j),tIc2(i,j)
ENDIF
IF ( tIc2(i,j).GT.0. _d 0) THEN
WRITE(stdUnit,'(A,I12,1PE14.6)')
& ' BBerr: Tice(2) > 0 ; it=', myIter, qIc2(i,j)
WRITE(stdUnit,'(A,4I5,2F11.4)')
& ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,tIc1(i,j),tIc2(i,j)
ENDIF
IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1010)
& 'ThSI_SOLVE4T: k, Ts, Tice=',0,tSrf1(i,j),tIc1(i,j),tIc2(i,j)
#endif
C-- Compute coefficients used in quadratic formula.
a10(i,j) = rhoi*cpIce *hIce(i,j)/(2. _d 0*dt) +
& k32*( 4. _d 0*dt*k32 + rhoi*cpIce *hIce(i,j) )
& / ( 6. _d 0*dt*k32 + rhoi*cpIce *hIce(i,j) )
b10(i,j) = -hIce(i,j)*
& ( rhoi*cpIce*tIc1(i,j) + rhoi*Lfresh*Tmlt1/tIc1(i,j) )
& /(2. _d 0*dt)
& - k32*( 4. _d 0*dt*k32*tFrz(i,j)
& +rhoi*cpIce*hIce(i,j)*tIc2(i,j) )
& / ( 6. _d 0*dt*k32 + rhoi*cpIce *hIce(i,j) )
& - fswint
c10(i,j) = rhoi*Lfresh*hIce(i,j)*Tmlt1 / (2. _d 0*dt)
ELSE
iceFlag(i,j) = 0. _d 0
ENDIF
ENDDO
ENDDO
#ifndef ALLOW_AUTODIFF
IF ( icount .gt. 0 ) THEN
WRITE(stdUnit,'(A,I5,A)')
& 'THSICE_SOLVE4TEMP: there are ',icount,
& ' case(s) where hIce WRITE(stdUnit,'(A,I3,A1,I3,A)')
& 'THSICE_SOLVE4TEMP: the last one was at (',ii,',',jj,
& ') with hIce = ', hIce(ii,jj)
WRITE( msgBuf, '(A)')
& 'THSICE_SOLVE4TEMP: should not enter if hIce CALL PRINT_ERROR( msgBuf , myThid )
STOP 'ABNORMAL END: S/R THSICE_SOLVE4TEMP'
ENDIF
#endif /* ALLOW_AUTODIFF */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE devdt(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE tsf(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE hsnow1(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE sh(:,:,:,:) = comlev1_bibj,key=ticekey,byte=isbyte
#endif
C-- Get surface fluxes over melting surface
#ifndef ALLOW_AUTODIFF
IF ( useBlkFlx .AND. iterate4Tsf ) THEN
#endif
DO j = jMin, jMax
DO i = iMin, iMax
Tsf(i,j) = 0.
ENDDO
ENDDO
#ifndef ALLOW_AUTODIFF
IF ( useEXF ) THEN
#endif
k = 1
CALL THSICE_GET_EXF( bi, bj, k,
I iMin, iMax, jMin, jMax,
I iceFlag, hSnow1, Tsf,
O flx0exSW, dFlxdT, evap00, dEvdT,
I myTime, myIter, myThid )
#ifndef ALLOW_AUTODIFF
C- could add this "ifdef" to hide THSICE_GET_BULKF from TAF
ELSEIF ( useBulkForce ) THEN
CALL THSICE_GET_BULKF( bi, bj,
I iMin, iMax, jMin, jMax,
I iceFlag, hSnow1, Tsf,
O flx0exSW, dFlxdT, evap00, dEvdT,
I myTime, myIter, myThid )
C--- end if: IF ( useEXF ) THEN
ENDIF
C--- end if: IF ( useBlkFlx .AND. iterate4Tsf ) THEN
ENDIF
#endif
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
C--- Compute new surface and internal temperatures; iterate until
C Tsfc converges.
DO j = jMin, jMax
DO i = iMin, iMax
Tsf(i,j) = tSrf1(i,j)
dTsrf1(i,j) = Terrmax
ENDDO
ENDDO
IF ( useBlkFlx ) THEN
iterMax = nitMaxTsf
ELSE
iterMax = 1
ENDIF
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE devdt(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE dflxdt(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE flx0exsw(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE flxtexsw(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE evap00(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE evapt(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
#endif
C ----- begin iteration -----
DO k = 1,iterMax
#ifndef ALLOW_AUTODIFF
IF ( iterate4Tsf ) THEN
iterate4Tsf = .FALSE.
C
IF ( useBlkFlx ) THEN
#else
kkey = (ticekey-1)*MaxTsf + k
CADJ STORE iceflag = comlev1_thsice_s4t, key=kkey, byte=isbyte
CADJ STORE hSnow1 = comlev1_thsice_s4t, key=kkey, byte=isbyte
CADJ STORE tsf = comlev1_thsice_s4t, key=kkey, byte=isbyte
CADJ STORE dEvdT = comlev1_thsice_s4t, key=kkey, byte=isbyte
#endif /* ALLOW_AUTODIFF */
C-- Compute top surface flux.
#ifndef ALLOW_AUTODIFF
IF ( useEXF ) THEN
#endif
CALL THSICE_GET_EXF( bi, bj, k+1,
I iMin, iMax, jMin, jMax,
I iceFlag, hSnow1, Tsf,
O flxTexSW, dFlxdT, evapT, dEvdT,
I myTime, myIter, myThid )
C- could add this "ifdef" to hide THSICE_GET_BULKF from TAF
#ifndef ALLOW_AUTODIFF
ELSEIF ( useBulkForce ) THEN
CALL THSICE_GET_BULKF( bi, bj,
I iMin, iMax, jMin, jMax,
I iceFlag, hSnow1, Tsf,
O flxTexSW, dFlxdT, evapT, dEvdT,
I myTime, myIter, myThid )
C--- end if: IF ( useEXF ) THEN
ENDIF
ELSE
DO j = jMin, jMax
DO i = iMin, iMax
IF ( iceFlag(i,j).GT.0. _d 0 ) THEN
flxTexSW(i,j) = flxExSW(i,j,1)
dFlxdT(i,j) = flxExSW(i,j,2)
ENDIF
ENDDO
ENDDO
C--- end if: IF ( useBlkFlx ) THEN
ENDIF
#endif /* ALLOW_AUTODIFF */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE devdt(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE dflxdt(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE flxtexsw(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE iceflag(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE tsf(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
#endif
C-- Compute new top layer and surface temperatures.
C If Tsfc is computed to be > 0 C, fix Tsfc = 0 and recompute T1
C with different coefficients.
DO j = jMin, jMax
DO i = iMin, iMax
IF ( iceFlag(i,j).GT.0. _d 0 ) THEN
flxNet = sHeat(i,j) + flxTexSW(i,j)
#ifdef ALLOW_DBUG_THSICE
IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1020)
& 'ThSI_SOLVE4T: flxNet,dFlxdT,k12,D=',
& flxNet, dFlxdT(i,j), k12(i,j), k12(i,j)-dFlxdT(i,j)
#endif
a1 = a10(i,j) - k12(i,j)*dFlxdT(i,j) / (k12(i,j)-dFlxdT(i,j))
b1 = b10(i,j) - k12(i,j)*(flxNet-dFlxdT(i,j)*Tsf(i,j))
& /(k12(i,j)-dFlxdT(i,j))
c1 = c10(i,j)
tIc1(i,j) = -(b1 + SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)
dTsrf1(i,j) = (flxNet + k12(i,j)*(tIc1(i,j)-Tsf(i,j)))
& /(k12(i,j)-dFlxdT(i,j))
Tsf(i,j) = Tsf(i,j) + dTsrf1(i,j)
C
IF ( Tsf(i,j) .GT. 0. _d 0 ) THEN
#ifdef ALLOW_DBUG_THSICE
IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1010)
& 'ThSI_SOLVE4T: k,ts,t1,dTs=',k,Tsf(i,j),tIc1(i,j),dTsrf1(i,j)
#endif
a1 = a10(i,j) + k12(i,j)
C note: b1 = b10 - k12*Tf0
b1 = b10(i,j)
tIc1(i,j) = (-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)
Tsf(i,j) = 0. _d 0
#ifndef ALLOW_AUTODIFF
IF ( useBlkFlx ) THEN
#endif /* ALLOW_AUTODIFF */
flxTexSW(i,j) = flx0exSW(i,j)
evapT(i,j) = evap00(i,j)
dTsrf1(i,j) = 0. _d 0
#ifndef ALLOW_AUTODIFF
ELSE
flxTexSW(i,j) = flxExSW(i,j,0)
dTsrf1(i,j) = 1000.
dFlxdT(i,j) = 0.
ENDIF
#endif /* ALLOW_AUTODIFF */
ENDIF
C-- Check for convergence. If no convergence, then repeat.
IF (ABS(dTsrf1(i,j)).GE.Terrmax) THEN
iceFlag(i,j) = 1. _d 0
ELSE
iceFlag(i,j) = 0. _d 0
ENDIF
iterate4Tsf = iterate4Tsf .OR. (iceFlag(i,j).GT.0. _d 0)
C Convergence test: Make sure Tsfc has converged, within prescribed error.
C (Energy conservation is guaranteed within machine roundoff, even
C if Tsfc has not converged.)
C If no convergence, then repeat.
#ifdef ALLOW_DBUG_THSICE
IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1010)
& 'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf(i,j),tIc1(i,j),dTsrf1(i,j)
IF ( useBlkFlx .AND. k.EQ.nitMaxTsf .AND.
& (iceFlag(i,j).GT.0. _d 0) ) THEN
WRITE(stdUnit,'(A,4I4,I12,F15.9)')
& ' BB: not converge: i,j,it,hi=',i,j,bi,bj,myIter,hIce(i,j)
WRITE(stdUnit,*)
& 'BB: not converge: Tsf, dTsf=', Tsf(i,j), dTsrf1(i,j)
WRITE(stdUnit,*)
& 'BB: not converge: flxNet,dFlxT=', flxNet, dFlxdT(i,j)
IF ( Tsf(i,j).LT.-70. _d 0 ) THEN
WRITE( msgBuf, '(A,2I4,2I3,I10,F12.3)')
& 'THSICE_SOLVE4TEMP: Too low Tsf in', i, j, bi, bj,
& myIter, Tsf(i,j)
CALL PRINT_ERROR( msgBuf , myThid )
STOP 'ABNORMAL END: S/R THSICE_SOLVE4TEMP'
ENDIF
ENDIF
#endif /* ALLOW_DBUG_THSICE */
ENDIF
ENDDO
ENDDO
#ifndef ALLOW_AUTODIFF
C--- end if: IF ( iterate4Tsf ) THEN
ENDIF
#endif /* ALLOW_AUTODIFF */
C--- end loop DO k = 1,iterMax
ENDDO
C ------ end iteration ------------
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
DO j = jMin, jMax
DO i = iMin, iMax
IF ( icMask(i,j).GT.0. _d 0 ) THEN
C-- Compute new bottom layer temperature.
k32 = 2. _d 0*kIce / hIce(i,j)
tIc2(i,j) = ( 2. _d 0*dt*k32*(tIc1(i,j)+2. _d 0*tFrz(i,j))
& + rhoi*cpIce*hIce(i,j)*tIc2(i,j))
& /(6. _d 0*dt*k32 + rhoi*cpIce*hIce(i,j))
#ifdef ALLOW_DBUG_THSICE
IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1010)
& 'ThSI_SOLVE4T: k, Ts, Tice=',k,Tsf(i,j),tIc1(i,j),tIc2(i,j)
netSW = flxAtm(i,j)
#endif
C-- Compute final flux values at surfaces.
tSrf1(i,j) = Tsf(i,j)
fct = k12(i,j)*(Tsf(i,j)-tIc1(i,j))
flxCnB(i,j) = 4. _d 0*kIce *(tIc2(i,j)-tFrz(i,j))/hIce(i,j)
flxNet = sHeat(i,j) + flxTexSW(i,j)
flxNet = flxNet + dFlxdT(i,j)*dTsrf1(i,j)
#ifndef ALLOW_AUTODIFF
IF ( useBlkFlx ) THEN
#endif
C- needs to update also Evap (Tsf changes) since Latent heat has been updated
evpAtm(i,j) = evapT(i,j) + dEvdT(i,j)*dTsrf1(i,j)
#ifndef ALLOW_AUTODIFF
ELSE
C- WARNING: Evap & +Evap*Lfresh are missing ! (but only affects Diagnostics)
evpAtm(i,j) = 0.
ENDIF
#endif
C- Update energy flux to Atmos with other than SW contributions;
C use latent heat = Lvap (energy=0 for liq. water at 0.oC)
flxAtm(i,j) = flxAtm(i,j) + flxTexSW(i,j)
& + dFlxdT(i,j)*dTsrf1(i,j) + evpAtm(i,j)*Lfresh
C- excess of energy @ surface (used for surface melting):
sHeat(i,j) = flxNet - fct
#ifdef ALLOW_DBUG_THSICE
IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1020)
& 'ThSI_SOLVE4T: flxNet,fct,Dif,flxCnB=',
& flxNet,fct,flxNet-fct,flxCnB(i,j)
#endif
C-- Compute new enthalpy for each layer.
qIc1(i,j) = -cpWater*Tmlt1 + cpIce *(Tmlt1-tIc1(i,j))
& + Lfresh*(1. _d 0-Tmlt1/tIc1(i,j))
qIc2(i,j) = -cpIce *tIc2(i,j) + Lfresh
#ifdef ALLOW_DBUG_THSICE
C-- Make sure internal ice temperatures do not exceed Tmlt.
C (This should not happen for reasonable values of i0swFrac)
IF (tIc1(i,j) .GE. Tmlt1) THEN
WRITE(stdUnit,'(A,2I4,2I3,1P2E14.6)')
& ' BBerr - Bug: IceT(1) > Tmlt',i,j,bi,bj,tIc1(i,j),Tmlt1
ENDIF
IF (tIc2(i,j) .GE. 0. _d 0) THEN
WRITE(stdUnit,'(A,2I4,2I3,1P2E14.6)')
& ' BBerr - Bug: IceT(2) > 0',i,j,bi,bj,tIc2(i,j)
ENDIF
IF ( dBug(i,j,bi,bj) ) THEN
WRITE(stdUnit,1020) 'ThSI_SOLV_4T: Tsf, Tice(1,2), dTsurf=',
& Tsf(i,j), tIc1(i,j), tIc2(i,j), dTsrf1(i,j)
WRITE(stdUnit,1020) 'ThSI_SOLV_4T: sHeat, flxCndBt, Qice =',
& sHeat(i,j), flxCnB(i,j), qIc1(i,j), qIc2(i,j)
WRITE(stdUnit,1020) 'ThSI_SOLV_4T: flxA, evpA, fxSW_bf,af=',
& flxAtm(i,j), evpAtm(i,j), netSW, flxSW(i,j)
ENDIF
#endif
ELSE
C-- ice-free grid point:
c tIc1 (i,j) = 0. _d 0
c tIc2 (i,j) = 0. _d 0
dTsrf1(i,j) = 0. _d 0
c sHeat (i,j) = 0. _d 0
c flxCnB(i,j) = 0. _d 0
c flxAtm(i,j) = 0. _d 0
c evpAtm(i,j) = 0. _d 0
ENDIF
ENDDO
ENDDO
#endif /* ALLOW_THSICE */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
RETURN
END