Skip to content

Instantly share code, notes, and snippets.

@JiaoyanHuang

JiaoyanHuang/ACAERO_EMIS.F

Last active March 2, 2016 03:14
Show Gist options
  • Save JiaoyanHuang/1717227d868d6ffc2b79 to your computer and use it in GitHub Desktop.
Save JiaoyanHuang/1717227d868d6ffc2b79 to your computer and use it in GitHub Desktop.
Download ZIP
aircraft aerosol emission
Raw
ACAERO_EMIS.F
!------------------------------------------------------------------------!
! 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