! #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_NH3 ! ! !DESCRIPTION: Perform NH3 emissions needed for TM5 CBM4 version. !\\ !\\ ! !INTERFACE: ! MODULE EMISSION_NH3 ! ! !USES: ! use GO, only : gol, goErr, goPr use tm5_distgrid, only : dgrid, get_distgrid, scatter, gather use dims, only : nregions, idate, okdebug use global_types, only : emis_data, d3_data use emission_read, only : used_providers, has_emis implicit none private ! ! !PUBLIC MEMBER FUNCTIONS: ! public :: Emission_NH3_init ! allocate public :: Emission_nh3_declare ! read monthly input public :: Emission_nh3_apply ! distribute & add emissions to tracer array public :: Emission_nh3_done ! deallocate ! ! !PRIVATE DATA MEMBERS: ! character(len=*), parameter :: mname = 'emission_nh3' type(emis_data), dimension(:,:), allocatable :: nh3_emis_2d type(d3_data), dimension(:,:), allocatable :: nh3_emis_3d integer :: nh3_3dsec, nh3_2dsec logical, allocatable :: has_data_3d(:), has_data_2d(:) ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - overhaul for AR5 ! 1 Dec 2011 - Narcisa Banda - added EDGAR 4 ! 25 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ CONTAINS !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_NH3_INIT ! ! !DESCRIPTION: Allocate space needed to handle the emissions !\\ !\\ ! !INTERFACE: ! SUBROUTINE EMISSION_NH3_INIT( status ) ! ! !USES: ! use dims, only : lm use emission_read, only : providers_def, numb_providers, ed42_nsect_nh3 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 ! 25 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Emission_NH3_Init' integer :: region integer :: lmr, lsec, lprov, i1, i2, j1, j2 ! --- begin -------------------------------------- status = 0 if(.not. has_emis) return ! nb of sectors nh3_2dsec = 0 nh3_3dsec = 0 do lprov = 1, numb_providers if (count(used_providers.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 nh3_2dsec = nh3_2dsec + n_ar5_ant_sec*count('NH3'.eq.ar5_cat_ant) + & n_ar5_shp_sec*count('NH3'.eq.ar5_cat_shp) if (LAR5BMB) nh3_2dsec = nh3_2dsec + n_ar5_bmb_sec*count('NH3'.eq.ar5_cat_bmb) nh3_3dsec = nh3_3dsec + n_ar5_air_sec*count('NH3'.eq.ar5_cat_air) elseif (trim(providers_def(lprov)%name) .eq. 'ED42') then nh3_2dsec = nh3_2dsec + ed42_nsect_nh3 ! no 3d sectors in EDGAR 4.2 else nh3_2dsec = nh3_2dsec + providers_def(lprov)%nsect2d nh3_3dsec = nh3_3dsec + providers_def(lprov)%nsect3d endif endif enddo ! basic check if (nh3_2dsec == 0) then write(gol,*) "WARNING - there is no 2D sectors for NH3 !"; call goPr end if if (nh3_3dsec == 0) then write(gol,*) "WARNING - there is no 3D sectors for NH3 !"; call goPr end if allocate( nh3_emis_2d( nregions, nh3_2dsec ) ) allocate( nh3_emis_3d( nregions, nh3_3dsec ) ) allocate( has_data_2d(nh3_2dsec)) ; has_data_2d=.false. allocate( has_data_3d(nh3_3dsec)) ; has_data_3d=.false. ! 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,nh3_2dsec allocate( nh3_emis_2d(region,lsec)%surf(i1:i2,j1:j2) ) end do do lsec=1,nh3_3dsec allocate( nh3_emis_3d(region,lsec)%d3(i1:i2,j1:j2,lmr) ) end do enddo ! ok status = 0 END SUBROUTINE EMISSION_NH3_INIT !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_NH3_DONE ! ! !DESCRIPTION: Free space after handling of the emissions !\\ !\\ ! !INTERFACE: ! SUBROUTINE EMISSION_NH3_DONE( status ) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - v0 ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Emission_NH3_Done' integer :: region, lsec ! --- begin --------------------------------- status = 0 if(.not. has_emis) return do region = 1, nregions do lsec=1,nh3_2dsec deallocate( nh3_emis_2d(region,lsec)%surf ) end do do lsec=1,nh3_3dsec deallocate( nh3_emis_3d(region,lsec)%d3 ) end do end do deallocate( nh3_emis_2d, nh3_emis_3d ) deallocate( has_data_2d, has_data_3d ) status = 0 END SUBROUTINE EMISSION_NH3_DONE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_NH3_DECLARE ! ! !DESCRIPTION: Opens, reads and evaluates input files (per month). ! Provides emissions on 2d/3d-arrays which are then added ! to tracers in routine *apply. !\\ !\\ ! !INTERFACE: ! SUBROUTINE EMISSION_NH3_DECLARE( status ) ! ! !USES: ! use toolbox, only : coarsen_emission use partools, only : isRoot, par_broadcast use dims, only : im, jm, lm, idate, sec_month, nlon360, nlat180, iglbsfc use chem_param, only : xmnh3 use emission_data, only : msg_emis, LAR5BMB ! ---------------- AR5 - EDGAR 4 - ETC. -------------------- use emission_data, only : emis_input_year 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_data, only : emis_input_dir_ed4 use emission_read, only : emission_ar5_regrid_aircraft use emission_read, only : emission_ar5_ReadSector use emission_read, only : emission_macc_ReadSector use emission_read, only : emission_ed4_ReadSector use emission_read, only : emission_gfed_ReadSector use emission_read, only : emission_retro_ReadSector use emission_read, only : sectors_def, numb_sectors use emission_read, only : ar5_dim_3ddata use emission_read, only : ed42_nh3_sectors ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - adapted for AR5 ! 1 Dec 2011 - Narcisa Banda - added EDGAR 4 ! 25 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/emission_nh3_declare' integer :: region, hasData integer,parameter :: add_field=0 integer,parameter :: amonth=2 integer :: imr, jmr, lmr, lsec ! AR5 real,dimension(:,:,:), allocatable :: field3d, field3d2 type(d3_data) :: emis3d, work(nregions) type(emis_data) :: wrk2D(nregions) integer :: seccount2d, seccount3d ! --- begin ----------------------------------------- status = 0 if(.not. has_emis) return write(gol,'(" EMISS-INFO ------------- read NH3 emissions -------------")'); call goPr ! reset arrays do region = 1, nregions do lsec=1,nh3_2dsec nh3_emis_2d(region,lsec)%surf = 0.0 end do do lsec=1,nh3_3dsec nh3_emis_3d(region,lsec)%d3 = 0.0 end do end do ! global arrays for coarsening do region = 1, nregions if (isRoot)then allocate(work(region)%d3(im(region),jm(region),lm(region))) else allocate(work(region)%d3(1,1,1)) end if enddo do region = 1, nregions wrk2D(region)%surf => work(region)%d3(:,:,1) end do ! -------------------------------- ! 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 else allocate( field3d( 1, 1, 1 ) ) end if sec : do lsec = 1, numb_sectors if (count(used_providers.eq.sectors_def(lsec)%prov).eq.0) cycle if ((trim(sectors_def(lsec)%prov).eq.'ED42') .and. (count(ed42_nh3_sectors.eq.sectors_def(lsec)%name) .eq. 0)) cycle if (associated(sectors_def(lsec)%species)) then ! AR5 checks if (count('NH3'.eq.sectors_def(lsec)%species).eq.0) cycle if ((trim(sectors_def(lsec)%catname) .eq. 'biomassburning').and.(.not.LAR5BMB)) cycle endif field3d = 0.0 if( sectors_def(lsec)%f3d ) then seccount3d = seccount3d + 1 else seccount2d = seccount2d + 1 end if if (isRoot) then ! READ select case( trim(sectors_def(lsec)%prov) ) case( 'AR5' ) ! Screen out solvent, waste and shipping sectors for NH3, ! because they are zero in the RCPs ! and not present in the historical files. if (trim(sectors_def(lsec)%name) .ne. 'emiss_slv' .and. & trim(sectors_def(lsec)%name) .ne. 'emiss_wst' .and. & trim(sectors_def(lsec)%name) .ne. 'emiss_shp') then call emission_ar5_ReadSector( 'NH3', emis_input_year, idate(2), lsec, field3d, status ) IF_NOTOK_RETURN(status=1) endif case( 'MACC' ) ! screen out 'nat', 'bio', 'air' sectors if ( ( .not. (trim(sectors_def(lsec)%name) .eq. 'emiss_nat')) .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, 'NH3', emis_input_year, idate(2), & '0.5x0.5_kg.nc', sectors_def(lsec)%name, 'kg / s', field3d, status ) IF_NOTOK_RETURN(status=1) endif case( 'ED41' ) select case(trim(sectors_def(lsec)%name)) case ('1A3b_c_e') call emission_ed4_ReadSector( emis_input_dir_ed4, 'NH3', 'nh3', 2005, idate(2), & lsec, trim(sectors_def(lsec)%prov), 'kg / s', field3d, status ) IF_NOTOK_RETURN(status=1) end select case( 'ED42' ) ! biomass burning (GFED/RETRO/AR5BMB) and transport (ED41) are excluded through ED42_NH3_SECTORS definition call emission_ed4_ReadSector( emis_input_dir_ed4, 'NH3', 'nh3', emis_input_year, idate(2), & lsec, trim(sectors_def(lsec)%prov), 'kg / s', field3d, status ) IF_NOTOK_RETURN(status=1) case('GFEDv3') call emission_gfed_ReadSector( emis_input_dir_gfed, 'nh3', emis_input_year, idate(2), & sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status ) IF_NOTOK_RETURN(status=1) case('RETRO') call emission_retro_ReadSector( emis_input_dir_retro, 'NH3', emis_input_year, idate(2), & sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status ) IF_NOTOK_RETURN(status=1) case('MEGAN') ! ! No biogenic NH3 emissions available in MEGAN ! case('DUMMY') case default write(gol,*) "Error in buidling list of providers USED_PROVIDERS"; 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 NH3 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 NH3 emissions for ",a," ",a," for month ",i2 )') & trim(sectors_def(lsec)%prov), trim(sectors_def(lsec)%name), idate(2) ; call goPr endif ! scale from kg/s to kg/month field3d = field3d * sec_month ! kg / month hasData=1 end if end if call Par_broadcast(hasData, status) IF_NOTOK_RETURN(status=1) if (hasData == 0) then cycle sec else if ( sectors_def(lsec)%f3d ) then has_data_3d(seccount3d)=.true. else has_data_2d(seccount2d)=.true. end if end if ! Distinguish 2d/3d sectors if( sectors_def(lsec)%f3d ) then write(gol,'("EMISS-ERROR - Unexpected 3D data: Uncomment code below ")'); call goErr status=1; TRACEBACK; return ! --------------------------- ! 3d data (AIRCRAFT) ! --------------------------- ! if (isRoot) then ! REGRID ! ! helper array for regridding ! allocate( field3d2( nlon360,nlat180,lm(1) ) ) ; field3d2 = 0.0 ! ! aircraft data: regrid vertically to model layers ! call emission_ar5_regrid_aircraft( iglbsfc, field3d, nlon360, nlat180, ar5_dim_3ddata, lm(1), field3d2, status ) ! IF_NOTOK_RETURN(status=1;deallocate(field3d,field3d2)) ! ! write some numbers ! call msg_emis( amonth, trim(sectors_def(lsec)%prov),sectors_def(lsec)%name, 'NH3', xmnh3, sum(field3d2) ) ! IF_NOTOK_RETURN(status=1;deallocate(field3d,field3d2)) ! ! distribute to nh3_emis in regions ! call Coarsen_Emission( 'NH3 '//trim(sectors_def(lsec)%name), nlon360, nlat180, lm(1), & ! field3d2, work, add_field, status ) ! IF_NOTOK_RETURN(status=1;deallocate(field3d,field3d2)) ! deallocate( field3d2 ) ! end if ! do region = 1, nregions ! call scatter(dgrid(region), nh3_emis_3d(region,seccount3d)%d3, work(region)%d3, 0, status) ! IF_NOTOK_RETURN(status=1) ! end do else ! ar5_sector is 2d ! --------------------------- ! 2d data (Anthropogenic, Ships, Biomassburning) ! --------------------------- if (isRoot) then ! print total & regrid call msg_emis( amonth, trim(sectors_def(lsec)%prov), sectors_def(lsec)%name, 'NH3', xmnh3, sum(field3d(:,:,1)) ) IF_NOTOK_RETURN(status=1) call coarsen_emission( 'NH3 '//sectors_def(lsec)%name, nlon360, nlat180, field3d(:,:,1), wrk2D, add_field, status) IF_NOTOK_RETURN(status=1) end if do region = 1, nregions call scatter(dgrid(region), nh3_emis_2d(region,seccount2d)%surf, work(region)%d3(:,:,1), 0, status) IF_NOTOK_RETURN(status=1) end do end if end do sec ! sectors deallocate( field3d ) do region = 1, nregions if (associated(wrk2D(region)%surf)) nullify(wrk2D(region)%surf) deallocate( work(region)%d3 ) end do ! check sectors found if( seccount2d /= nh3_2dsec ) then write(gol,'(80("-"))') ; call goPr write(gol,'("ERROR: 2d sectors do not equal total number:",i4," /= ",i4," !")') seccount2d, nh3_2dsec ; call goErr write(gol,'(80("-"))') ; call goPr status=1; return end if if( seccount3d /= nh3_3dsec ) then write(gol,'(80("-"))') ; call goPr write(gol,'("ERROR: 3d sectors do not equal total number:",i4," /= ",i4," !")') seccount3d, nh3_3dsec ; call goErr write(gol,'(80("-"))') ; call goPr status=1; return end if status = 0 END SUBROUTINE EMISSION_NH3_DECLARE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_NH3_APPLY ! ! !DESCRIPTION: Take monthly emissions, and ! - split them vertically ! - apply time splitting factors ! - add them up (add_3d) !\\ !\\ ! !INTERFACE: ! SUBROUTINE EMISSION_NH3_APPLY( region, status ) ! ! !USES: ! use dims, only : itaur, nsrce, tref use dims, only : im, jm, lm use datetime, only : tau2date use emission_data, only : emission_vdist_by_sector, LAR5BMB use emission_data, only : do_add_3d, do_add_3d_cycle, bb_cycle use emission_data, only : emis_bb_trop_cycle use chem_param, only : inh3, xmnh3, xmn use emission_read, only : sectors_def, numb_sectors use emission_read, only : ed42_nh3_sectors ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - adapted for new vertical distribution ! 25 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/emission_nh3_apply' integer, dimension(6) :: idater real :: dtime, fraction integer :: imr, jmr, lmr, lsec, i1, i2, j1, j2 integer :: seccount2d, seccount3d type(d3_data) :: emis3d ! --- begin ----------------------------------------- status = 0 if(.not. has_emis) return if( okdebug ) then write(gol,*) 'start of emission_nh3_apply'; call goPr end if call tau2date(itaur(region),idater) dtime=float(nsrce)/(2*tref(region)) !emissions are added in two steps...XYZECCEZYX. if(okdebug) then write(gol,*) 'emission_nh3_apply in region ',region,' at date: ',idater, ' with time step:', dtime ; call goPr endif ! get a working structure for 3d emissions call get_distgrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) allocate( emis3d%d3(i1:i2,j1:j2,lm(region)) ) ; emis3d%d3 = 0.0 ! count 2d and 3d sectors seccount2d = 0 seccount3d = 0 ! cycle over sectors do lsec = 1, numb_sectors if (count(used_providers.eq.sectors_def(lsec)%prov).eq.0) cycle if ((trim(sectors_def(lsec)%prov).eq.'ED42') .and. (count(ed42_nh3_sectors.eq.sectors_def(lsec)%name) .eq. 0)) cycle if (associated(sectors_def(lsec)%species)) then ! AR5 check if (count('NH3'.eq.sectors_def(lsec)%species).eq.0) cycle if ((trim(sectors_def(lsec)%catname) .eq. 'biomassburning').and.(.not.LAR5BMB)) cycle endif ! default: no additional splitting fraction = 1.0 ! ---------------------------------------------------------------------------------------- ! distinguish here between sectors and whether they should have additional splitting ! if( sectors_def(lsec)%catname == 'biomassburning' ) fraction = fraction * bb_frac etc... ! ---------------------------------------------------------------------------------------- ! distinguish between 2d/3d sectors if( sectors_def(lsec)%f3d ) then seccount3d = seccount3d + 1 if (.not.has_data_3d(seccount3d)) cycle emis3d%d3 = nh3_emis_3d(region,seccount3d)%d3 else seccount2d = seccount2d + 1 if (.not.has_data_2d(seccount2d)) cycle emis3d%d3 = 0.0 ! vertically distribute according to sector call emission_vdist_by_sector( sectors_def(lsec)%vdisttype, 'NH3', region, nh3_emis_2d(region,seccount2d), emis3d, status) IF_NOTOK_RETURN(status=1;deallocate(emis3d%d3)) end if ! add dataset according to sector and category if( emis_bb_trop_cycle .and. trim(sectors_def(lsec)%catname) == "biomassburning" ) then call do_add_3d_cycle( region, inh3, i1, j1, emis3d%d3, bb_cycle(region)%scalef, xmnh3, xmnh3, status, fraction ) IF_NOTOK_RETURN(status=1) else call do_add_3d( region, inh3, i1, j1, emis3d%d3, xmnh3, xmnh3, status, fraction ) IF_NOTOK_RETURN(status=1) end if end do deallocate( emis3d%d3 ) if(okdebug) then write(gol,*) 'end of emission_nh3_apply'; call goPr endif ! OK status = 0 end subroutine emission_nh3_apply !EOC END MODULE EMISSION_NH3