MODULE soil_evolution_mod


      IMPLICIT NONE
      
      CONTAINS    


     subroutine soil_pem_CN(ngrid,nsoil, &
               therm_i, timestep,tsurf,tsurf_prev,tsoil)


      use comsoil_h_PEM, only:   & fluxgeo,layer_PEM
      use comsoil_h,only: volcapa
      implicit none

!-----------------------------------------------------------------------
!  Author: LL
!  Purpose: Compute soil temperature using an implict 1st order scheme
!  
!  Note: depths of layers and mid-layers, soil thermal inertia and 
!        heat capacity are commons in comsoil_PEM.h
!        A convergence loop is added until equilibrium
!-----------------------------------------------------------------------

#include "dimensions.h"
!#include "dimphys.h"

!#include"comsoil.h"


!-----------------------------------------------------------------------
!  arguments
!  ---------
!  inputs:
      integer,intent(in) :: ngrid	! number of (horizontal) grid-points 
      integer,intent(in) :: nsoil	! number of soil layers  
      real,intent(in) :: therm_i(ngrid,nsoil) ! thermal inertia
      real,intent(in) :: timestep	    ! time step
      real,intent(in) :: tsurf(ngrid)   ! surface temperature
      real,intent(in) :: tsurf_prev(ngrid)   ! surface temperature at previous timestep

! outputs:
      real,intent(inout) :: tsoil(ngrid,nsoil) ! soil (mid-layer) temperature


      
! local variables:
      integer ig,ik,isoil
      REAL :: ice_inertia = 2120.
      real :: alph_PEM(nsoil)
      real :: beta_PEM(nsoil)
      real :: rhoc(nsoil)
      real :: volcapa_ice = 1.43e7

      real :: k_soil(nsoil)
      real :: d1(nsoil),d2(nsoil),d3(nsoil),re(nsoil)
      real :: Tcol(nsoil)


do ig = 1,ngrid

   Tcol(:) = tsoil(ig,:)
!0. Build thermal properties of the ground


  do isoil = 1,nsoil
    if(abs(ice_inertia-therm_i(ig,isoil)).lt.1e-6) then
      rhoc(isoil) = volcapa_ice
    else
      rhoc(isoil) = volcapa
    endif
    k_soil(isoil) = therm_i(ig,isoil)*therm_i(ig,isoil)/rhoc(isoil)
  enddo

!1. Build Crank-Nicholson scheme for heat equation, following Schorgofer derivations
!1.a all the points which are not impacted by boundary condition
  do isoil = 2,nsoil-1
    
    alph_pem(isoil) = timestep/((rhoc(isoil) + rhoc(isoil+1))/2.)*k_soil(isoil+1)/((layer_PEM(isoil+1) - layer_PEM(isoil))*(layer_PEM(isoil+1) - layer_PEM(isoil-1)))
    beta_pem(isoil) = timestep/((rhoc(isoil) + rhoc(isoil+1))/2.)*k_soil(isoil)/((layer_PEM(isoil) - layer_PEM(isoil-1))*(layer_PEM(isoil+1) - layer_PEM(isoil-1)))
  enddo

!1.b At the surface: T(z=0) = Tsurf.

alph_pem(1) = k_soil(2)*timestep/(rhoc(1)*layer_PEM(2)*(layer_PEM(2)-layer_PEM(1)))
beta_pem(1) = k_soil(1)*timestep/(rhoc(1)*layer_PEM(2)*layer_PEM(1))



!1.c At the bottom: dT/dz = -Fgeo/k.
beta_pem(nsoil) = timestep*k_soil(nsoil)/(2.*rhoc(nsoil)*(layer_PEM(isoil) - layer_PEM(isoil-1))**2)

!1 .d set the tridiagonal with diagonals d1,d2,d3, and right element re: 


  d1(:) = -beta_pem(:)   
  d2(:) = 1. + alph_pem(:) + beta_pem(:)
  d3(:) = -alph_pem(:)   
  d2(nsoil) = 1. + beta_pem(nsoil)
 
  re(1)= alph_pem(1)*Tcol(2) + (1.-alph_pem(1)-beta_pem(1))*Tcol(1) + beta_pem(1)*(tsurf(ig)+tsurf_prev(ig))
  do isoil=2,nsoil-1
     re(isoil) = beta_pem(isoil)*Tcol(isoil-1) + (1.-alph_pem(isoil)-beta_pem(isoil))*Tcol(isoil) + alph_pem(isoil)*Tcol(isoil+1)
  enddo
  re(nsoil) = beta_pem(nsoil)*Tcol(nsoil-1) + (1.-beta_pem(nsoil))*Tcol(nsoil) + &
  timestep/rhoc(nsoil)*fluxgeo/(layer_PEM(nsoil)-layer_PEM(nsoil-1)) 



!2. Update by tridiagonal inversion 
  call tridag(d1,d2,d3,re,Tcol,nsoil) 
  


      tsoil(ig,:) = Tcol(:)
  enddo


      END SUBROUTINE soil_pem_CN

!====================================

      SUBROUTINE TRIDAG(A,B,C,R,U,N)
! From Numerical Recipes in Fortran 77
      INTEGER, INTENT(IN)::  N
      INTEGER ,PARAMETER:: NMAX = 1000

      REAL, INTENT(IN) :: A(N),B(N),C(N),R(N)
      REAL,INTENT(INOUT) :: U(N)

      INTEGER j
      REAL bet,gam(NMAX)

      IF(B(1).EQ.0.)PAUSE
      BET=B(1)
      U(1)=R(1)/BET
      DO 11 J=2,N
        GAM(J)=C(J-1)/BET
        BET=B(J)-A(J)*GAM(J)
        IF(BET.EQ.0.)PAUSE
        U(J)=(R(J)-A(J)*U(J-1))/BET
11    CONTINUE
      DO 12 J=N-1,1,-1
        U(J)=U(J)-GAM(J+1)*U(J+1)
12    CONTINUE
      RETURN
      END  SUBROUTINE TRIDAG


      END MODULE soil_evolution_mod