pkg/thsice/thsice_solve4temp.F

C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_solve4temp.F,v 1.7 2006/03/13 03:55:39 jmc Exp $
C $Name:  $

#include "THSICE_OPTIONS.h"

CBOP
C     !ROUTINE: THSICE_SOLVE4TEMP
C     !INTERFACE:
      SUBROUTINE THSICE_SOLVE4TEMP(
     I                  bi, bj, siLo, siHi, sjLo, sjHi,
     I                  iMin,iMax, jMin,jMax, dBugFlag, useBlkFlx,
     I                  iceMask, hIce, hSnow, tFrz, flxExSW,
     U                  flxSW, tSrf, qIc1, qIc2,
     O                  tIc1, tIc2, dTsrf,
     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."  J. Geophys.
C..       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 "THSICE_SIZE.h"
#include "THSICE_PARAMS.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine Arguments ==
C     siLo,siHi   :: size of input/output array: 1rst dim. lower,higher bounds
C     sjLo,sjHi   :: size of input/output array: 2nd  dim. lower,higher bounds
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     useBlkFlx   :: use surf. fluxes from bulk-forcing external S/R
C---  Input:
C         iceMask :: sea-ice fractional mask [0-1]
C  hIce    (hi)   :: ice height [m]
C  hSnow   (hs)   :: snow height [m]
C  tFrz    (Tf)   :: 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  (netSW) :: net Short-Wave flux (+=down) [W/m2]: input= at surface
C           (=)   ::                 output= below sea-ice, into the ocean
C  tSrf    (Tsf)  :: surface (ice or snow) temperature
C  qIc1   (qicen) :: ice enthalpy (J/kg), 1rst level
C  qIc2   (qicen) :: ice enthalpy (J/kg), 2nd level
C---  Output
C  tIc1    (Tice) :: temperature of ice layer 1 [oC]
C  tIc2    (Tice) :: temperature of ice layer 2 [oC]
C  dTsrf   (dTsf) :: surf. temp adjusment: Ts^n+1 - Ts^n
C  sHeat(sHeating):: 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 siLo, siHi, sjLo, sjHi
      INTEGER bi,bj
      INTEGER iMin, iMax
      INTEGER jMin, jMax
      LOGICAL dBugFlag
      LOGICAL useBlkFlx
      _RL iceMask(siLo:siHi,sjLo:sjHi)
      _RL hIce   (siLo:siHi,sjLo:sjHi)
      _RL hSnow  (siLo:siHi,sjLo:sjHi)
      _RL tFrz   (siLo:siHi,sjLo:sjHi)
      _RL flxExSW(iMin:iMax,jMin:jMax,0:2)
      _RL flxSW  (siLo:siHi,sjLo:sjHi)
      _RL tSrf   (siLo:siHi,sjLo:sjHi)
      _RL qIc1   (siLo:siHi,sjLo:sjHi)
      _RL qIc2   (siLo:siHi,sjLo:sjHi)
      _RL tIc1   (siLo:siHi,sjLo:sjHi)
      _RL tIc2   (siLo:siHi,sjLo:sjHi)
c     _RL dTsrf  (siLo:siHi,sjLo:sjHi)
      _RL dTsrf  (iMin:iMax,jMin:jMax)
      _RL sHeat  (siLo:siHi,sjLo:sjHi)
      _RL flxCnB (siLo:siHi,sjLo:sjHi)
      _RL flxAtm (siLo:siHi,sjLo:sjHi)
      _RL evpAtm (siLo:siHi,sjLo:sjHi)
      _RL  myTime
      INTEGER myIter
      INTEGER myThid
CEOP

#ifdef ALLOW_THSICE
C     !LOCAL VARIABLES:
C---  local copy of input/output argument list variables (see description above)
c     _RL flxExcSw(0:2)
      _RL Tf
      _RL hi
      _RL hs
      _RL netSW
      _RL Tsf
      _RL qicen(nlyr)
      _RL Tice (nlyr)
c     _RL sHeating
c     _RL flxCnB
      _RL dTsf
c     _RL flxAtm
c     _RL evpAtm

C     == Local Variables ==
      INTEGER i,j
      INTEGER k, iterMax

      _RL  frsnow        ! fractional snow cover
      _RL  fswpen        ! SW penetrating beneath surface (W m-2)
      _RL  fswdn         ! SW absorbed at surface (W m-2)
      _RL  fswint        ! SW absorbed in ice (W m-2)
      _RL  fswocn        ! SW passed through ice to ocean (W m-2)
      _RL  flxExceptSw   ! net surface heat flux, except short-wave (W/m2)
C     evap           :: evaporation over snow/ice [kg/m2/s] (>0 if evaporate)
C     dEvdT          :: derivative of evap. with respect to Tsf [kg/m2/s/K]
      _RL  evap, dEvdT
      _RL  flx0          ! net surf heat flux, from Atmos. to sea-ice (W m-2)
      _RL  fct           ! heat conducted to top surface
      _RL  df0dT         ! deriv of flx0 wrt Tsf (W m-2 deg-1)
      _RL  k12, k32      ! thermal conductivity terms
      _RL  a10, b10      ! coefficients in quadratic eqn for T1
      _RL  a1, b1, c1    ! coefficients in quadratic eqn for T1
c     _RL  Tsf_start     ! old value of Tsf
      _RL  dt            ! timestep
      INTEGER iceornot

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-|--+----|

      dt  = thSIce_deltaT
      DO j = jMin, jMax
       DO i = iMin, iMax
        IF ( iceMask(i,j).GT.0. _d 0) THEN
          hi      = hIce(i,j)
          hs      = hSnow(i,j)
          Tf      = tFrz(i,j)
          netSW   = flxSW(i,j)
          Tsf     = tSrf(i,j)
          qicen(1)= qIc1(i,j)
          qicen(2)= qIc2(i,j)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_DBUG_THSICE
          IF ( dBug(i,j,bi,bj) ) WRITE(6,'(A,2I4,2I2)')
     &         'ThSI_SOLVE4T: i,j=',i,j,bi,bj
#endif
      IF ( hi.LT.himin ) THEN
C If hi < himin, melt the ice.
         STOP 'THSICE_SOLVE4TEMP: should not enter if hi<himin'
      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C fractional snow cover
      frsnow = 0. _d 0
      IF (hs .GT. 0. _d 0) frsnow = 1. _d 0

C Compute SW flux absorbed at surface and penetrating to layer 1.
      fswpen  = netSW * (1. _d 0 - frsnow) * i0
      fswocn = fswpen * exp(-ksolar*hi)
      fswint = fswpen - fswocn

      fswdn = netSW - fswpen

C Compute conductivity terms at layer interfaces.

      k12 = 4. _d 0*kice*ksnow / (ksnow*hi + 4. _d 0*kice*hs)
      k32 = 2. _d 0*kice  / hi

C compute ice temperatures
      a1 = cpice
      b1 = qicen(1) + (cpwater-cpice )*Tmlt1 - Lfresh
      c1 = Lfresh * Tmlt1
      Tice(1) = 0.5 _d 0 *(-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/a1
      Tice(2) = (Lfresh-qicen(2)) / cpice

      IF (Tice(1).GT.0. _d 0 .OR. Tice(2).GT.0. _d 0) THEN
          WRITE (6,*) 'BBerr Tice(1) > 0 = ',Tice(1)
          WRITE (6,*) 'BBerr Tice(2) > 0 = ',Tice(2)
      ENDIF
#ifdef ALLOW_DBUG_THSICE
      IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
     &  'ThSI_SOLVE4T: k, Ts, Tice=',0,Tsf,Tice
#endif

C Compute coefficients used in quadratic formula.

      a10 = rhoi*cpice *hi/(2. _d 0*dt) +
     &      k32 * (4. _d 0*dt*k32 + rhoi*cpice *hi)
     &       / (6. _d 0*dt*k32 + rhoi*cpice *hi)
      b10 = -hi*
     &      (rhoi*cpice*Tice(1)+rhoi*Lfresh*Tmlt1/Tice(1))
     &       /(2. _d 0*dt)
     &      - k32 * (4. _d 0*dt*k32*Tf+rhoi*cpice *hi*Tice(2))
     &       / (6. _d 0*dt*k32 + rhoi*cpice *hi) - fswint
      c1 = rhoi*Lfresh*hi*Tmlt1 / (2. _d 0*dt)

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C Compute new surface and internal temperatures; iterate until
C Tsfc converges.

      IF ( useBlkFlx ) THEN
        iterMax = nitMaxTsf
      ELSE
        iterMax = 1
      ENDIF
      dTsf = Terrmax

C ----- begin iteration  -----
      DO k = 1,iterMax
       IF ( ABS(dTsf).GE.Terrmax ) THEN

C Save temperatures at start of iteration.
c        Tsf_start = Tsf

        IF ( useBlkFlx ) THEN
C Compute top surface flux.
         IF (hs.GT.3. _d -1) THEN
              iceornot=2
         ELSE
              iceornot=1
         ENDIF
         CALL THSICE_GET_BULKF(
     I                        iceornot, Tsf,
     O                        flxExceptSw, df0dT, evap, dEvdT,
     I                        i,j,bi,bj,myThid )
        ELSE
         flxExceptSw = flxExSW(i,j,1)
         df0dT = flxExSW(i,j,2)
        ENDIF
        flx0 = fswdn + flxExceptSw
#ifdef ALLOW_DBUG_THSICE
      IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)
     &  'ThSI_SOLVE4T: flx0,df0dT,k12,D=', flx0,df0dT,k12,k12-df0dT
#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.

         a1 = a10 - k12*df0dT / (k12-df0dT)
         b1 = b10 - k12*(flx0-df0dT*Tsf) / (k12-df0dT)
         Tice(1) = -(b1 + SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)
         dTsf = (flx0 + k12*(Tice(1)-Tsf)) / (k12-df0dT)
         Tsf = Tsf + dTsf
         IF (Tsf .GT. 0. _d 0) THEN
#ifdef ALLOW_DBUG_THSICE
            IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
     &        'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf,Tice(1),dTsf
#endif
            a1 = a10 + k12
            b1 = b10          ! note b1 = b10 - k12*Tf0
            Tice(1) = (-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)
            Tsf = 0. _d 0
           IF ( useBlkFlx ) THEN
            IF (hs.GT.3. _d -1) THEN
                 iceornot=2
            ELSE
                 iceornot=1
            ENDIF
            CALL THSICE_GET_BULKF(
     I                        iceornot, Tsf,
     O                        flxExceptSw, df0dT, evap, dEvdT,
     I                        i,j,bi,bj,myThid )
            dTsf = 0. _d 0
           ELSE
            flxExceptSw = flxExSW(i,j,0)
            dTsf = 1000.
            df0dT = 0.
           ENDIF
           flx0 = fswdn + flxExceptSw
         ENDIF

C Check for convergence.  If no convergence, then repeat.
C
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(6,1010)
     &     'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf,Tice(1),dTsf
#endif
         IF ( useBlkFlx .AND. k.EQ.nitMaxTsf
     &                  .AND. ABS(dTsf).GE.Terrmax ) THEN
            WRITE (6,*) 'BB: thermw conv err ',i,j,bi,bj,dTsf
            WRITE (6,*) 'BB: thermw conv err, iceheight ', hi
            WRITE (6,*) 'BB: thermw conv err: Tsf, flx0', Tsf,flx0
            IF (Tsf.LT.-70. _d 0) STOP
         ENDIF

100   continue  ! surface temperature iteration
       ENDIF
      ENDDO
150   continue
C ------ end iteration ------------

C Compute new bottom layer temperature.

      Tice(2) = (2. _d 0*dt*k32*(Tice(1)+2. _d 0*Tf)
     &        + rhoi*cpice *hi*Tice(2))
     &         /(6. _d 0*dt*k32 + rhoi*cpice *hi)
#ifdef ALLOW_DBUG_THSICE
      IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
     &  'ThSI_SOLVE4T: k, Ts, Tice=',k,Tsf,Tice
#endif

C Compute final flux values at surfaces.

      fct    = k12*(Tsf-Tice(1))
      flxCnB(i,j) = 4. _d 0*kice *(Tice(2)-Tf)/hi
      flx0   = flx0 + df0dT*dTsf
      IF ( useBlkFlx ) THEN
C--   needs to update also Evap (Tsf changes) since Latent heat has been updated
        evpAtm(i,j) = evap + dEvdT*dTsf
      ELSE
C- WARNING: Evap & +Evap*Lfresh are missing ! (but only affects Diagnostics)
        evpAtm(i,j) = 0.
      ENDIF
C-    energy flux to Atmos: use net short-wave flux at surf. and
C     use latent heat = Lvap (energy=0 for liq. water at 0.oC)
      flxAtm(i,j) = netSW + flxExceptSw
     &                    + df0dT*dTsf + evpAtm(i,j)*Lfresh
C-    excess of energy @ surface (used for surface melting):
      sHeat(i,j) = flx0 - fct

C-    SW flux at sea-ice base left to the ocean
      flxSW(i,j) = fswocn

#ifdef ALLOW_DBUG_THSICE
      IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)
     &  'ThSI_SOLVE4T: flx0,fct,Dif,flxCnB=',
     &                 flx0,fct,flx0-fct,flxCnB(i,j)
