MODULE surface

#include "use_logging.h"

  IMPLICIT NONE
  PRIVATE
  SAVE

  REAL, PARAMETER :: pi=2.*ASIN(1.)

  ! common variables
  REAL, PUBLIC ::  I_mer,I_ter,Cd_mer,Cd_ter, &
       &           alb_mer,alb_ter,emi_mer,emi_ter

  ! precomputed variables
  REAL :: lambda
  REAL, ALLOCATABLE :: dz1(:),dz2(:)
  !$OMP THREADPRIVATE(dz1,dz2)
  REAL, ALLOCATABLE :: rnatur(:), albedo(:),emissiv(:), z0(:), inertie(:)
  !$OMP THREADPRIVATE( rnatur, albedo, emissiv, z0, inertie)

  ! internal state, written to / read from disk at checkpoint / restart
  REAL, ALLOCATABLE :: tsurf(:), tsoil(:,:)
  !$OMP THREADPRIVATE(tsurf, tsoil)
  ! variables below should be temporary arrays, not persistent
  REAL, ALLOCATABLE :: zc(:,:),zd(:,:), capcal(:), fluxgrd(:)
  !$OMP THREADPRIVATE(zc,zd, capcal, fluxgrd)

  PUBLIC :: init_soil, &
       soil, soil_new, soil_forward, soil_backward, &
       zc, zd, &
       rnatur, albedo, emissiv, z0, inertie, &
       tsurf, tsoil, capcal, fluxgrd

