MODULE surf_param_mod
!
  IMPLICIT NONE
!
CONTAINS
!
!-------------------------------------------------------------------------------
!
FUNCTION eff_surf_param(klon, nbtersurf, x, frac, hatype, Zref) RESULT(eff_param)
!
!-------------------------------------------------------------------------------
!
! Arguments:
  INTEGER, INTENT(IN)                          :: klon       ! grid point
  INTEGER, INTENT(IN)                          :: nbtersurf  ! number of land hetero. subsurfaces
  REAL, DIMENSION(klon, nbtersurf), INTENT(IN) :: x          ! variable or parameter to integrate
  REAL, DIMENSION(klon, nbtersurf), INTENT(IN) :: frac       ! fraction of each land hetero. subsurface
  CHARACTER(LEN=3), INTENT(IN)                 :: hatype     ! method to integrate the parameter
  REAL, OPTIONAL, DIMENSION(klon), INTENT(IN)  :: Zref       ! reference height for CDN averaging (m)
  REAL, DIMENSION(klon)                        :: eff_param  ! effective parameter
!-------------------------------------------------------------------------------
! Local variables:
  INTEGER ik, is
  REAL :: zrefd = 10.  ! default reference height for CDN averaging (m)
  REAL :: Cdref        ! reference height for CDN averaging (m)
!-------------------------------------------------------------------------------
!
  eff_param(:) = 0.
  DO ik = 1, klon
    DO is = 1, nbtersurf
      !
      ! arithmetic averaging
      IF (hatype == 'ARI') THEN
        eff_param(ik) = eff_param(ik) + frac(ik,is) * x(ik,is)
      !
      ! inverse averaging
      ELSEIF (hatype == 'INV') THEN
        IF (x(ik,is) .NE. 0.) THEN
          eff_param(ik) = eff_param(ik) + frac(ik,is) * 1./x(ik,is)
        ENDIF
      !
      ! inverse of square logarithm averaging
      ELSEIF (hatype == 'CDN') THEN
        IF (PRESENT(Zref)) THEN
          Cdref = Zref(ik)
        ELSE
          Cdref = zrefd
        ENDIF
      !
        IF (x(ik,is) .NE. 0.) THEN
          eff_param(ik) = eff_param(ik) + frac(ik,is) * 1./(LOG(Cdref/x(ik,is)))**2
        ENDIF
      !
      ELSE
        PRINT*, 'eff_surf_param: invalid averaging type: ', hatype
      ENDIF
    ENDDO
    !
    IF (hatype == 'CDN') THEN
      eff_param(ik) = Cdref * exp(-sqrt(1./eff_param(ik)))
    ENDIF
  !
  ENDDO
!
END FUNCTION eff_surf_param
!
!
!-------------------------------------------------------------------------------
!
FUNCTION average_surf_var(klon, nbtersurf, x, frac, hatype) RESULT(x_avg)
!
!-------------------------------------------------------------------------------
!
! Arguments:
  INTEGER, INTENT(IN)                          :: klon       ! grid point
  INTEGER, INTENT(IN)                          :: nbtersurf  ! number of land hetero. subsurfaces
  REAL, DIMENSION(klon, nbtersurf), INTENT(IN) :: x          ! variable or parameter to integrate
  REAL, DIMENSION(klon, nbtersurf), INTENT(IN) :: frac       ! fraction of each land hetero. subsurface
  CHARACTER(LEN=3), INTENT(IN)                 :: hatype     ! method to integrate the parameter
  REAL, DIMENSION(klon)                        :: x_avg      ! average variable
!
! Local variables:
  INTEGER ik, is
!
!-------------------------------------------------------------------------------
!
  x_avg(:) = 0.
  DO ik = 1, klon
    DO is = 1, nbtersurf
      !
      ! arithmetic averaging
      IF (hatype == 'ARI') THEN
        x_avg(ik) = x_avg(ik) + frac(ik,is) * x(ik,is)
      ELSE
        PRINT*, 'average_surf_var: invalid averaging type: ', hatype
      ENDIF
    ENDDO
  ENDDO
