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

  !   local saved variables:                                              
  !   ----------------------                                              
  REAL :: lambda 
  REAL,ALLOCATABLE :: dz1(:),dz2(:),zc(:,:),zd(:,:) 
  !$OMP THREADPRIVATE(dz1,dz2,zc,zd,lambda)                               

  PUBLIC :: soil

CONTAINS

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

    ! this is a function definition
    fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.) 
    
    !-----------------------------------------------------------------------
    !   ground levels                                                       
    !   grnd=z/l where l is the skin depth of the diurnal cycle:            
    !   --------------------------------------------------------            
    
    WRITELOG(*,*) 'nsoil,ngrid,firstcall=',nsoil,ngrid, .TRUE.

    ALLOCATE(dz1(nsoil),dz2(nsoil)) 
    ALLOCATE(zc(ngrid,nsoil),zd(ngrid,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')
  END SUBROUTINE init_soil
  
  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 za(ngrid),zb(ngrid) 
    REAL zdz2(nsoil),z1(ngrid) 
        
    IF (firstcall) THEN 
       CALL init_soil(ngrid, nsoil)
    ELSE
       !-----------------------------------------------------------------------
       !   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
       
    ENDIF

    !-----------------------------------------------------------------------
    !   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
    
  END SUBROUTINE soil
  
END MODULE surface