#if defined(L08_1A)
! *****************************COPYRIGHT*******************************

! (c) [University of Edinburgh] [2009]. All rights reserved.
! This routine has been licensed to the Met Office for use and
! distribution under the JULES collaboration agreement, subject
! to the terms and conditions set out therein.
! [Met Office Ref SC237]

! *****************************COPYRIGHT*******************************
!  SUBROUTINE SNOWPACK-------------------------------------------------

! Description:
!     Snow thermodynamics and hydrology

! Subroutine Interface:
SUBROUTINE snowpack ( land_pts,tile_pts,timestep,cansnowtile,     &
                      nsnow,tile_index,csnow,ei_tile,hcaps,hcons, &
                      ksnow,rho_snow_grnd,smcl1,snowfall,sthf1,   &
                      surf_htf_tile,tile_frac,v_sat1,ds,          &
                      melt_tile,sice,sliq,snomlt_sub_htf,         &
                      snowdepth,snowmass,tsnow,tsoil1,            &
                      snow_soil_htf,rho_snow,rho0,sice0,tsnow0 )

USE c_0_dg_c, ONLY :                                              &
!  imported scalar parameters
 tm              ! Temperature at which fresh water freezes
!                      ! and ice melts (K)

USE c_densty, ONLY :                                              &
!  imported scalar parameters
 rho_water       ! density of pure water (kg/m3)

USE c_lheat , ONLY :                                              &
!  imported scalar parameters
 lc                                                               &
                 ! latent heat of condensation of water
!                      ! at 0degC (J kg-1)
,lf              !  latent heat of fusion at 0degC (J kg-1)

USE c_perma, ONLY :                                               &
!  imported scalar parameters
 hcapi                                                            &
                 !  Specific heat capacity of ice (J/kg/K)
,hcapw           !  Specific heat capacity of water (J/kg/K)

USE soil_param, ONLY :                                            &
!  imported arrays with intent(in)
 dzsoil          !  Thicknesses of the soil layers (m)

USE ancil_info, ONLY :                                            &
!  imported scalars with intent(in)
 nsmax           ! Maximum possible number of snow layers

USE snow_param, ONLY :                                            &
 rho_snow_const                                                   &
                 !   constant density of lying snow (kg per m**3)
,snowliqcap      !   Liquid water holding capacity of lying snow
!                          as a fraction of snow mass.

USE parkind1, ONLY: jprb, jpim
USE yomhook, ONLY: lhook, dr_hook
IMPLICIT NONE

! Scalar parameters
REAL, PARAMETER :: GAMMA = 0.5
                 ! Implicit timestep weighting

! Scalar arguments with intent(in)
INTEGER, INTENT(IN) ::                                            &
 land_pts                                                         &
                 ! Total number of land points
,tile_pts        ! Number of tile points

REAL, INTENT(IN) ::                                               &
 timestep        ! Timestep (s)

LOGICAL, INTENT(IN) ::                                            &
 cansnowtile     ! Switch for canopy snow model

! Array arguments with intent(in)
INTEGER, INTENT(IN) ::                                            &
 nsnow(land_pts)                                                  &
                 ! Number of snow layers
,tile_index(land_pts)
                 ! Index of tile points

REAL, INTENT(IN) ::                                               &
 csnow(land_pts,nsmax)                                            &
                 ! Areal heat capacity of layers (J/K/m2)
,ei_tile(land_pts)                                                &
                 ! Sublimation of snow (kg/m2/s)
,hcaps(land_pts)                                                  &
                 ! Heat capacity of soil surface layer (J/K/m3)
,hcons(land_pts)                                                  &
                 ! Thermal conductivity of top soil layer,
!                      ! including water and ice (W/m/K)
,ksnow(land_pts,nsmax)                                            &
                 ! Thermal conductivity of layers (W/m/K)
,rho_snow_grnd(land_pts)                                          &
                 ! Snowpack bulk density (kg/m3)
,smcl1(land_pts)                                                  &
                 ! Moisture content of surface soil layer (kg/m2)
,snowfall(land_pts)                                               &
                 ! Snow reaching the ground (kg/m2)
,sthf1(land_pts)                                                  &
                 ! Frozen soil moisture content of surface layer
!                      ! as a fraction of saturation.
,surf_htf_tile(land_pts)                                          &
                 ! Snow surface heat flux (W/m2)
,tile_frac(land_pts)                                              &
                 ! Tile fractions
,v_sat1(land_pts)
                 ! Surface soil layer volumetric
!                      ! moisture concentration at saturation

! Array arguments with intent(inout)
REAL, INTENT(INOUT) ::                                            &
 ds(land_pts,nsmax)                                               &
                 ! Snow layer depths (m)
,melt_tile(land_pts)                                              &
                 ! Surface snowmelt rate (kg/m2/s)
,sice(land_pts,nsmax)                                             &
                 ! Ice content of snow layers (kg/m2)
,sliq(land_pts,nsmax)                                             &
                 ! Liquid content of snow layers (kg/m2)
,snomlt_sub_htf(land_pts)                                         &
                 ! Sub-canopy snowmelt heat flux (W/m2)
,snowdepth(land_pts)                                              &
                 ! Snow depth (m)
,snowmass(land_pts)                                               &
                 ! Snow mass on the ground (kg/m2)
,tsnow(land_pts,nsmax)                                            &
                 ! Snow layer temperatures (K)
,tsoil1(land_pts)
                 ! Soil surface layer temperature(K)

! Array arguments with intent(out)
REAL, INTENT(OUT) ::                                              &
 snow_soil_htf(land_pts)                                          &
                 ! Heat flux into the uppermost subsurface layer
!                      !        (W/m2)
!                      ! i.e. snow to ground, or into snow/soil
!                      !      composite layer
,rho_snow(land_pts,nsmax)                                         &
                 ! Density of snow layers (kg/m3)
,rho0(land_pts)                                                   &
                 ! Density of fresh snow (kg/m3)
,sice0(land_pts)                                                  &
                 ! Ice content of fresh snow (kg/m2)
!                      ! Where NSNOW=0, SICE0 is the mass of
!                      ! the snowpack.
,tsnow0(land_pts)
                 ! Temperature of fresh snow (K)

! Local scalars
INTEGER ::                                                        &
 i                                                                &
                 ! land point index
,k                                                                &
                 ! Tile point index
,n               ! Snow layer index

REAL ::                                                           &
 asoil                                                            &
                 ! 1 / (dz*hcap) for surface soil layer
,can_melt                                                         &
                 ! Melt of snow on the canopy (kg/m2/s)
,coldsnow                                                         &
                 ! layer cold content (J/m2)
,dsice                                                            &
                 ! Change in layer ice content (kg/m2)
,g_snow_surf                                                      &
                 ! Heat flux at the snow surface (W/m2)
,sliqmax                                                          &
                 ! Maximum liquid content for layer (kg/m2)
,submelt                                                          &
                 ! Melt of snow beneath canopy (kg/m2/s)
,smclf                                                            &
                 ! Frozen soil moisture content of
!                      ! surface layer (kg/m2)
,win             ! Water entering layer (kg/m2)

! Local arrays
REAL ::                                                           &
 asnow(nsmax)                                                     &
                 ! Effective thermal conductivity (W/m2/k)
,a(nsmax)                                                         &
                 ! Below-diagonal matrix elements
,b(nsmax)                                                         &
                 ! Diagonal matrix elements
,c(nsmax)                                                         &
                 ! Above-diagonal matrix elements
,dt(nsmax)                                                        &
                 ! Temperature increments (k)
,r(nsmax)        ! Matrix equation rhs

INTEGER(KIND=jpim), PARAMETER :: zhook_in  = 0
INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
REAL(KIND=jprb)               :: zhook_handle

!-----------------------------------------------------------------------
! External routines called:
!-----------------------------------------------------------------------
EXTERNAL tridag


!-----------------------------------------------------------------------

IF (lhook) CALL dr_hook('SNOWPACK',zhook_in,zhook_handle)
DO k=1,tile_pts
  i = tile_index(k)

  g_snow_surf = surf_htf_tile(i)
!-----------------------------------------------------------------------
! Add melt to snow surface heat flux, unless using the snow canopy model
!-----------------------------------------------------------------------
  IF ( .NOT. cansnowtile )                                        &
    g_snow_surf = g_snow_surf + lf*melt_tile(i)

  IF ( nsnow(i) == 0 ) THEN

! Add snowfall (including canopy unloading) to ground snowpack.
    snowmass(i) = snowmass(i) + snowfall(i)

    IF ( .NOT. cansnowtile ) THEN
!-----------------------------------------------------------------------
! Remove sublimation and melt from snowpack.
!-----------------------------------------------------------------------
      snowmass(i) = snowmass(i) -                                 &
                        ( ei_tile(i) + melt_tile(i) ) * timestep

    ELSE IF ( tsoil1(i) > tm ) THEN
