source: trunk/LMDZ.COMMON/libf/evolution/ini_soil_mod.F90 @ 2849

      MODULE ini_soil_mod

      IMPLICIT NONE
      
      CONTAINS    

     subroutine ini_icetable(timelen,ngrid,nsoil_PEM, &
               therm_i, timestep,tsurf_ave,tsoil_ave,tsurf_inst, tsoil_inst,q_co2,q_h2o,ps,ice_table)

#ifndef CPP_STD

      use vertical_layers_mod, only: ap,bp
      use comsoil_h_PEM, only:  fluxgeo,layer_PEM,inertiedat_PEM
      use comsoil_h,only: volcapa, nsoilmx
#endif

   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"

!-----------------------------------------------------------------------
!  arguments
!  ---------
!  inputs:
      integer,intent(in) :: timelen     ! Time length in for time-series data
      integer,intent(in) :: ngrid       ! number of (horizontal) grid-points 
      integer,intent(in) :: nsoil_PEM   ! number of soil layers  in the PEM

      real,intent(in) :: timestep           ! time step
      real,intent(in) :: tsurf_ave(ngrid)   ! surface temperature

      real,intent(in) :: q_co2(ngrid,timelen)                     ! MMR tracer co2 [kg/kg]
      real,intent(in) :: q_h2o(ngrid,timelen)                     ! MMR tracer h2o [kg/kg]
      real,intent(in) :: ps(ngrid,timelen)                        ! surface pressure [Pa]
      real,intent(in) :: tsurf_inst(ngrid,timelen) ! soil (mid-layer) temperature

! outputs:
      real,intent(inout) :: tsoil_ave(ngrid,nsoil_PEM) ! soil (mid-layer) temperature. 
      real,intent(inout) :: tsoil_inst(ngrid,nsoil_PEM,timelen) ! soil (mid-layer) temperature
      real,intent(out) :: ice_table(ngrid)                 ! ice table [m]
      real,intent(inout) :: therm_i(ngrid,nsoil_PEM) ! thermal inertia

! local variables:
      integer ig,isoil,it,k,iref
      REAL :: error_depth = 2.
      REAL :: tsoil_saved(nsoil_PEM)
      integer :: countmax = 20
      integer :: countloop 
      REAL :: tol_error = 0.1
      REAL :: ice_inertia = 2120.
      REAL :: alph_PEM(nsoil_PEM-1)
      REAL :: beta_PEM(nsoil_PEM-1)
      real :: rhoc(nsoil_PEM)
      real :: volcapa_ice = 1.43e7

      real :: k_soil(nsoil_PEM)
      real :: d1(nsoil_PEM),d2(nsoil_PEM),d3(nsoil_PEM),re(nsoil_PEM)
      real :: Tcol_saved(nsoil_PEM)
      real :: tsoil_prev(nsoil_PEM)
      real :: tsurf_prev
      real :: icedepth_prev

    real :: m_h2o = 18.01528E-3
    real :: m_co2 = 44.01E-3  
    real :: m_noco2 = 33.37E-3  
    real :: A,B,z1,z2
    real :: alpha = -6143.7 
    real :: beta = 29.9074
   
    real,allocatable :: mass_mean(:)                            ! mean mass above the surface
    real,allocatable :: zplev(:)                           ! pressure above the surface
    real,allocatable :: pvapor(:)                               ! partial pressure above the surface
   
    real,allocatable :: rhovapor(:)
    real :: rhovapor_avg                          ! mean vapor_density at the surface yearly averaged

    real :: psv_surf
    real,allocatable :: rho_soil(:)            ! water vapor in the soil
    real,allocatable :: rho_soil_avg(:)                ! water vapor in the soil yearly averaged

    real,allocatable :: diff_rho(:)                    ! difference of vapor content

#ifndef CPP_STD

    A =(1/m_co2 - 1/m_noco2)
    B=1/m_noco2

    ice_table(:) = 1.
    do ig = 1,ngrid
      error_depth = 1.
      countloop = 0
    
      do while(( error_depth.gt.tol_error).and.(countloop.lt.countmax).and.(ice_table(ig).gt.-1e-20))
    
        countloop = countloop +1
        Tcol_saved(:) = tsoil_ave(ig,:)

!2.  Compute ice table

! 2.1 Compute  water density at the surface, yearly averaged
        allocate(mass_mean(timelen))
!   1.1 Compute the partial pressure of vapor
! a. the molecular mass into the column
       mass_mean(:) = 1/(A*q_co2(ig,:) +B)
