source: trunk/LMDZ.PLUTO/libf/phypluto/soil.F90 @ 3978

      subroutine soil(ngrid,nsoil,firstcall,lastcall,   &
               therm_i, &
               timestep,tsurf,tsoil,    &
               capcal,fluxgrd)

      use comsoil_h, only: layer, mlayer, volcapa, inertiedat
      use comcstfi_mod, only: pi
      use time_phylmdz_mod, only: daysec
      use planete_mod, only: year_day
      use geometry_mod, only: longitude, latitude ! in radians
      use callkeys_mod, only: changeti,fluxgeo,fluxgeo2,patchflux,assymflux

      implicit none

!-----------------------------------------------------------------------
!  Author: Ehouarn Millour
!
!  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.h
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!  arguments
!  ---------
!  inputs:
      integer,intent(in) :: ngrid       ! number of (horizontal) grid-points
      integer,intent(in) :: nsoil       ! number of soil layers
      logical,intent(in) :: firstcall ! identifier for initialization call
      logical,intent(in) :: lastcall
      real,intent(in) :: therm_i(ngrid,nsoil) ! thermal inertia
      real,intent(in) :: timestep           ! time step
      real,intent(in) :: tsurf(ngrid)   ! surface temperature
! outputs:
      real,intent(out) :: tsoil(ngrid,nsoil) ! soil (mid-layer) temperature
      real,intent(out) :: capcal(ngrid) ! surface specific heat
      real,intent(out) :: fluxgrd(ngrid) ! surface diffusive heat flux

! local saved variables:
      real,dimension(:,:),save,allocatable :: mthermdiff ! mid-layer thermal diffusivity
      real,dimension(:,:),save,allocatable :: thermdiff ! inter-layer thermal diffusivity
      real,dimension(:),save,allocatable :: coefq ! q_{k+1/2} coefficients
      real,dimension(:,:),save,allocatable :: coefd ! d_k coefficients
      real,dimension(:,:),save,allocatable :: alph ! alpha_k coefficients
      real,dimension(:,:),save,allocatable :: beta ! beta_k coefficients
      real,save :: mu
!$OMP THREADPRIVATE(mthermdiff,thermdiff,coefq,coefd,alph,beta,mu)

! local variables:
      integer ig,ik

! 0. Initialisations and preprocessing step
      if (firstcall) then
      ! note: firstcall is set to .true. or .false. by the caller
      !       and not changed by soil.F

        ALLOCATE(mthermdiff(ngrid,0:nsoil-1)) ! mid-layer thermal diffusivity
        ALLOCATE(thermdiff(ngrid,nsoil-1))    ! inter-layer thermal diffusivity
        ALLOCATE(coefq(0:nsoil-1))              ! q_{k+1/2} coefficients
        ALLOCATE(coefd(ngrid,nsoil-1))        ! d_k coefficients
        ALLOCATE(alph(ngrid,nsoil-1))         ! alpha_k coefficients
        ALLOCATE(beta(ngrid,nsoil-1))         ! beta_k coefficients

        ! 0.05 Build coefficients mu, q_{k+1/2}, d_k, alpha_k and capcal
        ! mu
        mu=mlayer(0)/(mlayer(1)-mlayer(0))

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

      if (firstcall.or.changeti) then
! 0.1 Build mthermdiff(:), the mid-layer thermal diffusivities
      do ig=1,ngrid
        do ik=0,nsoil-1
          mthermdiff(ig,ik)=therm_i(ig,ik+1)*therm_i(ig,ik+1)/volcapa
          !write(*,*),'soil: ik: ',ik,' mthermdiff:',mthermdiff(ig,ik)
        enddo
      enddo

! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
      do ig=1,ngrid
        do ik=1,nsoil-1
      thermdiff(ig,ik)=((layer(ik)-mlayer(ik-1))*mthermdiff(ig,ik)  &
                     +(mlayer(ik)-layer(ik))*mthermdiff(ig,ik-1))   &
                         /(mlayer(ik)-mlayer(ik-1))
!       write(*,*),'soil: ik: ',ik,' thermdiff:',thermdiff(ig,ik)
        enddo
      enddo

! 0.3 Build coefficients d_k, alpha_k and capcal
      do ig=1,ngrid
        ! d_k
        do ik=1,nsoil-1
          coefd(ig,ik)=thermdiff(ig,ik)/(mlayer(ik)-mlayer(ik-1))
        enddo

        ! alph_{N-1}
        alph(ig,nsoil-1)=coefd(ig,nsoil-1)/ &
                       (coefq(nsoil-1)+coefd(ig,nsoil-1))
        ! alph_k
        do ik=nsoil-2,1,-1
          alph(ig,ik)=coefd(ig,ik)/(coefq(ik)+coefd(ig,ik+1)*   &
                                   (1.-alph(ig,ik+1))+coefd(ig,ik))
        enddo

        ! capcal
! Cstar
        capcal(ig)=volcapa*layer(1)+    &
                   (thermdiff(ig,1)/(mlayer(1)-mlayer(0)))* &
                   (timestep*(1.-alph(ig,1)))
! Cs
        capcal(ig)=capcal(ig)/(1.+mu*(1.0-alph(ig,1))*  &
                              thermdiff(ig,1)/mthermdiff(ig,0))
      !write(*,*)'soil: ig=',ig,' capcal(ig)=',capcal(ig)
      enddo ! of do ig=1,ngrid

      endif ! of if (firstcall)

      if (.not.firstcall) then

!  1. Compute soil temperatures
! First layer:
      do ig=1,ngrid
        tsoil(ig,1)=(tsurf(ig)+mu*beta(ig,1)*   &
                              thermdiff(ig,1)/mthermdiff(ig,0))/    &
                   (1.+mu*(1.0-alph(ig,1))* &
                    thermdiff(ig,1)/mthermdiff(ig,0))
      enddo
! Other layers:
      do ik=1,nsoil-1
        do ig=1,ngrid
          tsoil(ig,ik+1)=alph(ig,ik)*tsoil(ig,ik)+beta(ig,ik)
        enddo
      enddo

! Flux at the bottom
      do ig=1,ngrid

       if (assymflux) then
         if ( (latitude(ig)*180./pi.le.90.).and. &
                (latitude(ig)*180./pi.ge.0.)  ) then
            tsoil(ig,nsoil) = tsoil(ig,nsoil)  &
            + timestep*fluxgeo2/(volcapa*(layer(nsoil)-layer(nsoil-1)))
         else   
            tsoil(ig,nsoil) = tsoil(ig,nsoil) &
            + timestep*fluxgeo/(volcapa*(layer(nsoil)-layer(nsoil-1)))
         endif

       elseif(patchflux.eq.1) then
         if ( (longitude(ig)*180./pi.le.180.).and.(longitude(ig)*180./pi.ge.174.) &
         .and.(((latitude(ig)*180./pi.le.46.).and.(latitude(ig)*180./pi.ge.42.)) &
         .or.  ((latitude(ig)*180./pi.le.36.).and.(latitude(ig)*180./pi.ge.32.)) &
         .or.  ((latitude(ig)*180./pi.le.26.).and.(latitude(ig)*180./pi.ge.22.)) &
         .or.  ((latitude(ig)*180./pi.le.16.).and.(latitude(ig)*180./pi.ge.12.)) &
         .or.  ((latitude(ig)*180./pi.le.6.).and.(latitude(ig)*180./pi.ge.2.)) &
                ) ) then
            tsoil(ig,nsoil) = tsoil(ig,nsoil) &
            + timestep*fluxgeo2/(volcapa*(layer(nsoil)-layer(nsoil-1)))
         else   
            tsoil(ig,nsoil) = tsoil(ig,nsoil) &
             + timestep*fluxgeo/(volcapa*(layer(nsoil)-layer(nsoil-1)))
         endif

       else 
         tsoil(ig,nsoil) = tsoil(ig,nsoil) &
         + timestep*fluxgeo/(volcapa*(layer(nsoil)-layer(nsoil-1)))

       endif
      enddo

      endif! of if (firstcall)

!  2. Compute beta coefficients (preprocessing for next time step)
! Bottom layer, beta_{N-1}
      do ig=1,ngrid
        beta(ig,nsoil-1)=coefq(nsoil-1)*tsoil(ig,nsoil) &
                        /(coefq(nsoil-1)+coefd(ig,nsoil-1))
      enddo
! Other layers
      do ik=nsoil-2,1,-1
        do ig=1,ngrid
          beta(ig,ik)=(coefq(ik)*tsoil(ig,ik+1)+    &
                      coefd(ig,ik+1)*beta(ig,ik+1))/    &
                      (coefq(ik)+coefd(ig,ik+1)*(1.0-alph(ig,ik+1)) &
                       +coefd(ig,ik))
        enddo
      enddo


!  3. Compute surface diffusive flux & calorific capacity
      do ig=1,ngrid
! Cstar
!        capcal(ig)=volcapa(ig,1)*layer(ig,1)+
!     &              (thermdiff(ig,1)/(mlayer(ig,1)-mlayer(ig,0)))*
!     &              (timestep*(1.-alph(ig,1)))
! Fstar

!         print*,'this far in soil 1'
!         print*,'thermdiff=',thermdiff(ig,1)
!         print*,'mlayer=',mlayer
!         print*,'beta=',beta(ig,1)
!         print*,'alph=',alph(ig,1)
!         print*,'tsoil=',tsoil(ig,1)

        fluxgrd(ig)=(thermdiff(ig,1)/(mlayer(1)-mlayer(0)))*    &
                   (beta(ig,1)+(alph(ig,1)-1.0)*tsoil(ig,1))

!        mu=mlayer(ig,0)/(mlayer(ig,1)-mlayer(ig,0))
!        capcal(ig)=capcal(ig)/(1.+mu*(1.0-alph(ig,1))*
!     &                         thermdiff(ig,1)/mthermdiff(ig,0))
! Fs
        fluxgrd(ig)=fluxgrd(ig)+(capcal(ig)/timestep)*  &
                   (tsoil(ig,1)*(1.+mu*(1.0-alph(ig,1))*    &
                              thermdiff(ig,1)/mthermdiff(ig,0)) &
                    -tsurf(ig)-mu*beta(ig,1)*   &
                               thermdiff(ig,1)/mthermdiff(ig,0))
      enddo

      end