ch_aer_emission.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_aer_emission(PFLUX, PRHODREF, HSV, KSV_CHSBEG,  PFCO)
!!   ############################################################
!!
!!    PURPOSE
!!    -------
!!    Transforme les emissions  d'aérosol en masse kg.kg-1.m.s-1 en molecules.m-2.s-1 : flux du moment m3
!!    Calcule les flux des moments m0 et m6 à partir de sigma  et Rg (um)
!!
!!    REFERENCE
!!    ---------
!!    none
!!
!!    AUTHOR
!!    ------
!!    Pierre TULET (CNRM/GMEI)
!!
!!    MODIFICATIONS
!!    -------------
!!    none
!!
!!    EXTERNAL
!!    --------
!!    None
!!
USE modd_chs_aerosol
USE modd_dst_surf, ONLY : xdensity_dst
USE modi_abor1_sfx
!!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!!
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
!*      0.1    declarations of arguments
!
REAL,   DIMENSION(:,:),      INTENT(INOUT) :: PFLUX
REAL,   DIMENSION(:),        INTENT(IN) :: PRHODREF
INTEGER,                     INTENT(IN) :: KSV_CHSBEG
 CHARACTER(LEN=*), DIMENSION(:),  INTENT(IN)  :: HSV      ! name of chemical species
REAL,   DIMENSION(:),OPTIONAL, INTENT(IN)  :: PFCO   ! CO flux 

!
!
!*      0.2    declarations local variables
!
REAL :: ZDEN2MOL
        !  ZDEN2MOL = 6.0221367E+23 * 1E-6 / 28.9644E-3
        !  conversion factor density to mol/cm3
        !  n_molec (moelc./cm3):  M = 1E-6*RHO(kg/m3) * XAVOGADRO / XMD
REAL,DIMENSION(NSP+NCARB+NSOA) :: ZFAC, ZRHOI, ZMI
REAL,DIMENSION(SIZE(PFLUX,1),NSP+NCARB+NSOA,JPMODE) :: ZFCTOTA
REAL,DIMENSION(SIZE(PFLUX,1),JPIN) :: ZFM
REAL,DIMENSION(SIZE(PFLUX,1)) :: ZFCO
REAL,DIMENSION(SIZE(PFLUX,1)) :: ZCONVERSION
!
INTEGER :: JJ, JSV  ! loop counter
REAL   :: ZEMISRADIUSI, ZEMISRADIUSJ
REAL   :: ZVALBC, ZVALOC
INTEGER :: I_CH_M0i, I_CH_M0j, I_CH_M6i, I_CH_M6j, I_CH_H2Oi, I_CH_H2Oj,&
                  I_CH_SO4i,I_CH_SO4j, I_CH_NO3i, I_CH_NO3j, I_CH_NH3i, I_CH_NH3j,&
                  I_CH_OCi, I_CH_OCj, I_CH_BCi, I_CH_BCj  , I_CH_DSTi, I_CH_DSTj   
REAL(KIND=JPRB) :: ZHOOK_HANDLE

!
!-------------------------------------------------------------------------------
!
!*       1.     TRANSFER FROM GAS TO AEROSOL MODULE
!               ------------------------------------
!        1.1    initialisation 
!
IF (lhook) CALL dr_hook('CH_AER_EMISSION',0,zhook_handle)

i_ch_m0i=-999
i_ch_m0j=-999
i_ch_m6i=-999
i_ch_m6j=-999
i_ch_h2oi=-999
i_ch_h2oj=-999
i_ch_so4i=-999
i_ch_so4j=-999
i_ch_no3i=-999
i_ch_no3j=-999
i_ch_nh3i=-999
i_ch_nh3j=-999
i_ch_oci=-999
i_ch_ocj=-999
i_ch_bci=-999
i_ch_bcj=-999
i_ch_dsti=-999
i_ch_dstj=-999

