module snow_drvm 4
      implicit none
!@dbparam snow_cover_coef coefficient for topography variance in
!@+       snow cover parameterisation for albedo
      real*8 :: snow_cover_coef = .15d0
!@dbparam snow_cover_same_as_rad if > 0 use the same snow fraction
!@+  for the computation of albedo and inside the snow model
      integer :: snow_cover_same_as_rad = 0


      contains

      subroutine snow_drv( 1,6
ccc input
     &     fm, evap, snsh, srht, trht, canht,
     &     drips, dripw, htdrips, htdripw,
     &     devap_dt, dsnsh_dt, dts,
     &     tp_soil, dz_soil, nlsn, top_stdev,
ccc updated
     &       dzsn, wsn, hsn, nsn,
     &       fr_snow,
ccc output
     &     flmlt, fhsng, flmlt_scale, fhsng_scale, thrmsn,
     &     flux_snow
     &     )

      use snow_model, only: snow_adv, snow_fraction, snow_redistr
      implicit none

ccc input
!@var canht heat flux from canopy, i.e. stbo*t_can**4 (w/m**2)
      real*8 fm, evap, snsh, srht, trht, canht
      real*8 drips, dripw, htdrips, htdripw
      real*8 devap_dt, dsnsh_dt, dts
      real*8 tp_soil, dz_soil, top_stdev
      integer nlsn
ccc updated
      real*8 dzsn(nlsn+1), wsn(nlsn), hsn(nlsn), fr_snow
      integer nsn
ccc output
      real*8 flmlt, fhsng, flmlt_scale, fhsng_scale, thrmsn
      real*8 flux_snow(0:nlsn)

ccc  local vars:

      real*8 epotsn, srhtsn, trhtsn, snshsn
      real*8 prsn, htprsn
ccc      real*8 xkthsn(2),cthsn(2)

      real*8 devap_sn_dt, dsnsh_sn_dt
      real*8 fr_snow_old,epot_sn_old,tr_flux(0:nlsn)


      epotsn = fm*evap
      snshsn = fm*snsh
      srhtsn = fm*srht
      trhtsn = fm*trht + (1.d0-fm)*canht! i.e. thrm_can

ccc!!! xkthsn cthsn are actually not used at the moment - should include
c!!! should pass ground properties to snow_adv
c!      xkthsn(1)=xkh(1,1)
c!      xkthsn(2)=xkh(1,2)
c!      cthsn(1)=shc(1,1)
c!      cthsn(2)=shc(1,2)

      devap_sn_dt = fm*devap_dt
      dsnsh_sn_dt = fm*dsnsh_dt

    !  do ibv=i_bare,i_vege
      fr_snow_old = fr_snow
      epot_sn_old = epotsn

      if ( snow_cover_same_as_rad .ne. 0 ) then
        call snow_fraction1(dzsn, nsn, drips, dts, top_stdev,
     &       fr_snow_old, fr_snow )
      else
        call snow_fraction(dzsn, nsn, drips, dts,
     &       fr_snow_old, fr_snow )
      endif

      !print *, sum(dzsn(1:nsn)), top_stdev, fr_snow

      if ( fr_snow <= 0.d0 ) then ! no snow
        flmlt_scale = drips
        fhsng_scale = htdrips
        flmlt = 0.d0
        fhsng = 0.d0
        thrmsn = 0.d0           ! just in case, is * 0 anyway
        flux_snow(:) = 0.d0
        if ( fr_snow_old > 0.d0 ) then ! had snow on prev. step
          flmlt_scale  = flmlt_scale
     &         + sum( wsn(1:nsn) )*fr_snow_old/dts
          fhsng_scale = fhsng_scale
     &         + sum( hsn(1:nsn) )*fr_snow_old/dts
          wsn (1:nsn) = 0.d0
          hsn (1:nsn) = 0.d0
          dzsn(1:nsn) = 0.d0
          fr_snow = 0.d0
          nsn = 1
        endif
        return  ! nothing else to do ...
      endif
         
