Last active
March 2, 2016 03:14
-
-
Save JiaoyanHuang/1717227d868d6ffc2b79 to your computer and use it in GitHub Desktop.
aircraft aerosol emission
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
!------------------------------------------------------------------------! | |
! The Community Multiscale Air Quality (CMAQ) system software is in ! | |
! continuous development by various groups and is based on information ! | |
! from these groups: Federal Government employees, contractors working ! | |
! within a United States Government contract, and non-Federal sources ! | |
! including research institutions. These groups give the Government ! | |
! permission to use, prepare derivative works of, and distribute copies ! | |
! of their work in the CMAQ system to the public and to permit others ! | |
! to do so. The United States Environmental Protection Agency ! | |
! therefore grants similar permission to use the CMAQ system software, ! | |
! but users are requested to provide copies of derivative works or ! | |
! products designed to operate in the CMAQ system to the United States ! | |
! Government without restrictions as to use by others. Software ! | |
! that is used with the CMAQ system but distributed under the GNU ! | |
! General Public License or the GNU Lesser General Public License is ! | |
! subject to their copyright restrictions. ! | |
!------------------------------------------------------------------------! | |
C RCS file, release, date & time of last delta, author, state, [and locker] | |
C $Header: /project/yoj/arc/CCTM/src/emis/emis/ACAERO_EMIS.F,v 1.7 2012/01/19 15:21:13 yoj Exp $ | |
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: | |
module acaero_emis | |
C----------------------------------------------------------------------- | |
C Description: | |
C * Extracts aircraft emission from emission files | |
C | |
C 24 Nov 15 J. Huang created | |
C----------------------------------------------------------------------- | |
use aero_data, only: n_mode | |
implicit none | |
C number of chemical species in acaero aerosol | |
C integer, parameter :: nacaero_spc = 3 | |
integer, parameter :: ac_emlyr = 17 | |
integer :: IDX | |
C Aircraft Emissions Rates | |
real, allocatable, save :: ac_embuff ( :,:,:,: ) ! in all grid cells (C,R,L,S) | |
real, allocatable, save :: acaerooutm( :,:,:,:,: ) ! mass emission rates [ug/m**3/s](S,N,R,C,L) | |
real, allocatable, save :: acaerooutn( :,:,:,: ) ! number emission rates [1/m**3/s] (N,R,C,L) | |
real, allocatable, save :: acaeroouts( :,:,:,:) ! surface-area emisrates [m2/m**3/s] (N,R,C,L) | |
real, allocatable, save :: ac_em_rate ( : ) | |
real, allocatable, save :: ac_pm_em ( : ) | |
real, allocatable, save :: ac_em_part ( :, : ) | |
C 24 Nov 15 J. Huang | |
C emission split factor from Kinsey et al., 2009 (AEPX reprot 1-3) | |
C Kinsey et al., 2010 Atmos. Environ. 44 2147-2156 | |
C Petzold et al., 1998 JGR 104, D18,page 22, 171-22, 181 | |
C in our aircraft emissions, there are only PSO4, PEC, POC, three species of aerosols | |
real, parameter :: esplit_f( n_mode ) = (/0.5,0.5,0.000/) | |
C Factor for converting aerosol emissions from input units ... | |
real, allocatable, save :: ac_convem_pm (:) ! into [ug/m2/sec] | |
public acaerooutm, acaerooutn, acaeroouts, ! acaero_spc, nacaero_spc | |
& acaero_emis_init, get_acaero_emis | |
private | |
C Variables for converting emission rates into molar-mixing-ratio units | |
REAL, PARAMETER :: GPKG = 1.0E+03 ! g/kg | |
REAL, PARAMETER :: MGPG = 1.0E+06 ! ug/g | |
real(8) :: GSFAC | |
C diam`s from FRACMI,FRACMJ-weighted 2 2-bin averages of geom means | |
C Petzold et al., 1998 JGR 104, D18,page 22, 171-22, 181 | |
real, parameter :: dgvi = 0.030 ! geom mean diam of aitken mode [um] | |
real, parameter :: dgvj = 0.150 ! geom mean diam of accum mode [um] | |
real, parameter :: sigi = 1.55 ! geom std deviation of aitken mode flux | |
real, parameter :: sigj = 1.870 ! geom std deviation of accum mode flux | |
C Local Variables: | |
C Factors for converting 3rd moment emission rates into number and 2nd moment | |
C emission rates. (Diameters in [um] changed to [m] ) See Equations 7b and 7c | |
C of Binkowski & Roselle (2003) | |
real :: l2sgi ! [ln( sigi )] ** 2 | |
real :: l2sgj ! [ln( sigj )] ** 2 | |
real, save :: factnumi ! = exp( 4.5 * l2sgi ) / dgvi ** 3 * 1.0e18 | |
real, save :: factnumj ! = exp( 4.5 * l2sgj ) / dgvj ** 3 * 1.0e18 | |
real, save :: factm2i ! = exp( 0.5 * l2sgi ) / dgvi * 1.0e6 | |
real, save :: factm2j ! = exp( 0.5 * l2sgj ) / dgvj * 1.0e6 | |
real, save :: factsrfi ! = pi * factm2j | |
real, save :: factsrfj ! = pi * factm2k | |
C Domain decomposition info from emission and meteorology files | |
INTEGER, SAVE :: STARTCOL, ENDCOL, STARTROW, ENDROW | |
INTEGER, SAVE :: STRTCOLMC3, ENDCOLMC3, STRTROWMC3, ENDROWMC3 | |
C Setup log file | |
integer, save :: logdev | |
integer :: n | |
CONTAINS | |
C======================================================================= | |
function acaero_emis_init( jdate, jtime, tstep ) result( success ) | |
use hgrd_defn ! horizontal domain specifications | |
use aero_data ! aerosol species definitions | |
use utilio_defn | |
use grid_conf | |
C Arguments: | |
integer, intent( in ) :: jdate ! current model date, coded YYYYDDD | |
integer, intent( in ) :: jtime ! current model time, coded HHMMSS | |
integer, intent( in ) :: tstep ! output time step | |
logical success | |
C Includes: | |
include SUBST_FILES_ID ! file name parameters | |
include SUBST_CONST ! physical and mathematical constants | |
C External Functions: | |
integer, external :: setup_logdev | |
C Local variables: | |
character( 16 ) :: pname = 'ACAERO_EMIS_INIT' | |
character( 16 ) :: vname | |
character( 80 ) :: vardesc | |
character( 120 ) :: xmsg = ' ' | |
integer status | |
integer i, j, c, r, l, s | |
real(8) :: GSFAC | |
C Factor for converting aerosol emissions from input units ... | |
real :: ac_convem_pm_mass! into [ug/sec] | |
C Variables for the thickness and volume of each grid cell | |
REAL, ALLOCATABLE, SAVE :: cellhgt( : ) ! grid-cell height [sigma] | |
REAL, ALLOCATABLE, SAVE :: cellvol( : ) ! grid-cell volume [m2*sigma] | |
C Variables for calculating the volume of each grid cell | |
REAL DX1, DX2 ! grid-cell width and length [m] | |
REAL CELLAREA ! grid-cell area [m2] | |
C Domain decomposition info from emission and meteorology files | |
INTEGER GXOFF, GYOFF ! origin offset | |
logdev = setup_logdev() | |
success = .true. | |
allocate ( acaerooutm( n_aerospc,n_mode,ncols,nrows,ac_emlyr ), | |
& acaerooutn( n_mode,ncols,nrows,ac_emlyr ), | |
& acaeroouts( n_mode,ncols,nrows,ac_emlyr ), stat = status ) | |
if ( status .ne. 0 ) then | |
xmsg = '*** ACAEROOUTM, ACAEROOUTN or ACAEROOUTS memory allocation failed' | |
call m3warn ( pname, jdate, jtime, xmsg ) | |
success = .false.; return | |
end if | |
allocate ( ac_em_rate( n_emis_pm ), ac_pm_em( n_emis_pm ), stat = status ) | |
if ( status .ne. 0 ) then | |
xmsg = '*** ac_pm_em or ac_em_rate memory allocation failed' | |
call m3warn ( pname, jdate, jtime, xmsg ) | |
end if | |
allocate ( ac_em_part ( n_aerospc,n_mode ), stat = status ) | |
if ( status .ne. 0 ) then | |
xmsg = '*** ac_em_part memory allocation failed' | |
call m3warn ( pname, jdate, jtime, xmsg ) | |
end if | |
ALLOCATE ( ac_embuff( ncols,nrows,ac_emlyr,n_emis_pm ), | |
& STAT = STATUS ) | |
IF ( STATUS .NE. 0 ) THEN | |
XMSG = '*** ac_embuff or acaero_em_2 memory allocation failed' | |
CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) | |
END IF | |
ALLOCATE ( ac_convem_pm( ac_emlyr ), | |
& cellhgt( ac_emlyr ), | |
& cellvol( ac_emlyr ), STAT = STATUS ) | |
IF ( STATUS .NE. 0 ) THEN | |
XMSG = '*** CONVEM_PM, cellhgt or cellvol memory allocation failed' | |
CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) | |
END IF | |
C emission rates. (Diameters in [um] changed to [m] ) See Equations 7b and 7c | |
C of Binkowski & Roselle (2003) | |
l2sgi = log( sigi ) * log( sigi ) | |
l2sgj = log( sigj ) * log( sigj ) | |
factnumi = 1.0e18 * exp( 4.5 * l2sgi ) / dgvi ** 3 | |
factnumj = 1.0e18 * exp( 4.5 * l2sgj ) / dgvj ** 3 | |
factm2i = 1.0e06 * exp( 0.5 * l2sgi ) / dgvi | |
factm2j = 1.0e06 * exp( 0.5 * l2sgj ) / dgvj | |
factsrfi = pi * factm2i | |
factsrfj = pi * factm2j | |
#ifdef verbose | |
write( logdev,* ) ' ' | |
write( logdev,* ) ' l2sgi,l2sgj: ', l2sgi, l2sgj | |
write( logdev,* ) ' factnumi,factnumj: ', factnumi, factnumj | |
write( logdev,* ) ' factm2i,factm2j: ', factm2i, factm2j | |
write( logdev,* ) ' factsrfi,factsrfj: ', factsrfi, factsrfj | |
write( logdev,* ) ' modal avg dens(i/j): ', acaeromode_dens( 1 ), acaeromode_dens( 2 ) | |
write( logdev,* ) ' ' | |
#endif | |
C *** Get length and width of each grid cell | |
C note: crude estimate is made for LAT/LONG coordinate systems | |
IF ( GDTYP_GD .EQ. LATGRD3 ) THEN | |
DX1 = DG2M * XCELL_GD ! in m | |
DX2 = DG2M * YCELL_GD | |
& * COS( PI180*( YORIG_GD + YCELL_GD | |
& * FLOAT( GL_NROWS/2 ) ) ) ! in m | |
ELSE | |
DX1 = XCELL_GD ! in m | |
DX2 = YCELL_GD ! in m | |
END IF | |
C *** Get height of grid cell in each layer in sigma coordinates | |
C Multiply by grid area [m2] to obtain grid volume | |
CELLAREA = DX1 * DX2 | |
DO l = 1, ac_emlyr | |
cellhgt( l ) = X3FACE_GD( l ) - X3FACE_GD( l-1 ) | |
cellvol( l ) = cellhgt( l ) * CELLAREA | |
END DO | |
C *** Get scaling factor for converting aerosol emissions from | |
C their input units to [ug/s] and then to [ug/m2/s] using layer- | |
C specific grid-cell volume. Calling DESC3( EMIS_1 ) assumes EMIS_1 | |
C already opened. | |
IF ( PMEM_UNITS .EQ. 'G/S' .OR. | |
& PMEM_UNITS .EQ. 'g/s' ) THEN | |
ac_convem_pm_mass = MGPG ! (g/s) -> (ug/s) | |
ELSE IF ( PMEM_UNITS .EQ. 'KG/HR' .OR. | |
& PMEM_UNITS .EQ. 'kg/hr' ) THEN | |
ac_convem_pm_mass = GPKG * MGPG / 3600.0 ! (kg/hr) -> (ug/s) | |
ELSE | |
XMSG = 'Units incorrect on ' // EMIS_1 | |
CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT2 ) | |
END IF | |
do l = 1, ac_emlyr | |
ac_convem_pm( l ) = ac_convem_pm_mass / cellvol( l ) | |
end do | |
C *** Get domain window info from input files | |
CALL SUBHFILE ( ACAERO_EMIS_1, GXOFF, GYOFF, | |
& STARTCOL, ENDCOL, STARTROW, ENDROW ) | |
SUCCESS = .TRUE.; RETURN | |
end function acaero_emis_init | |
C======================================================================= | |
subroutine opacaero_emis ( jdate, jtime, tstep, nacaero_diag, wracaero_spc ) | |
C 27 Dec 10 J.Young: initial | |
use grid_conf ! horizontal & vertical domain specifications | |
use utilio_defn | |
implicit none | |
include SUBST_FILES_ID ! file name parameters | |
C Arguments: | |
integer jdate ! current model date, coded YYYYDDD | |
integer jtime ! current model time, coded HHMMSS | |
integer tstep ! output time step | |
integer nacaero_diag | |
character( 16 ) :: wracaero_spc( nacaero_diag ) | |
real(8) GSFAC | |
C Local variables: | |
character( 16 ) :: pname = 'OPACAERO_EMIS ' | |
character( 96 ) :: xmsg = ' ' | |
C CHARACTER( 16 ), SAVE :: ctm_acaero_emis_1 = 'ctm_acaero_emis_1' | |
C----------------------------------------------------------------------- | |
C Try to open existing file for update | |
if ( .not. open3( ACAERO_EMIS_1, fsrdwr3, pname ) ) then | |
xmsg = 'Could not open CTM_ACAERO_EMIS_1 for update - ' | |
& // 'try to open new' | |
call m3mesg( xmsg ) | |
end if | |
return | |
end subroutine opacaero_emis | |
C======================================================================= | |
subroutine get_acaero_emis( jdate, jtime, tstep, rjacm ) | |
use grid_conf ! horizontal & vertical domain specifications | |
use utilio_defn | |
use aero_data | |
C 8/18/11 D.Wong: incorporated twoway model implementation and change | |
C RC -> RCA and RN -> RNA and made it backward compatible | |
C Arguments: | |
integer, intent( in ) :: jdate ! current model date, coded YYYYDDD | |
integer, intent( in ) :: jtime ! current model time, coded HHMMSS | |
integer, intent( in ) :: tstep( 3 ) ! sync, output time step | |
real, intent( in ) :: rjacm( ncols,nrows,ac_emlyr ) ! reciprocal Jacobian [1/m] | |
C REAL (8) GSFAC | |
C Variables for the thickness and volume of each grid cell | |
C REAL, ALLOCATABLE, SAVE :: cellhgt( : ) ! grid-cell height [sigma] | |
C REAL, ALLOCATABLE, SAVE :: cellvol( : ) ! grid-cell volume [m2*sigma] | |
C Includes: | |
include SUBST_FILES_ID ! file name parameters | |
include SUBST_CONST ! constants (for PI) | |
! include SUBST_AE_EMIS ! aerosol emission surrogate names and map table | |
C External Functions: | |
C Parameters: | |
real, parameter :: f6dpi = 6.0 / pi | |
REAL, PARAMETER :: F6DPIM9 = 1.0E-9 * f6dpi ! 1.0E-9 = Kg/ug | |
real, parameter :: gpkg = 1.0e03 ! g/kg | |
character( 16 ) :: pname = 'GET_ACAERO_EMIS ' | |
character( 16 ) :: vname | |
character( 96 ) :: xmsg | |
integer c, r, j, m, v, s, l | |
integer, save :: wstep = 0 ! local write counter | |
integer :: mdate, mtime ! diagnostic file write date&time | |
C Third moment emissions rates [m3/m3/s] | |
REAL( 8 ) :: EMISM3( N_MODE ) | |
logical, save :: firstime = .true. | |
C----------------------------------------------------------------------- | |
if ( firstime ) then | |
firstime = .false. | |
if ( .not. desc3 ( ACAERO_EMIS_1 ) ) then | |
xmsg = 'Could not get ' // ACAERO_EMIS_1 //' file description' | |
call m3exit( pname, jdate, jtime, xmsg, xstat1 ) | |
end if | |
end if | |
ac_embuff = 0.0 | |
do s = 1, n_emis_pm | |
IDX = PMEM_MAP( s ) | |
IF ( .NOT. INTERPX( ACAERO_EMIS_1, AEROSPC( IDX )%EMIS, PNAME, | |
& STARTCOL,ENDCOL, STARTROW,ENDROW, 1,ac_emlyr, | |
& JDATE, JTIME, ac_embuff( 1,1,1,S ) ) ) THEN | |
XMSG = 'Could not read ' | |
& // TRIM( AEROSPC( IDX )%EMIS ) // ' from ' // ACAERO_EMIS_1 | |
CALL M3WARN ( PNAME, JDATE, JTIME, XMSG ) | |
END IF | |
end do | |
C all mass emission | |
acaerooutm = 0.0 | |
acaerooutn = 0.0 | |
acaeroouts = 0.0 | |
do l = 1, ac_emlyr | |
do r = 1, my_nrows | |
do c = 1, my_ncols | |
do s = 1, n_emis_pm | |
ac_pm_em( s ) = ac_embuff( c,r,l,s ) | |
end do ! species | |
C *** Calculate scaling factor for converting mass emissions into [ug/m3/s] | |
C note: RJACM converts grid heights from sigma coordinates to meters | |
C Also calculate scaling factors for converting to molar-mixing-ratio units | |
GSFAC = ac_convem_pm( l ) * rjacm( c,r,l ) | |
do s = 1, n_emis_pm | |
C *** Calculate speciated mass emission rates for fine aerosol [ug/m3/s] | |
ac_em_rate (s) = ac_pm_em(s) * GSFAC | |
end do ! species | |
C *** Calculate emissions rate for third moments [m3/m3/s] of each mode | |
C (excluding sea salt), as in Equation 7a of Binkowski & Roselle (2003). | |
ac_em_part = 0.0 ! array assignment | |
EMISM3 = 0.0 ! array assignment | |
do s = 1, n_emis_pm | |
IDX = PMEM_MAP( s ) | |
do n = 1, n_mode - 1 | |
ac_em_part( IDX,n ) = esplit_f( n )* ac_em_rate( s ) | |
EMISM3( n ) = EMISM3( n ) + ac_em_part( IDX,n ) | |
& * ( F6DPIM9 / AEROSPC( IDX )%DENSITY ) | |
end do ! mode | |
end do ! species | |
do n = 1, n_mode | |
do s = 1, n_aerospc | |
acaerooutm( s,n,c,r,l) = ac_em_part(s,n) | |
end do ! mode | |
end do ! species | |
C Mode-specific emission rates of particle number [1/m**3/s] | |
acaerooutn( 1,c,r,l ) = EMISM3(1) * factnumi | |
acaerooutn( 2,c,r,l ) = EMISM3(2) * factnumj | |
acaerooutn( 3,c,r,l ) = 0.0 | |
C Mode-specific dry surface area emission rates [m2/m**3/s]. | |
C 2nd moment multiplied by PI to obtain the surface area emissions rate. | |
acaeroouts( 1,c,r,l ) = EMISM3(1) * factsrfi | |
acaeroouts( 2,c,r,l ) = EMISM3(2) * factsrfj | |
acaeroouts( 3,c,r,l ) = 0.0 | |
end do ! col | |
end do ! row | |
end do ! layer | |
end subroutine get_acaero_emis | |
c----------------------------------------------------------------------- | |
end module acaero_emis |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment