source: dynamico_lmdz/simple_physics/phyparam/physics/surface.F90 @ 4241

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,zc,lambda)

  ! internal state variables needed for checkpoint/restart
  REAL, ALLOCATABLE :: zc(:,:),zd(:,:)
  !$OMP THREADPRIVATE(zc,zd)

  PUBLIC :: soil, zc, zd

CONTAINS

  SUBROUTINE init_soil(ngrid,nsoil)
    INTEGER, INTENT(IN) :: ngrid, 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,ngrid,firstcall=',nsoil,ngrid, .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

  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
       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