ccc have snow - process it

      tr_flux = 0
      call snow_redistr( dzsn, wsn, hsn,
     &     nsn, fr_snow_old/fr_snow, tr_flux, dts )

ccc for tracers
      flux_snow(:) = tr_flux(:)

      ! all snow falls on fr_snow, but all liquid H2O uniformly
      prsn = drips/fr_snow + dripw
      htprsn = htdrips/fr_snow + htdripw

!!! debug : check if  snow is redistributed correctly
      if ( dzsn(1) > 0.d0 .and.
     &     dzsn(1) + prsn*dts*100.d0/15.d0 < .099d0) then
        call stop_model("set_snow: error in dz",255)
      endif

      call snow_adv(dzsn, wsn, hsn, nsn,
     &     srhtsn, trhtsn, snshsn, htprsn,
     $     epotsn,
     $     prsn, dts,
     &     tp_soil, dz_soil,
     &     flmlt, fhsng,
     &     thrmsn, dsnsh_sn_dt, devap_sn_dt, ! fbfv(ibv),
     &     tr_flux )

      flux_snow(:) = flux_snow(:) + tr_flux(:)
      flmlt_scale = 0.d0
      fhsng_scale = 0.d0

ccc make sure we don''t get negative 
      if ( flmlt < -1.d-16 ) then
        print *,'snow_drv: melt water flux < 0'
        print *,'flmlt  = ',flmlt
        print *,'pr     = ',drips, dripw
        print *,'epotsn = ',epotsn,epot_sn_old
        print *,'frac   = ',fr_snow,fr_snow_old
        print *
        !call stop_model('flmlt(ibv) < 0.d0',255)
      endif
      flmlt = max( flmlt, 0.d0 )
 
ccc update fluxes that could have been changed by snow model
      if ( fm > 0.d0 ) then
        evap = epotsn/fm
        snsh = snshsn/fm
      endif

      return
      end subroutine snow_drv



      subroutine snow_fraction1(  1,3
     &     dz, nl, prsnow, dt,
     &     top_dev, fract_snow_old, fract_snow )
!@sum computes snow cover from snow water eq. and topography
!@var fract_snow snow cover fraction (0-1)
!@var top_dev standard deviation of the surface elevation
      use constant, only : teeny
      use snow_model, only : MIN_SNOW_THICKNESS, MIN_FRACT_COVER,
     &     rho_water, rho_fresh_snow
      real*8, intent(in) :: dz(:), prsnow, dt, fract_snow_old
      integer, intent(in) :: nl
      real*8, intent(out) :: fract_snow
      real*8, intent(in) :: top_dev
! locals
!@par rho_sn average snow density (kg/m^3)
      real*8, parameter :: rho_sn = 300.d0
      real*8 dz_aver, fresh_snow

      fresh_snow = rho_water/rho_fresh_snow * prsnow*dt
      dz_aver = max( sum(dz(1:nl))*fract_snow_old + fresh_snow, 0.d0 )

      fract_snow = min( 1.d0, dz_aver/MIN_SNOW_THICKNESS )

      ! using formula from the paper by A. Roesch et al
      ! (Climate Dynamics (2001), 17: 933-946)
!!      fract_snow = fract_snow *
ccc     $     .95d0 * tanh( 100.d0 * snow_water ) *
ccc                               currently using only topography part
!!     $     sqrt ( 1000.d0 * snow_water /
!!     $     (1000.d0 * snow_water + teeny + snow_cover_coef * top_dev) )
      fract_snow = fract_snow *
     &     sqrt ( rho_sn * dz_aver /
     &     (rho_sn * dz_aver + teeny + snow_cover_coef * top_dev) )

      if ( fract_snow < MIN_FRACT_COVER ) fract_snow = 0.d0

      if ( .not. fract_snow >= 0.d0 )
     &     call stop_model("NaN in snow_fraction1, 255")
      
      end subroutine snow_fraction1



      end module snow_drvm