ch_bvocemn.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 ch_bvocem_n (CHI, GB, I, &
                             psw_forbio,prhoa,psfts)
!!   ###############################
!!
!!***  *BVOCEM*
!! 
!!    PURPOSE
!!    -------
!!    Calculate the biogenic emission fluxes according to the 
!!    subgrid vegetation given by the soil interface
!!
!!    METHOD
!!    ------
!!
!!
!!    AUTHOR
!!    ------
!!    F. Solmon (LA) & V. Masson (CNRM)
!!    
!!    MODIFICATIONS
!!    -------------
!!    Original: 25/10/00
!!    P. Tulet  30/07/03 externalisation of biogenics fluxes (2D => 1D)
!!    R. Alkama 04/2012  add 7 new vegtype (19 instead 12)
!!
!!    EXTERNAL
!!    --------
!
USE modd_ch_isba_n, ONLY : ch_isba_t
USE modd_gr_biog_n, ONLY : gr_biog_t
USE modd_isba_n, ONLY : isba_t
!
USE modi_vegtype_to_patch
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
USE modd_bvoc_par
USE modd_csts,ONLY : xmd, xavogadro
USE modd_co2v_par
USE modd_surf_par,ONLY:xundef
USE modd_isba_par
USE modd_data_cover_par, ONLY : nvegtype, nvt_tebd, nvt_bone, nvt_trbe, &
                                nvt_trbd, nvt_tebe, nvt_tene, nvt_bobd, &
                                nvt_bond, nvt_shrb, nvt_bogr, nvt_gras, & 
                                nvt_trog, nvt_park, nvt_c3, nvt_c4,     &
                                nvt_irr 
!!
!!
!------------------------------------------------------------------------------
!
!*       0.   DECLARATIONS
!        -----------------
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE

!
TYPE(ch_isba_t), INTENT(INOUT) :: chi
TYPE(gr_biog_t), INTENT(INOUT) :: gb
TYPE(isba_t), INTENT(INOUT) :: i
!
REAL, DIMENSION(:,:), INTENT(IN)    :: psw_forbio
REAL, DIMENSION(:),   INTENT(IN)    :: prhoa
REAL, DIMENSION(:,:), INTENT(INOUT) :: psfts
!
!*       0.1  declaration of arguments
!
!*   0.1 Declaration of local variables
!
REAL, DIMENSION(SIZE(PSW_FORBIO,1)) :: zrad_par,  zlcor_rad
!                            PAR radiation in case of ISBA-STD use
!
REAL, DIMENSION(SIZE(PSW_FORBIO,1)) :: zfiso_for  , zfmono_for,   &
                                       zfiso_grass, zfmono_grass, &
                                       zfiso_crop , zfmono_crop     
!                                Fluxes coming from different landuse
REAL, DIMENSION(SIZE(PSW_FORBIO,1), NVEGTYPE) :: ztcor ,ztcorm
!
REAL, DIMENSION(SIZE(PSW_FORBIO,1),SIZE(I%XABC),NVEGTYPE) :: zbvocpar 
!                                PAR at gauss level in micromolphot/m2/s
!
REAL, DIMENSION(SIZE(PSW_FORBIO,1)) :: zisopot, zmonopot, zratio
!
INTEGER:: kngauss     
!                        nbre of gauss level in integration
!                        index of patch corresponding to forest(+ligneaous)
INTEGER:: jpatch, jsv
REAL(KIND=JPRB) :: zhook_handle
!
!------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('CH_BVOCEM_N',0,zhook_handle)
!
!* 1. Contribution of forest and ligneous vegetation 
!   from ISOPOT and MONOPOT maps 
!   ------------------------------------------------
!
!* 1.0 Preliminary : patch index corresponding to forest
!
!2.Contribution of other types of vegetation than forest, consider the vegtype fraction in the pixel 
!------------------------------------------------------------------------------------------
!
!* 2.0 Preliminary : patch index corresponding to grassland, crops (C3+C4)
!
!1.1.1 Using ISBA_Ags explicit light attenuation 
! number of g Gauss level for the integration 
IF (i%CPHOTO/='NON') THEN
  kngauss = SIZE(i%XABC)
ELSE
  !1.1.2 using isba std version 
  zrad_par(:)= 0.
  DO jpatch = 1,i%NPATCH
    zrad_par(:)= zrad_par(:) +(psw_forbio(:,jpatch)*i%XPATCH(:,jpatch) ) * xparcf * 4.7 
  END DO
  zlcor_rad(:) = zlcor_func(zrad_par(:))