! b. pressure level
        allocate(zplev(timelen))
        do it = 1,timelen
         zplev(it) = ap(1) + bp(1)*ps(ig,it)
        enddo

! c. Vapor pressure
        allocate(pvapor(timelen))
        pvapor(:) = mass_mean(:)/m_h2o*q_h2o(ig,:)*zplev(:)
        deallocate(zplev) 
        deallocate(mass_mean)
  
!  d!  Check if there is frost at the surface and then compute the density
!    at the surface
        allocate(rhovapor(timelen))

         do it = 1,timelen
           psv_surf = exp(alpha/tsurf_inst(ig,it) +beta)
           rhovapor(it) = min(psv_surf,pvapor(it))/tsurf_inst(ig,it)
         enddo
        deallocate(pvapor)
        rhovapor_avg =SUM(rhovapor(:),1)/timelen                     
        deallocate(rhovapor)

! 2.2 Compute  water density at the soil layer, yearly averaged

        allocate(rho_soil(timelen))
        allocate(rho_soil_avg(nsoil_PEM))

        do isoil = 1,nsoil_PEM
          do it = 1,timelen 
            rho_soil(it) = exp(alpha/tsoil_inst(ig,isoil,it) +beta)/tsoil_inst(ig,isoil,it)        
          enddo
         rho_soil_avg(isoil) = SUM(rho_soil(:),1)/timelen     
        enddo
        deallocate(rho_soil)
         
!2.3 Final: compute ice table
        icedepth_prev = ice_table(ig)
        ice_table(ig) = -1
        allocate(diff_rho(nsoil_PEM))
        do isoil = 1,nsoil_PEM
          diff_rho(isoil) = rhovapor_avg -  rho_soil_avg(isoil)
        enddo
        deallocate(rho_soil_avg)

        if(diff_rho(1) > 0) then
          ice_table(ig) = 0.
        else
          do isoil = 1,nsoil_PEM -1
            if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then
              z1 = (diff_rho(isoil) - diff_rho(isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1))
              z2 = -layer_PEM(isoil+1)*z1 +  diff_rho(isoil+1)
              ice_table(ig) = -z2/z1
              exit
            endif
          enddo
        endif
    
        deallocate(diff_rho)

!3. Update Soil Thermal Inertia

        if (ice_table(ig).gt. 0.) then
          if (ice_table(ig).lt. 1e-10) then
            do isoil = 1,nsoil_PEM
              therm_i(ig,isoil)=ice_inertia
            enddo
          else 
      ! 4.1 find the index of the mixed layer
          iref=0 ! initialize iref
          do k=1,nsoil_PEM ! loop on layers
            if (ice_table(ig).ge.layer_PEM(k)) then
              iref=k ! pure regolith layer up to here
            else
              ! correct iref was obtained in previous cycle
            exit
            endif
          enddo

      ! 4.2 Build the new ti
         do isoil=1,iref
           therm_i(ig,isoil) =inertiedat_PEM(ig,isoil) 
         enddo

         if (iref.lt.nsoil_PEM) then
           if (iref.ne.0) then
             ! mixed layer
             therm_i(ig,iref+1)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ &
              (((ice_table(ig)-layer_PEM(iref))/(inertiedat_PEM(ig,iref)**2))+ &
                      ((layer_PEM(iref+1)-ice_table(ig))/(ice_inertia**2))))
            
           else ! first layer is already a mixed layer
            ! (ie: take layer(iref=0)=0)
            therm_i(ig,1)=sqrt((layer_PEM(1))/ &
                          (((ice_table(ig))/(inertiedat_PEM(ig,1)**2))+ &
                           ((layer_PEM(1)-ice_table(ig))/(ice_inertia**2))))
           endif ! of if (iref.ne.0)        
           ! lower layers of pure ice
           do isoil=iref+2,nsoil_PEM
             therm_i(ig,isoil)=ice_inertia
           enddo
         endif ! of if (iref.lt.(nsoilmx))
       endif ! permanent glaciers

       call soil_pem_1D(nsoil_PEM,.true.,therm_i(ig,:),timestep,tsurf_ave(ig),tsoil_ave(ig,:),alph_PEM,beta_PEM) 
       call soil_pem_1D(nsoil_PEM,.false.,therm_i(ig,:),timestep,tsurf_ave(ig),tsoil_ave(ig,:),alph_PEM,beta_PEM) 

       do it = 1,timelen 
         tsoil_inst(ig,:,it) = tsoil_inst(ig,:,it) - (Tcol_saved(:) - tsoil_ave(ig,:))
       enddo

       error_depth = abs(icedepth_prev - ice_table(ig))      
       endif

     enddo

     error_depth = 1.
     countloop = 0
   enddo

