4.30 Fortran: Module Interface ice_tstm - energy conserving sea ice model thermodynamics (Source File: ice_tstm.F)

Energy-conserving sea ice model
Routines to solve heat-equation using linear method

See Bitz, C.M., and W.H. Lipscomb, 1999: An energy-conserving thermodynamic model of sea ice, J. Geophys. Res., 104, 15,669-15,677.

See Bitz, C.M., M.M. Holland, A.J. Weaver, M. Eby, 2001: Simulating the ice-thickness distribution in a climate model, J. Geophys. Res., 106, 2441-2464.

REVISION HISTORY:

   author: C. M. Bitz

INTERFACE:

       module ice_tstm

USES:

       use ice_constants
       use ice_itd

4.30.1 tstm - calculates surface and internal snow/ice temp

INTERFACE:

       subroutine tstm(   dtau,  tmz1d,  sal1d,     Tf1
      $                ,  area,     hi,     hs,    dswr
      $                , dswrv,  dswrn,   flwd,  dflwup
      $                ,  dflh,   dfsh,  asnow,      ni
      $                ,  tbot,     Ib,      F,   condb
      $                ,    ts,     ti,  flwup,     flh
      $                ,   fsh)

DESCRIPTION:

This routine calculates the evolution of the ice interior and surface temperature from the heat equation and surface energy balance

The albedo is fixed for this calculation

Solves the heat equation which is non-linear for saline ice by linearizing and then iterating a set of equations Scheme is backwards Euler giving a tridiagonal set of equations implicit in temperature
(Tried crank-nicholson but it behaves poorly when thin ice has a weird initial temperature profile.

NOTE there are two values for minimum snow thickness, hs-min used elsewhere in the code, and hsmin used here. CC says:
The solution to the heat equation ignores the insulating effect of snow if it less than 0.01 m thick, but I do not like to "kill" it when it is that thick because sometimes the snowfall rate is really small...
Must have a sufficient amount of snow to solve heat equation in snow hsmin is the minimum depth of snow in order to solve for ti(0) if snow thickness is less than hsmin then do not change ti(0)

The number of equations that must be solved by the tridiagonal solver depends on whether the surface is melting and whether there is snow. Four cases are possible:

1 = freezing w/ snow, 2 = freezing w/ no snow,
3 = melting w/ snow, and 4 = melting w/ no snow

REVISION HISTORY:

   author: C. M. Bitz

USES:

INPUT/OUTPUT PARAMETERS:

       integer (kind=int_kind), intent(in) :: ni ! number vertical layers in ice
       real (kind=dbl_kind), intent(in) ::
      &   dtau     ! timestep                                     (s)          
      &,  tmz1d(ni)! melting temp of each layer                   (C)
      &,  sal1d(ni+1) ! salinity of each layer                    (ppt)
      &,  area     ! area of the ice/snow
      &,  hi,hs    ! ice and snow thickness                       (m)
      &,  asnow    ! 1 - patchy snow frac
      &,  dswr     ! above srfc net dnwd shortwave, positive down (W/m**2)
      &,  dswrv    ! dswr in vis (wvlngth < 700nm)                (W/m**2)
      &,  dswrn    ! dswr in nir (wvlngth > 700nm)                (W/m**2)
      &,  flwd     ! dnwd longwave flux (always positive)         (W/m**2)
      &,  dfsh     ! derriv wrt ts of dnwd sensible flux          (W/m**2)
      &,  dflh     ! derriv wrt ts of dnwd latent flux            (W/m**2)
      &,  dflwup   ! derriv wrt ts of upwd longwave flux          (W/m**2)
      &,  Tf1      ! freezing temp of water below ice             (C)
 
       real (kind=dbl_kind), intent(out) ::
      &   tbot     ! bottom ice temp                              (C)
      &,  condb    ! conductive flux at bottom srfc               (W/m**2)
      &,  F        ! net flux at top srfc including conductive flux (W/m**2)
      &,  Ib       ! solar passing through the bottom ice srfc    (W/m**2)
 
       real (kind=dbl_kind), intent(inout) ::
      &   ts       ! srfc temp of snow or ice  (C)
      &,  ti(0:ni) ! snow(0) and ice(1:ni) interior temp (C)
      &,  fsh      ! dnwd sensible flux (W/m**2)
      &,  flh      ! dnwd latent flux (always negative) (W/m**2)
      &,  flwup    ! upwd longwave flux (always negative) (W/m**2)

4.30.2 getabc - computes elements of tridiagonal matrix

INTERFACE:

       subroutine getabc(a,b,c,r,ti,tbot,zeta,k,eta,ni,lfirst)

DESCRIPTION:

Compute elements of tridiagonal matrix

REVISION HISTORY:

   author: C. M. Bitz

USES:

INPUT/OUTPUT PARAMETERS:

       integer (kind=int_kind), intent(in) :: ni   ! number of layers
      $       ,lfirst        ! start with this layer 
 
       real (kind=dbl_kind), intent(in) ::
      &        ti   (0:ni)     ! temperature of ice-snow layers
      $       ,tbot            ! temperature of ice bottom srfc
      $       ,zeta (0:nmax)   ! 
      $       ,k    (0:nmax+1) ! ice-snow conductivitiy
      $       ,eta  (0:nmax)   !
 
       real (kind=dbl_kind), intent(out) ::
      &        a    (-1:nmax)  ! sub-diagonal elements
      $       ,b    (-1:nmax)  ! diagonal elements
      $       ,c    (-1:nmax)  ! super-diagonal elements
      $       ,r    (-1:nmax)  ! constants (indep. of ti)

4.30.3 tridag - solves tridiagonal equations

INTERFACE:

       subroutine tridag(a,b,c,r,u,ni)

DESCRIPTION:

Solves tridiagonal equations

REVISION HISTORY:

   author: C. M. Bitz

USES:

INPUT/OUTPUT PARAMETERS:

       integer (kind=int_kind), intent(in) :: ni  ! number of rows
       real (kind=dbl_kind), intent(in) ::
      &        a      (nmax)  ! sub-diagonal elements
      $       ,b      (nmax)  ! diagonal elements
      $       ,c      (nmax)  ! super-diagonal elements
      $       ,r      (nmax)  ! constants (indep. of ti)
 
       real (kind=dbl_kind), intent(out) ::
      &        u      (nmax)  ! solution

4.30.4 conductiv - calculates T,S dependent conductivity

INTERFACE:

       subroutine conductiv(sal1d,hsmin,ki,ti,tbot,hs,dhi,ni)

DESCRIPTION:

Calculates T,S dependent conductivity

REVISION HISTORY:

   author: C. M. Bitz

USES:

INPUT/OUTPUT PARAMETERS:

       integer (kind=int_kind), intent(in) :: ni    ! number of layers
       real (kind=dbl_kind), intent(in) ::
      &   hsmin         ! minimum allowable snow thickness for heat eq.
      &,  sal1d(ni+1)   ! salinity of each layer
      &,  hs,ti(0:ni),tbot,dhi
       real (kind=dbl_kind), intent(out) ::
      &   ki(0:nmax+1)