DO jsv=1, size(hsv)
  IF (trim(hsv(jsv)) == "M0I") i_ch_m0i=jsv-ksv_chsbeg+1
  IF (trim(hsv(jsv)) == "M0J") i_ch_m0j=jsv-ksv_chsbeg+1
  IF (trim(hsv(jsv)) == "M6I") i_ch_m6i=jsv-ksv_chsbeg+1
  IF (trim(hsv(jsv)) == "M6J") i_ch_m6j=jsv-ksv_chsbeg+1
  IF (trim(hsv(jsv)) == "H2OI") i_ch_h2oi=jsv-ksv_chsbeg+1
  IF (trim(hsv(jsv)) == "H2OJ") i_ch_h2oj=jsv-ksv_chsbeg+1
  IF (trim(hsv(jsv)) == "SO4I") i_ch_so4i=jsv-ksv_chsbeg+1
  IF (trim(hsv(jsv)) == "SO4J") i_ch_so4j=jsv-ksv_chsbeg+1
  IF (trim(hsv(jsv)) == "NO3I") i_ch_no3i=jsv-ksv_chsbeg+1
  IF (trim(hsv(jsv)) == "NO3J") i_ch_no3j=jsv-ksv_chsbeg+1
  IF (trim(hsv(jsv)) == "NH3I") i_ch_nh3i=jsv-ksv_chsbeg+1
  IF (trim(hsv(jsv)) == "NH3J") i_ch_nh3j=jsv-ksv_chsbeg+1
  IF (trim(hsv(jsv)) == "OCI") i_ch_oci=jsv-ksv_chsbeg+1
  IF (trim(hsv(jsv)) == "OCJ") i_ch_ocj=jsv-ksv_chsbeg+1
  IF (trim(hsv(jsv)) == "BCI") i_ch_bci=jsv-ksv_chsbeg+1
  IF (trim(hsv(jsv)) == "BCJ") i_ch_bcj=jsv-ksv_chsbeg+1
  IF (trim(hsv(jsv)) == "DSTI") i_ch_dsti=jsv-ksv_chsbeg+1
  IF (trim(hsv(jsv)) == "DSTJ") i_ch_dstj=jsv-ksv_chsbeg+1
END DO

IF (i_ch_m0i ==-999) CALL abor1_sfx ('WRONG VALUE FOR I_CH_M0i ')
IF (i_ch_m0j ==-999) CALL abor1_sfx ('WRONG VALUE FOR I_CH_M0j ')
IF (i_ch_m6i ==-999 .AND. lvarsigi) CALL abor1_sfx ('WRONG VALUE FOR I_CH_M6i ')
IF (i_ch_m6j ==-999 .AND. lvarsigj) CALL abor1_sfx ('WRONG VALUE FOR I_CH_M6j ')
IF (i_ch_h2oi==-999) CALL abor1_sfx ('WRONG VALUE FOR I_CH_H2Oi')
IF (i_ch_h2oj==-999) CALL abor1_sfx ('WRONG VALUE FOR I_CH_H2Oj')
IF (i_ch_so4i==-999) CALL abor1_sfx ('WRONG VALUE FOR I_CH_SO4i')
IF (i_ch_so4j==-999) CALL abor1_sfx ('WRONG VALUE FOR I_CH_SO4j')
IF (i_ch_no3i==-999) CALL abor1_sfx ('WRONG VALUE FOR I_CH_NO3i')
IF (i_ch_no3j==-999) CALL abor1_sfx ('WRONG VALUE FOR I_CH_NO3j')
IF (i_ch_nh3i==-999) CALL abor1_sfx ('WRONG VALUE FOR I_CH_NH3i')
IF (i_ch_nh3j==-999) CALL abor1_sfx ('WRONG VALUE FOR I_CH_NH3j')
IF (i_ch_oci ==-999) CALL abor1_sfx ('WRONG VALUE FOR I_CH_OCi ')
IF (i_ch_ocj ==-999) CALL abor1_sfx ('WRONG VALUE FOR I_CH_OCj ')
IF (i_ch_bci ==-999) CALL abor1_sfx ('WRONG VALUE FOR I_CH_BCi ')
IF (i_ch_bcj ==-999) CALL abor1_sfx ('WRONG VALUE FOR I_CH_BCj ')
IF (i_ch_dsti ==-999) CALL abor1_sfx ('WRONG VALUE FOR I_CH_DSTi ')
IF (i_ch_dstj ==-999) CALL abor1_sfx ('WRONG VALUE FOR I_CH_DSTj ')

zmi(:) = 250.
zmi(jp_aer_so4) = 98.
zmi(jp_aer_no3) = 63.
zmi(jp_aer_nh3) = 17.
zmi(jp_aer_h2o) = 18.
zmi(jp_aer_dst) = 100.

