init_isba_sbl.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 init_isba_sbl(HISBA, HCPSURF, KLVL, PTSTEP, PPA, PPS, PTA, PQA, PRHOA, PU, PV,   &
                               pdir_sw, psca_sw, psw_bands, prain, psnow,                       &
                               pzref, puref, ptg, ppatch, pwg, pwgi, pz0, psso_slope,           &
                               presa, pveg, plai, pwr, prgl, prsmin, pgamma, pwrmax_cf,         &
                               pz0_o_z0h, pwfc, pwsat, ptsnow, pz, pt, pq, pwind, ptke, pp)  
!     #################################################################################
!
!!****  *INIT_WATER_SBL* - inits water SBL profiles
!!
!!    PURPOSE
!!    -------
!
!!**  METHOD
!!    ------
!!
!!    REFERENCE
!!    ---------
!!      
!!
!!    AUTHOR
!!    ------
!!     S. Riette
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    03/2010
!!------------------------------------------------------------------
!
USE modd_type_snow
!
USE modd_csts,             ONLY : xcpd, xrd, xp00, xg, xlvtt
USE modd_surf_atm,         ONLY : lnosof
USE modd_canopy_turb,      ONLY : xalpsbl
!
USE modi_cls_tq
USE modi_isba_snow_frac
USE modi_wet_leaves_frac
USE modi_veg
USE modi_drag
USE modi_cls_wind
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
 CHARACTER(LEN=*)  , INTENT(IN)  :: hisba     ! type of ISBA version
 CHARACTER(LEN=*)  , INTENT(IN)  :: hcpsurf   ! specific heat at surface
