thermal_layers_conf.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     ######################################################################
      SUBROUTINE thermal_layers_conf(HTYPE,PD,PD_OUT,PHC,PHC_OUT,PTC,PTC_OUT)
!     ######################################################################
!
!!****  *THERMAL_LAYERS_CONF* 
!!
!!    PURPOSE
!!    -------
!     Adjust the thermal characteristics of the layers in road, wall, roof or
!     floor depending on the number of layers that the user wants to use during
!     the simulations.
!     Initial data are prescribed depending on user preference.
!     They have to be averaged on the layers use in the simulation
!  
!!
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!      V. Masson   *Meteo France*      
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    05/2012
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
USE modd_surf_par, ONLY : xundef
USE modi_abor1_sfx
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*       0.1   Declarations of arguments
!              -------------------------
!
 CHARACTER(LEN=5),     INTENT(IN)  :: HTYPE     ! type of surface
REAL, DIMENSION(:,:), INTENT(IN)  :: PD        ! input Layer Thickness
REAL, DIMENSION(:,:), INTENT(OUT) :: PD_OUT    ! output Layer Thickness
REAL, DIMENSION(:,:), INTENT(IN), OPTIONAL  :: PHC       ! input Heat Capacity
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PHC_OUT   ! output Heat Capacity
REAL, DIMENSION(:,:), INTENT(IN), OPTIONAL  :: PTC       ! input Thermal conductivity
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PTC_OUT   ! output Thermal conductivity
!
!*       0.2   Declarations of local variables
!
REAL :: ZD_TOT    ! Total depth
REAL :: ZD_HALF   ! Depth of the half of the total surface
!                                           ! (excluding central layer in case
!                                           ! of odd number of layers)
REAL  :: ZD_MID    ! Thickness of the layer in the middle
!                                           ! in case of odd number of layers
REAL, DIMENSION(0:SIZE(PD    ,2))::ZD_IN  ! Depth from the surface 
!                                          ! to the layer bottom
REAL, DIMENSION(0:SIZE(PD_OUT,2))::ZD_OUT ! Depth from the surface 
!                                                      ! to the layer bottom
REAL, DIMENSION(SIZE(PD,2))     :: ZW, ZHC     ! 1/TC
REAL, DIMENSION(SIZE(PD_OUT,2)) :: ZW_OUT, ZHC_OUT ! 1/TC
INTEGER                        :: IIN       ! Number of layer in input data
INTEGER                        :: IOUT      ! Number of layer in output fields
INTEGER                        :: JIN, JI       ! Loop counter on input layers
INTEGER                        :: JOUT      ! Loop counter on output layers
!
REAL, PARAMETER                :: ZD_G1 = 0.001  ! uppermost soil layer 
!                                                ! thickness/depth       ( m)
!                                                ! Can not be too thin as 
!                                                ! then definition of soil
!                                                ! properties (i.e. phyiscal
!                                                ! representation of) and 
!                                                ! accuarcy of
!                                                ! numerical solution come
!                                                ! into question. If it is too
!                                                ! thick, then resolution of
!                                                ! diurnal cycle not as valid.
REAL(KIND=JPRB) :: ZHOOK_HANDLE, ZHOOK_HANDLE_OMP
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('THERMAL_LAYER_CONF_1',0,zhook_handle)
!
IF (PRESENT(phc_out).AND..NOT.PRESENT(phc)) THEN
  CALL abor1_sfx("THERMAL_LAYERS_CONF:IF HC_OUT IS PRESENT, HC MUST BE PRESENT TOO.")
ELSEIF (PRESENT(ptc_out).AND..NOT.PRESENT(ptc)) THEN
  CALL abor1_sfx("THERMAL_LAYERS_CONF:IF TC_OUT IS PRESENT, TC MUST BE PRESENT TOO.")