ENDIF
!  
!
 CALL by_patch(nvt_tebd, ztcor(:,nvt_tebd), ztcorm(:,nvt_tebd))
 CALL by_patch(nvt_bone, ztcor(:,nvt_bone), ztcorm(:,nvt_bone))
 CALL by_patch(nvt_trbe, ztcor(:,nvt_trbe), ztcorm(:,nvt_trbe))
 CALL by_patch(nvt_trbd, ztcor(:,nvt_trbd), ztcorm(:,nvt_trbd))
 CALL by_patch(nvt_tebe, ztcor(:,nvt_tebe), ztcorm(:,nvt_tebe))
 CALL by_patch(nvt_tene, ztcor(:,nvt_tene), ztcorm(:,nvt_tene))
 CALL by_patch(nvt_bobd, ztcor(:,nvt_bobd), ztcorm(:,nvt_bobd))
 CALL by_patch(nvt_bond, ztcor(:,nvt_bond), ztcorm(:,nvt_bond))
 CALL by_patch(nvt_shrb, ztcor(:,nvt_shrb), ztcorm(:,nvt_shrb))
 CALL by_patch(nvt_bogr, ztcor(:,nvt_bogr), ztcorm(:,nvt_bogr))
 CALL by_patch(nvt_gras, ztcor(:,nvt_gras), ztcorm(:,nvt_gras))
 CALL by_patch(nvt_trog, ztcor(:,nvt_trog), ztcorm(:,nvt_trog))
 CALL by_patch(nvt_park, ztcor(:,nvt_park), ztcorm(:,nvt_park))
 CALL by_patch(nvt_c3  , ztcor(:,nvt_c3)  , ztcorm(:,nvt_c3)  )
 CALL by_patch(nvt_c4  , ztcor(:,nvt_c4)  , ztcorm(:,nvt_c4)  )
 CALL by_patch(nvt_irr , ztcor(:,nvt_irr) , ztcorm(:,nvt_irr) )
!
!
zratio(:) = i%XVEGTYPE(:,nvt_tebd) + i%XVEGTYPE(:,nvt_bone) + i%XVEGTYPE(:,nvt_trbe) + &
             i%XVEGTYPE(:,nvt_trbd) + i%XVEGTYPE(:,nvt_tebe) + i%XVEGTYPE(:,nvt_tene) + &
             i%XVEGTYPE(:,nvt_bobd) + i%XVEGTYPE(:,nvt_bond) + i%XVEGTYPE(:,nvt_shrb)
!
WHERE (zratio(:)/=0.)
  zisopot(:) = gb%XISOPOT (:) / zratio(:)
  zmonopot(:) = gb%XMONOPOT(:) / zratio(:) 
ELSEWHERE
  zisopot(:) = 0.
  zmonopot(:) = 0.
END WHERE
!                       
 CALL by_veg9(nvt_tebd, nvt_bone, nvt_trbe, nvt_trbd, nvt_tebe, nvt_tene, nvt_bobd, &
                    nvt_bond, nvt_shrb, zisopot, zmonopot, zfiso_for, zfmono_for)
!
zisopot(:) = xisopot_grass
zmonopot(:) = xmonopot_grass
 CALL by_veg4(nvt_gras, nvt_trog, nvt_park, nvt_bogr, zisopot, zmonopot, zfiso_grass, zfmono_grass)
!
zisopot(:) = xisopot_crop
zmonopot(:) = xmonopot_crop
 CALL by_veg3(nvt_c3, nvt_c4, nvt_irr, zisopot, zmonopot, zfiso_crop, zfmono_crop)
