mode_aer_surf.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.
!!   ########################
MODULE mode_aer_surf
!!   ########################
!!
!! MODULE DUST PSD (Particle Size Distribution)
!! Purpose: Contains subroutines to convert from transported variables (ppp)
!! to understandable aerosol variables, e.g. #/m3, kg/m3, sigma, R_{n}
!-------------------------------------------------------------------------------
!!    MODIFICATIONS
!!    -------------
!!
!!      J.Escobar     06/2013  for REAL4/8 add EPSILON management
!!    M.Leriche 2015 : masse molaire Black carbon à 12 g/mol
!!
!------------------------------------------------------------------------------- 
  USE modd_chs_aerosol
  USE modd_dst_surf, ONLY : xdensity_dst
!
  USE yomhook   ,ONLY : lhook,   dr_hook
  USE parkind1  ,ONLY : jprb
  USE modd_surf_par , ONLY : xsurf_tiny
!
  IMPLICIT NONE
!!
CONTAINS
!!
SUBROUTINE init_var(PSV,PFAC,PCTOTA)
!
IMPLICIT NONE
!
REAL,DIMENSION(:,:), INTENT(IN)    :: PSV      ! [aerosol concentration]
REAL,DIMENSION(:), INTENT(OUT)     :: PFAC     ! M3 / mass conversion factor
REAL,DIMENSION(:,:,:), INTENT(OUT) :: PCTOTA
!
REAL,DIMENSION(NSP+NCARB+NSOA)                     :: ZMI      ! [kg/mol] molar weight of aerosol
REAL,DIMENSION(NSP+NCARB+NSOA)                     :: ZRHOI    ! aerosol density
REAL                                               :: ZPI
INTEGER                                            :: JJ
REAL, PARAMETER                                    :: ZMOL = 6.0221367e+11
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (lhook) CALL dr_hook('MODE_AER_SURF:INIT_VAR',0,zhook_handle)
!
zrhoi(:) = 1.8e3
zrhoi(jp_aer_h2o) = 1.0e3   ! water
zrhoi(jp_aer_dst) = xdensity_dst
!
! Moments index
nm0(1) = 1 
nm3(1) = 2
nm6(1) = 3 
nm0(2) = 4 
nm3(2) = 5 
nm6(2) = 6 
!
!Set density of aerosol, here  (kg/m3)
! Aerosol Density
! Cf Ackermann (all to black carbon except water)
!
!Set molecular weightn g/mol 
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_bc)   = 12.
zmi(jp_aer_dst)  = 100.
IF (nsoa .EQ. 10) THEN
  zmi(jp_aer_soa1) = 88. 
  zmi(jp_aer_soa2) = 180.
  zmi(jp_aer_soa3) = 1.5374857e3
  zmi(jp_aer_soa4) = 1.9586780e3
  zmi(jp_aer_soa5) = 195.
  zmi(jp_aer_soa6) = 195.
  zmi(jp_aer_soa7) = 165.
  zmi(jp_aer_soa8) = 195.
  zmi(jp_aer_soa9) = 270.
  zmi(jp_aer_soa10) = 210.
END IF
!
zpi = 2.*asin(1.)
DO jj=1,nsp+ncarb+nsoa
  pfac(jj)=(4./3.)*zpi*zrhoi(jj)*1.e-9
ENDDO
!
!*       2    transfer aerosol mass from gas to aerosol variables
!               (and conversion of mol.cm-3 --> microgram/m3)
!
pctota(:,:,:) = 0.
! aerosol phase
pctota(:,jp_aer_so4,1) = psv(:,jp_ch_so4i)*zmi(jp_aer_so4)/zmol
pctota(:,jp_aer_so4,2) = psv(:,jp_ch_so4j)*zmi(jp_aer_so4)/zmol

pctota(:,jp_aer_no3,1) = psv(:,jp_ch_no3i)*zmi(jp_aer_no3)/zmol
pctota(:,jp_aer_no3,2) = psv(:,jp_ch_no3j)*zmi(jp_aer_no3)/zmol