#endif

C Compute new enthalpy for each layer.

      qicen(1) = -cpwater*Tmlt1 + cpice *(Tmlt1-Tice(1))
     &            + Lfresh*(1. _d 0-Tmlt1/Tice(1))
      qicen(2) = -cpice *Tice(2) + Lfresh

C Make sure internal ice temperatures do not exceed Tmlt.
C (This should not happen for reasonable values of i0.)

      IF (Tice(1) .GE. Tmlt1) THEN
        WRITE (6,'(A,2I4,2I3,1P2E14.6)')
     &   'BBerr - Bug: IceT(1) > Tmlt',i,j,bi,bj,Tice(1),Tmlt1
      ENDIF
      IF (Tice(2) .GE. 0. _d 0) THEN
       WRITE (6,'(A,2I4,2I3,1P2E14.6)')
     &   'BBerr - Bug: IceT(2) > 0',i,j,bi,bj,Tice(2)
      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C--    Update Sea-Ice state :
          tSrf(i,j)  = Tsf
          tIc1(i,j)  = Tice(1)
          tic2(i,j)  = Tice(2)
          qIc1(i,j)  = qicen(1)
          qIc2(i,j)  = qicen(2)
c         dTsrf(i,j) = dTsf
          IF ( .NOT.useBlkFlx ) dTsrf(i,j) = dTsf
