#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