ecearth3/sources/tm5mp/proj/cb05/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. This is used for diagnostics
!               purposes by some user_output_* routines.
!\\
!\\
! !INTERFACE: 
!
MODULE CALC_PM
  !
  ! !USES:
  !
  USE GO, only : gol, goErr, goPr

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  Public :: PMX_Integrate_3d,PMx_integrate_0d
  !
  ! !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 binas,       only : pi
    Use tracer_data, only : mass_dat

    ! I need to convert (/kg air) to (/m3 air)
    use chem_param,  only : iso4nus,isoanus
    use chem_param,  only : iso4ais, ibcais, ipomais, isoaais
    use chem_param,  only : iso4acs, ibcacs, ipomacs, isoaacs, issacs, iduacs
    use chem_param,  only : iso4cos, ibccos, ipomcos, isoacos, isscos, iducos
    use chem_param,  only : ibcaii, ipomaii, isoaaii, iduaci
    use chem_param,  only : inus_n, iais_n, iacs_n, icos_n, iaii_n, iaci_n, icoi_n
    use chem_param,  only : iducoi, ino3_a, inh4, imsa
    Use mo_aero_m7,  only : nmod, cmedr2mmedr, sigma, cmr2ram

    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 ! aerodynamic diameter in micrometers
    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
    !   21 Feb 2023 - T. van Noije - changed to cut-off based on aerodynamic diameter
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    Integer :: imod
    Real    :: z, rad, rad_limit
    Real    :: dry_mass, wet_mass, particle_number, mode_volume, mode_density
    Real    :: air_volume                ! For converting mass to density in kg/m3
    Real    :: modfrac

    Real, Dimension(nmod) :: lnsigma

    ! --- begin -------------------------------------

    ! PMx is defined as the mass below a threshold diameter x,
    ! where x is the aerodynamic diameter in micrometers.
    ! The aerodynamic diameter is defined as the diameter of a
    ! spherical particle with equal settling velocity but a material density of 1 g/cm^3.
    ! For particles larger than the mean free path of air,
    ! the aerodynamic diameter is approximated by the wet diameter times
    ! the square root of the particle density in g/cm^3.
    ! This approximation is often used for particles as small as 0.5 micron diameter 
    ! (e.g., US EPA 1996, Air Quality Criteria for Particulate Matter, Volume I, 1996)
    ! We currently do not correct for non-spherical shapes.
    ! That could be an additional correction.
    
    lnsigma = log(sigma)
    
    ! initialise target value
    pmx = 0.0

    Do imod = 1, nmod

       ! Calculate the dry and wet mass per mode. 
       ! The calculation is kept consistent with the earlier version of the code,
       ! where the mass of NO3_a and NH4 are added to the soluble accumulation mode.
       ! Consistent with the treatment in other parts of the code,
       ! we add the MSA mass also to the soluble accumulation code.
       ! Note that it would be possible to make the calculation more consistent with
       ! the calculation of the modal mass in M7.
       ! However, as I cannot exclude that this subroutine is called at initialization,
       ! this would require to add the dry (or wet) modal mass field to the restart file.
       ! Because of this extra complication, I haven't done that.

       select case (imod)
       case (1) 
          dry_mass = mass_dat(region)%rm(is,js,ls, iso4nus) + mass_dat(region)%rm(is,js,ls, isoanus)
          wet_mass = dry_mass + h2o_mode(region,1)%d3(is,js,ls)
          particle_number = mass_dat(region)%rm(is,js,ls, inus_n)
       case (2)
          dry_mass = mass_dat(region)%rm(is,js,ls, iso4ais) + mass_dat(region)%rm(is,js,ls, ibcais ) + &
                     mass_dat(region)%rm(is,js,ls, ipomais) + mass_dat(region)%rm(is,js,ls, isoaais)
          wet_mass = dry_mass + h2o_mode(region,2)%d3(is,js,ls)
          particle_number = mass_dat(region)%rm(is,js,ls, iais_n)
       case (3)
          dry_mass = mass_dat(region)%rm(is,js,ls, iso4acs) + mass_dat(region)%rm(is,js,ls, ibcacs ) + &
                     mass_dat(region)%rm(is,js,ls, ipomacs) + mass_dat(region)%rm(is,js,ls, isoaacs) + &
                     mass_dat(region)%rm(is,js,ls, issacs ) + mass_dat(region)%rm(is,js,ls, iduacs ) + &
                     mass_dat(region)%rm(is,js,ls, ino3_a ) + mass_dat(region)%rm(is,js,ls, inh4   ) + &
                     mass_dat(region)%rm(is,js,ls, imsa )
          wet_mass = dry_mass + h2o_mode(region,3)%d3(is,js,ls)
          particle_number = mass_dat(region)%rm(is,js,ls, iacs_n)
       case (4)
          dry_mass = mass_dat(region)%rm(is,js,ls, iso4cos) + mass_dat(region)%rm(is,js,ls, ibccos ) + &
                     mass_dat(region)%rm(is,js,ls, ipomcos) + mass_dat(region)%rm(is,js,ls, isoacos) + &
                     mass_dat(region)%rm(is,js,ls, isscos ) + mass_dat(region)%rm(is,js,ls, iducos )
          wet_mass = dry_mass + h2o_mode(region,4)%d3(is,js,ls)
          particle_number = mass_dat(region)%rm(is,js,ls, icos_n)
       case (5)
          dry_mass = mass_dat(region)%rm(is,js,ls, ibcaii ) + &
                     mass_dat(region)%rm(is,js,ls, ipomaii) + mass_dat(region)%rm(is,js,ls, isoaaii)
          wet_mass = dry_mass
          particle_number = mass_dat(region)%rm(is,js,ls, iaii_n)
       case (6)
          dry_mass = mass_dat(region)%rm(is,js,ls, iduaci )
          wet_mass = dry_mass
          particle_number = mass_dat(region)%rm(is,js,ls, iaci_n)
       case (7)
          dry_mass = mass_dat(region)%rm(is,js,ls, iducoi )
          wet_mass = dry_mass
          particle_number = mass_dat(region)%rm(is,js,ls, icoi_n)
       end select

       if ( dry_mass > 100*TINY(1.0) .and. particle_number > 100*TINY(1.0) ) then

          ! Calculate the median radius of the volume weighted distribution.
          rad = rw_mode(region,imod)%d3(is,js,ls) * cmedr2mmedr(imod) 
          
          ! In principe, as the dry mass in the considered mode is non-zero,
          ! the mode radius or volume cannot be zero at this point.
          ! However, as the mass calculated here is not fully consistent with that in M7,
          ! an additional check on the mode radius is included.
          ! The exact same condition is used as in the earlier version of the code,
          ! where PMx was calculated using for the threshold diameter x
          ! the geometric wet diameter instead of the aerodynamic diameter.
          if( rad < 100*TINY(1.0) ) then     
             ! Particles are so small that all mass should be included.
             modfrac  = 1.0
          else
             ! Calculate the mass density of ambient particles for each mode
             mode_volume = particle_number * ((rw_mode(region,imod)%d3(is,js,ls)*cmr2ram(imod))**3.0)*pi/0.75
             mode_density = wet_mass / mode_volume ! kg/m3

             ! Convert the threshold aerodynamic diameter in microns
             ! to a geometric radius in meters.
             rad_limit = isizelimit * ((1000./mode_density)**0.5) * 5.e-7

             z = ( log(rad_limit) - log(rad) ) / lnsigma(imod)
             modfrac = 0.5 + 0.5 * ERF(z*hr2)
          endif

          ! Sum contributions from the different modes to the dry mass below the threshold diameter
          pmx = pmx + dry_mass * modfrac
       endif
    end do
    
    ! Finally divide by the air volume to get PMx in kg/m3
    air_volume = ( gph_dat(region)%data(is,js,ls+1) - gph_dat(region)%data(is,js,ls) ) * &
         region_dat(region)%dxyp(js)

    pmx = pmx / air_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
             pmxl=0.0
             CALL pmx_integrate_0d( region, is, js, ls, isizelimit, pmxl, status ) 
             IF_NOTOK_RETURN( status=1 )
             !write(4444,*)pmxl,is,js,ls,ubound(pmx),lbound(pmx)

             pmx(is,js,ls) = pmxl
             !write(4444,*)pmx(is,js,ls),is,js,ls,ubound(pmx),lbound(pmx)
          END DO
       END DO
      ! write(5555,*)size(pmx),shape(pmx),i2,j2,lm(region),is,js,ls,isizelimit,pmxl
    END DO

    status = 0

  END SUBROUTINE PMX_Integrate_3d
  
