! #define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') 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: EMISSION_POM ! ! !DESCRIPTION: data and methods for Particulate Organic Matter (POM) emissions. ! ! POM is a sum of fossil and biofuel, vegetation fire, and SOA emissions. ! ! Dimension properties: ! 1) from fossil fuel = 0.015 um and sigma = 1.59 ! 2) from vegetation fires = 0.04 um sigma = 1.59 ! 3) POM from SOA has the characteristics of secondary ! products from monoterpine oxidation, which probably ! condense on pre-existing particles. Therefore only mass ! is calculated and added to the POM mass field if there ! are particles in the one or more soluble modes, if not ! they are formed in the Aitken soluble mode. ! The distribution of the total condensing SOA mass (M) in ! modes is ! ! Solubility properties (from Stier et al.ACP, 2005): ! 1) POM from fossil fuel emissions is considered 65% soluble ! 2) POM from vegetation fire emissions is considered 65% soluble ! 3) POM from SOA is considered 100% soluble (Kanakidou et al. ACP, 2004) ! !>>> TvN ! The latest revision allows to apply ! different sizes for the emissions from ! different sectors, as well as for ! biofuel emissions. ! Also, different solubility fractions can ! now be assumed for emissions from biomass burning, ! biofuel combustion and fossil fuel combustion. ! See settings in emission_data.F90. ! ! The soluble and insoluble SOA fractions are assumed to ! condense onto particles in the Aitken modes, ! as in the previous version by E. Vignati. ! However, this assumption is questionable ! in regions where biomass burning is important, ! since emissions from vegetation fires ! are now assumed to occur in the accumulation modes. ! (Stier et al. assume condensation ! in the soluble Aitken and accumulation modes, ! but it is not clear how the condensation is/should be distributed ! over the two soluble modes in that case.) ! ! The formation of new particles has been removed. ! Previously, a cut-off was applied based on emis_number ! (i.e. the number of particles emitted in a gridbox per month) ! from all sectors excluding the surrogate SOA emissions. ! This criterion is unrealistic ! in the sense that it does not describe any physical process, ! but was only used to select the cells without primary OA emissions. ! We have tested that the inclusion of such a cutoff for new particle formation ! has only very marginal impacts on the results, and can therefore be removed ! (in any case when starting from reasonable initial concentrations). ! ! Note that in the ECHAM version described by Zhang et al. ! fossil fuel and biofuel emissions are considered ! 100% and 65% soluble, respectively. ! Since these sources are not distinguished in the currently used ! emission data sets, we keep the solubility of both sources ! at 65%, as was done by Stier et al. !<<< TvN ! !\\ !\\ ! !INTERFACE: ! MODULE EMISSION_POM ! ! !USES: ! use GO, only : gol, goErr, goPr use tm5_distgrid, only : dgrid, get_distgrid, scatter, gather use dims, only : nregions, okdebug use global_types, only : emis_data, d3_data use emission_data, only : emis_input_dir_aerocom use emission_data, only : emis_input_year_oc use emission_read, only : used_providers_aer, has_aer_emis implicit none private ! ! !PUBLIC MEMBER FUNCTIONS: ! public :: emission_pom_init ! allocate public :: emission_pom_done ! deallocate public :: emission_pom_declare ! read data ! ! !PRIVATE DATA MEMBERS: ! character(len=*), parameter :: mname = 'emission_pom' type(emis_data), dimension(:,:), allocatable :: pom_emis_2d, pom_bf_emis_2d type(d3_data), dimension(:,:), allocatable :: pom_emis_3d integer :: pom_2dsec, pom_3dsec ! ! !REVISION HISTORY: ! ? ??? 2005 - Elisabetta Vignati - changed for coupling with M7 ! 1 Sep 2010 - Achim Strunk - introduced AR5 emissions ! - reorganised the array structures and ! vertical distribution facility ! - cleaning ! 26 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ CONTAINS !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_POM_INIT ! ! !DESCRIPTION: Allocate space needed to handle the emissions. !\\ !\\ ! !INTERFACE: ! SUBROUTINE EMISSION_POM_INIT( status ) ! ! !USES: ! use dims, only : lm use emission_read, only : providers_def, numb_providers use emission_data, only : LAR5BMB use emission_read, only : n_ar5_ant_sec, n_ar5_shp_sec, n_ar5_air_sec, n_ar5_bmb_sec use emission_read, only : ar5_cat_ant, ar5_cat_shp, ar5_cat_air, ar5_cat_bmb ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - v0 ! 26 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Emission_POM_Init' integer :: region, lsec integer :: lmr, lprov, i1, i2, j1, j2 ! --- begin -------------------------------------- status = 0 if(.not. has_aer_emis) return ! nb of sectors pom_2dsec = 0 pom_3dsec = 0 do lprov = 1, numb_providers if (count(used_providers_aer.eq.providers_def(lprov)%name)/=0) then if (trim(providers_def(lprov)%name) .eq. 'AR5') then ! nb of available sectors in AR5 depends on category pom_2dsec = pom_2dsec + n_ar5_ant_sec*count('OC'.eq.ar5_cat_ant) + & n_ar5_shp_sec*count('OC'.eq.ar5_cat_shp) if (LAR5BMB) pom_2dsec = pom_2dsec + n_ar5_bmb_sec*count('OC'.eq.ar5_cat_bmb) pom_3dsec = pom_3dsec + count('OC'.eq.ar5_cat_air) else pom_2dsec = pom_2dsec + providers_def(lprov)%nsect2d pom_3dsec = pom_3dsec + providers_def(lprov)%nsect3d endif endif enddo allocate( pom_emis_2d( nregions, pom_2dsec ) ) allocate( pom_bf_emis_2d( nregions, pom_2dsec ) ) allocate( pom_emis_3d( nregions, pom_3dsec ) ) ! allocate information arrays (2d and 3d) do region=1,nregions CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) lmr = lm(region) do lsec=1,pom_2dsec allocate( pom_emis_2d(region,lsec)%surf(i1:i2,j1:j2) ) allocate( pom_bf_emis_2d(region,lsec)%surf(i1:i2,j1:j2) ) end do do lsec=1,pom_3dsec allocate( pom_emis_3d(region,lsec)%d3(i1:i2,j1:j2,lmr) ) end do enddo ! ok status = 0 END SUBROUTINE EMISSION_POM_INIT !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_POM_DONE ! ! !DESCRIPTION: Free memory. !\\ !\\ ! !INTERFACE: ! SUBROUTINE EMISSION_POM_DONE( status ) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - v0 ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Emission_POM_Done' integer :: region, lsec ! --- begin -------------------------------------- status = 0 if(.not. has_aer_emis) return do region = 1, nregions do lsec=1,pom_2dsec deallocate( pom_emis_2d(region,lsec)%surf ) deallocate( pom_bf_emis_2d(region,lsec)%surf ) end do do lsec=1,pom_3dsec deallocate( pom_emis_3d(region,lsec)%d3 ) end do end do deallocate( pom_emis_2d ) deallocate( pom_bf_emis_2d ) deallocate( pom_emis_3d ) ! ok status = 0 END SUBROUTINE EMISSION_POM_DONE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_POM_DECLARE ! ! !DESCRIPTION: Opens, reads and evaluates input files (per month). ! Provides emissions on 2d/3d-arrays which are then given ! to emis_number and emis_mass, which are used in ! sedimentation. (no *_apply!) !\\ !\\ ! !INTERFACE: ! SUBROUTINE EMISSION_POM_DECLARE( status ) ! ! !USES: ! use binas, only : pi use partools, only : isRoot, par_broadcast use toolbox, only : coarsen_emission, distribute_emis2D use dims, only : im, jm, lm, idate, sec_month, nlon360, nlat180 use chem_param, only : xmc, sigma_lognormal, pom_density use chem_param, only : mode_aii, mode_ais, mode_acs use MDF, only : MDF_Open, MDF_NETCDF, MDF_Get_Var use MDF, only : MDF_Close, MDF_READ, MDF_Inq_VarID use emission_data, only : emis_mass, emis_number use emission_data, only : rad_emi_ff_sol, rad_emi_ff_insol use emission_data, only : rad_emi_ene_sol, rad_emi_ene_insol use emission_data, only : rad_emi_ind_sol, rad_emi_ind_insol use emission_data, only : rad_emi_tra_sol, rad_emi_tra_insol use emission_data, only : rad_emi_shp_sol, rad_emi_shp_insol use emission_data, only : rad_emi_air_sol, rad_emi_air_insol use emission_data, only : rad_emi_bf_sol, rad_emi_bf_insol use emission_data, only : rad_emi_bb_sol, rad_emi_bb_insol !use emission_data, only : rad_soa use emission_data, only : frac_pom_sol_ff, frac_pom_sol_bf, frac_pom_sol_bb, frac_soa_sol !use emission_data, only : oc2pom use emission_data, only : oc2pom_ff, oc2pom_bf, oc2pom_bb, oc2pom_soa use emission_data, only : msg_emis, emis_temp, emission_vdist_by_sector use emission_data, only : vd_class_name_len ! ---------------- CMIP6 - AR5 - GFED - RETRO - MACC -------------------- use emission_data, only : LAR5BMB use emission_data, only : emis_input_dir_mac use emission_data, only : emis_input_dir_retro use emission_data, only : emis_input_dir_gfed use emission_read, only : emission_ar5_regrid_aircraft use emission_read, only : emission_cmip6_ReadSector use emission_read, only : emission_cmip6bmb_ReadSector use emission_read, only : emission_ar5_ReadSector use emission_read, only : emission_macc_ReadSector use emission_read, only : emission_gfed_ReadSector use emission_read, only : emission_retro_ReadSector use emission_read, only : sector_name_len use emission_read, only : sectors_def, numb_sectors use emission_read, only : ar5_dim_3ddata !RM use mo_aero, only : nsoa ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - revamped for AR5 ! 26 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/emission_pom_declare' integer :: region, hasData integer, parameter :: add_field=0 integer, parameter :: amonth=2 !real :: rad_aver_mass real :: numbscale_exp integer :: imr, jmr, lmr, lsec, i1, i2, j1, j2, l real(kind=4), dimension(:,:,: ), allocatable :: hdfr3 integer :: fid, varid ! hdf related ! --------------------------------------------------------------- real, dimension(:,:), allocatable :: field2d real, dimension(:,:,:), allocatable :: field3d, field3d2 real, dimension(:,:), allocatable :: field2d_bf integer :: mode_sol real :: mass2numb_fact real :: mass2numb_ff_sol, mass2numb_ff_insol real :: mass2numb_ene_sol, mass2numb_ene_insol real :: mass2numb_ind_sol, mass2numb_ind_insol real :: mass2numb_tra_sol, mass2numb_tra_insol real :: mass2numb_shp_sol, mass2numb_shp_insol real :: mass2numb_air_sol, mass2numb_air_insol real :: mass2numb_bf_sol, mass2numb_bf_insol real :: mass2numb_bb_sol, mass2numb_bb_insol real :: mass2numb_nonbf_sol, mass2numb_nonbf_insol real :: oc2pom type(d3_data), dimension(nregions) :: emis3d, work, work_bf, work3d type(d3_data), dimension(nregions) :: frac_bf type(emis_data),dimension(nregions) :: wrk2D, wrk2D_bf, emis_glb integer :: seccount2d, seccount3d character(len=sector_name_len) :: tmpsector character(len=vd_class_name_len) :: tmpvsplit ! --- begin ----------------------------------------- status = 0 if(.not. has_aer_emis) return write(gol,'(" EMISS-INFO ------------- read POM emissions -------------")'); call goPr ! Reset emissions do region = 1, nregions do lsec=1,pom_2dsec pom_emis_2d(region,lsec)%surf = 0.0 pom_bf_emis_2d(region,lsec)%surf = 0.0 end do do lsec=1,pom_3dsec pom_emis_3d(region,lsec)%d3 = 0.0 end do end do ! Allocate work arrays do region = 1, nregions CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) lmr = lm(region) allocate( work3d(region)%d3 (i1:i2,j1:j2, ar5_dim_3ddata) ) ; work3d(region)%d3 = 0.0 allocate( emis3d(region)%d3 (i1:i2,j1:j2, lmr ) ) ; emis3d(region)%d3 = 0.0 allocate( frac_bf(region)%d3 (i1:i2,j1:j2, 1 ) ) end do ! Global arrays for coarsening do region = 1, nregions if (isRoot)then allocate(work(region)%d3(im(region),jm(region),ar5_dim_3ddata)) allocate(work_bf(region)%d3(im(region),jm(region),1)) else allocate(work(region)%d3(1,1,1)) allocate(work_bf(region)%d3(1,1,1)) end if enddo do region = 1, nregions wrk2D(region)%surf => work(region)%d3(:,:,1) wrk2D_bf(region)%surf => work_bf(region)%d3(:,:,1) end do ! ----------------------------- ! get emissions (ATTENTION: THIS IS Organic Carbon! ) write(gol,'(1x,80("-"))') ; call goPr write(gol,*) ' WARNING: OC emissions are used instead of POM !' ; call goPr !write(gol,*) ' masses are transformed using constant factor oc2pom' ; call goPr write(gol,*) ' masses are transformed using constant factors' ; call goPr write(gol,*) ' for vegetation fires and other sources' ; call goPr write(gol,'(1x,80("-"))') ; call goPr ! -------------------------------- ! do a loop over available sectors ! -------------------------------- ! count 2d and 3d sectors seccount2d = 0 seccount3d = 0 ! always allocate here 3d data set (for 2d sectors it will be filled in first layer only) if (isRoot) then allocate( field3d( nlon360, nlat180, ar5_dim_3ddata ) ) ; field3d = 0.0 allocate( field2d_bf( nlon360, nlat180 ) ) ; field2d_bf = 0.0 else allocate( field3d( 1, 1, 1 ) ) allocate( field2d_bf( 1, 1 ) ) end if ! mass to number conversion factors for the relevant modes numbscale_exp = EXP(1.5*(LOG(sigma_lognormal(mode_aii)))**2) mass2numb_fact = 3./(4.*pi*(numbscale_exp**3)*pom_density) mass2numb_ff_insol = mass2numb_fact/(rad_emi_ff_insol**3) mass2numb_ene_insol = mass2numb_fact/(rad_emi_ene_insol**3) mass2numb_ind_insol = mass2numb_fact/(rad_emi_ind_insol**3) mass2numb_tra_insol = mass2numb_fact/(rad_emi_tra_insol**3) mass2numb_shp_insol = mass2numb_fact/(rad_emi_shp_insol**3) mass2numb_air_insol = mass2numb_fact/(rad_emi_air_insol**3) mass2numb_bf_insol = mass2numb_fact/(rad_emi_bf_insol**3) mass2numb_bb_insol = mass2numb_fact/(rad_emi_bb_insol**3) numbscale_exp = EXP(1.5*(LOG(sigma_lognormal(mode_ais)))**2) mass2numb_fact = 3./(4.*pi*(numbscale_exp**3)*pom_density) mass2numb_ff_sol = mass2numb_fact/(rad_emi_ff_sol**3) mass2numb_ene_sol = mass2numb_fact/(rad_emi_ene_sol**3) mass2numb_ind_sol = mass2numb_fact/(rad_emi_ind_sol**3) mass2numb_tra_sol = mass2numb_fact/(rad_emi_tra_sol**3) mass2numb_shp_sol = mass2numb_fact/(rad_emi_shp_sol**3) mass2numb_air_sol = mass2numb_fact/(rad_emi_air_sol**3) numbscale_exp = EXP(1.5*(LOG(sigma_lognormal(mode_acs)))**2) mass2numb_fact = 3./(4.*pi*(numbscale_exp**3)*pom_density) mass2numb_bf_sol = mass2numb_fact/(rad_emi_bf_sol**3) mass2numb_bb_sol = mass2numb_fact/(rad_emi_bb_sol**3) sec : do lsec = 1, numb_sectors if (count(used_providers_aer.eq.sectors_def(lsec)%prov).eq.0) cycle if (associated(sectors_def(lsec)%species)) then if (count('OC'.eq.sectors_def(lsec)%species).eq.0) cycle if ((trim(sectors_def(lsec)%catname) .eq. 'biomassburning').and.(.not.LAR5BMB)) cycle endif field3d = 0.0 field2d_bf = 0.0 if( sectors_def(lsec)%f3d ) then seccount3d = seccount3d + 1 else seccount2d = seccount2d + 1 end if if( trim(sectors_def(lsec)%catname) == 'biomassburning' ) then oc2pom = oc2pom_bb else oc2pom = oc2pom_ff endif if (isRoot) then ! READ select case( trim(sectors_def(lsec)%prov) ) case( 'CMIP6' ) call emission_cmip6_ReadSector( 'OC', emis_input_year_oc, idate(2), lsec, field3d, status, field2d_bf ) IF_NOTOK_RETURN(status=1;deallocate(field3d,field2d_bf)) case( 'CMIP6BMB' ) call emission_cmip6bmb_ReadSector( 'OC', emis_input_year_oc, idate(2), lsec, field3d, status ) IF_NOTOK_RETURN(status=1;deallocate(field3d)) case( 'AR5' ) ! Screen out agricultural and solvent sectors for OC, ! because they are zero in the RCPs ! and not present in the historical files. if (trim(sectors_def(lsec)%name) .ne. 'emiss_agr' .and. & trim(sectors_def(lsec)%name) .ne. 'emiss_slv') then call emission_ar5_ReadSector( 'OC', emis_input_year_oc, idate(2), lsec, field3d, status ) IF_NOTOK_RETURN(status=1;deallocate(field3d)) endif case( 'MACC' ) ! Screen out 'soil', 'nat', 'oc', bio', and 'air' since they are not available for OC. if ( ( .not. (trim(sectors_def(lsec)%name) .eq. 'emiss_soil')) .and. & ( .not. (trim(sectors_def(lsec)%name) .eq. 'emiss_nat') ) .and. & ( .not. (trim(sectors_def(lsec)%name) .eq. 'emiss_oc') ) .and. & ( .not. (trim(sectors_def(lsec)%name) .eq. 'emiss_bio' ) ) .and. & ( .not. (trim(sectors_def(lsec)%name) .eq. 'emiss_air') ) ) then call emission_macc_ReadSector( emis_input_dir_mac, 'OC', emis_input_year_oc, idate(2), & '0.5x0.5_kg.nc', sectors_def(lsec)%name, 'kg / s', field3d, status ) IF_NOTOK_RETURN(status=1;deallocate(field3d)) end if case('GFEDv3') call emission_gfed_ReadSector( emis_input_dir_gfed, 'oc', emis_input_year_oc, idate(2), & 'GFED', 'kg / s', field3d(:,:,1), status ) IF_NOTOK_RETURN(status=1) case('RETRO') call emission_retro_ReadSector( emis_input_dir_retro, 'OC', emis_input_year_oc, idate(2), & sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status ) IF_NOTOK_RETURN(status=1) case default write(gol,*) "Error in list of providers for POM"; call goErr status=1; TRACEBACK; return end select ! nothing found??? if( sum(field3d) < 100.*TINY(1.0) ) then if (okdebug) then write(gol,'("EMISS-INFO - no POM emissions found for ",a," ",a," for month ",i2 )') & trim(sectors_def(lsec)%prov), trim(sectors_def(lsec)%name), idate(2) ; call goPr endif hasData=0 else if (okdebug) then write(gol,'("EMISS-INFO - found POM emissions for ",a," ",a," for month ",i2 )') & trim(sectors_def(lsec)%prov), trim(sectors_def(lsec)%name), idate(2) ; call goPr endif ! convert from OC to POM if( trim(sectors_def(lsec)%catname) == 'biomassburning' ) then field3d = field3d * oc2pom_bb else field3d = field3d * oc2pom_ff ! correct biofuel contribution in surface level field3d(:,:,1) = field3d(:,:,1) + field2d_bf(:,:) * (oc2pom_bf - oc2pom_ff) field2d_bf = field2d_bf * oc2pom_bf endif ! scale from kg/s to kg/month field3d = field3d * sec_month ! kg / month field2d_bf = field2d_bf * sec_month hasData=1 end if end if call Par_broadcast(hasData, status) IF_NOTOK_RETURN(status=1) if (hasData == 0) cycle sec ! distinguish 2d/3d sectors if( sectors_def(lsec)%f3d ) then ! --------------------------------------- ! 3d data (AIRCRAFT), available for CMIP6 ! --------------------------------------- if (isRoot) then ! write some numbers call msg_emis( amonth, trim(sectors_def(lsec)%prov), sectors_def(lsec)%name, 'POM', oc2pom_ff*xmc, sum(field3d) ) ! distribute to work arrays in regions call Coarsen_Emission( 'POM '//trim(sectors_def(lsec)%name), nlon360, nlat180, ar5_dim_3ddata, & field3d, work, add_field, status ) IF_NOTOK_RETURN(status=1) end if ! scatter, sum up on target array do region = 1, nregions call scatter(dgrid(region), work3d(region)%d3, work(region)%d3, 0, status) IF_NOTOK_RETURN( status=1 ) CALL GET_DISTGRID( dgrid(region), I_STRT=i1, J_STRT=j1) ! aircraft data: regrid vertically to model layers call emission_ar5_regrid_aircraft( region, i1, j1, work3d(region)%d3, emis3d(region)%d3, status ) IF_NOTOK_RETURN( status=1 ) pom_emis_3d(region,seccount3d)%d3 = pom_emis_3d(region,seccount3d)%d3 + emis3d(region)%d3 end do ! add to emis target arrays do region = 1, nregions emis_mass (region,mode_aii)%d4(:,:,:,2) = & emis_mass (region,mode_aii)%d4(:,:,:,2) + & (1.-frac_pom_sol_ff) * pom_emis_3d(region,seccount3d)%d3(:,:,:) emis_mass (region,mode_ais)%d4(:,:,:,3) = & emis_mass (region,mode_ais)%d4(:,:,:,3) + & ( frac_pom_sol_ff) * pom_emis_3d(region,seccount3d)%d3(:,:,:) emis_number(region,mode_aii)%d4(:,:,:,2) = & emis_number(region,mode_aii)%d4(:,:,:,2) + & (1.-frac_pom_sol_ff) * pom_emis_3d(region,seccount3d)%d3(:,:,:) * mass2numb_air_insol emis_number(region,mode_ais)%d4(:,:,:,3) = & emis_number(region,mode_ais)%d4(:,:,:,3) + & ( frac_pom_sol_ff) * pom_emis_3d(region,seccount3d)%d3(:,:,:) * mass2numb_air_sol end do else ! 2d sector ! --------------------------- ! 2d data (Anthropogenic, Ships, Biomassburning) ! --------------------------- if (isRoot) then ! print total & regrid ! in the following print statement ! a possible difference between oc2pom_bf and oc2pom_ff ! is not accounted for, ! but this has no further consequences call msg_emis( amonth, trim(sectors_def(lsec)%prov), sectors_def(lsec)%name, 'POM', oc2pom*xmc, & sum(field3d(:,:,1)) ) IF_NOTOK_RETURN(status=1) call coarsen_emission( 'POM '//sectors_def(lsec)%name, nlon360, nlat180, field3d(:,:,1), wrk2D, add_field, status) IF_NOTOK_RETURN(status=1) call coarsen_emission( 'POM solid biofuel '//sectors_def(lsec)%name, nlon360, nlat180, field2d_bf(:,:), wrk2D_bf, add_field, status) IF_NOTOK_RETURN(status=1) end if ! get temporary array for vertical distribution do region = 1, nregions call scatter(dgrid(region), pom_emis_2d(region,seccount2d)%surf, work(region)%d3(:,:,1), 0, status) IF_NOTOK_RETURN(status=1) call scatter(dgrid(region), pom_bf_emis_2d(region,seccount2d)%surf, work_bf(region)%d3(:,:,1), 0, status) IF_NOTOK_RETURN(status=1) emis3d(region)%d3 = 0.0 ! do vertical distribution call emission_vdist_by_sector( sectors_def(lsec)%vdisttype, 'POM', region, pom_emis_2d(region,seccount2d), & emis3d(region), status ) IF_NOTOK_RETURN(status=1) if( trim(sectors_def(lsec)%catname) == 'biomassburning' ) then ! add to emis target arrays emis_mass (region,mode_aii)%d4(:,:,:,2) = & emis_mass (region,mode_aii)%d4(:,:,:,2) + emis3d(region)%d3(:,:,:) * & (1.-frac_pom_sol_bb) emis_number(region,mode_aii)%d4(:,:,:,2) = & emis_number(region,mode_aii)%d4(:,:,:,2) + emis3d(region)%d3(:,:,:) * & (1.-frac_pom_sol_bb) * mass2numb_bb_insol emis_mass (region,mode_acs)%d4(:,:,:,3) = & emis_mass (region,mode_acs)%d4(:,:,:,3) + emis3d(region)%d3(:,:,:) * & frac_pom_sol_bb emis_number(region,mode_acs)%d4(:,:,:,3) = & emis_number(region,mode_acs)%d4(:,:,:,3) + emis3d(region)%d3(:,:,:) * & frac_pom_sol_bb * mass2numb_bb_sol else ! currenly only for CMIP6 sector names select case( trim(sectors_def(lsec)%name) ) case ('ENE') mass2numb_nonbf_sol = mass2numb_ene_sol mass2numb_nonbf_insol = mass2numb_ene_insol case ('IND') mass2numb_nonbf_sol = mass2numb_ind_sol mass2numb_nonbf_insol = mass2numb_ind_insol case ('TRA') mass2numb_nonbf_sol = mass2numb_tra_sol mass2numb_nonbf_insol = mass2numb_tra_insol case ('SHP') mass2numb_nonbf_sol = mass2numb_shp_sol mass2numb_nonbf_insol = mass2numb_shp_insol case default mass2numb_nonbf_sol = mass2numb_ff_sol mass2numb_nonbf_insol = mass2numb_ff_insol end select if ( trim(sectors_def(lsec)%prov) == 'CMIP6' ) then ! calculate mass fraction related to solid biofuel where ( pom_emis_2d(region,seccount2d)%surf(:,:) > tiny(0.0) ) frac_bf(region)%d3(:,:,1) = pom_bf_emis_2d(region,seccount2d)%surf(:,:) / & pom_emis_2d(region,seccount2d)%surf(:,:) elsewhere frac_bf(region)%d3(:,:,1) = 0.0 endwhere ! for safety, prevent fractions larger than unity. where (frac_bf(region)%d3(:,:,1) > 1.0 ) frac_bf(region)%d3(:,:,1) = 1.0 endwhere else frac_bf(region)%d3(:,:,1) = 0.0 endif do l = 1, lmr emis_mass (region,mode_aii)%d4(:,:,l,2) = & emis_mass (region,mode_aii)%d4(:,:,l,2) + emis3d(region)%d3(:,:,l) * & ( (1.-frac_bf(region)%d3(:,:,1)) * (1.-frac_pom_sol_ff) + & frac_bf(region)%d3(:,:,1) * (1.-frac_pom_sol_bf) ) emis_number(region,mode_aii)%d4(:,:,l,2) = & emis_number(region,mode_aii)%d4(:,:,l,2) + emis3d(region)%d3(:,:,l) * & ( (1.-frac_bf(region)%d3(:,:,1)) * (1.-frac_pom_sol_ff) * mass2numb_nonbf_insol + & frac_bf(region)%d3(:,:,1) * (1.-frac_pom_sol_bf) * mass2numb_bf_insol ) emis_mass (region,mode_ais)%d4(:,:,l,3) = & emis_mass (region,mode_ais)%d4(:,:,l,3) + emis3d(region)%d3(:,:,l) * & ( (1.-frac_bf(region)%d3(:,:,1)) * frac_pom_sol_ff ) emis_number(region,mode_ais)%d4(:,:,l,3) = & emis_number(region,mode_ais)%d4(:,:,l,3) + emis3d(region)%d3(:,:,l) * & ( (1.-frac_bf(region)%d3(:,:,1)) * frac_pom_sol_ff * mass2numb_nonbf_sol ) emis_mass (region,mode_acs)%d4(:,:,l,3) = & emis_mass (region,mode_acs)%d4(:,:,l,3) + emis3d(region)%d3(:,:,l) * & ( frac_bf(region)%d3(:,:,1) * frac_pom_sol_bf ) emis_number(region,mode_acs)%d4(:,:,l,3) = & emis_number(region,mode_acs)%d4(:,:,l,3) + emis3d(region)%d3(:,:,l) * & ( frac_bf(region)%d3(:,:,1) * frac_pom_sol_bf * mass2numb_bf_sol ) enddo endif end do end if end do sec ! sectors deallocate( field3d ) deallocate( field2d_bf ) do region = 1, nregions if (associated(wrk2D(region)%surf)) nullify(wrk2D(region)%surf) if (associated(wrk2D_bf(region)%surf)) nullify(wrk2D_bf(region)%surf) deallocate( work(region)%d3 ) deallocate( work_bf(region)%d3 ) deallocate( work3d(region)%d3 ) deallocate( frac_bf(region)%d3 ) end do ! check sectors found if( seccount2d /= pom_2dsec ) then write(gol,'(80("-"))') ; call goPr write(gol,'("ERROR: 2d sectors do not equal total number:",i4," /= ",i4," !")') seccount2d, pom_2dsec ; call goErr write(gol,'(80("-"))') ; call goPr status=1; return end if if( seccount3d /= pom_3dsec ) then write(gol,'(80("-"))') ; call goPr write(gol,'("ERROR: 3d sectors do not equal total number:",i4," /= ",i4," !")') seccount3d, pom_3dsec ; call goErr write(gol,'(80("-"))') ; call goPr status=1; return end if ! Old AEROCOM SOA emissions available only when nsoa == 0 ! Otherwise the production of SOA is calculated from isoprene and monoterpene ! through chemical reactions in ebischeme.f90. if (nsoa .EQ. 0) then ! --------------------------------------------------- ! SOA condenses onto existing particles in the Aitken soluble mode ! --------------------------------------------------- ! select case( emis_input_provider ) ! case( 'AR5','MACC' ) write(gol,'(80("-"))') ; call goPr write(gol,*) ' WARNING: SOA emissions are employed in boxes ' ; call goPr write(gol,*) ' with too few particles in mode_ais! ' ; call goPr write(gol,*) ' Masses of SOA are neglected !!! ' ; call goPr write(gol,'(80("-"))') ; call goPr ! end select tmpsector = 'anthropogenic' ! assume SOA emissions to be vertically split due to this class: tmpvsplit = 'nearsurface' ! work arrays if (isRoot) then do region = 1, nregions allocate(emis_glb(region)%surf(im(region), jm(region))) end do else do region = 1, nregions allocate(emis_glb(region)%surf(1,1)) end do end if ! Read if (isRoot) then allocate( hdfr3(nlon360,nlat180,1) ) ; hdfr3 = 0.0 allocate( field2d(nlon360,nlat180) ) ; field2d = 0.0 ! FIXME - this is untested, SOA.nc4 needs to be created - but not needed in EC-Earth call MDF_Open( trim(emis_input_dir_aerocom)//'/SOA.nc4', MDF_NETCDF, MDF_READ, fid, status) IF_NOTOK_RETURN(status=1) call MDF_Inq_VarID( fid, 'FIELD', varid, status) IF_NOTOK_RETURN(status=1) call MDF_Get_Var( fid, varid, hdfr3, status, start=(/1,1,idate(2)/), count=(/nlon360,nlat180,1/)) IF_NOTOK_RETURN(status=1) CALL MDF_Close( fid, status ) IF_NOTOK_RETURN(status=1) ! paste to 2d field field2d = hdfr3(:,:,1) !>>> TvN ! SOA emissions are provided in Tg POM, and are based on ! the assumption that 15% of the emitted terpene mass is converted to SOA ! (Kanakidou et al., ACP, 2005; Dentener et al., ACP, 2006). ! However, it seems that such an assumption ! is inconsistent with a POM to OC ratio of about 2.0 for freshly formed SOA ! (Aiken et al., Environ. Sci. Technol., 2008). ! The VOCs with the largest potential for SOA formation ! are mono-terpenes (C10H16), which account for 40-80% ! of the overall terpene emissions, and seqsquiterpenes (C15H24), ! which have a 100% SOA yield. ! The full weight to carbon weight ratios for these components ! are 1.133 and 1.16, respectively. ! Therefore, it seems reasonable to increase the SOA emissions ! as provided in the AeroCom input file ! by a factor of about ocpom_soa/1.15. ! Increasing the SOA emissions from AeroCom, ! which amount to 19.2 Tg/yr, is also consistent with AR5, ! where 20 Tg/yr is assumed to be a lower bound. field2d = field2d * oc2pom_soa/1.15 !<<< TvN deallocate(hdfr3) ! prompt some numbers call msg_emis( amonth, 'AEROCOM', 'SOA-'//trim(tmpsector), 'POM', oc2pom_soa*xmc, sum(field2d) ) IF_NOTOK_RETURN(status=1;deallocate(field2d)) ! initialise/reset emis_temp(regions) do region = 1, nregions emis_temp(region)%surf = 0.0 end do ! regrid from field2d to emis_glb call coarsen_emission( 'SOA '//trim(tmpsector), nlon360, nlat180, field2d, emis_glb, add_field, status) IF_NOTOK_RETURN(status=1;deallocate(field2d)) deallocate( field2d ) end if !>>> TvN ! EV total particle number emitted in a gridbox !rad_aver_mass = rad_soa*EXP(1.5*(LOG(sigma_lognormal(mode_ais)))**2) !mass_to_numb = 3./(4.*pi*(rad_aver_mass**3)*pom_density) ! Scatter, vertical distribution ! ------------------------------ do region = 1, nregions call scatter(dgrid(region), emis_temp(region)%surf, emis_glb(region)%surf, 0, status) IF_NOTOK_RETURN(status=1) emis3d(region)%d3 = 0.0 call emission_vdist_by_sector( tmpvsplit, 'SOA', region, emis_temp(region), emis3d(region), status ) IF_NOTOK_RETURN(status=1) ! -------------- ! Comment by EV: ! Let the SOA condense on the existing aerosols. Not emit the aerosol numbers. ! Only if there are insufficient aerosols (<1.0), we emit the aerosol numbers. ! add masses !emis_mass(region,mode_ais)%d4(:,:,:,3) = & ! emis_mass(region,mode_ais)%d4(:,:,:,3) + emis3d(region)%d3(:,:,:) ! TB use the new compounds available to track the emissions from SOA.hdf ! emis_mass(region,mode_ais)%D4(:,:,:,4) put mass as soa-compound ! emis_mass(region,mode_aii)%D4(:,:,:,3) put mass as soa-compound emis_mass(region,mode_ais)%d4(:,:,:,4) = & emis_mass(region,mode_ais)%d4(:,:,:,4) + & frac_soa_sol * emis3d(region)%d3(:,:,:) emis_mass(region,mode_aii)%d4(:,:,:,3) = & emis_mass(region,mode_aii)%d4(:,:,:,3) + & (1.-frac_soa_sol) * emis3d(region)%d3(:,:,:) ! add numbers where appropriate ! where(emis_number(region,mode_ais)%d4(:,:,:,3) < 1.) !emis_number(region,mode_ais)%d4(:,:,:,3) = & ! emis_number(region,mode_ais)%d4(:,:,:,3) + emis3d(region)%d3(:,:,:) * mass_to_numb ! emis_number(region,mode_ais)%d4(:,:,:,3) = & ! emis_number(region,mode_ais)%d4(:,:,:,3) + & ! frac_soa_sol * emis3d(region)%d3(:,:,:) * mass_to_numb ! endwhere ! where(emis_number(region,mode_aii)%d4(:,:,:,2) < 1.) ! emis_number(region,mode_aii)%d4(:,:,:,2) = & ! emis_number(region,mode_aii)%d4(:,:,:,2) + & ! (1.-frac_soa_sol) * emis3d(region)%d3(:,:,:) * mass_to_numb ! endwhere !<<< TvN enddo ! Done do region = 1, nregions deallocate(emis_glb(region)%surf) deallocate( emis3d(region)%d3 ) end do end if status = 0 END SUBROUTINE EMISSION_POM_DECLARE !EOC END MODULE EMISSION_POM