REAL,               INTENT(IN)   :: ptstep   ! timestep of the integration
INTEGER           , INTENT(IN)  :: klvl      ! number      of levels in canopy
REAL, DIMENSION(:), INTENT(IN)  :: ppa       ! pressure at forcing level             (Pa)
REAL, DIMENSION(:), INTENT(IN)  :: pps       ! pressure at atmospheric model surface (Pa)
REAL, DIMENSION(:), INTENT(IN)  :: pta       ! air temperature forcing               (K)
REAL, DIMENSION(:), INTENT(IN)  :: pqa       ! air humidity forcing                  (kg/m3)
REAL, DIMENSION(:), INTENT(IN)  :: prhoa     ! air density                           (kg/m3)
REAL, DIMENSION(:), INTENT(IN)  :: pu        ! zonal wind                            (m/s)
REAL, DIMENSION(:), INTENT(IN)  :: pv        ! meridian wind                         (m/s)
REAL, DIMENSION(:,:),INTENT(IN) :: pdir_sw   ! direct  solar radiation (on horizontal surf.)
!                                            !                                       (W/m2)
REAL, DIMENSION(:,:),INTENT(IN) :: psca_sw   ! diffuse solar radiation (on horizontal surf.)
!                                            !                                        (W/m2)
REAL, DIMENSION(:), INTENT(IN)  :: psw_bands ! mean wavelength of each shortwave band (m)
REAL, DIMENSION(:), INTENT(IN)  :: psnow     ! snow precipitation                    (kg/m2/s)
REAL, DIMENSION(:), INTENT(IN)  :: prain     ! liquid precipitation                  (kg/m2/s)
REAL, DIMENSION(:), INTENT(IN)  :: pzref     ! height of T,q forcing                 (m)
REAL, DIMENSION(:), INTENT(IN)  :: puref     ! height of wind forcing                (m)
REAL, DIMENSION(:,:), INTENT(IN):: ptg       ! surface and sub-surface soil temperature profile (K)   
REAL, DIMENSION(:,:), INTENT(IN):: ppatch    ! fraction of each tile/patch 
REAL, DIMENSION(:,:), INTENT(IN):: pwg       ! soil volumetric water content profile   (m3/m3)
REAL, DIMENSION(:,:), INTENT(IN):: pwgi      ! soil liquid water equivalent volumetric
REAL, DIMENSION(:,:), INTENT(IN):: pz0       ! roughness length
REAL, DIMENSION(:), INTENT(IN)  :: psso_slope! slope of S.S.O.                         (-)
REAL, DIMENSION(:,:), INTENT(IN):: presa     ! aerodynamic resistance                  (s/m)
REAL, DIMENSION(:,:), INTENT(IN):: pveg      ! vegetation cover fraction               (-)
REAL, DIMENSION(:,:), INTENT(IN):: plai      ! Leaf Area Index                         (m2/m2)
REAL, DIMENSION(:,:), INTENT(IN):: pwr       ! liquid water retained on the
!                                            ! foliage of the vegetation
!                                            ! canopy       
REAL, DIMENSION(:,:), INTENT(IN):: prgl      ! maximum solar radiation
!                                            ! usable in photosynthesis                (W/m2)
REAL, DIMENSION(:,:), INTENT(IN):: prsmin    ! minimum stomatal resistance             (s/m)
REAL, DIMENSION(:,:), INTENT(IN):: pgamma    ! coefficient for the calculation
!                                            ! of the surface stomatal
!                                            ! resistance
REAL, DIMENSION(:,:), INTENT(IN):: pwrmax_cf ! coefficient for maximum water 
!                                            ! interception 
!                                            ! storage capacity on the vegetation      (-)
REAL, DIMENSION(:,:), INTENT(IN):: pz0_o_z0h ! ratio of surface roughness lengths
!                                            ! (momentum to heat)                      (-)
REAL, DIMENSION(:,:), INTENT(IN):: pwfc      ! field capacity volumetric water content
!                                            ! profile                                 (m3/m3)
REAL, DIMENSION(:,:), INTENT(IN):: pwsat     ! porosity profile                        (m3/m3) 
TYPE(surf_snow)     , INTENT(IN):: ptsnow    ! snow state
REAL, DIMENSION(:,:), INTENT(IN):: pz        ! height of middle of each level grid   (m)
!
REAL, DIMENSION(:,:), INTENT(OUT) :: pt   ! temperature at each level in SBL      (m/s)
REAL, DIMENSION(:,:), INTENT(OUT) :: pq   ! humidity    at each level in SBL      (kg/m3)
REAL, DIMENSION(:,:), INTENT(OUT) :: pwind! wind        at each level in SBL      (m/s)
REAL, DIMENSION(:,:), INTENT(OUT) :: ptke ! Tke         at each level in SBL      (m2/s2)
REAL, DIMENSION(:,:), INTENT(OUT) :: pp   ! pressure    at each level in SBL      (kg/m3)
!
!*      0.2    declarations of local variables
!
!* forcing variables
!
REAL, DIMENSION(SIZE(PTA))   :: zwind    ! lowest atmospheric level wind speed           (m/s)
REAL, DIMENSION(SIZE(PTA))   :: zexna    ! Exner function at lowest SBL scheme level     (-)
REAL, DIMENSION(SIZE(PTA))   :: zqa      ! specific humidity                             (kg/m3)
!
! SBL turbulence scheme
!
REAL, DIMENSION(SIZE(PTA))   ::zri
REAL, DIMENSION(SIZE(PTA))   ::zcd
REAL, DIMENSION(SIZE(PTA))   ::zcdn
REAL, DIMENSION(SIZE(PTA))   ::zch
REAL, DIMENSION(SIZE(PTA))   ::ztnm
REAL, DIMENSION(SIZE(PTA))   ::zqnm
REAL, DIMENSION(SIZE(PTA))   ::zhunm
REAL, DIMENSION(SIZE(PTA))   ::zp_slope_cos
REAL, DIMENSION(SIZE(PTA))   ::zz0
REAL, DIMENSION(SIZE(PTA))   ::zz0h
REAL, DIMENSION(SIZE(PTA))   ::zexns
REAL, DIMENSION(SIZE(PTA))   ::zts
REAL, DIMENSION(SIZE(PTA))   ::zhu
REAL, DIMENSION(SIZE(PTA))   ::zqs
REAL, DIMENSION(SIZE(PTA))   ::zz0eff
REAL, DIMENSION(SIZE(PTA))   ::zwg
REAL, DIMENSION(SIZE(PTA))   ::zwgi
REAL, DIMENSION(SIZE(PTA))   ::zveg
REAL, DIMENSION(SIZE(PTA))   ::zresa
REAL, DIMENSION(SIZE(PTA))   ::zhug
REAL, DIMENSION(SIZE(PTA))   ::zhugi
REAL, DIMENSION(SIZE(PTA))   ::zhv
REAL, DIMENSION(SIZE(PTA))   ::zcps
REAL, DIMENSION(SIZE(PTA))   ::zwrmax_cf
REAL, DIMENSION(SIZE(PTA))   ::zwr
REAL, DIMENSION(SIZE(PTA))   ::zz0_with_snow
REAL, DIMENSION(SIZE(PTA))   ::zpsng
REAL, DIMENSION(SIZE(PTA))   ::zpsnv
REAL, DIMENSION(SIZE(PTA))   ::zpsnv_a
REAL, DIMENSION(SIZE(PTA))   ::zpsn
REAL, DIMENSION(SIZE(PTA))   ::zsnowalb
REAL, DIMENSION(SIZE(PTA),SIZE(PTSNOW%WSNOW,2)) ::zsnowswe
REAL, DIMENSION(SIZE(PTA),SIZE(PTSNOW%WSNOW,2)) ::zsnowrho
REAL, DIMENSION(SIZE(PTA))   ::zffg
REAL, DIMENSION(SIZE(PTA))   ::zffgnos
REAL, DIMENSION(SIZE(PTA))   ::zffv
REAL, DIMENSION(SIZE(PTA))   ::zffvnos
REAL, DIMENSION(SIZE(PTA))   ::zff
REAL, DIMENSION(SIZE(PTA))   ::zrs
REAL, DIMENSION(SIZE(PTA))   ::zp_global_sw
REAL, DIMENSION(SIZE(PTA))   ::zf2
REAL, DIMENSION(SIZE(PTA))   ::zf5
REAL, DIMENSION(SIZE(PTA))   ::zlai
REAL, DIMENSION(SIZE(PTA))   ::zgamma
REAL, DIMENSION(SIZE(PTA))   ::zrgl
REAL, DIMENSION(SIZE(PTA))   ::zrsmin
REAL, DIMENSION(SIZE(PTA))   ::zdelta
REAL, DIMENSION(SIZE(PTA))   ::zwrmax
REAL, DIMENSION(SIZE(PTA))   ::zcls_wind_zon
REAL, DIMENSION(SIZE(PTA))   ::zcls_wind_mer
REAL, DIMENSION(SIZE(PTA),SIZE(PTSNOW%WSNOW,2))   ::zsum_layer
REAL, DIMENSION(SIZE(PTA))   ::zsum
REAL, DIMENSION(SIZE(PTA))   :: zleg_delta  ! soil evaporation delta fn
REAL, DIMENSION(SIZE(PTA))   :: zlegi_delta ! soil sublimation delta fn
REAL, DIMENSION(SIZE(PTA))   :: zlvtt
!
INTEGER                     :: jswb
INTEGER                     :: jlayer
INTEGER                     :: jpatch
!
REAL, DIMENSION(SIZE(PTA),SIZE(PPATCH,2)) ::zwsnow
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('INIT_ISBA_SBL',0,zhook_handle)
!    
!Means over patches
zts     = sum(ptg(:,:)*ppatch(:,:) ,dim=2)
zwg     = sum(pwg(:,:)*ppatch(:,:) ,dim=2)
zwgi    = sum(pwgi(:,:)*ppatch(:,:),dim=2)
zz0     = sum(ppatch(:,:)*pz0(:,:) ,dim=2)
!
!We choose to set ZZ0EFF and ZZ0_WITH_SNOW equal to ZZ0
zz0eff        = zz0
zz0_with_snow = zz0
zz0h(:) = sum(ppatch(:,:) * pz0(:,:)/pz0_o_z0h(:,:),dim=2)
zveg(:) = sum(ppatch(:,:) * pveg(:,:)              ,dim=2)
!
zp_slope_cos(:) = 1./sqrt(1.+psso_slope(:)**2)
IF (lnosof) zp_slope_cos(:) = 1.0
!
zresa(:) = sum(ppatch(:,:)*presa(:,:),dim=2)
WHERE(zveg(:)>0)
  zlai(:)= sum(ppatch(:,:)*pveg(:,:)*plai(:,:)     ,dim=2,mask=pveg(:,:)>0) / zveg(:)
  zwrmax_cf(:)= sum(ppatch(:,:)*pveg(:,:)*pwrmax_cf(:,:),dim=2,mask=pveg(:,:)>0) / zveg(:)
  zwr(:)= sum(ppatch(:,:)*pveg(:,:)*pwr(:,:)      ,dim=2,mask=pveg(:,:)>0) / zveg(:)