!
!---------------------------------------------------------------------------------------
!
!3.Summation of different contribution for fluxes 
!------------------------------------------------
!
!isoprene in ppp.m.s-1
gb%XFISO (:)=(3.0012e-10/3600.) * ( zfiso_for(:) + zfiso_grass(:) + zfiso_crop(:) ) + 1e-17
!monoterpenes
gb%XFMONO(:)=(1.5006e-10/3600.) * ( zfmono_for(:) + zfmono_grass(:)+ zfmono_crop(:) ) + 1e-17
!
! conversion in molecules/m2/s
!
gb%XFISO(:)  = gb%XFISO(:)  * xavogadro * prhoa(:) / xmd
gb%XFMONO(:) = gb%XFMONO(:) * xavogadro * prhoa(:) / xmd
!
DO jsv=chi%SVI%NSV_CHSBEG,chi%SVI%NSV_CHSEND
  IF (chi%SVI%CSV(jsv) == "BIO") THEN
    ! RELACS CASE
    psfts(:,jsv) = psfts(:,jsv) + (gb%XFISO(:) + gb%XFMONO(:)) 
  ELSE IF (chi%SVI%CSV(jsv) == "ISO" .OR. chi%SVI%CSV(jsv) == "ISOP") THEN
    ! RACM CASE
    psfts(:,jsv) = psfts(:,jsv) + gb%XFISO(:)  
  ELSE IF (chi%SVI%CSV(jsv) == "API"  .OR. chi%SVI%CSV(jsv) == "LIM" .OR. &
           chi%SVI%CSV(jsv) == "BIOL" .OR. chi%SVI%CSV(jsv) == "BIOH" ) THEN
    ! RACM CASE
    ! CACM or RELACS 2 CASE     
    psfts(:,jsv) = psfts(:,jsv) + 0.5 * gb%XFMONO(:) 
  ENDIF
END DO
!
!**********************************************************************************
IF (lhook) CALL dr_hook('CH_BVOCEM_N',1,zhook_handle)
 CONTAINS
!
SUBROUTINE by_patch(NVT_VEGTYPE,PTCOR,PTCORM)
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: nvt_vegtype
REAL, DIMENSION(:), INTENT(OUT) :: ptcor
REAL, DIMENSION(:), INTENT(OUT) :: ptcorm
!
REAL, DIMENSION(SIZE(PSW_FORBIO,1)) :: zbvocsg
REAL, DIMENSION(SIZE(PSW_FORBIO,1),SIZE(I%XABC)) :: zbvocpar 
INTEGER:: ipatch, jlayer, it
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('CH_BVOCEM_N:BY_PATCH',0,zhook_handle)
!
ipatch = vegtype_to_patch(nvt_vegtype, i%NPATCH)
!
ptcor(:) = 0.
ptcorm(:) = 0.
DO it=1,SIZE(i%XTG,1)
  IF (i%XTG(it,1,ipatch).LE.1000.) THEN
    ptcorm(it)=ztcorm0_func(i%XTG(it,1,ipatch))
    ptcor(it)=ztcor0_func(i%XTG(it,1,ipatch))
  ENDIF
ENDDO
!
IF (i%CPHOTO/='NON') THEN
  !PAR over Forest canopies, in micro-molE.m-2.s-1 
  zbvocpar(:,:) = gb%XIACAN(:,:,ipatch)*4.7
  !Calculation of radiative attenuation effect in the canopy on correction factor
  zbvocsg(:) = 0.
  DO jlayer=1,kngauss
    zbvocsg(:) = zbvocsg(:) + i%XPOI(jlayer) * zlcor_func(zbvocpar(:,jlayer)) 
  ENDDO
  ptcor(:) = ptcor(:) * zbvocsg(:)
ELSE
  ptcor(:) = ptcor(:) * xcanfac * zlcor_rad(:)
ENDIF
!
IF (lhook) CALL dr_hook('CH_BVOCEM_N:BY_PATCH',1,zhook_handle)
!
END SUBROUTINE by_patch
!--------------------------------------------------------------------------
SUBROUTINE by_veg3(NVT_V1, NVT_V2, NVT_V3, &
                  pisopot, pmonopot, pfiso, pfmono)
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: nvt_v1
INTEGER, INTENT(IN) :: nvt_v2
INTEGER, INTENT(IN) :: nvt_v3
REAL, DIMENSION(:), INTENT(IN) :: pisopot
REAL, DIMENSION(:), INTENT(IN) :: pmonopot
REAL, DIMENSION(:), INTENT(OUT) :: pfiso
REAL, DIMENSION(:), INTENT(OUT) :: pfmono
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('CH_BVOCEM_N:BY_VEG3',0,zhook_handle)
!
!isoprene flux 
!!
!! warning, XISOPOT external map accounts for the total forest fraction
WHERE ( i%XVEGTYPE(:,nvt_v1) + i%XVEGTYPE(:,nvt_v2) + i%XVEGTYPE(:,nvt_v3) > 0. )
  !
  pfiso(:) = pisopot(:) *                   &
     ( ztcor(:,nvt_v1) * i%XVEGTYPE(:,nvt_v1) &
      +ztcor(:,nvt_v2) * i%XVEGTYPE(:,nvt_v2) &
      +ztcor(:,nvt_v3) * i%XVEGTYPE(:,nvt_v3) )
  !
  pfmono(:) = pmonopot(:) *                  &
     ( ztcorm(:,nvt_v1) * i%XVEGTYPE(:,nvt_v1) &
      +ztcorm(:,nvt_v2) * i%XVEGTYPE(:,nvt_v2) &
      +ztcorm(:,nvt_v3) * i%XVEGTYPE(:,nvt_v3) )          
  !