#endif
 
      END SUBROUTINE ini_icetable
      subroutine soil_pem_1D(nsoil,firstcall, &
               therm_i,                          &
               timestep,tsurf,tsoil,alph,beta)
#ifndef CPP_STD
      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,  &
                           mu_PEM,fluxgeo
      use comsoil_h,only: volcapa
#endif
      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) :: nsoil       ! number of soil layers  
      logical,intent(in) :: firstcall ! identifier for initialization call 
      real,intent(in) :: therm_i(nsoil) ! thermal inertia
      real,intent(in) :: timestep           ! time step
      real,intent(in) :: tsurf   ! surface temperature
 
! outputs:
      real,intent(inout) :: tsoil(nsoil) ! soil (mid-layer) temperature
      real,intent(inout) :: alph(nsoil-1)
      real,intent(inout) :: beta(nsoil-1)
      
! local variables:
      integer ik
      real :: thermdiff_PEM(nsoil-1)
      real :: mthermdiff_PEM(0:nsoil-1)
      real :: coefd_PEM(nsoil-1)
      real :: coefq_PEM(0:nsoil-1)

#ifndef CPP_STD

! 0. Initialisations and preprocessing step
 if (firstcall) then
    
! 0.1 Build mthermdiff_PEM(:), the mid-layer thermal diffusivities
        do ik=0,nsoil-1
          mthermdiff_PEM(ik)=therm_i(ik+1)*therm_i(ik+1)/volcapa
    
        enddo

! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
        do ik=1,nsoil-1
      thermdiff_PEM(ik)=((layer_PEM(ik)-mlayer_PEM(ik-1))*mthermdiff_PEM(ik) &
                     +(mlayer_PEM(ik)-layer_PEM(ik))*mthermdiff_PEM(ik-1))  &
                         /(mlayer_PEM(ik)-mlayer_PEM(ik-1))
        enddo

! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alpha_k and capcal
      ! mu_PEM
      mu_PEM=mlayer_PEM(0)/(mlayer_PEM(1)-mlayer_PEM(0))

      ! q_{1/2}
      coefq_PEM(0)=volcapa*layer_PEM(1)/timestep
        ! q_{k+1/2}
        do ik=1,nsoil-1
          coefq_PEM(ik)=volcapa*(layer_PEM(ik+1)-layer_PEM(ik))                  &
                      /timestep
        enddo

        ! d_k
        do ik=1,nsoil-1
          coefd_PEM(ik)=thermdiff_PEM(ik)/(mlayer_PEM(ik)-mlayer_PEM(ik-1))
        enddo
        
        ! alph_PEM_{N-1}
        alph(nsoil-1)=coefd_PEM(nsoil-1)/                        &
                       (coefq_PEM(nsoil-1)+coefd_PEM(nsoil-1))
        ! alph_PEM_k
        do ik=nsoil-2,1,-1
          alph(ik)=coefd_PEM(ik)/(coefq_PEM(ik)+coefd_PEM(ik+1)*    &
                                   (1.-alph(ik+1))+coefd_PEM(ik))
        enddo

endif ! of if (firstcall)

      IF (.not.firstcall) THEN
!  2. Compute soil temperatures
! First layer:
        tsoil(1)=(tsurf+mu_PEM*beta(1)* &  
                                  thermdiff_PEM(1)/mthermdiff_PEM(0))/ &
                   (1.+mu_PEM*(1.0-alph(1))*&
                    thermdiff_PEM(1)/mthermdiff_PEM(0))

! Other layers:
      do ik=1,nsoil-1
          tsoil(ik+1)=alph(ik)*tsoil(ik)+beta(ik)
       enddo
      
          ENDIF

!  2. Compute beta_PEM coefficients (preprocessing for next time step)
! Bottom layer, beta_PEM_{N-1}

        beta(nsoil-1)=coefq_PEM(nsoil-1)*tsoil(nsoil)            &
                        /(coefq_PEM(nsoil-1)+coefd_PEM(nsoil-1)) & 
                 +  fluxgeo/(coefq_PEM(nsoil-1)+coefd_PEM(nsoil-1))
! Other layers
      do ik=nsoil-2,1,-1
          beta(ik)=(coefq_PEM(ik)*tsoil(ik+1)+           &
                      coefd_PEM(ik+1)*beta(ik+1))/               &
                      (coefq_PEM(ik)+coefd_PEM(ik+1)*(1.0-alph(ik+1)) &
                       +coefd_PEM(ik))
      enddo
#endif
      end

      END MODULE ini_soil_mod