ELSEWHERE
  zlai(:) = plai(:,1)
  zwrmax_cf(:) = pwrmax_cf(:,1)
  zwr(:) = pwr(:,1)
ENDWHERE
!
zsum_layer(:,:) = 0.
zsum(:) = 0.
!
DO jlayer=1,SIZE(ptsnow%WSNOW,2)
  zsnowswe(:,jlayer) = sum(ppatch(:,:)*ptsnow%WSNOW(:,jlayer,:),dim=2)
  zsum_layer(:,jlayer) = sum(ppatch(:,:),dim=2,mask=ptsnow%WSNOW(:,jlayer,:)>0)
  WHERE(zsum_layer(:,jlayer)>0)      
    zsnowrho(:,jlayer)= sum( ppatch(:,:)*ptsnow%RHO(:,jlayer,:), dim=2, &
                             mask=ptsnow%WSNOW(:,jlayer,:)>0) / zsum_layer(:,jlayer)
  ELSEWHERE
    zsnowrho(:,jlayer)=ptsnow%RHO(:,jlayer,1)
  ENDWHERE
END DO
!
zsum(:)=sum(zsum_layer(:,:),dim=2)
!
zwsnow(:,:) = 0.
DO jpatch=1,SIZE(ptsnow%WSNOW,3)
  DO jlayer=1,SIZE(ptsnow%WSNOW,2)
    zwsnow(:,jpatch) = zwsnow(:,jpatch) + ptsnow%WSNOW(:,jlayer,jpatch)
  ENDDO