pctota(:,jp_aer_nh3,1) = psv(:,jp_ch_nh3i)*zmi(jp_aer_nh3)/zmol
pctota(:,jp_aer_nh3,2) = psv(:,jp_ch_nh3j)*zmi(jp_aer_nh3)/zmol
!
! water
pctota(:,jp_aer_h2o,1) = psv(:,jp_ch_h2oi)*zmi(jp_aer_h2o)/zmol
pctota(:,jp_aer_h2o,2) = psv(:,jp_ch_h2oj)*zmi(jp_aer_h2o)/zmol
!
! primary organic carbon
pctota(:,jp_aer_oc,1) = psv(:,jp_ch_oci)*zmi(jp_aer_oc)/zmol
pctota(:,jp_aer_oc,2) = psv(:,jp_ch_ocj)*zmi(jp_aer_oc)/zmol
!
! primary black carbon
pctota(:,jp_aer_bc,1) = psv(:,jp_ch_bci)*zmi(jp_aer_bc)/zmol
pctota(:,jp_aer_bc,2) = psv(:,jp_ch_bcj)*zmi(jp_aer_bc)/zmol
!
!dust
pctota(:,jp_aer_dst,1) = psv(:,jp_ch_dsti)*zmi(jp_aer_dst)/6.0221367e+11
pctota(:,jp_aer_dst,2) = psv(:,jp_ch_dstj)*zmi(jp_aer_dst)/6.0221367e+11
!
IF (nsoa .EQ. 10) THEN
  pctota(:,jp_aer_soa1,1) = psv(:,jp_ch_soa1i)*zmi(jp_aer_soa1)/zmol
  pctota(:,jp_aer_soa1,2) = psv(:,jp_ch_soa1j)*zmi(jp_aer_soa1)/zmol
  pctota(:,jp_aer_soa2,1) = psv(:,jp_ch_soa2i)*zmi(jp_aer_soa2)/zmol
  pctota(:,jp_aer_soa2,2) = psv(:,jp_ch_soa2j)*zmi(jp_aer_soa2)/zmol
  pctota(:,jp_aer_soa3,1) = psv(:,jp_ch_soa3i)*zmi(jp_aer_soa3)/zmol
  pctota(:,jp_aer_soa3,2) = psv(:,jp_ch_soa3j)*zmi(jp_aer_soa3)/zmol
  pctota(:,jp_aer_soa4,1) = psv(:,jp_ch_soa4i)*zmi(jp_aer_soa4)/zmol
  pctota(:,jp_aer_soa4,2) = psv(:,jp_ch_soa4j)*zmi(jp_aer_soa4)/zmol
  pctota(:,jp_aer_soa5,1) = psv(:,jp_ch_soa5i)*zmi(jp_aer_soa5)/zmol
  pctota(:,jp_aer_soa5,2) = psv(:,jp_ch_soa5j)*zmi(jp_aer_soa5)/zmol

  pctota(:,jp_aer_soa6,1) = psv(:,jp_ch_soa6i)*zmi(jp_aer_soa6)/zmol
  pctota(:,jp_aer_soa6,2) = psv(:,jp_ch_soa6j)*zmi(jp_aer_soa6)/zmol
  pctota(:,jp_aer_soa7,1) = psv(:,jp_ch_soa7i)*zmi(jp_aer_soa7)/zmol
  pctota(:,jp_aer_soa7,2) = psv(:,jp_ch_soa7j)*zmi(jp_aer_soa7)/zmol
  pctota(:,jp_aer_soa8,1) = psv(:,jp_ch_soa8i)*zmi(jp_aer_soa8)/zmol
  pctota(:,jp_aer_soa8,2) = psv(:,jp_ch_soa8j)*zmi(jp_aer_soa8)/zmol
  pctota(:,jp_aer_soa9,1) = psv(:,jp_ch_soa9i)*zmi(jp_aer_soa9)/zmol
  pctota(:,jp_aer_soa9,2) = psv(:,jp_ch_soa9j)*zmi(jp_aer_soa9)/zmol
  pctota(:,jp_aer_soa10,1) = psv(:,jp_ch_soa10i)*zmi(jp_aer_soa10)/zmol
  pctota(:,jp_aer_soa10,2) = psv(:,jp_ch_soa10j)*zmi(jp_aer_soa10)/zmol