! Aerosol Density
! Cf Ackermann (all to black carbon except water)
zrhoi(:) = 1.8e3
zrhoi(jp_aer_h2o) = 1.0e3   ! water
zrhoi(jp_aer_dst) = xdensity_dst

zden2mol = 1e-6 * xavogadro / xmd

IF (crgunit=="MASS") THEN
zemisradiusi = xemisradiusi * exp(-3.*(log(xemissigi))**2)
zemisradiusj = xemisradiusj * exp(-3.*(log(xemissigj))**2)
ELSE
zemisradiusi = xemisradiusi
zemisradiusj = xemisradiusj
END IF
!
DO jj=1,nsp+ncarb+nsoa
  zfac(jj) = (4./3.)*3.14292654*zrhoi(jj)*1.e-9
ENDDO

zfm(:,:)=0.
zfctota(:,:,:)=0.
pflux(:,:) = max(pflux(:,:),0.)
!
!*       0      conversion into  kg.kg-1.m.s-1 (due to PCONVERSION)

zvalbc  = 0.
zvaloc  = 0.
zfco(:) = 0.
IF ((lco2pm).AND.(PRESENT(pfco))) THEN
  !ZVALBC=2.748549E-02  ! CO / BC conversion factor
  !ZVALOC=4.947248E-02  ! CO / POM conversion factor
  zvalbc= 5.* 0.6e-10 / 0.4e-8  ! CO / BC conversion factor
  zvaloc= 5.* 0.3e-10 / 0.4e-8  ! CO / POM conversion factor
  zfco(:) = pfco(:)
END IF

! Initial aerosols fluxes have been transformed into molecu.m-2.s-1, 
! conversion into are in kg.kg-1.m.s-1 
!  conversion in kg.kg-1.m.s-1
!
zconversion(:) =  xavogadro * prhodref(:)
!
pflux(:,i_ch_so4i) = pflux(:,i_ch_so4i) / zconversion(:) * zmi(jp_aer_so4)*1e-3
pflux(:,i_ch_so4j) = pflux(:,i_ch_so4j) / zconversion(:) * zmi(jp_aer_so4)*1e-3
pflux(:,i_ch_no3i) = pflux(:,i_ch_no3i) / zconversion(:) * zmi(jp_aer_no3)*1e-3
pflux(:,i_ch_no3j) = pflux(:,i_ch_no3j) / zconversion(:) * zmi(jp_aer_no3)*1e-3
pflux(:,i_ch_nh3i) = pflux(:,i_ch_nh3i) / zconversion(:) * zmi(jp_aer_nh3)*1e-3
pflux(:,i_ch_nh3j) = pflux(:,i_ch_nh3j) / zconversion(:) * zmi(jp_aer_nh3)*1e-3
pflux(:,i_ch_h2oi) = pflux(:,i_ch_h2oi) / zconversion(:) * zmi(jp_aer_h2o)*1e-3
pflux(:,i_ch_h2oj) = pflux(:,i_ch_h2oj) / zconversion(:) * zmi(jp_aer_h2o)*1e-3
pflux(:,i_ch_oci)  = (pflux(:,i_ch_oci) + zfco(:) * zvaloc / 2.) / zconversion(:) * zmi(jp_aer_oc)*1e-3
pflux(:,i_ch_ocj)  = (pflux(:,i_ch_ocj) + zfco(:) * zvaloc ) / zconversion(:) * zmi(jp_aer_oc)*1e-3
pflux(:,i_ch_bci)  = (pflux(:,i_ch_bci) + zfco(:) * zvalbc / 2.) / zconversion(:) * zmi(jp_aer_bc)*1e-3
pflux(:,i_ch_bcj)  = (pflux(:,i_ch_bcj) + zfco(:) * zvalbc ) / zconversion(:) * zmi(jp_aer_bc)*1e-3
pflux(:,i_ch_dsti) = pflux(:,i_ch_dsti) / zconversion(:) * zmi(jp_aer_dst)*1e-3 
pflux(:,i_ch_dstj) = pflux(:,i_ch_dstj) / zconversion(:) * zmi(jp_aer_dst)*1e-3 
!
!*       1.0    transfer aerosol mass from gas to aerosol variables
!               (and conversion of kg.kg-1.m.s-1 --> microgram.m-2.s-1)
!
zfctota(:,jp_aer_so4,1) = pflux(:,i_ch_so4i) *1e+9 * prhodref(:) 
zfctota(:,jp_aer_so4,2) = pflux(:,i_ch_so4j) *1e+9 * prhodref(:)

