#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" #include "output.inc" ! !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !MODULE: EMISSION ! ! !DESCRIPTION: wrappers around various emissions (init/declare/apply/done) ! routines, needed for TM5 CBM4 version. ! Also hold emissions budget variables. ! !\\ !\\ ! !INTERFACE: ! MODULE EMISSION ! ! !USES: ! USE GO, ONLY : gol, goErr, goPr USE GO, ONLY : GO_Timer_Def, GO_Timer_End, GO_Timer_Start USE TM5_DISTGRID, ONLY : dgrid, Get_DistGrid, scatter, gather USE dims, ONLY : nregions, okdebug USE emission_data, ONLY : plandr, emis2D #ifdef with_budgets USE emission_data, ONLY : budemi_dat, budemi_data, sum_emission USE budget_global, ONLY : nbud_vg,nbudg USE chem_param, ONLY : ntracet #endif use emission_co2 , only : Emission_CO2_Init , Emission_CO2_Done , emission_co2_declare , emission_co2_apply IMPLICIT NONE PRIVATE ! ! !PUBLIC MEMBER FUNCTIONS: ! PUBLIC :: Emission_Init ! allocate/init budget var; call other emiss-related init PUBLIC :: Emission_Done ! gather/write final budget PUBLIC :: Declare_Emission ! allocate emiss var (new run), read emiss data (new month) PUBLIC :: Emission_Apply ! ! ! !PRIVATE DATA MEMBERS: ! ! These budgets should be communicated to chemistry #ifdef with_budgets REAL, DIMENSION(nbudg, nbud_vg, ntracet) :: budemig REAL, DIMENSION(nbudg, nbud_vg, ntracet), PUBLIC :: budemig_all ! for buggy MPI (see budget_global.F90) #endif integer :: itim_appl, itim_co2 CHARACTER(len=*), PARAMETER :: mname = 'emission' ! ! !REVISION HISTORY: ! 20 Aug 2010 - A. Strunk - Adapted to AR5 emissions + various other changes. ! 27 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! 14 May 2014 - T. van Noije- created from chemistry version ! production of CO2 from other C emissions is ! currently not included ! ! - NO MORE REVISION HISTORY. This is gathered from repository from now on. !EOP !---------------------------------------------------------------------- CONTAINS !---------------------------------------------------------------------- ! TM5 ! !---------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_INIT ! ! !DESCRIPTION: Initialise emission fields and parameters by reading ! the rc-file. Allocate and initialize budget variables. !\\ !\\ ! !INTERFACE: ! SUBROUTINE EMISSION_INIT( status ) ! ! !USES: ! use GO, only : TrcFile, Init, Done, ReadRc use meteodata, only : Set, oro_dat use dims, only : iglbsfc, nregions, lm use dims, only : idate, ndyn_max, tref use global_data, only : inputdir, rcfile use emission_data, only : emis_input_dir use emission_data, only : emis_input_year_co2 use emission_data, only : LCMIP6 use emission_data, only : LAR5 use emission_data, only : emis_input_dir_cmip6 use emission_data, only : emis_input_dir_ar5 use emission_read, only : emission_read_init ! ! !OUTPUT PARAMETERS: ! INTEGER, INTENT(out) :: status ! ! !REVISION HISTORY: ! 27 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------------ !BOC CHARACTER(len=*), PARAMETER :: rname = mname//'/Emission_Init' INTEGER :: region, i1, i2, j1, j2, ntim, lmr, imode TYPE(TrcFile) :: rcF REAL :: dtime integer :: emis_input_year ! ----------------------------------- ! read settings from rcfile ! ----------------------------------- call Init( rcF, rcfile, status ) IF_NOTOK_RETURN(status=1) if (okdebug) then write(gol,*) "EMISS-INFO - running year : ", idate(1) ; call goPr end if ! emission base year (assumption: no run overlapping more than one year) call ReadRc( rcF, 'input.emis.year', emis_input_year, status, default=idate(1) ) IF_ERROR_RETURN(status=1) write(gol,*) 'EMISS-INFO - Anthropogenic emissions'; call goPr write(gol,*) 'EMISS-INFO - requested for the following years:'; call goPr call ReadRc( rcF, 'input.emis.year.co2', emis_input_year_co2, status, default=emis_input_year ) IF_ERROR_RETURN(status=1) write(gol,*) 'EMISS-INFO - CO2: ', emis_input_year_co2; call goPr ! default directory for emissions data is "standard input files" dir emis_input_dir=trim(inputdir) ! directory of each data provider call ReadRc( rcF, 'input.emis.dir.CMIP6', emis_input_dir_cmip6, status, default=emis_input_dir ) IF_ERROR_RETURN(status=1) call ReadRc( rcF, 'input.emis.dir.AR5', emis_input_dir_ar5, status, default=emis_input_dir ) IF_ERROR_RETURN(status=1) ! Flags call ReadRc( rcF, 'use_cmip6', LCMIP6, status, default=.false. ) IF_ERROR_RETURN(status=1) call ReadRc( rcF, 'use_ar5', LAR5, status, default=.false. ) IF_ERROR_RETURN(status=1) ! very basic checks if (count((/ LCMIP6, LAR5 /)) > 1) then write(gol,*) 'ERROR: You use more than one ANTHROPOGENIC inventory'; call goErr status=1; TRACEBACK; return end if ! init providers info call emission_read_init( rcF, status ) ! used by vertical distribution: CALL Set( oro_dat(iglbsfc), status, used=.TRUE. ) ! Allocate data ! ------------- 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( plandr(region)%surf(i1:i2, j1:j2)) ALLOCATE( emis2D(region)%surf(i1:i2, j1:j2)) #ifdef with_budgets ALLOCATE( budemi_dat(region)%emi(i1:i2, j1:j2, nbud_vg, ntracet) ) budemi_dat(region)%emi = 0.0 sum_emission(region) = 0.0 #endif ENDDO ! Done ! ------------- call Done( rcF, status ) IF_NOTOK_RETURN(status=1) ! define timers: call GO_Timer_Def( itim_appl, 'emission appl', status ) IF_NOTOK_RETURN(status=1) call GO_Timer_Def( itim_co2, 'emission co2', status ) IF_NOTOK_RETURN(status=1) status = 0 END SUBROUTINE EMISSION_INIT !EOC !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: EMISSION_DONE ! ! !DESCRIPTION: calculate and write final budgets !\\ !\\ ! !INTERFACE: ! SUBROUTINE EMISSION_DONE( status ) ! ! !USES: ! USE dims, ONLY : nregions, im, jm #ifdef with_budgets USE chem_param, ONLY : ntracet, names USE budget_global, ONLY : budget_file_global, nbud_vg, budg_dat, nbudg, NHAB #ifdef with_hdf4 USE file_hdf, ONLY : THdfFile, TSds USE file_hdf, ONLY : Init, Done, WriteAttribute, WriteData, SetDim #endif USE Dims, ONLY : region_name USE partools, ONLY : isRoot, par_reduce, par_reduce_element #endif ! ! !OUTPUT PARAMETERS: ! INTEGER, INTENT(out) :: status ! ! !REVISION HISTORY: ! 16 Jul 2010 - A. Strunk - ! 27 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------ !BOC CHARACTER(len=*), PARAMETER :: rname = mname//'/Emission_Done' INTEGER :: region, i1, i2, j1, j2 #ifdef with_budgets #ifdef with_hdf4 TYPE(THdfFile) :: io TYPE(TSds) :: sds #endif REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: collect_emissions INTEGER :: nsend,j,i,n,nzone,nzone_v real, dimension(nregions) :: sum_emission_all #endif ! --- begin --------------------------------- #ifdef with_budgets ! add up contribution from all proc DO region = 1, nregions CALL PAR_REDUCE(sum_emission(region), 'SUM', sum_emission_all(region), status) IF_NOTOK_RETURN(status=1) END DO ! Write global budget of tracer #1 IF ( isRoot ) THEN write (gol,'("EMISS-INFO - ----------------------------------------------")'); call goPr write (gol,'("EMISS-INFO - Budget of tracer ",a," (kg) ")') trim(names(1)) ; call goPr write (gol,'("EMISS-INFO - ----------------------------------------------")'); call goPr do region = 1, nregions write (gol,'(A,E13.6)') 'EMISS-INFO - mass emitted : ',sum_emission_all(region); call goPr enddo #ifdef with_hdf4 CALL Init(io, budget_file_global, 'write', status) IF_NOTOK_RETURN(status=1) CALL WriteAttribute(io, 'sum_emission', sum_emission_all, status) IF_NOTOK_RETURN(status=1) #endif budemig = 0.0 END IF ! Gather budgets REG: DO region = 1, nregions CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) if (isRoot) then ALLOCATE(collect_emissions(im(region), jm(region), nbud_vg, ntracet)) else ALLOCATE(collect_emissions(1,1,1,1) ) end if CALL GATHER( dgrid(region), budemi_dat(region)%emi, collect_emissions, 0, status) IF_NOTOK_RETURN(status=1) #ifdef with_hdf4 ! write Not-Horizontally-Aggregated-Budgets IF (isRoot.and.NHAB) THEN CALL Init(Sds,io, 'budemi_dat_'//region_name(region),(/im(region),jm(region),nbud_vg,ntracet/), 'real(4)', status) CALL SetDim(Sds, 0, 'im_'//region_name(region),'longitude', (/(j, j=1,im(region))/), status) CALL SetDim(Sds, 1, 'jm_'//region_name(region),'latitude', (/(j, j=1,jm(region))/), status) CALL SetDim(Sds, 2, 'nbud_vg','vertical layer', (/(j, j=1,nbud_vg)/), status) CALL SetDim(Sds, 3, 'ntracet','tracer number', (/(j, j=1,ntracet)/), status) IF_NOTOK_RETURN(status=1) CALL WriteData(Sds,collect_emissions,status) IF_NOTOK_RETURN(status=1) CALL Done(Sds, status) IF_NOTOK_RETURN(status=1) ENDIF #endif ! horizontally aggregates budgets DO n=1,ntracet DO nzone_v=1,nbud_vg DO j=j1,j2 DO i=i1,i2 nzone = budg_dat(region)%nzong(i,j) budemig(nzone,nzone_v,n) = budemig(nzone,nzone_v,n) + budemi_dat(region)%emi(i,j,nzone_v,n) END DO END DO !j END DO !nzone_v END DO !nt DEALLOCATE( collect_emissions ) DEALLOCATE( budemi_dat(region)%emi ) ENDDO REG CALL PAR_REDUCE_ELEMENT( budemig, 'SUM', budemig_all, status) IF_NOTOK_RETURN(status=1) #ifdef with_hdf4 ! Write horizontally aggregated budget IF ( isRoot ) THEN CALL Init(Sds,io, 'budemi',(/nbudg,nbud_vg,ntracet/), 'real(8)', status) IF_NOTOK_RETURN(status=1) CALL SetDim(Sds, 0, 'nbudg','horizontal region', (/(j, j=1,nbudg)/), status) CALL SetDim(Sds, 1, 'nbud_vg','vertical layer', (/(j, j=1,nbud_vg)/), status) CALL SetDim(Sds, 2, 'ntracet','tracer number', (/(j, j=1,ntracet)/), status) IF_NOTOK_RETURN(status=1) CALL WriteData(Sds,budemig_all,status) IF_NOTOK_RETURN(status=1) CALL Done(Sds, status) IF_NOTOK_RETURN(status=1) CALL Done(io, status) IF_NOTOK_RETURN(status=1) ENDIF #endif #endif /* BUDGETS */ ! call other emission_*_done routines CALL FREE_EMISSION(status) IF_NOTOK_RETURN(status=1) ! ok status = 0 END SUBROUTINE EMISSION_DONE !EOC !--------------------------------------------------------------------------- ! TM5 ! !--------------------------------------------------------------------------- !BOP ! ! !IROUTINE: DECLARE_EMISSION ! ! !DESCRIPTION: Called at run start (init/allocate emiss data) and then at ! beginning of every month to just read in data. ! Called from SS_MONTHLY_UPDATE. !\\ !\\ ! !INTERFACE: ! SUBROUTINE DECLARE_EMISSION( status ) ! ! !USES: ! USE Grid, ONLY : FillGrid USE MDF, ONLY : MDF_Open, MDF_NETCDF, MDF_READ, MDF_Inq_VarID, MDF_Get_Var, MDF_Close USE dims, ONLY : im, jm, lm, newsrun USE dims, ONLY : nregions, iglbsfc, nlat180, nlon360 USE chem_param USE partools, ONLY : isRoot USE global_data, ONLY : emis_data USE meteodata, ONLY : global_lli ! AR5/EDGAR4 use emission_data, only : emis_input_dir ! ! !OUTPUT PARAMETERS: ! INTEGER, INTENT(out) :: status ! ! !REVISION HISTORY: ! 16 Jul 2010 - A. Strunk - Adapted to revised emission_*.F90 routines ! 27 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! - anything that is done only if it's a newrun, and that does not require meteo data, should go in INIT ! !EOP !------------------------------------------------------------------------------ !BOC CHARACTER(len=*), PARAMETER :: rname = mname//'/declare_emission' INTEGER :: region, imode, hid, varid REAL, DIMENSION(:,:), ALLOCATABLE :: pland type(emis_data), dimension(nregions) :: wrk ! --------------------------------------------------------------- ! ** land fraction ! --------------------------------------------------------------- IF( newsrun ) THEN if(isRoot)then ALLOCATE( pland(nlon360,nlat180) ) DO region=1,nregions allocate( wrk(region)%surf(im(region),jm(region)) ) end DO else ALLOCATE( pland(1,1) ) DO region=1,nregions allocate( wrk(region)%surf(1,1)) end DO end if if (isRoot)then CALL MDF_Open( TRIM(emis_input_dir)//'/land/landfraction.nc4', MDF_NETCDF, MDF_READ, hid, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Inq_VarID( hid, 'LANDFRACTION', varid, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Get_Var( hid, varid, pland, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Close( hid, status ) IF_NOTOK_RETURN(status=1) ! coarsen or distribute to zoom regions: DO region = 1, nregions ! convert grid: CALL FillGrid( global_lli(region), 'n', wrk(region)%surf, & global_lli(iglbsfc), 'n', pland, 'area-aver', status ) IF_NOTOK_RETURN(status=1) END DO end if DO region = 1, nregions call scatter( dgrid(region), plandr(region)%surf, wrk(region)%surf, 0, status) IF_NOTOK_RETURN(status=1) DEALLOCATE( wrk(region)%surf ) END DO DEALLOCATE( pland ) ENDIF ! --------------------------------------------------------------- ! ** init each constituent ! --------------------------------------------------------------- ! ** 1st time, initialise emissions IF ( newsrun ) THEN CALL Emission_CO2_Init( status ) IF_NOTOK_RETURN(status=1) END IF ! ** every month, read and re-grid CALL emission_co2_declare( status ) IF_NOTOK_RETURN(status=1) ! ok status = 0 END SUBROUTINE DECLARE_EMISSION !EOC !------------------------------------------------------------------- ! TM5 ! !------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_APPLY ! ! !DESCRIPTION: Call emission_apply methods of constituent modules. ! --> add current emissions to tracers array. !\\ !\\ ! !INTERFACE: ! SUBROUTINE EMISSION_APPLY( region, status ) ! ! !USES: ! USE chem_param ! ! !INPUT PARAMETERS: ! INTEGER, INTENT(in) :: region ! ! !OUTPUT PARAMETERS: ! INTEGER, INTENT(out) :: status ! ! !REVISION HISTORY: ! 16 Jul 2010 - A. Strunk - Adapted to revised emission_*.F90 routines ! 27 Mar 2012 - P. Le Sager - cleanup for lat-lon mpi decomposition ! !EOP !----------------------------------------------------------------- !BOC CHARACTER(len=*), PARAMETER :: rname = mname//'/emission_apply' ! --- begin -------------------------------------- ! start timing: call GO_Timer_Start( itim_appl, status ) IF_NOTOK_RETURN(status=1) IF (okdebug) then WRITE(gol,*) 'start of emission_apply for region:',region ; call goPr END IF ! CO2 emissions call GO_Timer_Start( itim_co2, status ) IF_NOTOK_RETURN(status=1) CALL emission_co2_apply( region, status ) IF_NOTOK_RETURN(status=1) call GO_Timer_End( itim_co2, status ) IF_NOTOK_RETURN(status=1) IF(okdebug) then WRITE(gol,*) 'End of adding emission '; call goPr END IF ! end timing: call GO_Timer_End( itim_appl, status ) IF_NOTOK_RETURN(status=1) ! ok status = 0 END SUBROUTINE EMISSION_APPLY !EOC !---------------------------------------------------------------------------- ! TM5 ! !---------------------------------------------------------------------------- !BOP ! ! !IROUTINE: FREE_EMISSION ! ! !DESCRIPTION: Deallocate space needed to handle the emissions by calling ! *done methods of constituents' modules. !\\ !\\ ! !INTERFACE: ! SUBROUTINE FREE_EMISSION( status ) ! ! !USES: ! USE dims, ONLY : nregions ! ! !OUTPUT PARAMETERS: ! INTEGER, INTENT(out) :: status ! ! !REVISION HISTORY: ! 16 Jul 2010 - A. Strunk - Adapted to revised emission_*.F90 routines ! 27 Mar 2012 - P. Le Sager - Adapted for lon-lat MPI domain decomposition ! !EOP !--------------------------------------------------------------------------- !BOC CHARACTER(len=*), PARAMETER :: rname = mname//'/free_emission' INTEGER :: region, imode ! --- begin ----------------------------------- DO region = 1, nregions DEALLOCATE(plandr(region)%surf) DEALLOCATE(emis2D(region)%surf) ENDDO CALL Emission_CO2_Done( status ) IF_NOTOK_RETURN(status=1) ! done status = 0 END SUBROUTINE FREE_EMISSION !EOC END MODULE EMISSION