ENDIF
!
iout= SIZE(pd_out,2)
!
iin = SIZE(pd,2)
!
IF (lhook) CALL dr_hook('THERMAL_LAYER_CONF_1',1,zhook_handle)
!
!$OMP PARALLEL PRIVATE(ZHOOK_HANDLE_OMP)
IF (lhook) CALL dr_hook('THERMAL_LAYER_CONF_2',0,zhook_handle_omp)
!$OMP DO PRIVATE(JI,ZD_IN,ZD_TOT,ZD_OUT,ZD_HALF,ZD_MID,ZW,ZHC,ZW_OUT,ZHC_OUT)
DO ji=1,SIZE(pd_out,1)
  !
  zd_in(0) = 0.
  DO jin=1,iin
    zd_in(jin) = zd_in(jin-1) + pd(ji,jin)
  END DO
  zd_tot = zd_in(iin)
  !
  !* Depths for the computational grid
  !
  !* surface like road or floor (thin grid at the surface, coarse at the bottom)
  !
  zd_out(0) = 0.
  !
  IF (htype=='ROAD ' .OR. htype=='FLOOR') THEN
    !
    CALL tebgrid(zd_tot,zd_out(1:),zd_g1)
    !
    DO jout=1,iout
      pd_out(ji,jout) = zd_out(jout) - zd_out(jout-1) ! Depths => Thickness of layer
    END DO
    !
  ELSE
    !
    !* surface like roof or wall (thin grid on both sides, coarse in the middle)
    !
    IF (mod(iout,2)==0) THEN   ! even number of output layers
      zd_half = zd_tot / 2.
    ELSE                       ! odd  number of output layers
      zd_mid = 2. * zd_tot / iout ! middle layer is arbitrarily fixed
      IF (iout==3) zd_mid = max(zd_mid,zd_tot-2.*zd_g1) ! to impose layers equal
                                                                   ! to ZD_G1 on both sides
      zd_half = (zd_tot-zd_mid) / 2.
      pd_out(ji,iout/2+1) = zd_mid
    END IF
    !
    CALL tebgrid(zd_half,zd_out(1:iout/2),zd_g1)
    !
    DO jout=1,iout
      !
      IF (jout<=iout/2) THEN
        pd_out(ji,jout) = zd_out(jout) - zd_out(jout-1) ! Depths => Thickness of layer
        pd_out(ji,iout+1-jout) = pd_out(ji,jout)
      ENDIF
      !
      !* recomputes Depths for further averagings
      IF (jout>1) zd_out(jout) = zd_out(jout-1) + pd_out(ji,jout)
      !
    END DO
    !
  ENDIF
  !
  IF (pd(ji,1)==xundef)  pd_out(ji,1:iout) = xundef
  !
  !-------------------------------------------------------------------------------
  !
  !* Averaging of the Heat Capacity and the Thermal conductivity
  !
  IF (PRESENT(ptc)) THEN
    zw(:) =1./ptc(ji,:)
  ELSE
    zw(:) =1.
  ENDIF
  IF (PRESENT(phc)) THEN
    zhc(:) = phc(ji,:)
  ELSE
    zhc(:) = 1.
  ENDIF
  !
  CALL av_thermal_data(iout,pd(ji,:),zd_in(1:iin),zd_out(1:iout),zhc,zw,zhc_out,zw_out)
  !
  IF (PRESENT(ptc_out)) THEN
    ptc_out(ji,:)=xundef
    WHERE (zw_out(:)/=xundef) ptc_out(ji,:)=1./zw_out(:)
  ENDIF
  IF (PRESENT(phc_out)) phc_out(ji,:) = zhc_out(:)
  !
ENDDO
!$OMP END DO
IF (lhook) CALL dr_hook('THERMAL_LAYER_CONF_2',1,zhook_handle_omp)
!$OMP END PARALLEL
!
!-------------------------------------------------------------------------------
CONTAINS
!-------------------------------------------------------------------------------
SUBROUTINE av_thermal_data(KOUT,PDD,PDD_IN,PDD_OUT,PF1,PF2,PF1_OUT,PF2_OUT)
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: KOUT
REAL, DIMENSION(:), INTENT(IN) :: PDD
REAL, DIMENSION(:), INTENT(IN) :: PDD_IN
REAL, DIMENSION(:), INTENT(IN) :: PDD_OUT
REAL, DIMENSION(:), INTENT(IN)  :: PF1
REAL, DIMENSION(:), INTENT(IN)  :: PF2
REAL, DIMENSION(:), INTENT(OUT) :: PF1_OUT
REAL, DIMENSION(:), INTENT(OUT) :: PF2_OUT
!
REAL    :: ZF1! ponderated field
REAL    :: ZF2! ponderated field
REAL    :: ZS ! sum of weights
REAL    :: ZC ! coefficient of ponderation
REAL    :: ZD_LIM ! limit of previous layer that has been treated
!
REAL    :: ZEPS=1.e-6
!
INTEGER :: JOUT, JIN
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!* total depth:
!
!
IF (pdd(1)==xundef) THEN
  pf1_out(:) = xundef
  pf2_out(:) = xundef
  !