CONTAINS

  SUBROUTINE init_soil(nsoil)
    INTEGER, INTENT(IN) :: nsoil
    REAL :: min_period,dalph_soil, rk,fz1,rk1,rk2
    INTEGER :: jk

    !-----------------------------------------------------------------------
    !   ground levels
    !   grnd=z/l where l is the skin depth of the diurnal cycle:
    !   --------------------------------------------------------

    WRITELOG(*,*) 'nsoil,firstcall=',nsoil, .TRUE.

    ALLOCATE(dz1(nsoil),dz2(nsoil))

    min_period=20000.
    dalph_soil=2.

    !   la premiere couche represente un dixieme de cycle diurne
    fz1=sqrt(min_period/pi)

    DO jk=1,nsoil
       rk1=jk
       rk2=jk-1
       dz2(jk)=fz(rk1)-fz(rk2)
    ENDDO
    DO jk=1,nsoil-1
       rk1=jk+.5
       rk2=jk-.5
       dz1(jk)=1./(fz(rk1)-fz(rk2))
    ENDDO
    lambda=fz(.5)*dz1(1)

    WRITELOG(*,*) 'full layers, intermediate layers (secoonds)'
    DO jk=1,nsoil
       rk=jk
       rk1=jk+.5
       rk2=jk-.5
       WRITELOG(*,*) fz(rk1)*fz(rk2)*pi,        &
            &        fz(rk)*fz(rk)*pi
    ENDDO
    LOG_INFO('init_soil')

  CONTAINS

    FUNCTION fz(rk) RESULT(val)
      REAL :: val, rk
      val = fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
    END FUNCTION fz

  END SUBROUTINE init_soil

  PURE SUBROUTINE soil_backward(ngrid,nsoil, zc,zd, ptsrf,ptsoil)
    INTEGER, INTENT(IN) :: ngrid, nsoil         ! number of columns, of soil layers
    REAL, INTENT(IN)    :: zc(ngrid, nsoil), zd(ngrid, nsoil) ! LU factorization
    REAL, INTENT(IN)    :: ptsrf(ngrid)         ! new surface temperature
    REAL, INTENT(INOUT) :: ptsoil(ngrid,nsoil)  ! soil temperature
    INTEGER :: ig, jk

    !-----------------------------------------------------------------------
    !  Computation of the soil temperatures using the Cgrd and Dgrd
    !  coefficient computed during the forward sweep
    !  -----------------------------------------------
    
    !  surface temperature => temperature in first soil layer
    DO ig=1,ngrid
       ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/                   &
            &      (lambda*(1.-zd(ig,1))+1.)
    ENDDO
    
    !   other temperatures
    DO jk=1,nsoil-1
       DO ig=1,ngrid
          ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk)
       ENDDO
    ENDDO
  END SUBROUTINE Soil_backward

  PURE SUBROUTINE soil_forward(ngrid, nsoil, ptimestep, ptherm_i, ptsrf, ptsoil, &
       &                       zc, zd, pcapcal, pfluxgrd)
    INTEGER, INTENT(IN) :: ngrid, nsoil         ! number of columns, of soil layers
    REAL, INTENT(IN)    :: ptimestep,         & ! time step
         &                 ptherm_i(ngrid),   & ! thermal inertia ??
         &                 ptsrf(ngrid),      & ! surface temperature before heat conduction
         &                 ptsoil(ngrid, nsoil) ! soil temperature before heat conduction
    REAL, INTENT(OUT)   :: zc(ngrid,nsoil),   &
         &                 zd(ngrid, nsoil),  & ! LU factorization for backward sweep
         &                 pcapcal(ngrid),    & ! effective calorific capacity
         &                 pfluxgrd(ngrid)      ! conductive heat flux at the ground
    REAL :: z1, zdz2(ngrid)
    INTEGER :: jk, ig

    !-----------------------------------------------------------------------
    !   Computation of the Cgrd and Dgrd coefficients the backward sweep :
    !   ---------------------------------------------------------------
    
    DO jk=1,nsoil
       zdz2(jk)=dz2(jk)/ptimestep
    ENDDO
    
    DO ig=1,ngrid
       z1=zdz2(nsoil)+dz1(nsoil-1)
       zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1
       zd(ig,nsoil-1)=dz1(nsoil-1)/z1
    ENDDO
    
    DO jk=nsoil-1,2,-1
       DO ig=1,ngrid
          z1=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk)))
          zc(ig,jk-1)=                                                &
               &      (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1
          zd(ig,jk-1)=dz1(jk-1)*z1
       ENDDO
    ENDDO
    
    !-----------------------------------------------------------------------
    !   computation of the surface diffusive flux from ground and
    !   calorific capacity of the ground:
    !   ---------------------------------
    
    DO ig=1,ngrid
       pfluxgrd(ig)=ptherm_i(ig)*dz1(1)*                              &
            &   (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1))
       z1=lambda*(1.-zd(ig,1))+1.
       pcapcal(ig)=ptherm_i(ig)*                                      &
            &   ptimestep*(zdz2(1)+(1.-zd(ig,1))*dz1(1))/z1
       pfluxgrd(ig)=pfluxgrd(ig)                                      &
            &   +pcapcal(ig)*(ptsoil(ig,1)*z1-lambda*zc(ig,1)-ptsrf(ig))   &
            &   /ptimestep
    ENDDO
  END SUBROUTINE soil_forward

  SUBROUTINE soil_new(ngrid,nsoil,ptimestep,ptherm_i, ptsrf,ptsoil, pcapcal,pfluxgrd)
    INTEGER, INTENT(IN) :: ngrid, nsoil         ! number of columns, of soil layers
    REAL, INTENT(IN)    :: ptimestep,         & ! time step
         &                 ptherm_i(ngrid)      ! thermal inertia ??
    REAL, INTENT(INOUT) :: ptsrf(ngrid),      & ! surface temperature
         &                 ptsoil(ngrid,nsoil)  ! soil temperature
    REAL, INTENT(OUT)   :: pcapcal(ngrid),    & ! effective calorific capacity
         &                 pfluxgrd(ngrid)      ! conductive heat flux at the ground
    CALL soil_backward(ngrid,nsoil, zc,zd, ptsrf,ptsoil)
    CALL soil_forward(ngrid, nsoil, ptimestep, ptherm_i, ptsrf, ptsoil, &
            &            zc, zd, pcapcal, pfluxgrd)
  END SUBROUTINE soil_new

  SUBROUTINE soil(ngrid,nsoil,firstcall,ptherm_i,          &
       &          ptimestep,ptsrf,ptsoil,                  &
       &          pcapcal,pfluxgrd)

    !=======================================================================
    !
    !   Auteur:  Frederic Hourdin     30/01/92
    !   -------
    !
    !   objet:  computation of : the soil temperature evolution
    !   ------                   the surfacic heat capacity "Capcal"
    !                            the surface conduction flux pcapcal
    !
    !
    !   Method: implicit time integration
    !   -------
    !   Consecutive ground temperatures are related by:
    !           T(k+1) = C(k) + D(k)*T(k)  (1)
    !   the coefficients C and D are computed at the t-dt time-step.
    !   Routine structure:
    !   1)new temperatures are computed  using (1)
    !   2)C and D coefficients are computed from the new temperature
    !     profile for the t+dt time-step
    !   3)the coefficients A and B are computed where the diffusive
    !     fluxes at the t+dt time-step is given by
    !            Fdiff = A + B Ts(t+dt)
    !     or     Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
    !            with F0 = A + B (Ts(t))
    !                 Capcal = B*dt
    !
    !   Interface:
    !   ----------
    !
    !   Arguments:
    !   ----------
    !   ngrid               number of grid-points
    !   ptimestep              physical timestep (s)
    !   pto(ngrid,nsoil)     temperature at time-step t (K)
    !   ptn(ngrid,nsoil)     temperature at time step t+dt (K)
    !   pcapcal(ngrid)      specific heat (W*m-2*s*K-1)
    !   pfluxgrd(ngrid)      surface diffusive flux from ground (Wm-2)
    !
    !=======================================================================
    !   declarations:
    !   -------------


    !-----------------------------------------------------------------------
    !  arguments
    !  ---------

    INTEGER ngrid,nsoil
    REAL ptimestep
    REAL ptsrf(ngrid),ptsoil(ngrid,nsoil),ptherm_i(ngrid)
    REAL pcapcal(ngrid),pfluxgrd(ngrid)
    LOGICAL firstcall


    !-----------------------------------------------------------------------
    !  local arrays
    !  ------------

    INTEGER ig,jk
    REAL zdz2(nsoil),z1(ngrid)

    IF (firstcall) THEN
       ! init_soil is now called by iniphyparam
       ! CALL init_soil(ngrid, nsoil)
    ELSE
       IF(.FALSE.) THEN
          !-----------------------------------------------------------------------
          !   Computation of the soil temperatures using the Cgrd and Dgrd
          !  coefficient computed at the previous time-step:
          !  -----------------------------------------------
          
          !    surface temperature
          DO ig=1,ngrid
             ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/                   &
                  &      (lambda*(1.-zd(ig,1))+1.)
          ENDDO
          
          !   other temperatures
          DO jk=1,nsoil-1
             DO ig=1,ngrid
                ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk)
             ENDDO
          ENDDO
       ELSE
          CALL soil_backward(ngrid,nsoil, zc,zd, ptsrf,ptsoil)
       END IF
       
    ENDIF

    IF(.FALSE.) THEN
       !-----------------------------------------------------------------------
       !   Computation of the Cgrd and Dgrd coefficient for the next step:
       !   ---------------------------------------------------------------
       
       DO jk=1,nsoil
          zdz2(jk)=dz2(jk)/ptimestep
       ENDDO
       
       DO ig=1,ngrid
          z1(ig)=zdz2(nsoil)+dz1(nsoil-1)
          zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1(ig)
          zd(ig,nsoil-1)=dz1(nsoil-1)/z1(ig)
       ENDDO
       
       DO jk=nsoil-1,2,-1
          DO ig=1,ngrid
             z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk)))
             zc(ig,jk-1)=                                                &
                  &      (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1(ig)
             zd(ig,jk-1)=dz1(jk-1)*z1(ig)
          ENDDO
       ENDDO
       
       !-----------------------------------------------------------------------
       !   computation of the surface diffusive flux from ground and
       !   calorific capacity of the ground:
       !   ---------------------------------
       
       DO ig=1,ngrid
          pfluxgrd(ig)=ptherm_i(ig)*dz1(1)*                              &
               &   (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1))
          z1(ig)=lambda*(1.-zd(ig,1))+1.
          pcapcal(ig)=ptherm_i(ig)*                                      &
               &   ptimestep*(zdz2(1)+(1.-zd(ig,1))*dz1(1))/z1(ig)
          pfluxgrd(ig)=pfluxgrd(ig)                                      &
               &   +pcapcal(ig)*(ptsoil(ig,1)*z1(ig)-lambda*zc(ig,1)-ptsrf(ig))   &
               &   /ptimestep
       ENDDO
    ELSE
       CALL soil_forward(ngrid, nsoil, ptimestep, ptherm_i, ptsrf, ptsoil, &
            &            zc, zd, pcapcal, pfluxgrd)
    END IF
  END SUBROUTINE soil
     
   END MODULE surface