zfctota(:,jp_aer_nh3,1) = pflux(:,i_ch_nh3i) *1e+9 * prhodref(:)
zfctota(:,jp_aer_nh3,2) = pflux(:,i_ch_nh3j) *1e+9 * prhodref(:)

zfctota(:,jp_aer_no3,1) = pflux(:,i_ch_no3i) *1e+9 * prhodref(:)
zfctota(:,jp_aer_no3,2) = pflux(:,i_ch_no3j) *1e+9 * prhodref(:)

zfctota(:,jp_aer_h2o,1) = pflux(:,i_ch_h2oi) *1e+9 * prhodref(:)
zfctota(:,jp_aer_h2o,2) = pflux(:,i_ch_h2oj) *1e+9 * prhodref(:)

zfctota(:,jp_aer_oc,1)  = pflux(:,i_ch_oci)  *1e+9 * prhodref(:)
zfctota(:,jp_aer_oc,2)  = pflux(:,i_ch_ocj)  *1e+9 * prhodref(:)

zfctota(:,jp_aer_bc,1)  = pflux(:,i_ch_bci)  *1e+9 * prhodref(:)
zfctota(:,jp_aer_bc,2)  = pflux(:,i_ch_bcj)  *1e+9 * prhodref(:)

zfctota(:,jp_aer_dst,1)  = pflux(:,i_ch_dsti)  *1e+9 * prhodref(:)
zfctota(:,jp_aer_dst,2)  = pflux(:,i_ch_dstj)  *1e+9 * prhodref(:)
!
!*       1.1    calculate moment 3 flux from total aerosol mass
!
zfm(:,2) = 0.
zfm(:,5) = 0.
DO jj = 1,nsp+ncarb+nsoa
  zfm(:,2) = zfm(:,2)+zfctota(:,jj,1)/zfac(jj)
  zfm(:,5) = zfm(:,5)+zfctota(:,jj,2)/zfac(jj)
