#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