!NOT-USED  SUBROUTINE PMX_Integrate_2d( region, isizelimit, pmx, status )
!NOT-USED    !
!NOT-USED    ! !USES:
!NOT-USED    !
!NOT-USED    USE DIMS,          ONLY : im, jm, lm
!NOT-USED    USE TM5_DISTGRID,  ONLY : dgrid, Get_DistGrid, gather, scatter_i_band
!NOT-USED    !
!NOT-USED    ! !INPUT PARAMETERS:
!NOT-USED    !
!NOT-USED    INTEGER,                INTENT(in)    :: region
!NOT-USED    REAL,                   INTENT(in)    :: isizelimit
!NOT-USED    !
!NOT-USED    ! !OUTPUT PARAMETERS:
!NOT-USED    !
!NOT-USED    INTEGER,                INTENT(out)   :: status
!NOT-USED    REAL, DIMENSION(:,:), INTENT(out)   :: pmx
!NOT-USED    !
!NOT-USED    ! !REVISION HISTORY: 
!NOT-USED    !    7 Oct 2010 - Achim Strunk - v0
!NOT-USED    !   19 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
!NOT-USED    !
!NOT-USED    ! !REMARKS:
!NOT-USED    !
!NOT-USED    !EOP
!NOT-USED    !------------------------------------------------------------------------
!NOT-USED    !BOC
!NOT-USED 
!NOT-USED    CHARACTER(len=*), PARAMETER :: rname = mname//'/PMX_Integrate_3d'
!NOT-USED    INTEGER                     :: is, js, ls, i1, i2, j1, j2
!NOT-USED    REAL                        :: pmxl 
!NOT-USED 
!NOT-USED    ! reset value
!NOT-USED    pmx = 0.0
!NOT-USED 
!NOT-USED    CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
!NOT-USED 
!NOT-USED    
!NOT-USED    DO js = j1, j2
!NOT-USED    DO is = i1, i2
!NOT-USED 
!NOT-USED             CALL pmx_integrate_0d( region, is, js, ls, isizelimit, pmxl, status ) 
!NOT-USED             IF_NOTOK_RETURN( status=1 )
!NOT-USED 
!NOT-USED             pmx(is,js) = pmxl
!NOT-USED    END DO
!NOT-USED    END DO
!NOT-USED    
!NOT-USED 
!NOT-USED    status = 0
!NOT-USED 
!NOT-USED  END SUBROUTINE PMX_Integrate_2d
!NOT-USED  !EOC

END MODULE CALC_PM