ENDDO
!
!*       1.2    calculate moment 0 flux from dispersion and mean radius Rg
!
zfm(:,1)= zfm(:,2) / ((zemisradiusi**3)*exp(4.5 * (log(xemissigi))**2)) 
!
zfm(:,4)= zfm(:,5) / ((zemisradiusj**3)*exp(4.5 * (log(xemissigj))**2)) 
!
!*       1.3    calculate moment 6 flux from dispersion and mean diameter
!
zfm(:,3) = zfm(:,1) * (zemisradiusi**6) *exp(18 *(log(xemissigi))**2)
!
zfm(:,6) = zfm(:,4) * (zemisradiusj**6) *exp(18 *(log(xemissigj))**2)
!
!*       1.4    conversion en ppp.m.s-1
!
! conversion in atmospheric unit only for moments 0 and 6 
pflux(:,i_ch_m0i) = zfm(:,1) * 1e-6 / (zden2mol * prhodref(:))
pflux(:,i_ch_m0j) = zfm(:,4) * 1e-6 / (zden2mol * prhodref(:))
!
IF (lvarsigi) pflux(:,i_ch_m6i) = zfm(:,3) / (zden2mol * prhodref(:))
IF (lvarsigj) pflux(:,i_ch_m6j) = zfm(:,6) / (zden2mol * prhodref(:))
!
! aerosol phase conversion kg/kg.m.s-1 into ppp.m.s-1
pflux(:,i_ch_so4i) = pflux(:,i_ch_so4i) * xmd / (zmi(jp_aer_so4)*1e-3)
pflux(:,i_ch_so4j) = pflux(:,i_ch_so4j) * xmd / (zmi(jp_aer_so4)*1e-3)
pflux(:,i_ch_no3i) = pflux(:,i_ch_no3i) * xmd / (zmi(jp_aer_no3)*1e-3)
pflux(:,i_ch_no3j) = pflux(:,i_ch_no3j) * xmd / (zmi(jp_aer_no3)*1e-3)
pflux(:,i_ch_nh3i) = pflux(:,i_ch_nh3i) * xmd / (zmi(jp_aer_nh3)*1e-3)
pflux(:,i_ch_nh3j) = pflux(:,i_ch_nh3j) * xmd / (zmi(jp_aer_nh3)*1e-3)
pflux(:,i_ch_h2oi) = pflux(:,i_ch_h2oi) * xmd / (zmi(jp_aer_h2o)*1e-3)
pflux(:,i_ch_h2oj) = pflux(:,i_ch_h2oj) * xmd / (zmi(jp_aer_h2o)*1e-3)
!
pflux(:,i_ch_oci) = pflux(:,i_ch_oci) * xmd / (zmi(jp_aer_oc)*1e-3)
pflux(:,i_ch_ocj) = pflux(:,i_ch_ocj) * xmd / (zmi(jp_aer_oc)*1e-3)
pflux(:,i_ch_bci) = pflux(:,i_ch_bci) * xmd / (zmi(jp_aer_bc)*1e-3)
pflux(:,i_ch_bcj) = pflux(:,i_ch_bcj) * xmd / (zmi(jp_aer_bc)*1e-3)
pflux(:,i_ch_dsti) = pflux(:,i_ch_dsti) * xmd / (zmi(jp_aer_dst)*1e-3)
pflux(:,i_ch_dstj) = pflux(:,i_ch_dstj) * xmd / (zmi(jp_aer_dst)*1e-3)
!
! then conversion in  molecules.m-2.s-1
pflux(:,i_ch_m0i) = pflux(:,i_ch_m0i) * zconversion(:)
pflux(:,i_ch_m0j) = pflux(:,i_ch_m0j) * zconversion(:)
!
IF (lvarsigi) pflux(:,i_ch_m6i) = pflux(:,i_ch_m6i) * zconversion(:) 
IF (lvarsigj) pflux(:,i_ch_m6j) = pflux(:,i_ch_m6j) * zconversion(:)
!
pflux(:,i_ch_so4i) = pflux(:,i_ch_so4i) * zconversion(:) / (zmi(jp_aer_so4)*1e-3)
pflux(:,i_ch_so4j) = pflux(:,i_ch_so4j) * zconversion(:) / (zmi(jp_aer_so4)*1e-3)
pflux(:,i_ch_no3i) = pflux(:,i_ch_no3i) * zconversion(:) / (zmi(jp_aer_no3)*1e-3)
pflux(:,i_ch_no3j) = pflux(:,i_ch_no3j) * zconversion(:) / (zmi(jp_aer_no3)*1e-3)
pflux(:,i_ch_nh3i) = pflux(:,i_ch_nh3i) * zconversion(:) / (zmi(jp_aer_nh3)*1e-3)
pflux(:,i_ch_nh3j) = pflux(:,i_ch_nh3j) * zconversion(:) / (zmi(jp_aer_nh3)*1e-3)
pflux(:,i_ch_h2oi) = pflux(:,i_ch_h2oi) * zconversion(:) / (zmi(jp_aer_h2o)*1e-3)
pflux(:,i_ch_h2oj) = pflux(:,i_ch_h2oj) * zconversion(:) / (zmi(jp_aer_h2o)*1e-3)
pflux(:,i_ch_oci) = pflux(:,i_ch_oci)   * zconversion(:) / (zmi(jp_aer_oc)*1e-3)
pflux(:,i_ch_ocj) = pflux(:,i_ch_ocj)   * zconversion(:) / (zmi(jp_aer_oc)*1e-3)
pflux(:,i_ch_bci) = pflux(:,i_ch_bci)   * zconversion(:) / (zmi(jp_aer_bc)*1e-3)
pflux(:,i_ch_bcj) = pflux(:,i_ch_bcj)   * zconversion(:) / (zmi(jp_aer_bc)*1e-3)
pflux(:,i_ch_dsti) = pflux(:,i_ch_dsti) * zconversion(:) / (zmi(jp_aer_dst)*1e-3)
pflux(:,i_ch_dstj) = pflux(:,i_ch_dstj) * zconversion(:) / (zmi(jp_aer_dst)*1e-3)
!
IF (lhook) CALL dr_hook('CH_AER_EMISSION',1,zhook_handle)
!
END SUBROUTINE ch_aer_emission