ELSE
  !
  zf1 = 0.
  zf2 = 0.
  zs  = 0.
  jin = 1
  jout= 1
  zd_lim = 0.
  !
  DO
    !
    IF (jout>kout) EXIT
    !
    IF (pdd_in(jin)< pdd_out(jout)-zeps) THEN
      !
      !ZC = PDD_IN(JIN) - MAX(PDD_IN(JIN-1),PDD_OUT(JOUT-1))
      zc = pdd_in(jin) - zd_lim
      zf1 = zf1 + zc * pf1(jin)
      zf2 = zf2 + zc * pf2(jin)
      zs = zs + zc
      zd_lim = pdd_in(jin)
      !
      jin=jin+1
      !
    ELSE
      !
      !ZC = PDD_OUT(JOUT) - MAX(PDD_IN(JIN-1),PDD_OUT(JOUT-1))
      zc = pdd_out(jout) - zd_lim
      zf1 = zf1 + zc * pf1(jin)
      zf2 = zf2 + zc * pf2(jin)
      zs = zs + zc
      pf1_out(jout) = zf1/zs
      pf2_out(jout) = zf2/zs
      zd_lim = pdd_out(jout)
      !
      jout = jout+1
      zf1 = 0.
      zf2 = 0.
      zs  = 0.
      !
    END IF
    !
  END DO
  !
ENDIF
!
END SUBROUTINE av_thermal_data
!
!     #########
      SUBROUTINE tebgrid( PSOILDEPTH, PD_G, PD_G1 )

!     ##########################################################################
!
!!****  *TEBGRID*  
!!
!!    PURPOSE
!!    -------
!
!     Calculates the soil grid configuration using a simple
!     geometric relation for all sub-surface layers.
!     This algorithm assumes the total soil depth > 0 m
!         
!     
!!**  METHOD
!!    ------
!
!     Direct calculation
!
!!    EXTERNAL
!!    --------
!
!     None
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!    Noilhan and Planton (1989)
!!    Belair (1995)
!!    Boone (2000)
!!    Boone et al. (2000)
!!    Habets et al. (2003)
!!      
!!    AUTHOR
!!    ------
!!      A. Boone           * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    12/04/03
!!      B. Decharme    12/10 uppermost soil layer set to 1cm
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
REAL,  INTENT(IN)  :: PSOILDEPTH  ! total soil depth            (m)
!                                   
REAL, DIMENSION(:), INTENT(OUT) :: PD_G        ! depth of base of soil layers (m)
REAL, OPTIONAL,       INTENT(IN)  :: PD_G1       ! depth of first layer
!
!
!*      0.2    declarations of local variables
!
INTEGER                           :: JJ, JI, JNLVL
!
!
REAL, PARAMETER                   :: ZGRIDFACTOR = 3.0 ! soil depth factor
!                                                      ! of increase with depth
!                                                      ! for all *sub-surface* 
!                                                      ! layers. Note, uppermost
!                                                      ! layer fixed by other
!                                                      ! constraints.          (-)
!
REAL                              :: ZD_G1 = 0.01      ! uppermost soil layer 
!                                                      ! thickness/depth       (m)
!                                                      ! Can not be too thin as 
!                                                      ! then definition of soil
!                                                      ! properties (i.e. phyiscal
!                                                      ! representation of) and 
!                                                      ! accuarcy of
!                                                      ! numerical solution come
!                                                      ! into question. If it is too
!                                                      ! thick, then resolution of
!                                                      ! diurnal cycle not as valid.
!                                                      ! Also chosen to comply with
!                                                      ! remotely sensed soil moisture.
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
!        0.     Initialization
!               --------------
!
jnlvl = SIZE(pd_g)
!
IF (PRESENT(pd_g1)) zd_g1 = pd_g1
!
IF (psoildepth < jnlvl*zd_g1) THEN
  !
  !*       3.     In the LIMIT For extremely thin soils
  !               ------------------------------------------
  !               This should be a RARE occurance, but 
  !               accounted for none-the-less ...:
  !               hold the ratio between all layer 
  !               thicknesses constant. 
  DO jj = 1,jnlvl
    pd_g(jj) = jj*psoildepth/jnlvl
  ENDDO
  !
ELSE
  !
  pd_g(1)     = zd_g1
  pd_g(jnlvl) = psoildepth
  !
  DO jj=jnlvl-1,2,-1
    !*       1.     Assign soil layer depths
    !               ------------------------
    !               using a geometric relation
    !               for layers 2...N
    !               This is GENERAL rule.
    !               Note that the first soil layer
    !               is FIXED except for VERY thin
    !               soils (see #3 below).
    pd_g(jj) = pd_g(jj+1)/zgridfactor
    !*       2.     When the soil is sufficiently thin
    !               ------------------------------------------
    !               We recalculate layer depths such
    !               that all layer thicknesses are >= ZD_G1
    !               We favor keeping a minimum grid thickness
    !               OVER maintaining geometric relation
    !               for increasingly thin soils. This means
    !               that uppermost soil moisture is readily
    !               comparable (i.e. for same layer thickness)
    !               EVERYWHERE except for most thin soils (below).
    pd_g(jj) = max(pd_g(jj), jj*zd_g1)
    !
  ENDDO     
  !
ENDIF
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE tebgrid
!
END SUBROUTINE thermal_layers_conf