!-----------------------------------------------------------------------
! For canopy model, calculate melt of snow on ground underneath canopy.
!-----------------------------------------------------------------------
      smclf = rho_water * dzsoil(1) * v_sat1(i) * sthf1(i)
      asoil = 1./ ( dzsoil(1) * hcaps(i) +                        &
                      hcapw * (smcl1(i) - smclf) + hcapi*smclf )
      submelt = MIN( snowmass(i) / timestep,                      &
                     (tsoil1(i) - tm) / (lf * asoil * timestep) )
      snowmass(i) = snowmass(i) - submelt * timestep
      tsoil1(i) = tsoil1(i) -                                     &
                   tile_frac(i) * asoil * timestep * lf * submelt
      melt_tile(i) = melt_tile(i) + submelt
      snomlt_sub_htf(i) = snomlt_sub_htf(i) +                     &
                             tile_frac(i) * lf * submelt

    END IF

! Set flux into uppermost snow/soil layer (after melting).
    snow_soil_htf(i) = surf_htf_tile(i)

! Diagnose snow depth.
    snowdepth(i) = snowmass(i) / rho_snow_grnd(i)

    rho0(i) = rho_snow_const
    sice0(i) = snowmass(i)
    tsnow0(i) = MIN( tsoil1(i), tm )

  ELSE

!-----------------------------------------------------------------------
! There is at least one snow layer. Calculate heat conduction between
!   layers and temperature increments.
!-----------------------------------------------------------------------

! Save rate of melting of snow on the canopy.
    IF ( cansnowtile ) THEN
      can_melt = melt_tile(i)
    ELSE
      can_melt = 0.0
    END IF

    IF ( nsnow(i) == 1 ) THEN

! Single layer of snow.
      asnow(1) = 2.0 /                                            &
                 ( snowdepth(i)/ksnow(i,1) + dzsoil(1)/hcons(i) )
      snow_soil_htf(i) = asnow(1) * ( tsnow(i,1) - tsoil1(i) )
      dt(1) = ( g_snow_surf - snow_soil_htf(i) ) * timestep /     &
                     ( csnow(i,1) + GAMMA * asnow(1) * timestep )
      snow_soil_htf(i) = asnow(1) *                                 &
                         ( tsnow(i,1) + GAMMA*dt(1) - tsoil1(i) )
      tsnow(i,1) = tsnow(i,1) + dt(1)

    ELSE

! Multiple snow layers.
      DO n=1,nsnow(i)-1
        asnow(n) = 2.0 /                                          &
                   ( ds(i,n)/ksnow(i,n) + ds(i,n+1)/ksnow(i,n+1) )
      END DO
      n = nsnow(i)
      asnow(n) = 2.0 / ( ds(i,n)/ksnow(i,n) + dzsoil(1)/hcons(i) )
      a(1) = 0.
      b(1) = csnow(i,1) + GAMMA*asnow(1)*timestep
      c(1) = -GAMMA * asnow(1) * timestep
      r(1) = ( g_snow_surf - asnow(1)*(tsnow(i,1)-tsnow(i,2)) )   &
                                                      * timestep
      DO n=2,nsnow(i)-1
        a(n) = -GAMMA * asnow(n-1) * timestep
        b(n) = csnow(i,n) + GAMMA * ( asnow(n-1) + asnow(n) )     &
                                                       * timestep
        c(n) = -GAMMA * asnow(n) * timestep
        r(n) = asnow(n-1)*(tsnow(i,n-1) - tsnow(i,n) ) * timestep &
               +  asnow(n)*(tsnow(i,n+1) - tsnow(i,n)) * timestep
      END DO
      n = nsnow(i)
      a(n) = -GAMMA * asnow(n-1) * timestep
      b(n) = csnow(i,n) + GAMMA * (asnow(n-1)+asnow(n)) * timestep
      c(n) = 0.
      r(n) = asnow(n-1)*( tsnow(i,n-1) - tsnow(i,n) ) * timestep  &
              + asnow(n) * ( tsoil1(i) - tsnow(i,n) ) * timestep