ELSEWHERE
  !
  pfiso(:) = 0.
  pfmono(:) = 0.
  !
END WHERE
!
IF (lhook) CALL dr_hook('CH_BVOCEM_N:BY_VEG3',1,zhook_handle)
!
END SUBROUTINE by_veg3
!--------------------------------------------------------------------------
SUBROUTINE by_veg4(NVT_V1, NVT_V2, NVT_V3, NVT_V4,&
                  pisopot, pmonopot, pfiso, pfmono)
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: nvt_v1
INTEGER, INTENT(IN) :: nvt_v2
INTEGER, INTENT(IN) :: nvt_v3
INTEGER, INTENT(IN) :: nvt_v4
REAL, DIMENSION(:), INTENT(IN) :: pisopot
REAL, DIMENSION(:), INTENT(IN) :: pmonopot
REAL, DIMENSION(:), INTENT(OUT) :: pfiso
REAL, DIMENSION(:), INTENT(OUT) :: pfmono
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('CH_BVOCEM_N:BY_VEG4',0,zhook_handle)
!
!isoprene flux 
!!
!! warning, XISOPOT external map accounts for the total forest fraction
WHERE ( i%XVEGTYPE(:,nvt_v1) + i%XVEGTYPE(:,nvt_v2) + i%XVEGTYPE(:,nvt_v3) &
       +i%XVEGTYPE(:,nvt_v4) > 0. )
  !
  pfiso(:) = pisopot(:) *                   &
     ( ztcor(:,nvt_v1) * i%XVEGTYPE(:,nvt_v1) &
      +ztcor(:,nvt_v2) * i%XVEGTYPE(:,nvt_v2) &
      +ztcor(:,nvt_v3) * i%XVEGTYPE(:,nvt_v3) &
      +ztcor(:,nvt_v4) * i%XVEGTYPE(:,nvt_v4) )
  !
  pfmono(:) = pmonopot(:) *                  &
     ( ztcorm(:,nvt_v1) * i%XVEGTYPE(:,nvt_v1) &
      +ztcorm(:,nvt_v2) * i%XVEGTYPE(:,nvt_v2) &
      +ztcorm(:,nvt_v3) * i%XVEGTYPE(:,nvt_v3) &
      +ztcorm(:,nvt_v4) * i%XVEGTYPE(:,nvt_v4) )            
  !
ELSEWHERE
  !
  pfiso(:) = 0.
  pfmono(:) = 0.
  !
