mode_dslt_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_dslt_surf
!!   ########################
!!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
 CONTAINS
!!
!!   ############################################################
SUBROUTINE massflux2momentflux(     &
          pflux,                    & ! [kg/m2/s] for M3, zero for other]
          prhodref,                 & ! [kg/m3] air density
          pemisradius,              & ! [um] emitted radius for the different modes
          pemissig,                 & ! [-] emitted sigma for the different modes
          kmde,                     &
          pconvertfacm0,            &
          pconvertfacm6,            &
          pconvertfacm3,            &
          ovarsig,                  &
          orgfix                    &
          )  
!!   ############################################################
!!
!!    PURPOSE
!!    -------
!!    Transform emissions in mass (kg/m2/sec) to emissions of moments which have
!!    a bit strange units
!!    MESONH carries the following units during transport:
!!    M0=#/molec_{air}
!!    M3=molec_{dst}/molec_{air}
!!    M6=um6/molec_{air}*1.d6
!!    The surface model should have (for dust)
!!    M0=#/m3*[kg_{dst}/mole_{dst}/XAVOGADRO]
!!    M3=kg/m3
!!    M6=um6/m3
!!
!!    REFERENCE
!!    ---------
!!    Tulet et al, ORILAM manuscript for transformation of modal parameters
!!    J. Geophys. Res., 110, D18201, doi:10.1029/2004JD005716
!!
!!    AUTHOR
!!    ------
!!    Alf Grini and Pierre TULET (CNRM/GMEI)
!!
!!    MODIFICATIONS
!!    -------------
!!      J.Escobar     06/2013  for REAL4/8 add EPSILON management
!!    
!!
!!    EXTERNAL
!!    --------
!!    None
!!
IMPLICIT NONE
!!
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
!*      0.1    declarations of arguments
!
REAL,   DIMENSION(:,:),      INTENT(INOUT) :: pflux     !In; kg/m2/s (index #2, #5, #8 etc)
                                                        !Out: mole particles per mole air m/s *(MWdst/MWair*rhoair)(index #1)
                                                        !Out: kg/m2/s (index #2)
                                                        !Out: moles m6/moles air m/s *(MWdst/MWair*rhoair)(index #3)
REAL,   DIMENSION(:),        INTENT(IN)    :: prhodref  !I [kg/m3] density of air
REAL,   DIMENSION(:),        INTENT(IN)    :: pemisradius !I [um] emitted radius
REAL,   DIMENSION(:),        INTENT(IN)    :: pemissig    !I [-] emitted sigma for the modes
INTEGER,        INTENT(IN)    :: kmde
REAL,           INTENT(IN)    :: pconvertfacm0
REAL,           INTENT(IN)    :: pconvertfacm6
REAL,           INTENT(IN)    :: pconvertfacm3
LOGICAL,        INTENT(IN)    :: ovarsig
LOGICAL,        INTENT(IN)    :: orgfix
!
!*      0.2    declarations local variables
!
REAL,DIMENSION(SIZE(PFLUX,1),3) :: zfm               !Intermediate variable to get moments
!
INTEGER   :: jmode  ! Counter for dust modes
INTEGER   :: jsv_idx ! Counter for dust scalar variables
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------
!
!MESONH carries the following units during transport:
!M0=#/molec_{air}
!M3=molec_{dst}/molec_{air}
!M6=um6/molec_{air}*1.d6

!The surface should get the following units
!M0=#/m3*MW_DST/XAVOGADRO
!M3=kg/m3
!M6=um6/m3*1.d6 MW_DST/XAVOGADRO

!Emissions of dust are in kg/m2/sec for mode 3 at this point

!Factor which is needed so that all gains normal units when leaving ground paramn
IF (lhook) CALL dr_hook('MODE_DSLT_SURF:MASSFLUX2MOMENTFLUX',0,zhook_handle)

!Initialize initermediate moments
zfm(:,:)=0.

DO jmode=1,kmde

  !Make index which is 0 for first mode, 3 for second, 6 for third etc
  IF (ovarsig) THEN
    jsv_idx = (jmode-1)*3
  ELSE IF (orgfix) THEN
    jsv_idx = jmode-2
  ELSE
    jsv_idx = (jmode-1)*2
  END IF

  !IN THIS VERSION, MASS FLUX (kg/m2/sec) IS SENT IN INDEX #2, #5, #8 if 3 moments per mode
  !IF TWO MOMENTS PER MODE, MASS FLUX IS SENT AS INDEX #2, #4, #6

  !Get flux of number in #/m2/sec from flux of mass in kg/m2/sec
  zfm(:,1) = pflux(:,jsv_idx+2)                     & ! kg_{dst}/m2/sec
             / pemisradius(jmode)**3                & ! *um^{-3} ==> #/m2/sec*(m3/um3)
             * exp(-4.5*(log(pemissig(jmode)))**2)  & ! Take into account size distribution  
             / pconvertfacm3           ! /(kg_{dst}/m^{3}_{dst)} ==> m^3_{dst}/m2/sec  
 
   
  ! Get flux of moment 6 consistent with the other moments
  zfm(:,3) = zfm(:,1)                              & ! [#/m3]
             * (pemisradius(jmode)**6)             & ! *um6 ==> um6/m2/sec 
             * exp(18. *(log(pemissig(jmode)))**2)   ! Take into account size distribution  

  !Get flux of Moment 0 in transport units
  IF (.NOT.orgfix) THEN
    pflux(:,jsv_idx+1) = zfm(:,1)            & ! particles/m^2/sec
                         * pconvertfacm0       ! ==> particles/m2/sec * kg_dst/m3_{air}  
  END IF
   
  ! Flux moment 6
  IF (ovarsig) THEN
    pflux(:,jsv_idx+3) = zfm(:,3)          & ! um^6/m^2/sec
                         * pconvertfacm6     ! ==>   
  ENDIF

  !Multiply with molecular weights so that you get back the units described above when
  !when multiply with the opposite variable in ground_paramn.f90
  !PFLUX(:,JSV_IDX+1) = PFLUX(:,JSV_IDX+1) * 100.E-3 * PRHODREF(:) / ZMD   !#_{aer}/molec_{air} m/s * kg_{aer}/m^3_{air}
  !IF (LVARSIG) PFLUX(:,JSV_IDX+3) = PFLUX(:,JSV_IDX+3) * 100.E-3 * PRHODREF(:) / ZMD   !um^6_{aer}/molec_{air}*cm^3/m^3 m/s kg_{aer}/m^3_{air}
  !
  !
ENDDO !Loop on modes
!
IF (lhook) CALL dr_hook('MODE_DSLT_SURF:MASSFLUX2MOMENTFLUX',1,zhook_handle)
!
END SUBROUTINE massflux2momentflux

!**********************************************************************
!**********************************************************************
!**********************************************************************

SUBROUTINE dsltmoment2size(       &
          psvt,                   & !I [XX/m3] input scalar variables (moment of distribution)
          prhodref,               & !I [kg/m3] density of air       
          pemisradius,            & ![um] emitted radius for the different modes
          pemissig,               & ![-] emitted sigma for the different modes   
          km0,                    &
          km3,                    &
          km6,                    &
          pconvertfacm0,          &
          pconvertfacm6,          &
          pconvertfacm3,          &
          ovarsig,                &
          orgfix,                 &   
          psig1d,                 & !O [-] standard deviation of aerosol distribution
          prg1d,                  & !O [um] number median diameter of aerosol distribution
          pn1d,                   & !O [#/m3] number concentration of aerosols
          pmass1d,                & !O [kg/m3] mass concentration of aerosol
          pm1d                    & !O aerosols 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)
!!    At this point, M3 is in kg/m3, M0 in #/m3*(kg_{dst}/mole), M6 in um6/m3*1.d6*(kg_{dst}/mole)
!!
!!    All the moments have been transformed in MESONH (atmospheric model) so that the surface gets
!!    M0 [#/m3] *XMOLARWEIGHT_DST/XAVOGADRO
!!    M3 [kg/m3]
!!    M6 [um6/m3*1.d6] *XMOLARWEIGHT_DST/XAVOGADRO
!!   
!!    REFERENCE
!!    ---------
!!    Tulet et al, ORILAM manuscript for transformation of modal parameters
!!    J. Geophys. Res., 110, D18201, doi:10.1029/2004JD005716
!!
!!    AUTHOR
!!    ------
!!    Pierre TULET (LA)
!!
!!    MODIFICATIONS
!!    -------------
!!    Alf Grini (CNRM)
!!
!!    EXTERNAL
!!    --------
!!    
USE modd_surf_par , ONLY : xsurf_tiny
!!
IMPLICIT NONE
!!
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
!*      0.1    declarations of arguments
!
!INPUT
REAL,       DIMENSION(:,:),  INTENT(IN)     :: psvt      !I [ppp] moments in surface units
REAL,       DIMENSION(:),    INTENT(IN)     :: prhodref  !I [kg/m3] density of air
REAL,       DIMENSION(:),    INTENT(IN)     :: pemissig  
REAL,       DIMENSION(:),    INTENT(IN)     :: pemisradius
INTEGER,DIMENSION(:),    INTENT(IN) :: km0             ! [idx] index for Mode 0 in passed variables
INTEGER,DIMENSION(:),    INTENT(IN) :: km3             ! [idx] indexes for Mode 3 in passed variables
INTEGER,DIMENSION(:),    INTENT(IN) :: km6             ! [idx] indexes for Mode 6 in passed variables
REAL,           INTENT(IN)    :: pconvertfacm0
REAL,           INTENT(IN)    :: pconvertfacm6
REAL,           INTENT(IN)    :: pconvertfacm3
LOGICAL,        INTENT(IN)    :: ovarsig
LOGICAL,        INTENT(IN)    :: orgfix
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)     :: pmass1d  !O [kg_{aer}/m3] mass concentration
REAL,       DIMENSION(:,:),  OPTIONAL, INTENT(OUT)     :: pm1d     !O aerosols moments (MESONH units)
!
!*      0.2    declarations local variables
!
REAL,DIMENSION(SIZE(PSVT,1), SIZE(PSVT,2)) :: zsv    ! [dusts moment concentration]
REAL,DIMENSION(SIZE(PSVT,1), SIZE(KM0)*3)  :: zm     ! [moments] local array for moments  
REAL,DIMENSION(SIZE(PSVT,1))               :: zsigma ! [-] standard deviation
REAL,DIMENSION(SIZE(PSVT,1))               :: zrg    ! [um] number median diameter
INTEGER                   :: jn, j0, j3, j6          ! [idx] loop counters
REAL(KIND=JPRB) :: zhook_handle
!
!        1.1    initialisation 
! 
!Get the conversion factors
IF (lhook) CALL dr_hook('MODE_DSLT_SURF:DSLTMOMENT2SIZE',0,zhook_handle)
!  
!Get scalar variable indexes
!
!Save the moments in a local array
zsv(:,:) = max(psvt(:,:),  xsurf_tiny)
!
DO jn=1,SIZE(km0)

  j0 = 1 + (jn-1)*3 
  j3 = 2 + (jn-1)*3
  j6 = 3 + (jn-1)*3

  !calculate moment 3 from total aerosol mass in kg/m3 ==> um3/m3
  zm(:,j3) = zsv(:,km3(jn))          & ! kg_{aer}/m3_{air}
             / pconvertfacm3           ! ==> m3_{dst}/m3_{air}

  IF (ovarsig) THEN ! give M6 (case of variable standard deviation)
              
    !Get number concentration (#/molec_{air}==>#/m3)
    zm(:,j0) = zsv(:,km0(jn))          & ! #/m3air*M_{dst}/avogadro 
               / pconvertfacm0           ! ==> #/m3  
                 
    !Calculate moment 6 from the sent value
    zm(:,j6) = zsv(:,km6(jn))          & ! um6/m3_{air}*(cm3/m3)*M_{dst}/Avogadro
               / pconvertfacm6           ! ==> um6/m3  

    !Get sigma (only if sigma is allowed to vary)
    !Get intermediate values for sigma M3^2/(M0*M6) (ORILAM paper, eqn 8)
    zsigma(:) = zm(:,j3)**2 / (zm(:,j0)*zm(:,j6))
    !Limit the intermediate value, can not be larger than 1
    zsigma(:) = min(1-1e-10,zsigma(:))
    !Limit the value for intermediate, can not be smaller than 0
    zsigma(:) = max(1e-10,zsigma(:))
    !Calculate log(sigma)
    zsigma(:) = log(zsigma(:))
    !Finally get the real sigma the negative sign is because of 
    !The way the equation is written (M3^2/(M0*M6)) instead of (M0*M6)/M3^3
    zsigma(:) = exp(1./3.*sqrt(-zsigma(:)))
        
  ELSE IF (orgfix) THEN ! compute M6 from M3, Rg and SIGMA    

    !Get the emitted sigma for this mode
    zsigma(:) = pemissig(jn)

    zm(:,j0) = zm(:,j3) /               &
              ((pemisradius(km3(jn))**3)*exp(4.5 * log(zsigma(:))**2))  

  ELSE ! compute M6 from M0, M3 and SIGMA
          
    !Get the emitted sigma for this mode
    zsigma(:) = pemissig(jn)

    !Get number concentration (#/molec_{air}==>#/m3)
    zm(:,j0) = zsv(:,km0(jn))        & ! #/m3air*M_{dst}/avogadro 
               / pconvertfacm0         ! ==> #/m3  

  END IF

  !Calculate moment 6 from this emitted sigma
  zm(:,j6) = zm(:,j0) * ((zm(:,j3)/zm(:,j0))**(1./3.) &
            * exp(-(3./2.)*log(zsigma(:))**2))**6     &
            * exp(18.*log(zsigma(:))**2)  

  !Get number median radius (eqn. 7 in Orilam manuscript)
  zrg(:) = ((zm(:,j3)**4) / (zm(:,j6)*zm(:,j0)**3)) ** (1./6.)     
  
  !Give the sigma-values to the passed array
  IF(present(psig1d)) psig1d(:,jn) = zsigma(:)

  !Set the number concentrations in the passed array
  IF(present(pn1d)) pn1d(:,jn) = zm(:,j0)

  !Get the number median radius
  IF(present(prg1d)) prg1d(:,jn)= zrg(:)
    
  IF(present(pmass1d))THEN
    pmass1d(:,jn)=  zm(:,j0)          &!#/m^3_{air}
                    * pconvertfacm3   &
                    * zrg(:)**3 * exp(4.5*(log(zsigma(:)))**2)
  ENDIF

END DO  !Loop on modes

IF(present(pm1d)) pm1d(:,:) = zm(:,:)

IF (lhook) CALL dr_hook('MODE_DST_SURF:DUSTMOMENT2SIZE',1,zhook_handle)
!
!
END SUBROUTINE dsltmoment2size


END MODULE mode_dslt_surf