!-----------------------------------------------------------------------
! Call the tridiagonal solver.
!-----------------------------------------------------------------------
! DEPENDS ON: tridag
      CALL tridag( nsnow(i),nsmax,a,b,c,r,dt )

      n = nsnow(i)
      snow_soil_htf(i) = asnow(n) *                               &
                         ( tsnow(i,n) + GAMMA*dt(n) - tsoil1(i) )
      DO n=1,nsnow(i)
        tsnow(i,n) = tsnow(i,n) + dt(n)
      END DO
    END IF  !  NSNOW

!-----------------------------------------------------------------------
! Melt snow in layers with temperature exceeding melting point
!-----------------------------------------------------------------------
    DO n=1,nsnow(i)
      coldsnow = csnow(i,n)*(tm - tsnow(i,n))
      IF ( coldsnow < 0 ) THEN
        tsnow(i,n) = tm
        dsice = -coldsnow / lf
        IF ( dsice > sice(i,n) ) dsice = sice(i,n)
        ds(i,n) = ( 1.0 - dsice/sice(i,n) ) * ds(i,n)
        sice(i,n) = sice(i,n) - dsice
        sliq(i,n) = sliq(i,n) + dsice
      END IF
    END DO
! Melt still > 0? - no snow left

!-----------------------------------------------------------------------
! Remove snow by sublimation unless snow is beneath canopy
!-----------------------------------------------------------------------
    IF ( .NOT. cansnowtile ) THEN
      dsice = MAX( ei_tile(i), 0. ) * timestep
      IF ( dsice > 0.0 ) THEN
        DO n=1,nsnow(i)
          IF ( dsice > sice(i,n) ) THEN
! Layer sublimates completely
            dsice = dsice - sice(i,n)
            sice(i,n) = 0.
            ds(i,n) = 0.
          ELSE
! Layer sublimates partially
            ds(i,n) = (1.0 - dsice/sice(i,n))*ds(i,n)
            sice(i,n) = sice(i,n) - dsice
            EXIT    !   sublimation exhausted
          END IF
        END DO
      END IF  !  DSICE>0
    END IF  !  CANSNOWTILE

!-----------------------------------------------------------------------
! Move liquid water in excess of holding capacity downwards or refreeze.
!-----------------------------------------------------------------------
    win = 0.0
    DO n=1,nsnow(i)
      sliq(i,n) = sliq(i,n) + win
      win = 0.0
      sliqmax = snowliqcap * rho_water * ds(i,n)
      IF (sliq(i,n) > sliqmax) THEN
! Liquid capacity exceeded
        win = sliq(i,n) - sliqmax
        sliq(i,n) = sliqmax
      END IF
      coldsnow = csnow(i,n)*(tm - tsnow(i,n))
      IF (coldsnow > 0) THEN
! Liquid can freeze
        dsice = MIN(sliq(i,n), coldsnow / lf)
        sliq(i,n) = sliq(i,n) - dsice
        sice(i,n) = sice(i,n) + dsice
        tsnow(i,n) = tsnow(i,n) + lf*dsice/csnow(i,n)
      END IF
    END DO

!-----------------------------------------------------------------------
! The remaining liquid water flux is melt.
! Include any separate canopy melt in this diagnostic.
!-----------------------------------------------------------------------
    melt_tile(i) = ( win / timestep ) + can_melt

!-----------------------------------------------------------------------
! Diagnose layer densities
!-----------------------------------------------------------------------
    DO n=1,nsnow(i)
      rho_snow(i,n) = 0.
      IF ( ds(i,n) > EPSILON(ds) )                                &
        rho_snow(i,n) = (sice(i,n) + sliq(i,n)) / ds(i,n)
    END DO

!-----------------------------------------------------------------------
! Add snowfall and frost as layer 0.
!-----------------------------------------------------------------------
    sice0(i) = snowfall(i)
    IF ( .NOT. cansnowtile )                                      &
      sice0(i) = snowfall(i) - MIN(ei_tile(i), 0.) * timestep
    tsnow0(i) = tsnow(i,1)
    rho0(i) = rho_snow_const

!-----------------------------------------------------------------------
! Diagnose total snow depth and mass
!-----------------------------------------------------------------------
    snowdepth(i) = sice0(i) / rho0(i)
    snowmass(i) = sice0(i)
    DO n=1,nsnow(i)
      snowdepth(i) = snowdepth(i) + ds(i,n)
      snowmass(i) = snowmass(i) + sice(i,n) + sliq(i,n)
    END DO

  END IF    !  NSNOW

END DO  !  k (points)
IF (lhook) CALL dr_hook('SNOWPACK',zhook_out,zhook_handle)
RETURN

END SUBROUTINE snowpack
#endif