END IF
!
IF (lhook) CALL dr_hook('MODE_AER_SURF:INIT_VAR',1,zhook_handle)
!
END SUBROUTINE init_var
!
!!   ############################################################
SUBROUTINE ppp2aero_surf(             &
          PSVT,                       &!I [ppp] input scalar variables (moment of distribution)
          prhodref,                   &!I [kg/m3] density of air       
          psig1d,                     &!O [-] standard deviation of aerosol distribution
          prg1d,                      &!O [um] number median diameter of aerosol distribution
          pn1d,                       &!O [#/m3] number concentration of aerosols
          pctota,                     &!O [ug/m3] mass of each aerosol compounds
          pm1d                        &!moments 0, 3 and 6
         )  
!!   ############################################################
!
!!
!!    PURPOSE
!!    -------
!!    Translate the three moments M0, M3 and M6 given in ppp into
!!    Values which can be understood more easily (R, sigma, N, M)
!! 
!!    CALLING STRUCTURE NOTE: OPTIONAL VARIABLES
!!    -------
!!    CALL PPP2AERO_SURFS(PSVT, PRHODREF, PSIG1D=SIGVAR,  &
!!       PRG1D=RVAR, PN1D=NVAR, PM1D=ZM)
!!
!!    REFERENCE
!!    ---------
!!    none
!!
!!    AUTHOR
!!    ------
!!    Pierre TULET (LA)
!!
!!    MODIFICATIONS
!!    -------------
!!    Alf Grini (CNRM)
!!
!!    EXTERNAL
!!    --------
!!
    IMPLICIT NONE
!!
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
!*      0.1    declarations of arguments
!
!INPUT
REAL,       DIMENSION(:,:),  INTENT(IN)     :: PSVT      !I [#/molec_{air}] first moment
                                                         !I [molec_{aer}/molec_{air} 3rd moment
                                                         !I [um6/molec_{air}*(cm3/m3)] 6th moment
REAL,       DIMENSION(:),    INTENT(IN)      :: PRHODREF !I [kg/m3] density of air

!OUTPUT
REAL,       DIMENSION(:,:),  OPTIONAL, INTENT(OUT)     :: PSIG1D   !O [-] standard deviation
REAL,       DIMENSION(:,:),  OPTIONAL, INTENT(OUT)     :: PRG1D    !O [um] number median diameter
REAL,       DIMENSION(:,:),  OPTIONAL, INTENT(OUT)     :: PN1D     !O [#/m3] number concentration
REAL,       DIMENSION(:,:,:),OPTIONAL, INTENT(OUT)     :: PCTOTA   !O [ug/m3] mass of each component
REAL,       DIMENSION(:,:),  OPTIONAL, INTENT(OUT)     :: PM1D     !O moments 0,3 and 6 
!
!*      0.2    declarations local variables
!    
REAL,DIMENSION(SIZE(PSVT,1), SIZE(PSVT,2))         :: ZSV      ! [aerosol concentration]
REAL,DIMENSION(SIZE(PSVT,1))                       :: ZSIGMA   ! [-] standard deviation
REAL,DIMENSION(SIZE(PSVT,1),NSP+NCARB+NSOA,JPMODE) :: ZCTOTA
REAL,DIMENSION(SIZE(PSVT,1),JPMODE*3)              :: ZM
REAL,DIMENSION(JPMODE*3)                           :: ZMMIN
!
REAL,DIMENSION(NSP+NCARB+NSOA)                     :: ZFAC     ! M3 / mass conversion factor
REAL, PARAMETER                                    :: ZDEN2MOL = 1e-6 * 6.0221367e+23  / 28.9644e-3
INTEGER                                            :: JJ, JN   ! [idx] loop counters
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_AER_SURF:PPP2AERO_SURF',0,zhook_handle)
!
!        1.    initialisation
!
DO jj=1, SIZE(psvt,2)
  zsv(:,jj) =  psvt(:,jj) * zden2mol * prhodref(:)
  zsv(:,jj) = max(zsv(:,jj),1e-40 * zden2mol * prhodref(:))
ENDDO
!
 CALL init_var(zsv,zfac,zctota)
!
!*       2    calculate moment 3 from total aerosol mass
!
zm(:,2) = 0.
zm(:,5) = 0.
DO jj = 1,nsp+ncarb+nsoa
  zm(:,2) = zm(:,2)+zctota(:,jj,1)/zfac(jj) !==>um3_{aer}/m3_{air} (volume ==> 3rd moment)
  zm(:,5) = zm(:,5)+zctota(:,jj,2)/zfac(jj) !==>um3_{aer}/m3_{air} (volume ==> 3rd moment)
ENDDO
!
!
!*       3    set  moment 0 
!
zm(:,1)=   max(zsv(:,jp_ch_m0i) * 1e+6, xsurf_tiny) ! molec_{aer}/m3_{air}
zm(:,4)=   max(zsv(:,jp_ch_m0j) * 1e+6, xsurf_tiny) ! molec_{aer}/m3_{air}
WHERE ((zm(:,1) .LT. zmmin(1))) !.OR.(ZM(:,2) .LT. ZMMIN(2)))
  zm(:,1)= zmmin(1)
!  ZM(:,2)= ZMMIN(2)
!
!  ZCTOTA(:,JP_AER_H2O,1) = 0.
!  ZCTOTA(:,JP_AER_NH3,1) = 0.
!  ZCTOTA(:,JP_AER_SO4,1) = 0.
!  ZCTOTA(:,JP_AER_NO3,1) = 0.
!  ZCTOTA(:,JP_AER_BC,1) = 0.5 * ZM(:,2) * ZFAC(JP_AER_BC)
!  ZCTOTA(:,JP_AER_OC,1) = 0.5 * ZM(:,2) * ZFAC(JP_AER_OC)
END WHERE
!!
WHERE ((zm(:,4) .LT. zmmin(4))) !.OR.(ZM(:,5) .LT. ZMMIN(5)))  
  zm(:,4)= zmmin(4)
!  ZM(:,5)= ZMMIN(5)
!
!  ZCTOTA(:,JP_AER_H2O,2) = 0.
!  ZCTOTA(:,JP_AER_NH3,2) = 0.
!  ZCTOTA(:,JP_AER_SO4,2) = 0.
!  ZCTOTA(:,JP_AER_NO3,2) = 0.
!  ZCTOTA(:,JP_AER_BC,2) = 0.5 * ZM(:,5) * ZFAC(JP_AER_BC)
!  ZCTOTA(:,JP_AER_OC,2) = 0.5 * ZM(:,5) * ZFAC(JP_AER_OC)
END WHERE
!
!*       4    set moment 6  ==> um6_{aer}/m3_{air}
!
IF (lvarsigi) THEN ! set M6 variable standard deviation
  zm(:,3) = max(zsv(:,jp_ch_m6i), xsurf_tiny)

  zsigma(:)=zm(:,2)**2/(zm(:,1)*zm(:,3))
  zsigma(:)=min(1-1e-10,zsigma(:))
  zsigma(:)=max(1e-10,zsigma(:))
  zsigma(:)= log(zsigma(:))
  zsigma(:)= exp(1./3.*sqrt(-zsigma(:)))
  zm(:,3) = zm(:,1) &
            * ( (zm(:,2)/zm(:,1))**(1./3.)  &
            * exp(-(3./2.)*log(zsigma(:))**2))**6 &
            * exp(18.*log(zsigma(:))**2)  

  IF(PRESENT(psig1d)) psig1d(:,1) = zsigma(:)

ELSE ! fixed standard deviation
  zm(:,3) = zm(:,1) &
            * ( (zm(:,2)/zm(:,1))**(1./3.)  &
            * exp(-(3./2.)*log(xemissigi)**2))**6 &
            * exp(18.*log(xemissigi)**2) 

  IF(PRESENT(psig1d)) psig1d(:,1) = xemissigi
END IF

IF (lvarsigj) THEN ! set M6 variable standard deviation
  zm(:,6) = max(zsv(:,jp_ch_m6j), xsurf_tiny)

  zsigma(:)=zm(:,5)**2/(zm(:,4)*zm(:,6))
  zsigma(:)=min(1-1e-10,zsigma(:))
  zsigma(:)=max(1e-10,zsigma(:))
  zsigma(:)= log(zsigma(:))
  zsigma(:)= exp(1./3.*sqrt(-zsigma(:)))

  zm(:,6) = zm(:,4) &
            * ( (zm(:,5)/zm(:,4))**(1./3.)  &
            * exp(-(3./2.)*log(zsigma(:))**2))**6 &
            * exp(18.*log(zsigma(:))**2)  

  IF(PRESENT(psig1d)) psig1d(:,2) = zsigma(:)

ELSE ! fixed standard deviation
  zm(:,6) = zm(:,4) &
            * ( (zm(:,5)/zm(:,4))**(1./3.)  &
            * exp(-(3./2.)*log(xemissigj)**2))**6 &
            * exp(18.*log(xemissigj)**2)  

  IF(PRESENT(psig1d)) psig1d(:,2) = xemissigj
END IF


!*       6    calculate modal parameters from moments
DO jn=1,jpmode
  IF(PRESENT(pn1d)) pn1d(:,jn) = zm(:,nm0(jn))

  IF(PRESENT(prg1d)) prg1d(:,jn)=(zm(:,nm3(jn))**4. &
                                 / (zm(:,nm6(jn))*zm(:,nm0(jn))**3.))**(1./6.)  
ENDDO

IF(PRESENT(pctota)) pctota(:,:,:) = zctota(:,:,:)
IF(PRESENT(pm1d)) pm1d(:,:) = zm(:,:)
IF (lhook) CALL dr_hook('MODE_AER_SURF:PPP2AERO_SURF',1,zhook_handle)
!
END SUBROUTINE ppp2aero_surf
!!   ############################################################
SUBROUTINE aero2ppp_surf(    &
            PSVT,            &!IO [ppp] input scalar variables (moment of distribution)
            prhodref,        &!I [kg/m3] density of air       
            psig1d,          &!I [-] standard deviation of aerosol distribution
            prg1d            &!I [um] number median diameter of aerosol distribution
            )  
!!   ############################################################
!
!!
!!    PURPOSE
!!    -------
!!    Translate the aerosol Mass, RG and SIGMA in the  three moments M0, M3 and M6 given in ppp 
!!
!!    REFERENCE
!!    ---------
!!    none
!!
!!    AUTHOR
!!    ------
!!    Pierre TULET (LA)
!!
!!    MODIFICATIONS
!!    -------------
!!    Alf Grini (CNRM)
!!
!!    EXTERNAL
!!    --------
!!
IMPLICIT NONE
!!
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
!*      0.1    declarations of arguments
!
!INPUT
REAL,       DIMENSION(:,:),  INTENT(INOUT)  :: PSVT      !I [#/molec_{air}] first moment
                                                         !I [molec_{aer}/molec_{air} 3rd moment
                                                         !I [um6/molec_{air}*(cm3/m3)] 6th moment
REAL,       DIMENSION(:),  INTENT(IN)       :: PRHODREF  !I [kg/m3] density of air

!OUTPUT
REAL,       DIMENSION(:,:),  INTENT(IN)     :: PSIG1D   !O [-] standard deviation
REAL,       DIMENSION(:,:),  INTENT(IN)     :: PRG1D    !O [um] number median diameter
!
!*      0.2    declarations local variables
!
REAL,DIMENSION(SIZE(PSVT,1),NSP+NCARB+NSOA,JPMODE):: ZCTOTA
REAL,DIMENSION(SIZE(PSVT,1), JPMODE*3) :: ZM                  ! [aerosol units] local array which goes to output later
!
REAL,DIMENSION(NSP+NCARB+NSOA)       :: ZFAC                ! M3 / mass conversion factor
REAL, PARAMETER                      :: ZDEN2MOL = 1e-6 * 6.0221367e+23  / 28.9644e-3
INTEGER                              :: JJ                  ! [idx] loop counters
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_AER_SURF:AERO2PPP_SURF',0,zhook_handle)
!
!        1.    initialisation 
!
DO jj=1, SIZE(psvt, 2)
  psvt(:,jj) =  psvt(:,jj) * zden2mol * prhodref(:)
ENDDO
!
 CALL init_var(psvt,zfac,zctota)
!
!*       3    calculate moment 3 from total aerosol mass
!
zm(:,2) = 0.
zm(:,5) = 0.
DO jj = 1,nsp+ncarb+nsoa
  zm(:,2) = zm(:,2)+zctota(:,jj,1)/zfac(jj)
  zm(:,5) = zm(:,5)+zctota(:,jj,2)/zfac(jj)
ENDDO
!
!
!*       4    calculate moment 0 from dispersion and mean radius
!
zm(:,1)= zm(:,2)/ ( (prg1d(:,1)**3)*exp(4.5 * log(psig1d(:,1))**2) )  
zm(:,4)= zm(:,5)/ ( (prg1d(:,2)**3)*exp(4.5 * log(psig1d(:,2))**2) )  
!
!*       5    calculate moment 6 from dispersion and mean radius
!
zm(:,3) = zm(:,1)*(prg1d(:,1)**6) * exp(18 *(log(psig1d(:,1)))**2)  
zm(:,6) = zm(:,4)*(prg1d(:,2)**6) * exp(18 *(log(psig1d(:,2)))**2)  
!
!*       6    return to ppp
!
psvt(:,jp_ch_m0i) = zm(:,1) * 1e-6 
psvt(:,jp_ch_m0j) = zm(:,4) * 1e-6
!
IF (lvarsigi) psvt(:,jp_ch_m6i) = zm(:,3) 
IF (lvarsigj) psvt(:,jp_ch_m6j) = zm(:,6)
!
DO jj=1,SIZE(psvt,2)
  psvt(:,jj) =  psvt(:,jj) / (zden2mol * prhodref(:))
ENDDO
!
IF (lhook) CALL dr_hook('MODE_AER_SURF:AERO2PPP_SURF',1,zhook_handle)
!
END SUBROUTINE aero2ppp_surf

END MODULE mode_aer_surf