ecearth3/sources/tm5mp/proj/cbm4/calc_pm.F90

#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') rname, __FILE__, __LINE__ ; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:      CALC_PM
!
! !DESCRIPTION: PM module to calculate particulate matter concentrations for
!               user defined aerodynamic diameters
!\\
!\\
! !INTERFACE: 
!
MODULE CALC_PM
  !
  ! !USES:
  !
  USE GO, only : gol, goErr, goPr

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  Public :: PMX_Integrate_3d
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter  ::  mname = 'calc_pm'
  Real,             Parameter  ::  hr2 = 0.5*sqrt(2.0)
  !
  ! !REVISION HISTORY: 
  !    1 Oct 2010 - Achim Strunk - Revised and made suitable for 3D and stations
  !   19 Jun 2012 -  P. Le Sager - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !   (1) module used only if with_m7 is used
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PMX_Integrate_0d
  !
  ! !DESCRIPTION: This routine generates a pmx value for a given grid cell
  !               Arbitrary radii may be specified. 
  !               NO3 and NH4 come from EQSAM and are considered to be in 
  !               the accumulation mode.
  !\\
  !\\
  ! !INTERFACE:
  !
  Subroutine PMX_Integrate_0d(region, is, js, ls, isizelimit, pmx, status)
    !
    ! !USES:
    !
    Use tracer_data, only : mass_dat

    ! I need to convert (/kg air) to (/m3 air)
    use chem_param,  only : iso4nus
    use chem_param,  only : iso4ais, ibcais, ipomais
    use chem_param,  only : iso4acs, ibcacs, ipomacs, issacs, iduacs
    use chem_param,  only : iso4cos, ibccos, ipomcos, isscos, iducos
    use chem_param,  only : ibcaii, ipomaii, iduaci
    use chem_param,  only : iducoi, ino3_a, inh4
    Use mo_aero_m7,  only : nmod, cmedr2mmedr, sigma

    Use m7_data,     only : h2o_mode, rw_mode

    ! To rewrite ug / cell to ug/m3
    Use global_data, only : region_dat
    Use Meteodata,   only : gph_dat
    !
    ! !INPUT PARAMETERS:
    !
    Integer,                Intent(In)    :: region
    real,                   intent(in)    :: isizelimit
    integer,                Intent(in)    :: is, js, ls
    !
    ! !OUTPUT PARAMETERS:
    !
    Integer,                Intent(out)   :: status
    real,                   Intent(out)   :: pmx   ! units = 'kg/m3'
    !
    ! !REVISION HISTORY: 
    !    7 Oct 2010 - Achim Strunk - v0
    !   19 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter    :: rname = mname//'/PMX_Integrate_0d'

    Integer :: imod
    Real    :: z, rad, sizelimit
    Real    :: volume                ! For converting mass to mass/volume
    Real, Dimension(nmod) :: modfrac ! Length: Number of M7 modes (7)

    Real, Dimension(nmod) :: lnsigma

    ! --- begin -------------------------------------
    
    lnsigma = log(sigma)
    
    ! We can convert micrometers diameter to meters radius. One micrometer diameter is 5.e-7 meter radius.
    sizelimit = isizelimit * 5.e-7

    ! initialise target value
    pmx = 0.0

    Do imod = 1, nmod
       ! Calculate the median radius of the volume weighted distribution.
       rad = rw_mode(region,imod)%d3(is,js,ls) * cmedr2mmedr(imod)

       !!             if( rad == 0.0 ) then ! changed to a precision depending value
       if( rad < 100*TINY(1.0) ) then
          modfrac(imod) = 1.0 ! Should not matter, because if the radius is zero, 
          ! there are no aerosols. But in principle, it would 
          ! mean that all aerosols are infinitely small, so 
          ! they would fit in any PM-class.
       else
          z = ( log(sizelimit) - log(rad) ) / lnsigma(imod)
          modfrac(imod) = 0.5 + 0.5 * ERF(z*hr2)
       end if
    end do
    
    ! And the volume. This is the nicest and cleanest line of code of this whole module.
    volume = ( gph_dat(region)%data(is,js,ls+1) - gph_dat(region)%data(is,js,ls) ) * &
         region_dat(region)%dxyp(js)

    pmx = pmx + mass_dat(region)%rm(is,js,ls, iso4nus) * modfrac(1) / volume 
    pmx = pmx + mass_dat(region)%rm(is,js,ls, iso4ais) * modfrac(2) / volume 
    pmx = pmx + mass_dat(region)%rm(is,js,ls, iso4acs) * modfrac(3) / volume 
    pmx = pmx + mass_dat(region)%rm(is,js,ls, iso4cos) * modfrac(4) / volume 
    pmx = pmx + mass_dat(region)%rm(is,js,ls, ibcais ) * modfrac(2) / volume 
    pmx = pmx + mass_dat(region)%rm(is,js,ls, ibcacs ) * modfrac(3) / volume 
    pmx = pmx + mass_dat(region)%rm(is,js,ls, ibccos ) * modfrac(4) / volume 
    pmx = pmx + mass_dat(region)%rm(is,js,ls, ibcaii ) * modfrac(5) / volume 
    pmx = pmx + mass_dat(region)%rm(is,js,ls, ipomais) * modfrac(2) / volume 
    pmx = pmx + mass_dat(region)%rm(is,js,ls, ipomacs) * modfrac(3) / volume 
    pmx = pmx + mass_dat(region)%rm(is,js,ls, ipomcos) * modfrac(4) / volume 
    pmx = pmx + mass_dat(region)%rm(is,js,ls, ipomaii) * modfrac(5) / volume 
    pmx = pmx + mass_dat(region)%rm(is,js,ls, issacs ) * modfrac(3) / volume 
    pmx = pmx + mass_dat(region)%rm(is,js,ls, isscos ) * modfrac(4) / volume 
    pmx = pmx + mass_dat(region)%rm(is,js,ls, iduacs ) * modfrac(3) / volume 
    pmx = pmx + mass_dat(region)%rm(is,js,ls, iducos ) * modfrac(4) / volume 
    pmx = pmx + mass_dat(region)%rm(is,js,ls, iduaci ) * modfrac(6) / volume 
    pmx = pmx + mass_dat(region)%rm(is,js,ls, iducoi ) * modfrac(7) / volume 
    pmx = pmx + mass_dat(region)%rm(is,js,ls, ino3_a ) * modfrac(3) / volume 
    pmx = pmx + mass_dat(region)%rm(is,js,ls, inh4   ) * modfrac(3) / volume 

    ! Count the water
    pmx = pmx + h2o_mode(region,1)%d3(is,js,ls)*modfrac(1) / volume 
    pmx = pmx + h2o_mode(region,2)%d3(is,js,ls)*modfrac(2) / volume 
    pmx = pmx + h2o_mode(region,3)%d3(is,js,ls)*modfrac(3) / volume 
    pmx = pmx + h2o_mode(region,4)%d3(is,js,ls)*modfrac(4) / volume 

    status = 0

  End Subroutine PMX_INTEGRATE_0D
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  PMX_Integrate_3d
  !
  ! !DESCRIPTION: This routine uses PMX_Integrate_0d for generating a 3d 
  !               array (over im, jm, lm) of pmx. 
  !               Arbitrary radii may be specified. 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PMX_Integrate_3d( region, isizelimit, pmx, status )
    !
    ! !USES:
    !
    USE DIMS,          ONLY : im, jm, lm
    USE TM5_DISTGRID,  ONLY : dgrid, Get_DistGrid, gather, scatter_i_band
    !
    ! !INPUT PARAMETERS:
    !
    INTEGER,                INTENT(in)    :: region
    REAL,                   INTENT(in)    :: isizelimit
    !
    ! !OUTPUT PARAMETERS:
    !
    INTEGER,                INTENT(out)   :: status
    REAL, DIMENSION(:,:,:), INTENT(out)   :: pmx
    !
    ! !REVISION HISTORY: 
    !    7 Oct 2010 - Achim Strunk - v0
    !   19 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    CHARACTER(len=*), PARAMETER :: rname = mname//'/PMX_Integrate_3d'
    INTEGER                     :: is, js, ls, i1, i2, j1, j2
    REAL                        :: pmxl 

    ! reset value
    pmx = 0.0

    CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

    DO ls =  1, lm(region)
    DO js = j1, j2
    DO is = i1, i2

             CALL pmx_integrate_0d( region, is, js, ls, isizelimit, pmxl, status ) 
             IF_NOTOK_RETURN( status=1 )

             pmx(is,js,ls) = pmxl
    END DO
    END DO
    END DO

    status = 0

  END SUBROUTINE PMX_Integrate_3d
  !EOC

END MODULE CALC_PM