!
END FUNCTION average_surf_var
!
!-------------------------------------------------------------------------------
!
FUNCTION interpol_tsoil(klon, nbtersurf, nsoilmx, nbtsoildepths, alpha, period, inertie, hcond, tsoil_depth, tsurf, tsoil_) RESULT(tsoil)
!
!-------------------------------------------------------------------------------
!
! Arguments:
  INTEGER, INTENT(IN)                                         :: klon          ! grid point
  INTEGER, INTENT(IN)                                         :: nbtersurf     ! number of land hetero. subsurfaces
  INTEGER, INTENT(IN)                                         :: nsoilmx       ! number of soil layers in the model grid
  INTEGER, INTENT(IN)                                         :: nbtsoildepths ! number of soil depths for soil temperature initialization
  REAL, INTENT(IN)                                            :: alpha         ! parameter for soil discretization
  REAL, INTENT(IN)                                            :: period        ! parameter for soil discretization
  REAL, DIMENSION(klon, nbtersurf), INTENT(IN)                :: inertie       ! soil thermal inertia
  REAL, DIMENSION(klon, nbtersurf), INTENT(IN)                :: hcond         ! soil heat conductivity
  REAL, DIMENSION(klon, nbtsoildepths, nbtersurf), INTENT(IN) :: tsoil_depth   ! soil depth at which temperature is given (m)
  REAL, DIMENSION(klon, nbtersurf), INTENT(IN)                :: tsurf         ! surface temperature
  REAL, DIMENSION(klon, nbtsoildepths, nbtersurf), INTENT(IN) :: tsoil_        ! soil temperature given at tsoil_depths
  REAL, DIMENSION(klon, nsoilmx, nbtersurf)                   :: tsoil         ! soil temperature interpolated in the model grid
!
! Local variables:
  INTEGER ik, is, iq, it
  REAL pi, slope, inter
  REAL, DIMENSION(klon, nbtersurf)          :: z1, hcap  ! first layer depth and soil heat capacity
  REAL, DIMENSION(klon, nsoilmx, nbtersurf) :: z         ! depth of the middle of the soil layer in the model grid (m)
!
!-------------------------------------------------------------------------------
!
  pi = ACOS(-1.)
  !
  DO ik = 1, klon
    DO is = 1, nbtersurf
      !
      hcap(ik,is) = inertie(ik,is)*inertie(ik,is)/hcond(ik,is)
      z1(ik,is) = SQRT(period*hcond(ik,is)/(pi*hcap(ik,is)))
      !
      DO iq = 1, nsoilmx
        ! compute depth of middle soil layer (in m)
        z(ik,iq,is) = z1(ik,is) * (alpha**(iq-0.5) - 1.) / (alpha - 1.)
        ! if z is between the surface and first tsoil_depth
        IF ((z(ik,iq,is) .GT. 0.) .AND. (z(ik,iq,is) .LE. tsoil_depth(ik,1,is))) THEN
          slope = (tsoil_(ik,1,is) - tsurf(ik,is)) / (tsoil_depth(ik,1,is) - 0.)
          inter = tsurf(ik,is)
          tsoil(ik,iq,is) = slope * z(ik,iq,is) + inter
        ENDIF
        ! other levels
        DO it = 1, nbtsoildepths-1
          IF ((z(ik,iq,is) .GT. tsoil_depth(ik,it,is)) .AND. (z(ik,iq,is) .LE. tsoil_depth(ik,it+1,is))) THEN
            slope = (tsoil_(ik,it+1,is) - tsoil_(ik,it,is)) / (tsoil_depth(ik,it+1,is) - tsoil_depth(ik,it,is))
            inter = tsoil_(ik,it,is) - slope * tsoil_depth(ik,it,is)
            tsoil(ik,iq,is) = slope * z(ik,iq,is) + inter
          ENDIF
        ENDDO
        ! for layers below
        IF (z(ik,iq,is) .GT. tsoil_depth(ik,nbtsoildepths,is)) THEN
          tsoil(ik,iq,is) = tsoil_(ik,nbtsoildepths,is)
        ENDIF
      ENDDO
      !
    ENDDO
  ENDDO
!
END FUNCTION interpol_tsoil
!
!-------------------------------------------------------------------------------
!
END MODULE surf_param_mod