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,PHC,PTC,PD,PHC_OUT,PTC_OUT,PD_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_tebgrid
!
IMPLICIT NONE
!
!*       0.1   Declarations of arguments
!              -------------------------
!
 CHARACTER(LEN=5),     INTENT(IN)  :: htype     ! type of surface
REAL, DIMENSION(:,:), INTENT(IN)  :: phc       ! input Heat Capacity
REAL, DIMENSION(:,:), INTENT(IN)  :: ptc       ! input Thermal conductivity
REAL, DIMENSION(:,:), INTENT(IN)  :: pd        ! input Layer Thickness
REAL, DIMENSION(:,:), INTENT(OUT) :: phc_out   ! output Heat Capacity
REAL, DIMENSION(:,:), INTENT(OUT) :: ptc_out   ! output Thermal conductivity
REAL, DIMENSION(:,:), INTENT(OUT) :: pd_out    ! output Layer Thickness
!
!*       0.2   Declarations of local variables
!
REAL, DIMENSION(SIZE(PHC,1))   :: zd_tot    ! Total depth
REAL, DIMENSION(SIZE(PHC,1))   :: zd_half   ! Depth of the half of the total surface
!                                           ! (excluding central layer in case
!                                           ! of odd number of layers)
REAL, DIMENSION(SIZE(PHC,1))   :: zd_mid    ! Thickness of the layer in the middle
!                                           ! in case of odd number of layers
REAL, DIMENSION(SIZE(PHC,1),0:SIZE(PHC    ,2))::zd_in  ! Depth from the surface 
!                                                      ! to the layer bottom
REAL, DIMENSION(SIZE(PHC,1),0:SIZE(PHC_OUT,2))::zd_out ! Depth from the surface 
!                                                      ! to the layer bottom
REAL, DIMENSION(SIZE(PHC,1),SIZE(PHC,2))     :: zw     ! 1/TC
REAL, DIMENSION(SIZE(PHC,1),SIZE(PHC_OUT,2)) :: zw_out ! 1/TC
INTEGER                        :: iin       ! Number of layer in input data
INTEGER                        :: iout      ! Number of layer in output fields
INTEGER                        :: jin       ! 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.
!-------------------------------------------------------------------------------
!
iin = SIZE(phc,2)
iout= SIZE(phc_out,2)
!
!-------------------------------------------------------------------------------
!
!* Depths for the computational grid
!
!* total depth:
!
!
zd_in(:,0) = 0.
DO jin=1,iin
  zd_in(:,jin) = zd_in(:,jin-1) + pd(:,jin)
END DO
zd_tot(:) = zd_in(:,iin)
!
!* surface like road or floor (thin grid at the surface, coarse at the bottom)
!
IF (htype=='ROAD ' .OR. htype=='FLOOR') THEN
  zd_out(:,0) = 0.
  CALL tebgrid(zd_tot,zd_out(:,1:),zd_g1)
  pd_out(:,1) = zd_out(:,1)
  DO jout=2,iout
    pd_out(:,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(:,iout/2+1) = zd_mid(:)
  END IF
  zd_out(:,0) = 0.
  CALL tebgrid(zd_half,zd_out(:,1:iout/2),zd_g1)
  pd_out(:,1) = zd_out(:,1)
  DO jout=2,iout/2
    pd_out(:,jout) = zd_out(:,jout) - zd_out(:,jout-1) ! Depths => Thickness of layer
  END DO
  DO jout=1,iout/2
    pd_out(:,iout+1-jout) = pd_out(:,jout)
  END DO
  !* recomputes Depths for further averagings
  DO jout=2,iout
    zd_out(:,jout) = zd_out(:,jout-1) + pd_out(:,jout)
  END DO

END IF
!
DO jout=1,iout
  WHERE (pd(:,1)==xundef)  pd_out(:,jout) = xundef
END DO
!-------------------------------------------------------------------------------
!
!* Averaging of the Heat Capacity and the Thermal conductivity
!
zw=1./ptc(:,:)
 CALL av_thermal_data(phc,zw,phc_out,zw_out)
ptc_out=xundef
WHERE (zw_out/=xundef) ptc_out=1./zw_out
!
!-------------------------------------------------------------------------------
 CONTAINS
!-------------------------------------------------------------------------------
SUBROUTINE av_thermal_data(PF1,PF2,PF1_OUT,PF2_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
!
INTEGER :: jl ! loop counter on spatial points
REAL    :: zeps=1.e-6
!
DO jl=1,SIZE(pf1,1)
 IF (pd(jl,1)==xundef) THEN
   pf1_out(jl,:) = xundef
   pf2_out(jl,:) = xundef
   cycle
 END IF
 !
 zf1 = 0.
 zf2 = 0.
 zs  = 0.
 jin = 1
 jout= 1
 zd_lim = 0.
 DO
  IF (jout>iout) EXIT
  !
  IF (zd_in(jl,jin)< zd_out(jl,jout)-zeps) THEN
!    ZC = ZD_IN(JL,JIN) - MAX(ZD_IN(JL,JIN-1),ZD_OUT(JL,JOUT-1))
    zc = zd_in(jl,jin) - zd_lim
    zf1 = zf1 + zc * pf1(jl,jin)
    zf2 = zf2 + zc * pf2(jl,jin)
    zs = zs + zc
    zd_lim = zd_in(jl,jin)
    !
    jin=jin+1
  ELSE
!    ZC = ZD_OUT(JL,JOUT) - MAX(ZD_IN(JL,JIN-1),ZD_OUT(JL,JOUT-1))
    zc = zd_out(jl,jout) - zd_lim
    zf1 = zf1 + zc * pf1(jl,jin)
    zf2 = zf2 + zc * pf2(jl,jin)
    zs = zs + zc
    pf1_out(jl,jout) = zf1/zs
    pf2_out(jl,jout) = zf2/zs
    zd_lim = zd_out(jl,jout)
    !
    jout = jout+1
    zf1 = 0.
    zf2 = 0.
    zs  = 0.
  END IF
 END DO
END DO
!
END SUBROUTINE
!
END SUBROUTINE thermal_layers_conf