END WHERE
!
IF (lhook) CALL dr_hook('CH_BVOCEM_N:BY_VEG4',1,zhook_handle)
!
END SUBROUTINE by_veg4
!--------------------------------------------------------------------------
SUBROUTINE by_veg9(NVT_V1, NVT_V2, NVT_V3, NVT_V4, NVT_V5, NVT_V6,     &
               nvt_v7, nvt_v8, nvt_v9, pisopot, pmonopot, pfiso, pfmono)
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: nvt_v1
INTEGER, INTENT(IN) :: nvt_v2
INTEGER, INTENT(IN) :: nvt_v3
INTEGER, INTENT(IN) :: nvt_v4
INTEGER, INTENT(IN) :: nvt_v5
INTEGER, INTENT(IN) :: nvt_v6
INTEGER, INTENT(IN) :: nvt_v7
INTEGER, INTENT(IN) :: nvt_v8
INTEGER, INTENT(IN) :: nvt_v9
REAL, DIMENSION(:), INTENT(IN) :: pisopot
REAL, DIMENSION(:), INTENT(IN) :: pmonopot
REAL, DIMENSION(:), INTENT(OUT) :: pfiso
REAL, DIMENSION(:), INTENT(OUT) :: pfmono
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('CH_BVOCEM_N:BY_VEG9',0,zhook_handle)
!
!isoprene flux 
!!
!! warning, XISOPOT external map accounts for the total forest fraction
WHERE ( i%XVEGTYPE(:,nvt_v1) + i%XVEGTYPE(:,nvt_v2) + i%XVEGTYPE(:,nvt_v3) &
       +i%XVEGTYPE(:,nvt_v4) + i%XVEGTYPE(:,nvt_v5) + i%XVEGTYPE(:,nvt_v6) &
       +i%XVEGTYPE(:,nvt_v7) + i%XVEGTYPE(:,nvt_v8) + i%XVEGTYPE(:,nvt_v9) > 0. )
  !
  pfiso(:) = pisopot(:) *                   &
     ( ztcor(:,nvt_v1) * i%XVEGTYPE(:,nvt_v1) &
      +ztcor(:,nvt_v2) * i%XVEGTYPE(:,nvt_v2) &
      +ztcor(:,nvt_v3) * i%XVEGTYPE(:,nvt_v3) &
      +ztcor(:,nvt_v4) * i%XVEGTYPE(:,nvt_v4) &
      +ztcor(:,nvt_v5) * i%XVEGTYPE(:,nvt_v5) &
      +ztcor(:,nvt_v6) * i%XVEGTYPE(:,nvt_v6) &
      +ztcor(:,nvt_v7) * i%XVEGTYPE(:,nvt_v7) &
      +ztcor(:,nvt_v8) * i%XVEGTYPE(:,nvt_v8) &
      +ztcor(:,nvt_v9) * i%XVEGTYPE(:,nvt_v9) )
  !
  pfmono(:) = pmonopot(:) *                  &
     ( ztcorm(:,nvt_v1) * i%XVEGTYPE(:,nvt_v1) &
      +ztcorm(:,nvt_v2) * i%XVEGTYPE(:,nvt_v2) &
      +ztcorm(:,nvt_v3) * i%XVEGTYPE(:,nvt_v3) &
      +ztcorm(:,nvt_v4) * i%XVEGTYPE(:,nvt_v4) & 
      +ztcorm(:,nvt_v5) * i%XVEGTYPE(:,nvt_v5) & 
      +ztcorm(:,nvt_v6) * i%XVEGTYPE(:,nvt_v6) & 
      +ztcorm(:,nvt_v7) * i%XVEGTYPE(:,nvt_v7) & 
      +ztcorm(:,nvt_v8) * i%XVEGTYPE(:,nvt_v8) &
      +ztcorm(:,nvt_v9) * i%XVEGTYPE(:,nvt_v9) ) 
  !
ELSEWHERE
  !
  pfiso(:) = 0.
  pfmono(:) = 0.
  !
END WHERE
!
IF (lhook) CALL dr_hook('CH_BVOCEM_N:BY_VEG9',1,zhook_handle)
!
END SUBROUTINE by_veg9
!--------------------------------------------------------------------------
FUNCTION zlcor_func(ZX)

REAL, DIMENSION(:)          :: zx
REAL, DIMENSION(SIZE(ZX))   :: zlcor_func
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('CH_BVOCEM_N:ZLCOR_FUNC',0,zhook_handle)
zlcor_func(:)=0.
zlcor_func(:) = zx(:)*xiso_cl*xiso_alf/(1+(xiso_alf**2)*(zx(:)**2))**0.5
IF (lhook) CALL dr_hook('CH_BVOCEM_N:ZLCOR_FUNC',1,zhook_handle)
!
END FUNCTION zlcor_func
!---------------------------------------------------------------------------
FUNCTION ztcor0_func(ZX)

REAL, PARAMETER             :: r   = 8.314
REAL          :: zx
REAL   :: ztcor0_func
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('CH_BVOCEM_N:ZTCOR0_FUNC',0,zhook_handle)
!
ztcor0_func=0.
ztcor0_func = exp(xiso_ct1*(zx-xiso_bts)/(r*xiso_bts*zx))     &
          /(1+exp(xiso_ct2*(zx-xiso_btm)/(r*xiso_bts*zx)))
       !
IF (lhook) CALL dr_hook('CH_BVOCEM_N:ZTCOR0_FUNC',1,zhook_handle)
END FUNCTION ztcor0_func
!---------------------------------------------------------------------------
FUNCTION ztcorm0_func(ZX)

REAL           :: zx
REAL  :: ztcorm0_func
REAL(KIND=JPRB) :: zhook_handle

!      
IF (lhook) CALL dr_hook('CH_BVOCEM_N:ZTCORM0_FUNC',0,zhook_handle)
ztcorm0_func= 0.
ztcorm0_func = exp(xmono_beta*(zx-xmono_t3))
IF (lhook) CALL dr_hook('CH_BVOCEM_N:ZTCORM0_FUNC',1,zhook_handle)
!
END FUNCTION ztcorm0_func
!
!---------------------------------------------------------------------------
!
END SUBROUTINE ch_bvocem_n