c         sHeat(i,j) = sHeating
c         flxCnB(i,j)= flxCnB
c         flxAtm(i,j)= flxAtm
c         evpAtm(i,j)= evpAtm
#ifdef ALLOW_DBUG_THSICE
          IF ( dBug(i,j,bi,bj) ) THEN
           WRITE(6,1020) 'ThSI_SOLV_4T: Tsf, Tice(1,2), dTsurf=',
     &                              Tsf, Tice, dTsf
           WRITE(6,1020) 'ThSI_SOLV_4T: sHeat, flxCndBt, Qice =',
     &           sHeat(i,j), flxCnB(i,j), qicen
           WRITE(6,1020) 'ThSI_SOLV_4T: flxA, evpA, fxSW_bf,af=',
     &           flxAtm(i,j), evpAtm(i,j), netSW, flxSW(i,j)
          ENDIF
#endif
        ELSE
          IF ( .NOT.useBlkFlx ) dTsrf(i,j) = 0. _d 0
        ENDIF
       ENDDO
      ENDDO
#endif  /* ALLOW_THSICE */

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      RETURN
      END
1	jmc	1.8	C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_solve4temp.F,v 1.7 2006/03/13 03:55:39 jmc Exp $
2	jmc	1.1	C $Name: $
3
4			#include "THSICE_OPTIONS.h"
5
6			CBOP
7			C !ROUTINE: THSICE_SOLVE4TEMP
8			C !INTERFACE:
9			SUBROUTINE THSICE_SOLVE4TEMP(
10	jmc	1.8	I bi, bj, siLo, siHi, sjLo, sjHi,
11			I iMin,iMax, jMin,jMax, dBugFlag, useBlkFlx,
12			I iceMask, hIce, hSnow, tFrz, flxExSW,
13			U flxSW, tSrf, qIc1, qIc2,
14			O tIc1, tIc2, dTsrf,
15			O sHeat, flxCnB, flxAtm, evpAtm,
16			I myTime, myIter, myThid )
17	jmc	1.1	C !DESCRIPTION: \bv
18			C ==========================================================
19			C \| S/R THSICE_SOLVE4TEMP
20			C ==========================================================
21			C \| Solve (implicitly) for sea-ice and surface temperature
22			C ==========================================================
23			C \ev
24
25	jmc	1.8	C ADAPTED FROM:
26			C LANL CICE.v2.0.2
27			C-----------------------------------------------------------------------
28			C.. thermodynamics (vertical physics) based on M. Winton 3-layer model
29			C.. See Bitz, C. M. and W. H. Lipscomb, 1999: "An energy-conserving
30			C.. thermodynamic sea ice model for climate study." J. Geophys.
31			C.. Res., 104, 15669 - 15677.
32			C.. Winton, M., 1999: "A reformulated three-layer sea ice model."
33			C.. Submitted to J. Atmos. Ocean. Technol.
34			C.. authors Elizabeth C. Hunke and William Lipscomb
35			C.. Fluid Dynamics Group, Los Alamos National Laboratory
36			C-----------------------------------------------------------------------
37			Cc****subroutine thermo_winton(n,fice,fsnow,dqice,dTsfc)
38			C.. Compute temperature change using Winton model with 2 ice layers, of
39			C.. which only the top layer has a variable heat capacity.
40
41	jmc	1.1	C !USES:
42			IMPLICIT NONE
43
44			C == Global variables ===
45	jmc	1.5	#include "EEPARAMS.h"
46	jmc	1.1	#include "THSICE_SIZE.h"
47			#include "THSICE_PARAMS.h"
48
49			C !INPUT/OUTPUT PARAMETERS:
50			C == Routine Arguments ==
51	jmc	1.8	C siLo,siHi :: size of input/output array: 1rst dim. lower,higher bounds
52			C sjLo,sjHi :: size of input/output array: 2nd dim. lower,higher bounds
53			C bi,bj :: tile indices
54			C iMin,iMax :: computation domain: 1rst index range
55			C jMin,jMax :: computation domain: 2nd index range
56			C dBugFlag :: allow to print debugging stuff (e.g. on 1 grid point).
57			C useBlkFlx :: use surf. fluxes from bulk-forcing external S/R
58			C--- Input:
59			C iceMask :: sea-ice fractional mask [0-1]
60			C hIce (hi) :: ice height [m]
61			C hSnow (hs) :: snow height [m]
62			C tFrz (Tf) :: sea-water freezing temperature [oC] (function of S)
63			C flxExSW (=) :: surf. heat flux (+=down) except SW, function of surf. temp Ts:
64			C 0: Flx(Ts=0) ; 1: Flx(Ts=Ts^n) ; 2: d.Flx/dTs(Ts=Ts^n)
65			C--- Modified (input&output):
66			C flxSW (netSW) :: net Short-Wave flux (+=down) [W/m2]: input= at surface
67			C (=) :: output= below sea-ice, into the ocean
68			C tSrf (Tsf) :: surface (ice or snow) temperature
69			C qIc1 (qicen) :: ice enthalpy (J/kg), 1rst level
70			C qIc2 (qicen) :: ice enthalpy (J/kg), 2nd level
71			C--- Output
72			C tIc1 (Tice) :: temperature of ice layer 1 [oC]
73			C tIc2 (Tice) :: temperature of ice layer 2 [oC]
74			C dTsrf (dTsf) :: surf. temp adjusment: Ts^n+1 - Ts^n
75			C sHeat(sHeating):: surf heating flux left to melt snow or ice (= Atmos-conduction)
76			C flxCnB (=) :: heat flux conducted through the ice to bottom surface
77			C flxAtm (=) :: net flux of energy from the atmosphere [W/m2] (+=down)
78			C without snow precip. (energy=0 for liquid water at 0.oC)
79			C evpAtm (=) :: evaporation to the atmosphere [kg/m2/s] (>0 if evaporate)
80			C--- Input:
81			C myTime :: current Time of simulation [s]
82			C myIter :: current Iteration number in simulation
83			C myThid :: my Thread Id number
84			INTEGER siLo, siHi, sjLo, sjHi
85			INTEGER bi,bj
86			INTEGER iMin, iMax
87			INTEGER jMin, jMax
88			LOGICAL dBugFlag
89	jmc	1.1	LOGICAL useBlkFlx
90	jmc	1.8	_RL iceMask(siLo:siHi,sjLo:sjHi)
91			_RL hIce (siLo:siHi,sjLo:sjHi)
92			_RL hSnow (siLo:siHi,sjLo:sjHi)
93			_RL tFrz (siLo:siHi,sjLo:sjHi)
94			_RL flxExSW(iMin:iMax,jMin:jMax,0:2)
95			_RL flxSW (siLo:siHi,sjLo:sjHi)
96			_RL tSrf (siLo:siHi,sjLo:sjHi)
97			_RL qIc1 (siLo:siHi,sjLo:sjHi)
98			_RL qIc2 (siLo:siHi,sjLo:sjHi)
99			_RL tIc1 (siLo:siHi,sjLo:sjHi)
100			_RL tIc2 (siLo:siHi,sjLo:sjHi)
101			c _RL dTsrf (siLo:siHi,sjLo:sjHi)
102			_RL dTsrf (iMin:iMax,jMin:jMax)
103			_RL sHeat (siLo:siHi,sjLo:sjHi)
104			_RL flxCnB (siLo:siHi,sjLo:sjHi)
105			_RL flxAtm (siLo:siHi,sjLo:sjHi)
106			_RL evpAtm (siLo:siHi,sjLo:sjHi)
107			_RL myTime
108			INTEGER myIter
109			INTEGER myThid
110			CEOP
111
112			#ifdef ALLOW_THSICE
113			C !LOCAL VARIABLES:
114			C--- local copy of input/output argument list variables (see description above)
115			c _RL flxExcSw(0:2)
116	jmc	1.1	_RL Tf
117			_RL hi
118			_RL hs
119	jmc	1.8	_RL netSW
120	jmc	1.1	_RL Tsf
121			_RL qicen(nlyr)
122			_RL Tice (nlyr)
123	jmc	1.8	c _RL sHeating
124			c _RL flxCnB
125	jmc	1.1	_RL dTsf
126	jmc	1.8	c _RL flxAtm
127			c _RL evpAtm
128	jmc	1.1
129			C == Local Variables ==
130	jmc	1.8	INTEGER i,j
131	jmc	1.6	INTEGER k, iterMax
132	jmc	1.1
133			_RL frsnow ! fractional snow cover
134			_RL fswpen ! SW penetrating beneath surface (W m-2)
135			_RL fswdn ! SW absorbed at surface (W m-2)
136			_RL fswint ! SW absorbed in ice (W m-2)
137			_RL fswocn ! SW passed through ice to ocean (W m-2)
138			_RL flxExceptSw ! net surface heat flux, except short-wave (W/m2)
139			C evap :: evaporation over snow/ice [kg/m2/s] (>0 if evaporate)
140			C dEvdT :: derivative of evap. with respect to Tsf [kg/m2/s/K]
141			_RL evap, dEvdT
142			_RL flx0 ! net surf heat flux, from Atmos. to sea-ice (W m-2)
143			_RL fct ! heat conducted to top surface
144			_RL df0dT ! deriv of flx0 wrt Tsf (W m-2 deg-1)
145			_RL k12, k32 ! thermal conductivity terms
146			_RL a10, b10 ! coefficients in quadratic eqn for T1
147			_RL a1, b1, c1 ! coefficients in quadratic eqn for T1
148			c _RL Tsf_start ! old value of Tsf
149			_RL dt ! timestep
150	jmc	1.8	INTEGER iceornot
151	jmc	1.1
152	jmc	1.8	C- define grid-point location where to print debugging values
153			#include "THSICE_DEBUG.h"
154	jmc	1.1
155	jmc	1.7	1010 FORMAT(A,I3,3F11.6)
156			1020 FORMAT(A,1P4E14.6)
157	jmc	1.1
158	jmc	1.8	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
159
160	jmc	1.1	dt = thSIce_deltaT
161	jmc	1.8	DO j = jMin, jMax
162			DO i = iMin, iMax
163			IF ( iceMask(i,j).GT.0. _d 0) THEN
164			hi = hIce(i,j)
165			hs = hSnow(i,j)
166			Tf = tFrz(i,j)
167			netSW = flxSW(i,j)
168			Tsf = tSrf(i,j)
169			qicen(1)= qIc1(i,j)
170			qicen(2)= qIc2(i,j)
171			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
172			#ifdef ALLOW_DBUG_THSICE
173			IF ( dBug(i,j,bi,bj) ) WRITE(6,'(A,2I4,2I2)')
174			& 'ThSI_SOLVE4T: i,j=',i,j,bi,bj
175			#endif
176	jmc	1.6	IF ( hi.LT.himin ) THEN
177	jmc	1.1	C If hi < himin, melt the ice.
178	jmc	1.3	STOP 'THSICE_SOLVE4TEMP: should not enter if hi<himin'
179	jmc	1.6	ENDIF
180	jmc	1.1
181			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
182
183			C fractional snow cover
184			frsnow = 0. _d 0
185	jmc	1.6	IF (hs .GT. 0. _d 0) frsnow = 1. _d 0
186	jmc	1.1
187			C Compute SW flux absorbed at surface and penetrating to layer 1.
188	jmc	1.8	fswpen = netSW * (1. _d 0 - frsnow) * i0
189	jmc	1.1	fswocn = fswpen * exp(-ksolar*hi)
190			fswint = fswpen - fswocn
191
192	jmc	1.8	fswdn = netSW - fswpen
193	jmc	1.1
194			C Compute conductivity terms at layer interfaces.
195
196			k12 = 4. _d 0kiceksnow / (ksnowhi + 4. _d 0kice*hs)
197			k32 = 2. _d 0*kice / hi
198
199			C compute ice temperatures
200			a1 = cpice
201			b1 = qicen(1) + (cpwater-cpice )*Tmlt1 - Lfresh
202			c1 = Lfresh * Tmlt1
203	jmc	1.6	Tice(1) = 0.5 _d 0 (-b1 - SQRT(b1b1-4. _d 0a1c1))/a1
204	jmc	1.1	Tice(2) = (Lfresh-qicen(2)) / cpice
205
206	jmc	1.6	IF (Tice(1).GT.0. _d 0 .OR. Tice(2).GT.0. _d 0) THEN
207			WRITE (6,*) 'BBerr Tice(1) > 0 = ',Tice(1)
208			WRITE (6,*) 'BBerr Tice(2) > 0 = ',Tice(2)
209			ENDIF
210	jmc	1.8	#ifdef ALLOW_DBUG_THSICE
211			IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
212			& 'ThSI_SOLVE4T: k, Ts, Tice=',0,Tsf,Tice
213			#endif
214	jmc	1.1
215			C Compute coefficients used in quadratic formula.
216
217			a10 = rhoicpice hi/(2. _d 0*dt) +
218			& k32 * (4. _d 0dtk32 + rhoicpice hi)
219			& / (6. _d 0dtk32 + rhoicpice hi)
220			b10 = -hi*
221			& (rhoicpiceTice(1)+rhoiLfreshTmlt1/Tice(1))
222			& /(2. _d 0*dt)
223			& - k32 * (4. _d 0dtk32Tf+rhoicpice hiTice(2))
224			& / (6. _d 0dtk32 + rhoicpice hi) - fswint
225			c1 = rhoiLfreshhiTmlt1 / (2. _d 0dt)
226
227	jmc	1.4	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
228	jmc	1.1	C Compute new surface and internal temperatures; iterate until
229			C Tsfc converges.
230
231	jmc	1.6	IF ( useBlkFlx ) THEN
232			iterMax = nitMaxTsf
233			ELSE
234			iterMax = 1
235			ENDIF
236			dTsf = Terrmax
237
238	jmc	1.1	C ----- begin iteration -----
239	jmc	1.6	DO k = 1,iterMax
240			IF ( ABS(dTsf).GE.Terrmax ) THEN
241	jmc	1.1
242			C Save temperatures at start of iteration.
243			c Tsf_start = Tsf
244
245			IF ( useBlkFlx ) THEN
246			C Compute top surface flux.
247	jmc	1.6	IF (hs.GT.3. _d -1) THEN
248	jmc	1.1	iceornot=2
249	jmc	1.6	ELSE
250	jmc	1.1	iceornot=1
251	jmc	1.6	ENDIF
252	jmc	1.1	CALL THSICE_GET_BULKF(
253			I iceornot, Tsf,
254			O flxExceptSw, df0dT, evap, dEvdT,
255			I i,j,bi,bj,myThid )
256			ELSE
257	jmc	1.8	flxExceptSw = flxExSW(i,j,1)
258			df0dT = flxExSW(i,j,2)
259	jmc	1.1	ENDIF
260	jmc	1.7	flx0 = fswdn + flxExceptSw
261	jmc	1.8	#ifdef ALLOW_DBUG_THSICE
262			IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)
263			& 'ThSI_SOLVE4T: flx0,df0dT,k12,D=', flx0,df0dT,k12,k12-df0dT
264			#endif
265	jmc	1.1
266			C Compute new top layer and surface temperatures.
267			C If Tsfc is computed to be > 0 C, fix Tsfc = 0 and recompute T1
268	jmc	1.7	C with different coefficients.
269	jmc	1.1
270			a1 = a10 - k12*df0dT / (k12-df0dT)
271			b1 = b10 - k12(flx0-df0dTTsf) / (k12-df0dT)
272	jmc	1.6	Tice(1) = -(b1 + SQRT(b1b1-4. _d 0a1c1))/(2. _d 0a1)
273	jmc	1.1	dTsf = (flx0 + k12*(Tice(1)-Tsf)) / (k12-df0dT)
274			Tsf = Tsf + dTsf
275	jmc	1.6	IF (Tsf .GT. 0. _d 0) THEN
276	jmc	1.8	#ifdef ALLOW_DBUG_THSICE
277			IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
278			& 'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf,Tice(1),dTsf
279			#endif
280	jmc	1.1	a1 = a10 + k12
281			b1 = b10 ! note b1 = b10 - k12*Tf0
282	jmc	1.6	Tice(1) = (-b1 - SQRT(b1b1-4. _d 0a1c1))/(2. _d 0a1)
283	jmc	1.1	Tsf = 0. _d 0
284			IF ( useBlkFlx ) THEN
285	jmc	1.6	IF (hs.GT.3. _d -1) THEN
286	jmc	1.1	iceornot=2
287	jmc	1.6	ELSE
288	jmc	1.1	iceornot=1
289	jmc	1.6	ENDIF
290	jmc	1.1	CALL THSICE_GET_BULKF(
291			I iceornot, Tsf,
292			O flxExceptSw, df0dT, evap, dEvdT,
293			I i,j,bi,bj,myThid )
294			dTsf = 0. _d 0
295			ELSE
296	jmc	1.8	flxExceptSw = flxExSW(i,j,0)
297	jmc	1.1	dTsf = 1000.
298			df0dT = 0.
299			ENDIF
300	jmc	1.7	flx0 = fswdn + flxExceptSw
301	jmc	1.6	ENDIF
302	jmc	1.1
303			C Check for convergence. If no convergence, then repeat.
304			C
305	jmc	1.7	C Convergence test: Make sure Tsfc has converged, within prescribed error.
306	jmc	1.1	C (Energy conservation is guaranteed within machine roundoff, even
307			C if Tsfc has not converged.)
308			C If no convergence, then repeat.
309
310	jmc	1.8	#ifdef ALLOW_DBUG_THSICE
311			IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
312			& 'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf,Tice(1),dTsf
313			#endif
314	jmc	1.6	IF ( useBlkFlx .AND. k.EQ.nitMaxTsf
315			& .AND. ABS(dTsf).GE.Terrmax ) THEN
316			WRITE (6,*) 'BB: thermw conv err ',i,j,bi,bj,dTsf
317			WRITE (6,*) 'BB: thermw conv err, iceheight ', hi
318			WRITE (6,*) 'BB: thermw conv err: Tsf, flx0', Tsf,flx0
319			IF (Tsf.LT.-70. _d 0) STOP
320	jmc	1.1	ENDIF
321
322			100 continue ! surface temperature iteration
323	jmc	1.6	ENDIF
324			ENDDO
325	jmc	1.1	150 continue
326			C ------ end iteration ------------
327
328			C Compute new bottom layer temperature.
329
330			Tice(2) = (2. _d 0dtk32(Tice(1)+2. _d 0Tf)
331			& + rhoicpice hi*Tice(2))
332			& /(6. _d 0dtk32 + rhoicpice hi)
333	jmc	1.8	#ifdef ALLOW_DBUG_THSICE
334			IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
335			& 'ThSI_SOLVE4T: k, Ts, Tice=',k,Tsf,Tice
336			#endif
337	jmc	1.1
338			C Compute final flux values at surfaces.
339
340			fct = k12*(Tsf-Tice(1))
341	jmc	1.8	flxCnB(i,j) = 4. _d 0kice (Tice(2)-Tf)/hi
342	jmc	1.1	flx0 = flx0 + df0dT*dTsf
343			IF ( useBlkFlx ) THEN
344			C-- needs to update also Evap (Tsf changes) since Latent heat has been updated
345	jmc	1.8	evpAtm(i,j) = evap + dEvdT*dTsf
346	jmc	1.1	ELSE
347	jmc	1.7	C- WARNING: Evap & +Evap*Lfresh are missing ! (but only affects Diagnostics)
348	jmc	1.8	evpAtm(i,j) = 0.
349	jmc	1.1	ENDIF
350	jmc	1.7	C- energy flux to Atmos: use net short-wave flux at surf. and
351			C use latent heat = Lvap (energy=0 for liq. water at 0.oC)
352	jmc	1.8	flxAtm(i,j) = netSW + flxExceptSw
353			& + df0dTdTsf + evpAtm(i,j)Lfresh
354	jmc	1.7	C- excess of energy @ surface (used for surface melting):
355	jmc	1.8	sHeat(i,j) = flx0 - fct
356	jmc	1.1
357			C- SW flux at sea-ice base left to the ocean
358	jmc	1.8	flxSW(i,j) = fswocn
359	jmc	1.1
360	jmc	1.8	#ifdef ALLOW_DBUG_THSICE
361			IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)
362			& 'ThSI_SOLVE4T: flx0,fct,Dif,flxCnB=',
363			& flx0,fct,flx0-fct,flxCnB(i,j)
364			#endif
365	jmc	1.1
366			C Compute new enthalpy for each layer.
367
368	jmc	1.7	qicen(1) = -cpwaterTmlt1 + cpice (Tmlt1-Tice(1))
369			& + Lfresh*(1. _d 0-Tmlt1/Tice(1))
370	jmc	1.1	qicen(2) = -cpice *Tice(2) + Lfresh
371
372			C Make sure internal ice temperatures do not exceed Tmlt.
373			C (This should not happen for reasonable values of i0.)
374
375	jmc	1.7	IF (Tice(1) .GE. Tmlt1) THEN
376	jmc	1.6	WRITE (6,'(A,2I4,2I3,1P2E14.6)')
377	jmc	1.1	& 'BBerr - Bug: IceT(1) > Tmlt',i,j,bi,bj,Tice(1),Tmlt1
378	jmc	1.6	ENDIF
379			IF (Tice(2) .GE. 0. _d 0) THEN
380			WRITE (6,'(A,2I4,2I3,1P2E14.6)')
381	jmc	1.1	& 'BBerr - Bug: IceT(2) > 0',i,j,bi,bj,Tice(2)
382	jmc	1.6	ENDIF
383	jmc	1.1
384	jmc	1.8	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
385			C-- Update Sea-Ice state :
386			tSrf(i,j) = Tsf
387			tIc1(i,j) = Tice(1)
388			tic2(i,j) = Tice(2)
389			qIc1(i,j) = qicen(1)
390			qIc2(i,j) = qicen(2)
391			c dTsrf(i,j) = dTsf
392			IF ( .NOT.useBlkFlx ) dTsrf(i,j) = dTsf
393			c sHeat(i,j) = sHeating
394			c flxCnB(i,j)= flxCnB
395			c flxAtm(i,j)= flxAtm
396			c evpAtm(i,j)= evpAtm
397			#ifdef ALLOW_DBUG_THSICE
398			IF ( dBug(i,j,bi,bj) ) THEN
399			WRITE(6,1020) 'ThSI_SOLV_4T: Tsf, Tice(1,2), dTsurf=',
400			& Tsf, Tice, dTsf
401			WRITE(6,1020) 'ThSI_SOLV_4T: sHeat, flxCndBt, Qice =',
402			& sHeat(i,j), flxCnB(i,j), qicen
403			WRITE(6,1020) 'ThSI_SOLV_4T: flxA, evpA, fxSW_bf,af=',
404			& flxAtm(i,j), evpAtm(i,j), netSW, flxSW(i,j)
405			ENDIF
406			#endif
407			ELSE
408			IF ( .NOT.useBlkFlx ) dTsrf(i,j) = 0. _d 0
409			ENDIF
410			ENDDO
411			ENDDO
412	jmc	1.1	#endif /* ALLOW_THSICE */
413
414			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
415
416			RETURN
417			END