ENDDO    
!
WHERE(zsum(:)>0)         
  zsnowalb(:) = sum(ppatch(:,:)*ptsnow%ALB(:,:),dim=2,mask=zwsnow(:,:)>0) / zsum(:)      
ELSEWHERE
  zsnowalb(:) = ptsnow%ALB(:,1)
ENDWHERE
!
zrgl(:) = sum(ppatch(:,:) * prgl(:,:),dim=2)
zrsmin(:) = sum(ppatch(:,:) * prsmin(:,:),dim=2)
zgamma(:) = sum(ppatch(:,:) * pgamma(:,:),dim=2)
!
zexna(:) = (ppa(:)/xp00)**(xrd/xcpd)
zexns(:) = (pps(:)/xp00)**(xrd/xcpd)
zqa(:) = pqa(:) / prhoa(:)
zwind(:) = sqrt(pu**2+pv**2)
!
!We compute the snow fractions
 CALL isba_snow_frac(ptsnow%SCHEME,                      &
                    zsnowswe, zsnowrho, zsnowalb,       &
                    zveg, zlai, zz0,                    &
                    zpsn, zpsnv_a, zpsng, zpsnv         )  
!
!We compute total shortwave incoming radiation needed by veg
zp_global_sw(:) = 0.
DO jswb=1,SIZE(psw_bands)
  zp_global_sw(:)   = zp_global_sw(:) + (pdir_sw(:,jswb) + psca_sw(:,jswb))
END DO
!
!We choose the case HPHOTO=='NON' and a humid soil (ZF2=1) to compute ZRS
zf2(:)=1.0
 CALL veg(zp_global_sw, pta, zqa, pps, zrgl, zlai, zrsmin, zgamma, zf2, zrs)
!Calculation of ZDELTA
 CALL wet_leaves_frac(zwr, zveg, zwrmax_cf, zz0_with_snow, zlai, zwrmax, zdelta)
!
!We choose the case LFLOOD=false
zffg(:) = 0.0
zffgnos(:) = 0.0
zffv(:) = 0.0
zffvnos(:) = 0.0
zff(:) = 0.0
!
zf5(:) = 1.0
zlvtt(:) = xlvtt
!We compute ZCD, ZCH and ZRI
 CALL drag(hisba, ptsnow%SCHEME, hcpsurf,  ptstep,                            &
          zts, zwg, zwgi, zexns, zexna, pta,                                  &
          zwind, zqa, prain, psnow, pps, zrs,                                 &
          zveg, zz0, zz0eff, zz0h, pwfc(:,1), pwsat(:,1),                     &
          zpsng, zpsnv, pzref, puref, zp_slope_cos, zdelta, zf5,              &
          zresa, zch, zcd, zcdn, zri, zhug, zhugi, zhv, zhu, zcps,            &
          zqs, zffg, zffv, zff, zffgnos, zffvnos, zleg_delta, zlegi_delta,    &
          zwr, prhoa, zlvtt                                                   )  
!
!Initialisation of T, Q, Wind and TKE on all canopy levels
DO jlayer=1,klvl
  !
  CALL cls_tq(pta, zqa, ppa, pps, pzref, zcd, zch, zri, zts, zhu, zz0h, &
              pz(:,jlayer), ztnm, zqnm, zhunm                           ) 
  ! 
  pt(:,jlayer)=ztnm
  pq(:,jlayer)=zqnm
  !
  CALL cls_wind(pu, pv, puref, zcd, zcdn, zri, pz(:,jlayer), &
                zcls_wind_zon, zcls_wind_mer                 )
  !
  pwind(:,jlayer) = sqrt( zcls_wind_zon(:)**2 + zcls_wind_mer(:)**2 )
  ptke(:,jlayer) = xalpsbl * zcd(:) * ( pu(:)**2 + pv(:)**2 )
  pp(:,jlayer) = ppa(:) + xg * prhoa(:) * (pz(:,klvl) - pz(:,jlayer))
  !
ENDDO
!
IF (lhook) CALL dr_hook('INIT_ISBA_SBL',1,zhook_handle) 
!
END SUBROUTINE init_isba_sbl