LCOV - code coverage report
Current view: top level - columnphysics - icepack_therm_bl99.F90 (source / functions) Hit Total Coverage
Test: 200704-015006:b1e41d9f12:3:base,travis,quick Lines: 339 543 62.43 %
Date: 2020-07-03 20:11:40 Functions: 5 6 83.33 %

          Line data    Source code
       1             : !=========================================================================
       2             : !
       3             : ! Update ice and snow internal temperatures
       4             : ! using Bitz and Lipscomb 1999 thermodynamics
       5             : !
       6             : ! authors: William H. Lipscomb, LANL
       7             : !          C. M. Bitz, UW
       8             : !          Elizabeth C. Hunke, LANL
       9             : !
      10             : ! 2012: Split from icepack_therm_vertical.F90
      11             : 
      12             :       module icepack_therm_bl99
      13             : 
      14             :       use icepack_kinds
      15             :       use icepack_parameters, only: c0, c1, c2, p1, p5, puny
      16             : #ifdef CESMCOUPLED
      17             :       use icepack_parameters, only: p01
      18             : #endif
      19             :       use icepack_parameters, only: rhoi, rhos, hs_min, cp_ice, cp_ocn, depressT, Lfresh, ksno, kice
      20             :       use icepack_parameters, only: conduct, calc_Tsfc, solve_zsal
      21             :       use icepack_warnings, only: warnstr, icepack_warnings_add
      22             :       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
      23             : 
      24             :       use icepack_therm_shared, only: ferrmax, l_brine
      25             :       use icepack_therm_shared, only: surface_heat_flux, dsurface_heat_flux_dTsf
      26             : 
      27             :       implicit none
      28             : 
      29             :       private
      30             :       public :: surface_fluxes, temperature_changes
      31             : 
      32             :       real (kind=dbl_kind), parameter :: &
      33             :          betak   = 0.13_dbl_kind, & ! constant in formula for k (W m-1 ppt-1)
      34             :          kimin   = 0.10_dbl_kind    ! min conductivity of saline ice (W m-1 deg-1)
      35             : 
      36             : !=======================================================================
      37             : 
      38             :       contains
      39             : 
      40             : !=======================================================================
      41             : !
      42             : ! Compute new surface temperature and internal ice and snow
      43             : ! temperatures.  Include effects of salinity on sea ice heat
      44             : ! capacity in a way that conserves energy (Bitz and Lipscomb, 1999).
      45             : !
      46             : ! New temperatures are computed iteratively by solving a tridiagonal
      47             : ! system of equations; heat capacity is updated with each iteration.
      48             : ! Finite differencing is backward implicit.
      49             : !
      50             : ! See Bitz, C.M., and W.H. Lipscomb, 1999:
      51             : ! An energy-conserving thermodynamic model of sea ice,
      52             : ! J. Geophys. Res., 104, 15,669-15,677.
      53             : !
      54             : ! authors William H. Lipscomb, LANL
      55             : !         C. M. Bitz, UW
      56             : 
      57      884989 :       subroutine temperature_changes (dt,                 & 
      58             :                                       nilyr,    nslyr,    &
      59             :                                       rhoa,     flw,      &
      60             :                                       potT,     Qa,       &
      61             :                                       shcoef,   lhcoef,   &
      62             :                                       fswsfc,   fswint,   &
      63      884989 :                                       Sswabs,   Iswabs,   &
      64             :                                       hilyr,    hslyr,    &
      65     1769978 :                                       zqin,     zTin,     &
      66     1769978 :                                       zqsn,     zTsn,     &
      67      884989 :                                       zSin,               & 
      68             :                                       Tsf,      Tbot,     &
      69             :                                       fsensn,   flatn,    &
      70             :                                       flwoutn,  fsurfn,   &
      71             :                                       fcondtopn,fcondbot, &
      72             :                                       einit               )
      73             : 
      74             :       integer (kind=int_kind), intent(in) :: &
      75             :          nilyr , & ! number of ice layers
      76             :          nslyr     ! number of snow layers
      77             : 
      78             :       real (kind=dbl_kind), intent(in) :: &
      79             :          dt              ! time step
      80             : 
      81             :       real (kind=dbl_kind), &
      82             :          intent(in) :: &
      83             :          rhoa        , & ! air density (kg/m^3)
      84             :          flw         , & ! incoming longwave radiation (W/m^2)
      85             :          potT        , & ! air potential temperature  (K)
      86             :          Qa          , & ! specific humidity (kg/kg)
      87             :          shcoef      , & ! transfer coefficient for sensible heat
      88             :          lhcoef      , & ! transfer coefficient for latent heat
      89             :          Tbot            ! ice bottom surface temperature (deg C)
      90             : 
      91             :       real (kind=dbl_kind), &
      92             :          intent(inout) :: &
      93             :          fswsfc      , & ! SW absorbed at ice/snow surface (W m-2)
      94             :          fswint          ! SW absorbed in ice interior below surface (W m-2)
      95             : 
      96             :       real (kind=dbl_kind), intent(in) :: &
      97             :          hilyr       , & ! ice layer thickness (m)
      98             :          hslyr       , & ! snow layer thickness (m)
      99             :          einit           ! initial energy of melting (J m-2)
     100             : 
     101             :       real (kind=dbl_kind), dimension (nslyr), &
     102             :          intent(inout) :: &
     103             :          Sswabs          ! SW radiation absorbed in snow layers (W m-2)
     104             : 
     105             :       real (kind=dbl_kind), dimension (nilyr), &
     106             :          intent(inout) :: &
     107             :          Iswabs          ! SW radiation absorbed in ice layers (W m-2)
     108             : 
     109             :       real (kind=dbl_kind), intent(inout):: &
     110             :          fsurfn      , & ! net flux to top surface, excluding fcondtopn
     111             :          fcondtopn   , & ! downward cond flux at top surface (W m-2)
     112             :          fsensn      , & ! surface downward sensible heat (W m-2)
     113             :          flatn       , & ! surface downward latent heat (W m-2)
     114             :          flwoutn         ! upward LW at surface (W m-2)
     115             : 
     116             :       real (kind=dbl_kind), intent(out):: &
     117             :          fcondbot        ! downward cond flux at bottom surface (W m-2)
     118             : 
     119             :       real (kind=dbl_kind), &
     120             :          intent(inout):: &
     121             :          Tsf             ! ice/snow surface temperature, Tsfcn
     122             : 
     123             :       real (kind=dbl_kind), dimension (nilyr), &
     124             :          intent(inout) :: &
     125             :          zqin        , & ! ice layer enthalpy (J m-3)
     126             :          zTin            ! internal ice layer temperatures
     127             : 
     128             :       real (kind=dbl_kind), dimension (nilyr), &
     129             :          intent(in) :: &
     130             :          zSin            ! internal ice layer salinities
     131             : 
     132             :       real (kind=dbl_kind), dimension (nslyr), &
     133             :          intent(inout) :: &
     134             :          zqsn        , & ! snow layer enthalpy (J m-3)
     135             :          zTsn            ! internal snow layer temperatures
     136             : 
     137             :      ! local variables
     138             : 
     139             :       integer (kind=int_kind), parameter :: &
     140             :          nitermax = 100  ! max number of iterations in temperature solver
     141             : 
     142             :       real (kind=dbl_kind), parameter :: &
     143             :          Tsf_errmax = 5.e-4_dbl_kind ! max allowed error in Tsf
     144             :                                      ! recommend Tsf_errmax < 0.01 K
     145             : 
     146             :       integer (kind=int_kind) :: &
     147             :          k           , & ! ice layer index
     148             :          niter       , & ! iteration counter in temperature solver
     149             :          nmat            ! matrix dimension
     150             : 
     151             :       logical (kind=log_kind) :: &
     152             :          l_snow      , & ! true if snow temperatures are computed
     153             :          l_cold          ! true if surface temperature is computed
     154             : 
     155             :       real (kind=dbl_kind) :: &
     156      315220 :          Tsf_start   , & ! Tsf at start of iteration
     157      315220 :          dTsf        , & ! Tsf - Tsf_start
     158      315220 :          dTi1        , & ! Ti1(1) - Tin_start(1)
     159      315220 :          dfsurf_dT   , & ! derivative of fsurf wrt Tsf
     160      315220 :          avg_Tsi     , & ! = 1. if new snow/ice temps averaged w/starting temps
     161      315220 :          enew            ! new energy of melting after temp change (J m-2)
     162             : 
     163             :       real (kind=dbl_kind) :: &
     164      315220 :          dTsf_prev   , & ! dTsf from previous iteration
     165      315220 :          dTi1_prev   , & ! dTi1 from previous iteration
     166      315220 :          dfsens_dT   , & ! deriv of fsens wrt Tsf (W m-2 deg-1)
     167      315220 :          dflat_dT    , & ! deriv of flat wrt Tsf (W m-2 deg-1)
     168      315220 :          dflwout_dT  , & ! deriv of flwout wrt Tsf (W m-2 deg-1)
     169      315220 :          dt_rhoi_hlyr, & ! dt/(rhoi*hilyr)
     170      315220 :          einex       , & ! excess energy from dqmat to ocean
     171      315220 :          ferr            ! energy conservation error (W m-2)
     172             : 
     173             :       real (kind=dbl_kind), dimension (nilyr) :: &
     174     3661298 :          Tin_init    , & ! zTin at beginning of time step
     175     3661298 :          Tin_start   , & ! zTin at start of iteration
     176     3661298 :          dTmat       , & ! zTin - matrix solution before limiting
     177     3661298 :          dqmat       , & ! associated enthalpy difference
     178     3661298 :          Tmlts           ! melting temp, -depressT * salinity
     179             : 
     180             :       real (kind=dbl_kind), dimension (nslyr) :: &
     181     2092194 :          Tsn_init    , & ! zTsn at beginning of time step
     182     2092194 :          Tsn_start   , & ! zTsn at start of iteration
     183     2092194 :          etas            ! dt / (rho * cp * h) for snow layers
     184             : 
     185             :       real (kind=dbl_kind), dimension (nilyr+nslyr+1) :: &
     186     4613954 :          sbdiag      , & ! sub-diagonal matrix elements
     187     5244394 :          diag        , & ! diagonal matrix elements
     188     4613954 :          spdiag      , & ! super-diagonal matrix elements
     189     4613954 :          rhs         , & ! rhs of tri-diagonal matrix equation
     190     4613954 :          Tmat            ! matrix output temperatures
     191             : 
     192             :       real (kind=dbl_kind), dimension (nilyr) :: &
     193     3661298 :          etai            ! dt / (rho * cp * h) for ice layers
     194             : 
     195             :       real (kind=dbl_kind), dimension(nilyr+nslyr+1):: &
     196     4613954 :          kh              ! effective conductivity at interfaces (W m-2 deg-1)
     197             : 
     198             :       real (kind=dbl_kind) :: &
     199      315220 :          ci          , & ! specific heat of sea ice (J kg-1 deg-1)
     200      315220 :          avg_Tsf     , & ! = 1. if Tsf averaged w/Tsf_start, else = 0.
     201      315220 :          Iswabs_tmp  , & ! energy to melt through fraction frac of layer
     202      315220 :          Sswabs_tmp  , & ! same for snow
     203      315220 :          dswabs      , & ! difference in swabs and swabs_tmp
     204      315220 :          frac        , & ! fraction of layer that can be melted through
     205      315220 :          dTemp           ! minimum temperature difference for absorption
     206             : 
     207             :       logical (kind=log_kind) :: &
     208             :          converged       ! = true when local solution has converged
     209             : 
     210             :       logical (kind=log_kind) , dimension (nilyr) :: &
     211     1200209 :          reduce_kh       ! reduce conductivity when T exceeds Tmlt
     212             : 
     213             :       character(len=*),parameter :: subname='(temperature_changes)'
     214             : 
     215             :       !-----------------------------------------------------------------
     216             :       ! Initialize
     217             :       !-----------------------------------------------------------------
     218             : 
     219      884989 :       converged  = .false.
     220      884989 :       l_snow     = .false.
     221      884989 :       l_cold     = .true.
     222      884989 :       fcondbot   = c0
     223      884989 :       dTsf_prev  = c0
     224      884989 :       dTi1_prev  = c0
     225      884989 :       dfsens_dT  = c0
     226      884989 :       dflat_dT   = c0
     227      884989 :       dflwout_dT = c0  
     228      884989 :       einex      = c0
     229      884989 :       dt_rhoi_hlyr = dt / (rhoi*hilyr)  ! hilyr > 0
     230      884989 :       if (hslyr > hs_min/real(nslyr,kind=dbl_kind)) &
     231      490719 :            l_snow = .true.
     232             : 
     233     2678890 :       do k = 1, nslyr
     234     1793901 :          Tsn_init (k) = zTsn(k) ! beginning of time step
     235     1793901 :          Tsn_start(k) = zTsn(k) ! beginning of iteration
     236     2678890 :          if (l_snow) then
     237      618559 :             etas(k) = dt/(rhos*cp_ice*hslyr)
     238             :          else
     239     1175342 :             etas(k) = c0
     240             :          endif
     241             :       enddo                     ! k
     242             : 
     243     7079912 :       do k = 1, nilyr
     244     6194923 :          Tin_init (k) =  zTin(k)   ! beginning of time step
     245     6194923 :          Tin_start(k) =  zTin(k)   ! beginning of iteration
     246     7079912 :          Tmlts    (k) = -zSin(k) * depressT
     247             :       enddo
     248             : 
     249             :       !-----------------------------------------------------------------
     250             :       ! Compute thermal conductivity at interfaces (held fixed during
     251             :       !  subsequent iterations).
     252             :       ! Ice and snow interfaces are combined into one array (kh) to
     253             :       !  simplify the logic.
     254             :       !-----------------------------------------------------------------
     255             : 
     256             :       call conductivity (l_snow,                    &
     257             :                          nilyr,    nslyr,           &
     258             :                          hilyr,    hslyr,           &
     259      884989 :                          zTin,     kh,      zSin)    
     260      884989 :       if (icepack_warnings_aborted(subname)) return
     261             : 
     262             :       !-----------------------------------------------------------------
     263             :       ! Check for excessive absorbed solar radiation that may result in
     264             :       ! temperature overshoots. Convergence is particularly difficult
     265             :       ! if the starting temperature is already very close to the melting 
     266             :       ! temperature and extra energy is added.   In that case, or if the
     267             :       ! amount of energy absorbed is greater than the amount needed to
     268             :       ! melt through a given fraction of a layer, we put the extra 
     269             :       ! energy into the surface.
     270             :       ! NOTE: This option is not available if the atmosphere model
     271             :       !       has already computed fsurf.  (Unless we adjust fsurf here)
     272             :       !-----------------------------------------------------------------
     273             : !mclaren: Should there be an if calc_Tsfc statement here then?? 
     274             : 
     275             : #ifdef CESMCOUPLED
     276             :       frac = c1
     277             :       dTemp = p01
     278             : #else
     279      884989 :       frac = 0.9_dbl_kind
     280      884989 :       dTemp = 0.02_dbl_kind
     281             : #endif
     282      884989 :       if (solve_zsal) dTemp = p1  ! lower tolerance with dynamic salinity
     283     7079912 :       do k = 1, nilyr
     284             : 
     285     6194923 :          Iswabs_tmp = c0 ! all Iswabs is moved into fswsfc
     286     6194923 :          if (Tin_init(k) <= Tmlts(k) - dTemp) then
     287     6027544 :             if (l_brine) then
     288     6027544 :                ci = cp_ice - Lfresh * Tmlts(k) / (Tin_init(k)**2)
     289     2150264 :                Iswabs_tmp = min(Iswabs(k), &
     290     6027544 :                                 frac*(Tmlts(k)-Tin_init(k))*ci/dt_rhoi_hlyr)
     291             :             else
     292           0 :                ci = cp_ice
     293           0 :                Iswabs_tmp = min(Iswabs(k), &
     294           0 :                                 frac*(-Tin_init(k))*ci/dt_rhoi_hlyr)
     295             :             endif
     296             :          endif
     297     6194923 :          if (Iswabs_tmp < puny) Iswabs_tmp = c0
     298             : 
     299     6194923 :          dswabs = min(Iswabs(k) - Iswabs_tmp, fswint)
     300             : 
     301     6194923 :          fswsfc   = fswsfc + dswabs
     302     6194923 :          fswint   = fswint - dswabs
     303     7079912 :          Iswabs(k) = Iswabs_tmp
     304             : 
     305             :       enddo
     306             : 
     307             : #ifdef CESMCOUPLED
     308             :       frac = 0.9_dbl_kind
     309             : #endif
     310     2678890 :       do k = 1, nslyr
     311     2678890 :          if (l_snow) then
     312             : 
     313      618559 :             Sswabs_tmp = c0
     314      618559 :             if (Tsn_init(k) <= -dTemp) then
     315      210013 :                Sswabs_tmp = min(Sswabs(k), &
     316      563024 :                                 -frac*Tsn_init(k)/etas(k))
     317             :             endif
     318      618559 :             if (Sswabs_tmp < puny) Sswabs_tmp = c0
     319             : 
     320      618559 :             dswabs = min(Sswabs(k) - Sswabs_tmp, fswint)
     321             : 
     322      618559 :             fswsfc   = fswsfc + dswabs
     323      618559 :             fswint   = fswint - dswabs
     324      618559 :             Sswabs(k) = Sswabs_tmp
     325             : 
     326             :          endif
     327             :       enddo
     328             : 
     329             :       !-----------------------------------------------------------------
     330             :       ! Solve for new temperatures.
     331             :       ! Iterate until temperatures converge with minimal energy error.
     332             :       !-----------------------------------------------------------------
     333      884989 :       converged  = .false.
     334             : 
     335    89383889 :       do niter = 1, nitermax
     336             : 
     337             :       !-----------------------------------------------------------------
     338             :       ! Identify cells, if any, where calculation has not converged.
     339             :       !-----------------------------------------------------------------
     340             : 
     341    89383889 :          if (.not.converged) then
     342             : 
     343             :       !-----------------------------------------------------------------
     344             :       ! Allocate and initialize
     345             :       !-----------------------------------------------------------------
     346             : 
     347     2009203 :             converged = .true.
     348     2009203 :             dfsurf_dT = c0
     349     2009203 :             avg_Tsi   = c0
     350     2009203 :             enew      = c0
     351     2009203 :             einex     = c0
     352             : 
     353             :       !-----------------------------------------------------------------
     354             :       ! Update specific heat of ice layers.
     355             :       ! To ensure energy conservation, the specific heat is a function of
     356             :       ! both the starting temperature and the (latest guess for) the
     357             :       ! final temperature.
     358             :       !-----------------------------------------------------------------
     359             : 
     360    16073624 :             do k = 1, nilyr
     361             : 
     362    14064421 :                if (l_brine) then
     363     4993681 :                   ci = cp_ice - Lfresh*Tmlts(k) /  &
     364    14064421 :                                 (zTin(k)*Tin_init(k))
     365             :                else
     366           0 :                   ci = cp_ice
     367             :                endif
     368    16073624 :                etai(k) = dt_rhoi_hlyr / ci
     369             : 
     370             :             enddo
     371             : 
     372     2009203 :             if (calc_Tsfc) then
     373             : 
     374             :       !-----------------------------------------------------------------
     375             :       ! Update radiative and turbulent fluxes and their derivatives
     376             :       ! with respect to Tsf.
     377             :       !-----------------------------------------------------------------
     378             : 
     379             :                ! surface heat flux
     380             :                call surface_heat_flux(Tsf    , fswsfc, &
     381             :                                       rhoa   , flw   , &
     382             :                                       potT   , Qa    , &
     383             :                                       shcoef , lhcoef, &
     384             :                                       flwoutn, fsensn, &
     385     2009203 :                                       flatn  , fsurfn)
     386     2009203 :                if (icepack_warnings_aborted(subname)) return
     387             : 
     388             :                ! derivative of heat flux with respect to surface temperature
     389             :                call dsurface_heat_flux_dTsf(Tsf      , rhoa      , & 
     390             :                                             shcoef   , lhcoef    , &
     391             :                                             dfsurf_dT, dflwout_dT, &
     392     2009203 :                                             dfsens_dT, dflat_dT  )
     393     2009203 :                if (icepack_warnings_aborted(subname)) return
     394             : 
     395             :       !-----------------------------------------------------------------
     396             :       ! Compute conductive flux at top surface, fcondtopn.
     397             :       ! If fsurfn < fcondtopn and Tsf = 0, then reset Tsf to slightly less
     398             :       !  than zero (but not less than -puny).
     399             :       !-----------------------------------------------------------------
     400             :                
     401     2009203 :                if (l_snow) then
     402     1124504 :                   fcondtopn = kh(1) * (Tsf - zTsn(1))
     403             :                else
     404      884699 :                   fcondtopn = kh(1+nslyr) * (Tsf - zTin(1))
     405             :                endif
     406             : 
     407     2009203 :                if (Tsf >= c0 .and. fsurfn < fcondtopn) &
     408        9231 :                     Tsf = -puny
     409             : 
     410             :       !-----------------------------------------------------------------
     411             :       ! Save surface temperature at start of iteration
     412             :       !-----------------------------------------------------------------
     413             :                
     414     2009203 :                Tsf_start = Tsf
     415             : 
     416     2009203 :                if (Tsf < c0) then
     417     1499661 :                   l_cold = .true.
     418             :                else
     419      509542 :                   l_cold = .false.
     420             :                endif
     421             : 
     422             :       !-----------------------------------------------------------------
     423             :       ! Compute elements of tridiagonal matrix.
     424             :       !-----------------------------------------------------------------
     425             :                
     426             :                call get_matrix_elements_calc_Tsfc (nilyr, nslyr, &
     427             :                                    l_snow,      l_cold,      &
     428             :                                    Tsf,         Tbot,        &
     429             :                                    fsurfn,      dfsurf_dT,   &
     430           0 :                                    Tin_init,    Tsn_init,    &
     431           0 :                                    kh,          Sswabs,      &
     432           0 :                                    Iswabs,                   &
     433           0 :                                    etai,        etas,        &
     434             :                                    sbdiag,      diag,        &
     435     2009203 :                                    spdiag,      rhs)   
     436     2009203 :                if (icepack_warnings_aborted(subname)) return
     437             : 
     438             :             else
     439             :                
     440             :                call get_matrix_elements_know_Tsfc (nilyr, nslyr, &
     441             :                                    l_snow,      Tbot,        &
     442           0 :                                    Tin_init,    Tsn_init,    &
     443           0 :                                    kh,          Sswabs,      &
     444           0 :                                    Iswabs,                   &
     445           0 :                                    etai,        etas,        &
     446             :                                    sbdiag,      diag,        &
     447             :                                    spdiag,      rhs,         &
     448           0 :                                    fcondtopn)
     449           0 :                if (icepack_warnings_aborted(subname)) return
     450             : 
     451             :             endif  ! calc_Tsfc
     452             : 
     453             :       !-----------------------------------------------------------------
     454             :       ! Solve tridiagonal matrix to obtain the new temperatures.
     455             :       !-----------------------------------------------------------------
     456             : 
     457     2009203 :             nmat = nslyr + nilyr + 1  ! matrix dimension
     458             : 
     459           0 :             call tridiag_solver (nmat,    sbdiag(:),   &
     460           0 :                                  diag(:), spdiag(:),   &
     461     2009203 :                                  rhs(:),  Tmat(:))
     462     2009203 :             if (icepack_warnings_aborted(subname)) return
     463             : 
     464             :       !-----------------------------------------------------------------
     465             :       ! Determine whether the computation has converged to an acceptable
     466             :       ! solution.  Five conditions must be satisfied:
     467             :       !
     468             :       !    (1) Tsf <= 0 C.
     469             :       !    (2) Tsf is not oscillating; i.e., if both dTsf(niter) and
     470             :       !        dTsf(niter-1) have magnitudes greater than puny, then
     471             :       !        dTsf(niter)/dTsf(niter-1) cannot be a negative number
     472             :       !        with magnitude greater than 0.5.  
     473             :       !    (3) abs(dTsf) < Tsf_errmax
     474             :       !    (4) If Tsf = 0 C, then the downward turbulent/radiative 
     475             :       !        flux, fsurfn, must be greater than or equal to the downward
     476             :       !        conductive flux, fcondtopn.
     477             :       !    (5) The net energy added to the ice per unit time must equal 
     478             :       !        the net change in internal ice energy per unit time,
     479             :       !        withinic the prescribed error ferrmax.
     480             :       !
     481             :       ! For briny ice (the standard case), zTsn and zTin are limited
     482             :       !  to prevent them from exceeding their melting temperatures.
     483             :       !  (Note that the specific heat formula for briny ice assumes
     484             :       !  that T < Tmlt.)  
     485             :       ! For fresh ice there is no limiting, since there are cases
     486             :       !  when the only convergent solution has zTsn > 0 and/or zTin > 0.
     487             :       !  Above-zero temperatures are then reset to zero (with melting 
     488             :       !  to conserve energy) in the thickness_changes subroutine.
     489             :       !-----------------------------------------------------------------
     490             : 
     491     2009203 :             if (calc_Tsfc) then
     492             : 
     493             :       !-----------------------------------------------------------------
     494             :       ! Reload Tsf from matrix solution
     495             :       !-----------------------------------------------------------------
     496             : 
     497     2009203 :                if (l_cold) then
     498     1499661 :                   if (l_snow) then
     499     1081554 :                      Tsf = Tmat(1)
     500             :                   else
     501      418107 :                      Tsf = Tmat(1+nslyr)
     502             :                   endif
     503             :                else                ! melting surface
     504      509542 :                   Tsf = c0
     505             :                endif
     506             : 
     507             :       !-----------------------------------------------------------------
     508             :       ! Initialize convergence flag (true until proven false), dTsf,
     509             :       !  and temperature-averaging coefficients.
     510             :       ! Average only if test 1 or 2 fails.
     511             :       ! Initialize energy.
     512             :       !-----------------------------------------------------------------
     513             :                
     514     2009203 :                dTsf = Tsf - Tsf_start
     515     2009203 :                avg_Tsf  = c0
     516             : 
     517             :       !-----------------------------------------------------------------
     518             :       ! Condition 1: check for Tsf > 0
     519             :       ! If Tsf > 0, set Tsf = 0, then average zTsn and zTin to force
     520             :       ! internal temps below their melting temps.
     521             :       !-----------------------------------------------------------------
     522             : 
     523     2009203 :                if (Tsf > puny) then
     524      123153 :                   Tsf = c0
     525      123153 :                   dTsf = -Tsf_start
     526      123153 :                   if (l_brine) avg_Tsi = c1   ! avg with starting temp
     527      123153 :                   converged = .false.
     528             : 
     529             :       !-----------------------------------------------------------------
     530             :       ! Condition 2: check for oscillating Tsf
     531             :       ! If oscillating, average all temps to increase rate of convergence.
     532             :       !-----------------------------------------------------------------
     533             : 
     534             :                elseif (niter > 1 &                ! condition (2)
     535             :                     .and. Tsf_start <= -puny &
     536             :                     .and. abs(dTsf) > puny &
     537             :                     .and. abs(dTsf_prev) > puny &
     538     1886050 :                     .and. -dTsf/(dTsf_prev+puny*puny) > p5) then
     539             : 
     540        1835 :                   if (l_brine) then ! average with starting temp
     541        1835 :                      avg_Tsf  = c1    
     542        1835 :                      avg_Tsi = c1
     543             :                   endif
     544        1835 :                   dTsf = p5 * dTsf
     545        1835 :                   converged = .false.
     546             :                endif
     547             : 
     548             : !!!            dTsf_prev = dTsf
     549             : 
     550             :       !-----------------------------------------------------------------
     551             :       ! If condition 2 failed, average new surface temperature with
     552             :       !  starting value.
     553             :       !-----------------------------------------------------------------
     554             :                Tsf  = Tsf &
     555     2009203 :                     + avg_Tsf * p5 * (Tsf_start - Tsf)
     556             : 
     557             :             endif   ! calc_Tsfc
     558             : 
     559     5924974 :             do k = 1, nslyr
     560             : 
     561             :       !-----------------------------------------------------------------
     562             :       ! Reload zTsn from matrix solution
     563             :       !-----------------------------------------------------------------
     564             : 
     565     3915771 :                if (l_snow) then
     566     1398680 :                   zTsn(k) = Tmat(k+1)
     567             :                else
     568     2517091 :                   zTsn(k) = c0
     569             :                endif
     570     3915771 :                if (l_brine) zTsn(k) = min(zTsn(k), c0)
     571             : 
     572             :       !-----------------------------------------------------------------
     573             :       ! If condition 1 or 2 failed, average new snow layer
     574             :       !  temperatures with their starting values.
     575             :       !-----------------------------------------------------------------
     576     1378807 :                zTsn(k) = zTsn(k) &
     577     3915771 :                        + avg_Tsi*p5*(Tsn_start(k)-zTsn(k))
     578             : 
     579             :       !-----------------------------------------------------------------
     580             :       ! Compute zqsn and increment new energy.
     581             :       !-----------------------------------------------------------------
     582     3915771 :                zqsn(k) = -rhos * (Lfresh - cp_ice*zTsn(k))
     583     3915771 :                enew  = enew + hslyr * zqsn(k)
     584             : 
     585     5924974 :                Tsn_start(k) = zTsn(k) ! for next iteration
     586             : 
     587             :             enddo                  ! nslyr
     588             : 
     589    16073624 :             dTmat(:) = c0
     590    16073624 :             dqmat(:) = c0
     591    16073624 :             reduce_kh(:) = .false.
     592    16073624 :             do k = 1, nilyr
     593             : 
     594             :       !-----------------------------------------------------------------
     595             :       ! Reload zTin from matrix solution
     596             :       !-----------------------------------------------------------------
     597             : 
     598    14064421 :                zTin(k) = Tmat(k+1+nslyr)
     599             : 
     600    14064421 :                if (l_brine .and. zTin(k) > Tmlts(k) - puny) then
     601           3 :                   dTmat(k) = zTin(k) - Tmlts(k)
     602           1 :                   dqmat(k) = rhoi * dTmat(k) &
     603           3 :                            * (cp_ice - Lfresh * Tmlts(k)/zTin(k)**2)
     604             : ! use this for the case that Tmlt changes by an amount dTmlt=Tmltnew-Tmlt(k)
     605             : !                             + rhoi * dTmlt &
     606             : !                             * (cp_ocn - cp_ice + Lfresh/zTin(k))
     607           3 :                   zTin(k) = Tmlts(k)
     608           3 :                   reduce_kh(k) = .true.
     609             :                endif
     610             : 
     611             :       !-----------------------------------------------------------------
     612             :       ! Condition 2b: check for oscillating zTin(1)
     613             :       ! If oscillating, average all ice temps to increase rate of convergence.
     614             :       !-----------------------------------------------------------------
     615             : 
     616    14064421 :                if (k==1 .and. .not.calc_Tsfc) then
     617           0 :                   dTi1 = zTin(k) - Tin_start(k)
     618             : 
     619             :                   if (niter > 1 &                    ! condition 2b    
     620             :                       .and. abs(dTi1) > puny &
     621             :                       .and. abs(dTi1_prev) > puny &
     622           0 :                       .and. -dTi1/(dTi1_prev+puny*puny) > p5) then
     623             : 
     624           0 :                      if (l_brine) avg_Tsi = c1
     625           0 :                      dTi1 = p5 * dTi1
     626           0 :                      converged = .false.
     627             :                   endif
     628           0 :                   dTi1_prev = dTi1
     629             :                endif   ! k = 1 .and. calc_Tsfc = F
     630             : 
     631             :       !-----------------------------------------------------------------
     632             :       ! If condition 1 or 2 failed, average new ice layer
     633             :       !  temperatures with their starting values.
     634             :       !-----------------------------------------------------------------
     635     4993681 :                zTin(k) = zTin(k) &
     636    14064421 :                        + avg_Tsi*p5*(Tin_start(k)-zTin(k))
     637             : 
     638             :       !-----------------------------------------------------------------
     639             :       ! Compute zqin and increment new energy.
     640             :       !-----------------------------------------------------------------
     641    14064421 :                if (l_brine) then
     642     4993681 :                   zqin(k) = -rhoi * (cp_ice*(Tmlts(k)-zTin(k)) &
     643     4993681 :                                    + Lfresh*(c1-Tmlts(k)/zTin(k)) &
     644    14064421 :                                    - cp_ocn*Tmlts(k))
     645             :                else
     646           0 :                   zqin(k) = -rhoi * (-cp_ice*zTin(k) + Lfresh)
     647             :                endif
     648    14064421 :                enew = enew + hilyr * zqin(k)
     649    14064421 :                einex = einex + hilyr * dqmat(k)
     650             : 
     651    16073624 :                Tin_start(k) = zTin(k) ! for next iteration
     652             : 
     653             :             enddo                  ! nilyr
     654             : 
     655     2009203 :             if (calc_Tsfc) then
     656             : 
     657             :       !-----------------------------------------------------------------
     658             :       ! Condition 3: check for large change in Tsf
     659             :       !-----------------------------------------------------------------
     660             :                
     661     2009203 :                if (abs(dTsf) > Tsf_errmax) then
     662      697711 :                   converged = .false.
     663             :                endif
     664             : 
     665             :       !-----------------------------------------------------------------
     666             :       ! Condition 4: check for fsurfn < fcondtopn with Tsf >= 0
     667             :       !-----------------------------------------------------------------
     668             :                
     669     2009203 :                fsurfn = fsurfn + dTsf*dfsurf_dT
     670     2009203 :                if (l_snow) then
     671     1124504 :                   fcondtopn = kh(1) * (Tsf-zTsn(1))
     672             :                else
     673      884699 :                   fcondtopn = kh(1+nslyr) * (Tsf-zTin(1))
     674             :                endif
     675             : 
     676     2009203 :                if (Tsf >= c0 .and. fsurfn < fcondtopn) then
     677        2722 :                   converged = .false.
     678             :                endif
     679             : 
     680     2009203 :                dTsf_prev = dTsf
     681             : 
     682             :             endif                   ! calc_Tsfc
     683             : 
     684             :       !-----------------------------------------------------------------
     685             :       ! Condition 5: check for energy conservation error
     686             :       ! Change in internal ice energy should equal net energy input.
     687             :       !-----------------------------------------------------------------
     688             : 
     689      713383 :             fcondbot = kh(1+nslyr+nilyr) * &
     690     2722586 :                        (zTin(nilyr) - Tbot)
     691             : 
     692             :             ! Flux extra energy out of the ice
     693     2009203 :             fcondbot = fcondbot + einex/dt 
     694             : 
     695             :             ferr = abs( (enew-einit)/dt &
     696     2009203 :                  - (fcondtopn - fcondbot + fswint) )
     697             : 
     698             :             ! factor of 0.9 allows for roundoff errors later
     699     2009203 :             if (ferr > 0.9_dbl_kind*ferrmax) then         ! condition (5)
     700             : 
     701      894629 :                converged = .false.
     702             : 
     703             :                ! reduce conductivity for next iteration
     704     7157032 :                do k = 1, nilyr
     705     7157032 :                   if (reduce_kh(k) .and. dqmat(k) > c0) then
     706           3 :                      frac = max(0.5*(c1-ferr/abs(fcondtopn-fcondbot)),p1)
     707             : !                     frac = p1
     708           3 :                      kh(k+nslyr+1) = kh(k+nslyr+1) * frac
     709           3 :                      kh(k+nslyr)   = kh(k+nslyr+1)
     710             :                   endif
     711             :                enddo
     712             : 
     713             :             endif               ! ferr 
     714             : 
     715             :          endif ! convergence
     716             : 
     717             :       enddo                     ! temperature iteration niter
     718             : 
     719             :       !-----------------------------------------------------------------
     720             :       ! Check for convergence failures.
     721             :       !-----------------------------------------------------------------
     722      884989 :       if (.not.converged) then
     723           0 :          write(warnstr,*) subname, 'Thermo iteration does not converge,'
     724           0 :          call icepack_warnings_add(warnstr)
     725           0 :          write(warnstr,*) subname, 'Ice thickness:',  hilyr*nilyr
     726           0 :          call icepack_warnings_add(warnstr)
     727           0 :          write(warnstr,*) subname, 'Snow thickness:', hslyr*nslyr
     728           0 :          call icepack_warnings_add(warnstr)
     729           0 :          write(warnstr,*) subname, 'dTsf, Tsf_errmax:',dTsf_prev, &
     730           0 :               Tsf_errmax
     731           0 :          call icepack_warnings_add(warnstr)
     732           0 :          write(warnstr,*) subname, 'Tsf:', Tsf
     733           0 :          call icepack_warnings_add(warnstr)
     734           0 :          write(warnstr,*) subname, 'fsurf:', fsurfn
     735           0 :          call icepack_warnings_add(warnstr)
     736           0 :          write(warnstr,*) subname, 'fcondtop, fcondbot, fswint', &
     737           0 :               fcondtopn, fcondbot, fswint
     738           0 :          call icepack_warnings_add(warnstr)
     739           0 :          write(warnstr,*) subname, 'fswsfc', fswsfc
     740           0 :          call icepack_warnings_add(warnstr)
     741           0 :          write(warnstr,*) subname, 'Iswabs',(Iswabs(k),k=1,nilyr)
     742           0 :          call icepack_warnings_add(warnstr)
     743           0 :          write(warnstr,*) subname, 'Flux conservation error =', ferr
     744           0 :          call icepack_warnings_add(warnstr)
     745           0 :          write(warnstr,*) subname, 'Initial snow temperatures:'
     746           0 :          call icepack_warnings_add(warnstr)
     747           0 :          write(warnstr,*) subname, (Tsn_init(k),k=1,nslyr)
     748           0 :          call icepack_warnings_add(warnstr)
     749           0 :          write(warnstr,*) subname, 'Initial ice temperatures:'
     750           0 :          call icepack_warnings_add(warnstr)
     751           0 :          write(warnstr,*) subname, (Tin_init(k),k=1,nilyr)
     752           0 :          call icepack_warnings_add(warnstr)
     753           0 :          write(warnstr,*) subname, 'Matrix ice temperature diff:'
     754           0 :          call icepack_warnings_add(warnstr)
     755           0 :          write(warnstr,*) subname, (dTmat(k),k=1,nilyr)
     756           0 :          call icepack_warnings_add(warnstr)
     757           0 :          write(warnstr,*) subname, 'dqmat*hilyr/dt:'
     758           0 :          call icepack_warnings_add(warnstr)
     759           0 :          write(warnstr,*) subname, (hilyr*dqmat(k)/dt,k=1,nilyr)
     760           0 :          call icepack_warnings_add(warnstr)
     761           0 :          write(warnstr,*) subname, 'Final snow temperatures:'
     762           0 :          call icepack_warnings_add(warnstr)
     763           0 :          write(warnstr,*) subname, (zTsn(k),k=1,nslyr)
     764           0 :          call icepack_warnings_add(warnstr)
     765           0 :          write(warnstr,*) subname, 'Matrix ice temperature diff:'
     766           0 :          call icepack_warnings_add(warnstr)
     767           0 :          write(warnstr,*) subname, (dTmat(k),k=1,nilyr)
     768           0 :          call icepack_warnings_add(warnstr)
     769           0 :          write(warnstr,*) subname, 'dqmat*hilyr/dt:'
     770           0 :          call icepack_warnings_add(warnstr)
     771           0 :          write(warnstr,*) subname, (hilyr*dqmat(k)/dt,k=1,nilyr)
     772           0 :          call icepack_warnings_add(warnstr)
     773           0 :          write(warnstr,*) subname, 'Final ice temperatures:'
     774           0 :          call icepack_warnings_add(warnstr)
     775           0 :          write(warnstr,*) subname, (zTin(k),k=1,nilyr)
     776           0 :          call icepack_warnings_add(warnstr)
     777           0 :          write(warnstr,*) subname, 'Ice melting temperatures:'
     778           0 :          call icepack_warnings_add(warnstr)
     779           0 :          write(warnstr,*) subname, (Tmlts(k),k=1,nilyr)
     780           0 :          call icepack_warnings_add(warnstr)
     781           0 :          write(warnstr,*) subname, 'Ice bottom temperature:', Tbot
     782           0 :          call icepack_warnings_add(warnstr)
     783           0 :          write(warnstr,*) subname, 'dT initial:'
     784           0 :          call icepack_warnings_add(warnstr)
     785           0 :          write(warnstr,*) subname, (Tmlts(k)-Tin_init(k),k=1,nilyr)
     786           0 :          call icepack_warnings_add(warnstr)
     787           0 :          write(warnstr,*) subname, 'dT final:'
     788           0 :          call icepack_warnings_add(warnstr)
     789           0 :          write(warnstr,*) subname, (Tmlts(k)-zTin(k),k=1,nilyr)
     790           0 :          call icepack_warnings_add(warnstr)
     791           0 :          write(warnstr,*) subname, 'zSin'
     792           0 :          call icepack_warnings_add(warnstr)
     793           0 :          write(warnstr,*) subname, (zSin(k),k=1,nilyr)
     794           0 :          call icepack_warnings_add(warnstr)
     795           0 :          call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     796           0 :          call icepack_warnings_add(subname//" temperature_changes: Thermo iteration does not converge" ) 
     797           0 :          return
     798             :       endif
     799             : 
     800      884989 :       if (calc_Tsfc) then
     801             : 
     802             :          ! update fluxes that depend on Tsf
     803      884989 :          flwoutn = flwoutn + dTsf_prev * dflwout_dT
     804      884989 :          fsensn  = fsensn  + dTsf_prev * dfsens_dT
     805      884989 :          flatn   = flatn   + dTsf_prev * dflat_dT
     806             : 
     807             :       endif                        ! calc_Tsfc
     808             : 
     809             :       end subroutine temperature_changes
     810             : 
     811             : !=======================================================================
     812             : !
     813             : ! Compute thermal conductivity at interfaces (held fixed during
     814             : !  the subsequent iteration).
     815             : !
     816             : ! NOTE: Ice conductivity must be >= kimin
     817             : !
     818             : ! authors William H. Lipscomb, LANL
     819             : !         C. M. Bitz, UW
     820             : 
     821      884989 :       subroutine conductivity (l_snow,                  &
     822             :                                nilyr,    nslyr,         &
     823             :                                hilyr,    hslyr,         &
     824      884989 :                                zTin,     kh,       zSin)
     825             : 
     826             :       logical (kind=log_kind), intent(in) :: &
     827             :          l_snow          ! true if snow temperatures are computed
     828             : 
     829             :       integer (kind=int_kind), intent(in) :: & 
     830             :          nilyr , & ! number of ice layers
     831             :          nslyr     ! number of snow layers
     832             : 
     833             :       real (kind=dbl_kind), intent(in) :: &
     834             :          hilyr       , & ! ice layer thickness (same for all ice layers)
     835             :          hslyr           ! snow layer thickness (same for all snow layers)
     836             : 
     837             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
     838             :          zTin         , & ! internal ice layer temperatures
     839             :          zSin             ! internal ice layer salinities
     840             : 
     841             :       real (kind=dbl_kind), dimension (nilyr+nslyr+1), &
     842             :          intent(out) :: &
     843             :          kh              ! effective conductivity at interfaces (W m-2 deg-1)
     844             : 
     845             :       ! local variables
     846             : 
     847             :       integer (kind=int_kind) :: &
     848             :          k               ! vertical index
     849             : 
     850             :       real (kind=dbl_kind), dimension (nilyr) :: &
     851     3661298 :          kilyr           ! thermal cond at ice layer midpoints (W m-1 deg-1)
     852             : 
     853             :       real (kind=dbl_kind), dimension (nslyr) :: &
     854     2092194 :          kslyr           ! thermal cond at snow layer midpoints (W m-1 deg-1)
     855             : 
     856             :       character(len=*),parameter :: subname='(conductivity)'
     857             : 
     858             :       ! interior snow layers (simple for now, but may be fancier later)
     859     2678890 :       do k = 1, nslyr
     860     2678890 :          kslyr(k) = ksno
     861             :       enddo                     ! nslyr
     862             : 
     863             :       ! interior ice layers
     864      884989 :       if (conduct == 'MU71') then
     865             :          ! Maykut and Untersteiner 1971 form (with Wettlaufer 1991 constants)
     866     2812016 :          do k = 1, nilyr
     867     2460514 :             kilyr(k) = kice + betak*zSin(k)/min(-puny,zTin(k))
     868     2812016 :             kilyr(k) = max (kilyr(k), kimin)
     869             :          enddo                     ! nilyr
     870             :       else
     871             :          ! Pringle et al JGR 2007 'bubbly brine'
     872     4267896 :          do k = 1, nilyr
     873     1321901 :             kilyr(k) = (2.11_dbl_kind - 0.011_dbl_kind*zTin(k) &
     874     1321901 :                      + 0.09_dbl_kind*zSin(k)/min(-puny,zTin(k))) &
     875     3734409 :                      * rhoi / 917._dbl_kind
     876     4267896 :             kilyr(k) = max (kilyr(k), kimin)
     877             :          enddo                     ! nilyr
     878             :       endif ! conductivity
     879             : 
     880             :       ! top snow interface, top and bottom ice interfaces
     881             :          ! top of snow layer; top surface of top ice layer
     882      884989 :       if (l_snow) then
     883      490719 :          kh(1)       = c2 * kslyr(1) / hslyr
     884      347398 :          kh(1+nslyr) = c2 * kslyr(nslyr) * kilyr(1) / &
     885      173699 :                        ( kslyr(nslyr)*hilyr +  &
     886     1011816 :                          kilyr(1    )*hslyr )
     887             :       else
     888      394270 :          kh(1)       = c0
     889      394270 :          kh(1+nslyr) = c2 * kilyr(1) / hilyr
     890             :       endif
     891             : 
     892             :       ! bottom surface of bottom ice layer
     893      884989 :       kh(1+nslyr+nilyr) = c2 * kilyr(nilyr) / hilyr
     894             : 
     895             :       ! interior snow interfaces
     896             : 
     897      884989 :       if (nslyr > 1) then
     898     1136140 :          do k = 2, nslyr
     899     1136140 :             if (l_snow) then
     900      127840 :                kh(k) = c2 * kslyr(k-1) * kslyr(k) / &
     901      191760 :                          ((kslyr(k-1) + kslyr(k))*hslyr)
     902             :             else
     903      781072 :                kh(k) = c0
     904             :             endif
     905             :          enddo                     ! nilyr
     906             :       endif ! nslyr > 1
     907             : 
     908             :       ! interior ice interfaces
     909     6194923 :       do k = 2, nilyr
     910     3782640 :          kh(k+nslyr) = c2 * kilyr(k-1) * kilyr(k) / &
     911     9977563 :                          ((kilyr(k-1) + kilyr(k))*hilyr)
     912             :       enddo                     ! nilyr
     913             : 
     914      884989 :       end subroutine conductivity
     915             : 
     916             : !=======================================================================
     917             : !
     918             : ! Compute radiative and turbulent fluxes and their derivatives
     919             : ! with respect to Tsf.
     920             : !
     921             : ! authors William H. Lipscomb, LANL
     922             : !         C. M. Bitz, UW
     923             : 
     924      417295 :       subroutine surface_fluxes (Tsf,        fswsfc,            &
     925             :                                  rhoa,       flw,               &
     926             :                                  potT,       Qa,                &
     927             :                                  shcoef,     lhcoef,            &
     928             :                                  flwoutn,    fsensn,            &
     929             :                                  flatn,      fsurfn,            &
     930             :                                  dflwout_dT, dfsens_dT,         &
     931             :                                  dflat_dT,   dfsurf_dT)
     932             : 
     933             :       real (kind=dbl_kind), intent(in) :: &
     934             :          Tsf             ! ice/snow surface temperature, Tsfcn
     935             : 
     936             :       real (kind=dbl_kind), intent(in) :: &
     937             :          fswsfc      , & ! SW absorbed at ice/snow surface (W m-2)
     938             :          rhoa        , & ! air density (kg/m^3)
     939             :          flw         , & ! incoming longwave radiation (W/m^2)
     940             :          potT        , & ! air potential temperature  (K)
     941             :          Qa          , & ! specific humidity (kg/kg)
     942             :          shcoef      , & ! transfer coefficient for sensible heat
     943             :          lhcoef          ! transfer coefficient for latent heat
     944             : 
     945             :       real (kind=dbl_kind), &
     946             :          intent(inout) :: &
     947             :          fsensn      , & ! surface downward sensible heat (W m-2)
     948             :          flatn       , & ! surface downward latent heat (W m-2)
     949             :          flwoutn     , & ! upward LW at surface (W m-2)
     950             :          fsurfn          ! net flux to top surface, excluding fcondtopn
     951             : 
     952             :       real (kind=dbl_kind), &
     953             :          intent(inout) :: &
     954             :          dfsens_dT   , & ! deriv of fsens wrt Tsf (W m-2 deg-1)
     955             :          dflat_dT    , & ! deriv of flat wrt Tsf (W m-2 deg-1)
     956             :          dflwout_dT      ! deriv of flwout wrt Tsf (W m-2 deg-1)
     957             : 
     958             :       real (kind=dbl_kind), &
     959             :          intent(inout) :: &
     960             :          dfsurf_dT       ! derivative of fsurfn wrt Tsf
     961             : 
     962             :       character(len=*),parameter :: subname='(surface_fluxes)'
     963             : 
     964             :       ! surface heat flux
     965             :       call surface_heat_flux(Tsf,     fswsfc, &
     966             :                              rhoa,    flw,    &
     967             :                              potT,    Qa,     &
     968             :                              shcoef,  lhcoef, &
     969             :                              flwoutn, fsensn, &
     970      417295 :                              flatn,   fsurfn)
     971      417295 :       if (icepack_warnings_aborted(subname)) return
     972             : 
     973             :       ! derivative of heat flux with respect to surface temperature
     974             :       call dsurface_heat_flux_dTsf(Tsf,       rhoa,       &
     975             :                                    shcoef,    lhcoef,     &
     976             :                                    dfsurf_dT, dflwout_dT, &
     977      417295 :                                    dfsens_dT, dflat_dT)
     978      417295 :       if (icepack_warnings_aborted(subname)) return
     979             : 
     980             :       end subroutine surface_fluxes
     981             : 
     982             : !=======================================================================
     983             : !
     984             : ! Compute terms in tridiagonal matrix that will be solved to find
     985             : !  the new vertical temperature profile
     986             : ! This routine is for the case in which Tsfc is being computed.
     987             : !
     988             : ! authors William H. Lipscomb, LANL
     989             : !         C. M. Bitz, UW
     990             : !
     991             : ! March 2004 by William H. Lipscomb for multiple snow layers
     992             : ! April 2008 by E. C. Hunke, divided into two routines based on calc_Tsfc 
     993             : 
     994     2009203 :       subroutine get_matrix_elements_calc_Tsfc (nilyr, nslyr, &
     995             :                                       l_snow,   l_cold,           &
     996             :                                       Tsf,      Tbot,             &
     997             :                                       fsurfn,   dfsurf_dT,        &
     998     4018406 :                                       Tin_init, Tsn_init,         &
     999     2009203 :                                       kh,       Sswabs,           &
    1000     2009203 :                                       Iswabs,                     &
    1001     4018406 :                                       etai,     etas,             &
    1002     2009203 :                                       sbdiag,   diag,             &
    1003     2009203 :                                       spdiag,   rhs)
    1004             : 
    1005             :       integer (kind=int_kind), intent(in) :: & 
    1006             :          nilyr , & ! number of ice layers
    1007             :          nslyr     ! number of snow layers
    1008             : 
    1009             :       logical (kind=log_kind), &
    1010             :          intent(in) :: &
    1011             :          l_snow      , & ! true if snow temperatures are computed
    1012             :          l_cold          ! true if surface temperature is computed
    1013             : 
    1014             :       real (kind=dbl_kind), intent(in) :: &
    1015             :          Tsf             ! ice/snow top surface temp (deg C)
    1016             : 
    1017             :       real (kind=dbl_kind), intent(in) :: &
    1018             :          fsurfn      , & ! net flux to top surface, excluding fcondtopn (W/m^2)
    1019             :          Tbot            ! ice bottom surface temperature (deg C)
    1020             : 
    1021             :       real (kind=dbl_kind), intent(in) :: &
    1022             :          dfsurf_dT       ! derivative of fsurf wrt Tsf
    1023             : 
    1024             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
    1025             :          etai        , & ! dt / (rho*cp*h) for ice layers
    1026             :          Tin_init    , & ! ice temp at beginning of time step
    1027             :          Sswabs      , & ! SW radiation absorbed in snow layers (W m-2)
    1028             :          Iswabs      , & ! absorbed SW flux in ice layers
    1029             :          etas        , & ! dt / (rho*cp*h) for snow layers
    1030             :          Tsn_init        ! snow temp at beginning of time step
    1031             :                          ! Note: no absorbed SW in snow layers
    1032             : 
    1033             :       real (kind=dbl_kind), dimension (nslyr+nilyr+1), &
    1034             :          intent(in) :: &
    1035             :          kh              ! effective conductivity at layer interfaces
    1036             : 
    1037             :       real (kind=dbl_kind), dimension (nslyr+nilyr+1), &
    1038             :          intent(inout) :: &
    1039             :          sbdiag      , & ! sub-diagonal matrix elements
    1040             :          diag        , & ! diagonal matrix elements
    1041             :          spdiag      , & ! super-diagonal matrix elements
    1042             :          rhs             ! rhs of tri-diagonal matrix eqn.
    1043             : 
    1044             :       ! local variables
    1045             : 
    1046             :       integer (kind=int_kind) :: &
    1047             :          k, ki, kr       ! vertical indices and row counters
    1048             : 
    1049             :       character(len=*),parameter :: subname='(get_matrix_elements_calc_Tsrf)'
    1050             : 
    1051             :       !-----------------------------------------------------------------
    1052             :       ! Initialize matrix elements.
    1053             :       ! Note: When we do not need to solve for the surface or snow
    1054             :       !       temperature, we solve dummy equations with solution T = 0.
    1055             :       !       Ice layers are fully initialized below.
    1056             :       !-----------------------------------------------------------------
    1057             : 
    1058     7934177 :       do k = 1, nslyr+1
    1059     5924974 :          sbdiag(k) = c0
    1060     5924974 :          diag  (k) = c1
    1061     5924974 :          spdiag(k) = c0
    1062     7934177 :          rhs   (k) = c0
    1063             :       enddo
    1064             :             
    1065             :       !-----------------------------------------------------------------
    1066             :       ! Compute matrix elements
    1067             :       !
    1068             :       ! Four possible cases to solve:
    1069             :       !   (1) Cold surface (Tsf < 0), snow present
    1070             :       !   (2) Melting surface (Tsf = 0), snow present
    1071             :       !   (3) Cold surface (Tsf < 0), no snow
    1072             :       !   (4) Melting surface (Tsf = 0), no snow
    1073             :       !-----------------------------------------------------------------
    1074             : 
    1075             :       !-----------------------------------------------------------------
    1076             :       ! Tsf equation for case of cold surface (with or without snow)
    1077             :       !-----------------------------------------------------------------
    1078             : 
    1079     2009203 :       if (l_cold) then
    1080     1499661 :          if (l_snow) then
    1081     1081554 :             k = 1
    1082             :          else                ! no snow
    1083      418107 :             k = 1 + nslyr
    1084             :          endif
    1085     1499661 :          kr = k
    1086             : 
    1087     1499661 :          sbdiag(kr) = c0
    1088     1499661 :          diag  (kr) = dfsurf_dT - kh(k)
    1089     1499661 :          spdiag(kr) = kh(k)
    1090     1499661 :          rhs   (kr) = dfsurf_dT*Tsf - fsurfn
    1091             :       endif                  ! l_cold
    1092             : 
    1093             :       !-----------------------------------------------------------------
    1094             :       ! top snow layer
    1095             :       !-----------------------------------------------------------------
    1096             : !           k = 1
    1097             : !           kr = 2
    1098             : 
    1099     2009203 :       if (l_snow) then
    1100     1124504 :          if (l_cold) then
    1101     1081554 :             sbdiag(2) = -etas(1) * kh(1)
    1102     1081554 :             spdiag(2) = -etas(1) * kh(2)
    1103      376739 :             diag  (2) = c1 &
    1104     1081554 :                       + etas(1) * (kh(1) + kh(2))
    1105      376739 :             rhs   (2) = Tsn_init(1) &
    1106     1081554 :                       + etas(1) * Sswabs(1)
    1107             :          else                ! melting surface
    1108       42950 :             sbdiag(2) = c0
    1109       42950 :             spdiag(2) = -etas(1) * kh(2)
    1110       21163 :             diag  (2) = c1 &
    1111       42950 :                       + etas(1) * (kh(1) + kh(2))
    1112       21163 :             rhs   (2) = Tsn_init(1) &
    1113       21163 :                       + etas(1)*kh(1)*Tsf &
    1114       42950 :                       + etas(1) * Sswabs(1)
    1115             :          endif               ! l_cold
    1116             :       endif                  ! l_snow
    1117             : 
    1118             :       !-----------------------------------------------------------------
    1119             :       ! remaining snow layers
    1120             :       !-----------------------------------------------------------------
    1121             : 
    1122     2009203 :       if (nslyr > 1) then
    1123             : 
    1124     2383210 :          do k = 2, nslyr
    1125     1906568 :             kr = k + 1
    1126             : 
    1127     2383210 :             if (l_snow) then
    1128      274176 :                sbdiag(kr) = -etas(k) * kh(k)
    1129      274176 :                spdiag(kr) = -etas(k) * kh(k+1)
    1130      137088 :                diag  (kr) = c1 &
    1131      274176 :                           + etas(k) * (kh(k) + kh(k+1))
    1132      137088 :                rhs   (kr) = Tsn_init(k) &
    1133      274176 :                           + etas(k) * Sswabs(k)
    1134             :             endif
    1135             :          enddo                  ! nslyr
    1136             : 
    1137             :       endif                     ! nslyr > 1
    1138             : 
    1139     2009203 :       if (nilyr > 1) then
    1140             : 
    1141             :       !-----------------------------------------------------------------
    1142             :       ! top ice layer
    1143             :       !-----------------------------------------------------------------
    1144             : 
    1145     2009203 :          ki = 1
    1146     2009203 :          k  = ki + nslyr
    1147     2009203 :          kr = k + 1
    1148             : 
    1149     2009203 :          if (l_snow .or. l_cold) then
    1150     1542611 :             sbdiag(kr) = -etai(ki) * kh(k)
    1151     1542611 :             spdiag(kr) = -etai(ki) * kh(k+1)
    1152      532831 :             diag  (kr) = c1 &
    1153     1542611 :                        + etai(ki) * (kh(k) + kh(k+1))
    1154      532831 :             rhs   (kr) = Tin_init(ki) &
    1155     1542611 :                        + etai(ki)*Iswabs(ki)
    1156             :          else    ! no snow, warm surface
    1157      466592 :             sbdiag(kr) = c0
    1158      466592 :             spdiag(kr) = -etai(ki) * kh(k+1)
    1159      180552 :             diag  (kr) = c1 &
    1160      466592 :                        + etai(ki) * (kh(k) + kh(k+1))
    1161      180552 :             rhs   (kr) = Tin_init(ki) &
    1162      180552 :                        + etai(ki)*Iswabs(ki) &
    1163      466592 :                        + etai(ki)*kh(k)*Tsf
    1164             :          endif
    1165             : 
    1166             :       !-----------------------------------------------------------------
    1167             :       ! bottom ice layer
    1168             :       !-----------------------------------------------------------------
    1169             : 
    1170     2009203 :          ki = nilyr
    1171     2009203 :          k  = ki + nslyr
    1172     2009203 :          kr = k + 1
    1173             :  
    1174     2009203 :          sbdiag(kr) = -etai(ki) * kh(k)
    1175     2009203 :          spdiag(kr) = c0
    1176      713383 :          diag  (kr) = c1  &
    1177     2009203 :                     + etai(ki) * (kh(k) + kh(k+1))
    1178      713383 :          rhs   (kr) = Tin_init(ki) &
    1179      713383 :                     + etai(ki)*Iswabs(ki) &
    1180     2009203 :                     + etai(ki)*kh(k+1)*Tbot
    1181             :       
    1182             :       else         ! nilyr = 1
    1183             : 
    1184             :       !-----------------------------------------------------------------
    1185             :       ! single ice layer
    1186             :       !-----------------------------------------------------------------
    1187             : 
    1188           0 :          ki = 1
    1189           0 :          k  = ki + nslyr
    1190           0 :          kr = k + 1
    1191             : 
    1192           0 :          if (l_snow .or. l_cold) then
    1193           0 :             sbdiag(kr) = -etai(ki) * kh(k)
    1194           0 :             spdiag(kr) = c0
    1195           0 :             diag  (kr) = c1 &
    1196           0 :                        + etai(ki) * (kh(k) + kh(k+1))
    1197           0 :             rhs   (kr) = Tin_init(ki) &
    1198           0 :                        + etai(ki) * Iswabs(ki) &
    1199           0 :                        + etai(ki) * kh(k+1)*Tbot
    1200             :          else   ! no snow, warm surface
    1201           0 :             sbdiag(kr) = c0
    1202           0 :             spdiag(kr) = c0
    1203           0 :             diag  (kr) = c1 &
    1204           0 :                        + etai(ki) * (kh(k) + kh(k+1))
    1205           0 :             rhs   (kr) = Tin_init(ki) &
    1206           0 :                        + etai(ki) * Iswabs(ki) &
    1207           0 :                        + etai(ki) * kh(k)*Tsf &
    1208           0 :                        + etai(ki) * kh(k+1)*Tbot
    1209             :          endif
    1210             :  
    1211             :       endif        ! nilyr > 1
    1212             : 
    1213             :       !-----------------------------------------------------------------
    1214             :       ! interior ice layers
    1215             :       !-----------------------------------------------------------------
    1216             : 
    1217    12055218 :       do ki = 2, nilyr-1
    1218             :            
    1219    10046015 :          k  = ki + nslyr
    1220    10046015 :          kr = k + 1
    1221             : 
    1222    10046015 :             sbdiag(kr) = -etai(ki) * kh(k)
    1223    10046015 :             spdiag(kr) = -etai(ki) * kh(k+1)
    1224     3566915 :             diag  (kr) = c1 &
    1225    10046015 :                        + etai(ki) * (kh(k) + kh(k+1))
    1226     3566915 :             rhs   (kr) = Tin_init(ki) &
    1227    12055218 :                        + etai(ki)*Iswabs(ki)
    1228             :       enddo                     ! nilyr
    1229             : 
    1230     2009203 :       end subroutine get_matrix_elements_calc_Tsfc
    1231             : 
    1232             : !=======================================================================
    1233             : !
    1234             : ! Compute terms in tridiagonal matrix that will be solved to find
    1235             : !  the new vertical temperature profile
    1236             : ! This routine is for the case in which Tsfc is already known.
    1237             : !
    1238             : ! authors William H. Lipscomb, LANL
    1239             : !         C. M. Bitz, UW
    1240             : !
    1241             : ! March 2004 by William H. Lipscomb for multiple snow layers
    1242             : ! April 2008 by E. C. Hunke, divided into two routines based on calc_Tsfc 
    1243             : 
    1244           0 :       subroutine get_matrix_elements_know_Tsfc (nilyr, nslyr, &
    1245             :                                       l_snow,   Tbot,             &
    1246           0 :                                       Tin_init, Tsn_init,         &
    1247           0 :                                       kh,       Sswabs,           &
    1248           0 :                                       Iswabs,                     &
    1249           0 :                                       etai,     etas,             &
    1250           0 :                                       sbdiag,   diag,             &
    1251           0 :                                       spdiag,   rhs,              &
    1252             :                                       fcondtopn)
    1253             : 
    1254             :       integer (kind=int_kind), intent(in) :: & 
    1255             :          nilyr , & ! number of ice layers
    1256             :          nslyr     ! number of snow layers
    1257             : 
    1258             :       logical (kind=log_kind), &
    1259             :          intent(in) :: &
    1260             :          l_snow          ! true if snow temperatures are computed
    1261             : 
    1262             :       real (kind=dbl_kind), intent(in) :: &
    1263             :          Tbot            ! ice bottom surface temperature (deg C)
    1264             : 
    1265             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
    1266             :          etai        , & ! dt / (rho*cp*h) for ice layers
    1267             :          Tin_init    , & ! ice temp at beginning of time step
    1268             :          Sswabs      , & ! SW radiation absorbed in snow layers (W m-2)
    1269             :          Iswabs      , & ! absorbed SW flux in ice layers
    1270             :          etas        , & ! dt / (rho*cp*h) for snow layers
    1271             :          Tsn_init        ! snow temp at beginning of time step
    1272             :                          ! Note: no absorbed SW in snow layers
    1273             : 
    1274             :       real (kind=dbl_kind), dimension (nslyr+nilyr+1), &
    1275             :          intent(in) :: &
    1276             :          kh              ! effective conductivity at layer interfaces
    1277             : 
    1278             :       real (kind=dbl_kind), dimension (nslyr+nilyr+1), &
    1279             :          intent(inout) :: &
    1280             :          sbdiag      , & ! sub-diagonal matrix elements
    1281             :          diag        , & ! diagonal matrix elements
    1282             :          spdiag      , & ! super-diagonal matrix elements
    1283             :          rhs             ! rhs of tri-diagonal matrix eqn.
    1284             : 
    1285             :       real (kind=dbl_kind), intent(in),  &
    1286             :          optional :: &
    1287             :          fcondtopn       ! conductive flux at top sfc, positive down (W/m^2)
    1288             : 
    1289             :       ! local variables
    1290             : 
    1291             :       integer (kind=int_kind) :: &
    1292             :          k, ki, kr       ! vertical indices and row counters
    1293             : 
    1294             :       character(len=*),parameter :: subname='(get_matrix_elements_know_Tsrf)'
    1295             : 
    1296             :       !-----------------------------------------------------------------
    1297             :       ! Initialize matrix elements.
    1298             :       ! Note: When we do not need to solve for the surface or snow
    1299             :       !       temperature, we solve dummy equations with solution T = 0.
    1300             :       !       Ice layers are fully initialized below.
    1301             :       !-----------------------------------------------------------------
    1302             : 
    1303           0 :       do k = 1, nslyr+1
    1304           0 :          sbdiag(k) = c0
    1305           0 :          diag  (k) = c1
    1306           0 :          spdiag(k) = c0
    1307           0 :          rhs   (k) = c0
    1308             :       enddo
    1309             :             
    1310             :       !-----------------------------------------------------------------
    1311             :       ! Compute matrix elements
    1312             :       !
    1313             :       ! Four possible cases to solve:
    1314             :       !   (1) Cold surface (Tsf < 0), snow present
    1315             :       !   (2) Melting surface (Tsf = 0), snow present
    1316             :       !   (3) Cold surface (Tsf < 0), no snow
    1317             :       !   (4) Melting surface (Tsf = 0), no snow
    1318             :       !-----------------------------------------------------------------
    1319             : 
    1320             :       !-----------------------------------------------------------------
    1321             :       ! top snow layer
    1322             :       !-----------------------------------------------------------------
    1323             : !        k = 1
    1324             : !        kr = 2
    1325             : 
    1326           0 :       if (l_snow) then
    1327           0 :          sbdiag(2) = c0
    1328           0 :          spdiag(2) = -etas(1) * kh(2)
    1329           0 :          diag  (2) = c1 &
    1330           0 :                    + etas(1) * kh(2)
    1331           0 :          rhs   (2) = Tsn_init(1) &
    1332           0 :                    + etas(1) * Sswabs(1) &
    1333           0 :                    + etas(1) * fcondtopn
    1334             :       endif   ! l_snow
    1335             : 
    1336             :       !-----------------------------------------------------------------
    1337             :       ! remaining snow layers
    1338             :       !-----------------------------------------------------------------
    1339             : 
    1340           0 :       if (nslyr > 1) then
    1341             : 
    1342           0 :          do k = 2, nslyr
    1343           0 :             kr = k + 1
    1344             : 
    1345           0 :             if (l_snow) then
    1346           0 :                sbdiag(kr) = -etas(k) * kh(k)
    1347           0 :                spdiag(kr) = -etas(k) * kh(k+1)
    1348           0 :                diag  (kr) = c1 &
    1349           0 :                           + etas(k) * (kh(k) + kh(k+1))
    1350           0 :                rhs   (kr) = Tsn_init(k) &
    1351           0 :                           + etas(k) * Sswabs(k)
    1352             :             endif
    1353             : 
    1354             :          enddo                  ! nslyr
    1355             : 
    1356             :       endif                     ! nslyr > 1
    1357             : 
    1358           0 :       if (nilyr > 1) then
    1359             : 
    1360             :       !-----------------------------------------------------------------
    1361             :       ! top ice layer
    1362             :       !-----------------------------------------------------------------
    1363             : 
    1364           0 :          ki = 1
    1365           0 :          k  = ki + nslyr
    1366           0 :          kr = k + 1
    1367             : 
    1368           0 :          if (l_snow) then
    1369             : 
    1370           0 :             sbdiag(kr) = -etai(ki) * kh(k)
    1371           0 :             spdiag(kr) = -etai(ki) * kh(k+1)
    1372           0 :             diag  (kr) = c1 &
    1373           0 :                        + etai(ki) * (kh(k) + kh(k+1))
    1374           0 :             rhs   (kr) = Tin_init(ki) &
    1375           0 :                        + etai(ki) * Iswabs(ki)
    1376             :          else                  
    1377           0 :             sbdiag(kr) = c0
    1378           0 :             spdiag(kr) = -etai(ki) * kh(k+1)
    1379           0 :             diag  (kr) = c1 &
    1380           0 :                        + etai(ki) * kh(k+1)
    1381           0 :             rhs   (kr) = Tin_init(ki) &
    1382           0 :                        + etai(ki) * Iswabs(ki) &
    1383           0 :                        + etai(ki) * fcondtopn
    1384             :          endif  ! l_snow
    1385             : 
    1386             :       !-----------------------------------------------------------------
    1387             :       ! bottom ice layer
    1388             :       !-----------------------------------------------------------------
    1389             : 
    1390           0 :          ki = nilyr
    1391           0 :          k  = ki + nslyr
    1392           0 :          kr = k + 1
    1393             :       
    1394           0 :          sbdiag(kr) = -etai(ki) * kh(k)
    1395           0 :          spdiag(kr) = c0
    1396           0 :          diag  (kr) = c1  &
    1397           0 :                     + etai(ki) * (kh(k) + kh(k+1))
    1398           0 :          rhs   (kr) = Tin_init(ki) &
    1399           0 :                     + etai(ki)*Iswabs(ki) &
    1400           0 :                     + etai(ki)*kh(k+1)*Tbot
    1401             :       
    1402             :       else         ! nilyr = 1
    1403             : 
    1404             :       !-----------------------------------------------------------------
    1405             :       ! single ice layer
    1406             :       !-----------------------------------------------------------------
    1407             : 
    1408           0 :          ki = 1
    1409           0 :          k  = ki + nslyr
    1410           0 :          kr = k + 1
    1411             : 
    1412           0 :          if (l_snow) then
    1413           0 :             sbdiag(kr) = -etai(ki) * kh(k)
    1414           0 :             spdiag(kr) = c0
    1415           0 :             diag  (kr) = c1 &
    1416           0 :                        + etai(ki) * (kh(k) + kh(k+1))
    1417           0 :             rhs   (kr) = Tin_init(ki) &
    1418           0 :                        + etai(ki) * Iswabs(ki) &
    1419           0 :                        + etai(ki) * kh(k+1)*Tbot
    1420             :          else
    1421           0 :             sbdiag(kr) = c0
    1422           0 :             spdiag(kr) = c0
    1423           0 :             diag  (kr) = c1 &
    1424           0 :                        + etai(ki) * kh(k+1)
    1425           0 :             rhs   (kr) = Tin_init(ki) &
    1426           0 :                        + etai(ki) * Iswabs(ki) &
    1427           0 :                        + etai(ki) * fcondtopn &
    1428           0 :                        + etai(ki) * kh(k+1)*Tbot
    1429             :          endif
    1430             : 
    1431             :       endif        ! nilyr > 1
    1432             : 
    1433             :       !-----------------------------------------------------------------
    1434             :       ! interior ice layers
    1435             :       !-----------------------------------------------------------------
    1436             : 
    1437           0 :       do ki = 2, nilyr-1
    1438             :            
    1439           0 :          k  = ki + nslyr
    1440           0 :          kr = k + 1
    1441             : 
    1442           0 :          sbdiag(kr) = -etai(ki) * kh(k)
    1443           0 :          spdiag(kr) = -etai(ki) * kh(k+1)
    1444           0 :          diag  (kr) = c1 &
    1445           0 :                     + etai(ki) * (kh(k) + kh(k+1))
    1446           0 :          rhs   (kr) = Tin_init(ki) &
    1447           0 :                     + etai(ki)*Iswabs(ki)
    1448             : 
    1449             :       enddo                     ! nilyr
    1450             : 
    1451           0 :       end subroutine get_matrix_elements_know_Tsfc
    1452             : 
    1453             : !=======================================================================
    1454             : !
    1455             : ! Tridiagonal matrix solver--used to solve the implicit vertical heat
    1456             : ! equation in ice and snow
    1457             : !
    1458             : ! authors William H. Lipscomb, LANL
    1459             : !         C. M. Bitz, UW
    1460             : 
    1461     2009203 :       subroutine tridiag_solver (nmat,     sbdiag,   &
    1462     4018406 :                                  diag,     spdiag,   &
    1463     2009203 :                                  rhs,      xout)
    1464             : 
    1465             :       integer (kind=int_kind), intent(in) :: &
    1466             :          nmat            ! matrix dimension
    1467             : 
    1468             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
    1469             :          sbdiag      , & ! sub-diagonal matrix elements
    1470             :          diag        , & ! diagonal matrix elements
    1471             :          spdiag      , & ! super-diagonal matrix elements
    1472             :          rhs             ! rhs of tri-diagonal matrix eqn.
    1473             : 
    1474             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
    1475             :          xout            ! solution vector
    1476             : 
    1477             :       ! local variables
    1478             : 
    1479             :       integer (kind=int_kind) :: &
    1480             :          k               ! row counter
    1481             : 
    1482             :       real (kind=dbl_kind) :: &
    1483      713383 :          wbeta           ! temporary matrix variable
    1484             : 
    1485             :       real (kind=dbl_kind), dimension(nmat) :: &
    1486    10390894 :          wgamma          ! temporary matrix variable
    1487             : 
    1488             :       character(len=*),parameter :: subname='(tridiag_solver)'
    1489             : 
    1490     2009203 :       wbeta = diag(1)
    1491     2009203 :       xout(1) = rhs(1) / wbeta
    1492             : 
    1493    19989395 :       do k = 2, nmat
    1494    17980192 :          wgamma(k) = spdiag(k-1) / wbeta
    1495    17980192 :          wbeta = diag(k) - sbdiag(k)*wgamma(k)
    1496    12744976 :          xout(k) = (rhs(k) - sbdiag(k)*xout(k-1)) &
    1497    26361883 :                     / wbeta
    1498             :       enddo                     ! k
    1499             : 
    1500    19989395 :       do k = nmat-1, 1, -1
    1501    19989395 :          xout(k) = xout(k) - wgamma(k+1)*xout(k+1)
    1502             :       enddo                     ! k
    1503             : 
    1504     2009203 :       end subroutine tridiag_solver
    1505             : 
    1506             : !=======================================================================
    1507             : 
    1508             :       end module icepack_therm_bl99
    1509             : 
    1510             : !=======================================================================

Generated by: LCOV version 1.14-6-g40580cd