! #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_CH4 ! ! !DESCRIPTION: Perform CH4 emissions needed for TM5 CBM4 version. ! For AR5 runs, fix the CH4 concentration and keep track ! of the 3D field where CH4 is added ! (with_ch4_emissions). !\\ !\\ ! !INTERFACE: ! MODULE EMISSION_CH4 ! ! !USES: ! use GO, only : gol, goPr, goErr use tm5_distgrid, only : dgrid, get_distgrid, scatter, scatter_i_band use partools, only : isRoot, par_broadcast use dims, only : nregions, idate, okdebug use global_types, only : emis_data, d3_data, d2_data, d23_data use emission_data, only : LCMIP6_CH4 use emission_data, only : emis_ch4_single, emis_ch4_fix3d use emission_data, only : emis_ch4_fixed_ppb use emission_read, only : used_providers_ch4, has_ch4_emis IMPLICIT NONE PRIVATE ! ! !PUBLIC MEMBER FUNCTIONS: ! public :: Emission_CH4_Init public :: Emission_CH4_Done public :: emission_ch4_declare public :: emission_ch4_apply ! ! !PRIVATE DATA MEMBERS: ! character(len=*), parameter :: mname = 'emission_ch4' #ifdef with_ch4_emis type(emis_data), dimension(:,:), allocatable :: ch4_emis_2d type(d3_data), dimension(:,:), allocatable :: ch4_emis_3d integer :: ch4_2dsec, ch4_3dsec #endif logical, allocatable :: has_data_3d(:), has_data_2d(:) type(d2_data), target :: zch4(nregions), tmp2d(nregions) type(d23_data), target :: wrk_c(nregions) type(d2_data), target :: zch4_p(nregions) ! previous month type(d23_data), target :: wrk_p(nregions) type(d2_data), target :: zch4_n(nregions) ! next month type(d23_data), target :: wrk_n(nregions) ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - overhaul for AR5 ! 28 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition ! - made zch4 of type d2_data ! 20 Nov 2012 - Ph. Le Sager - fix runs longer than a month ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ CONTAINS !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_CH4_INIT ! ! !DESCRIPTION: Allocate memory. !\\ !\\ ! !INTERFACE: ! subroutine Emission_CH4_Init( status ) ! ! !USES: ! use dims, only : jm, lm use emission_read, only : providers_def, numb_providers, ed42_nsect_ch4 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 - adapted for AR5 ! 28 Jun 2012 - Ph. Le Sager - always allocate zch4 ! - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Emission_CH4_Init' integer :: region, i1, i2, j1, j2 integer :: jmr, lmr, lsec, lprov ! --- begin -------------------------------------- #ifdef with_ch4_emis ! --- count sectors if(has_ch4_emis) then ch4_2dsec = 0 ch4_3dsec = 0 do lprov = 1, numb_providers if (count(used_providers_ch4.eq.providers_def(lprov)%name)/=0) then if (trim(providers_def(lprov)%name) .eq. 'CMIP6') then ! Historical emissions of CH4 from aircraft are zero everywhere ! and it is assumed that this will also be the case for the future scenarios ch4_2dsec = ch4_2dsec + providers_def(lprov)%nsect2d elseif (trim(providers_def(lprov)%name) .eq. 'AR5') then ! nb of available sectors in AR5 depends on category ch4_2dsec = ch4_2dsec + n_ar5_ant_sec*count('CH4'.eq.ar5_cat_ant) + & n_ar5_shp_sec*count('CH4'.eq.ar5_cat_shp) if (LAR5BMB) ch4_2dsec = ch4_2dsec + n_ar5_bmb_sec*count('CH4'.eq.ar5_cat_bmb) ch4_3dsec = ch4_3dsec + n_ar5_air_sec*count('CH4'.eq.ar5_cat_air) elseif (trim(providers_def(lprov)%name) .eq. 'ED42') then ch4_2dsec = ch4_2dsec + ed42_nsect_ch4 ! no 3d sectors in EDGAR 4.2 else ch4_2dsec = ch4_2dsec + providers_def(lprov)%nsect2d ch4_3dsec = ch4_3dsec + providers_def(lprov)%nsect3d endif endif enddo allocate( ch4_emis_2d( nregions, ch4_2dsec ) ) allocate( ch4_emis_3d( nregions, ch4_3dsec ) ) allocate( has_data_2d(ch4_2dsec)) ; has_data_2d=.false. allocate( has_data_3d(ch4_3dsec)) ; has_data_3d=.false. end if #endif REG: do region=1,nregions lmr = lm(region) ; jmr = jm(region) CALL GET_DISTGRID( dgrid(region), i_strt=i1, j_strt=j1, i_stop=i2, j_stop=j2) #ifdef with_ch4_emis if(has_ch4_emis) then ! --- allocate information arrays (2d and 3d) do lsec=1,ch4_2dsec allocate( ch4_emis_2d(region,lsec)%surf(i1:i2, j1:j2) ) end do do lsec=1,ch4_3dsec allocate( ch4_emis_3d(region,lsec)%d3(i1:i2, j1:j2, lmr) ) end do end if #endif ! 1d-constraining surface array allocate( zch4(region)%d2(j1:j2) ) allocate( zch4_p(region)%d2(j1:j2) ) allocate( zch4_n(region)%d2(j1:j2) ) allocate( tmp2d(region)%d2(j2-j1+1)) zch4(region)%d2 = 0.0 zch4_p(region)%d2 = 0.0 zch4_n(region)%d2 = 0.0 ! work arrays if(isRoot) then allocate( wrk_c(region)%d23(1,jmr)) allocate( wrk_p(region)%d23(1,jmr)) allocate( wrk_n(region)%d23(1,jmr)) else allocate( wrk_c(region)%d23(1,1)) allocate( wrk_p(region)%d23(1,1)) allocate( wrk_n(region)%d23(1,1)) end if enddo REG status = 0 END SUBROUTINE EMISSION_CH4_INIT !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_CH4_DONE ! ! !DESCRIPTION: Free memory !\\ !\\ ! !INTERFACE: ! subroutine Emission_CH4_Done( status ) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - adapted to new structures ! 28 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Emission_CH4_Done' integer :: region, lsec ! --- begin -------------------------------------- reg: do region = 1, nregions deallocate( zch4(region)%d2 ) deallocate( zch4_n(region)%d2 ) deallocate( zch4_p(region)%d2 ) deallocate( tmp2d(region)%d2 ) deallocate( wrk_c(region)%d23 ) deallocate( wrk_p(region)%d23 ) deallocate( wrk_n(region)%d23 ) #ifdef with_ch4_emis if (has_ch4_emis) then do lsec=1,ch4_2dsec deallocate( ch4_emis_2d(region,lsec)%surf ) end do do lsec=1,ch4_3dsec deallocate( ch4_emis_3d(region,lsec)%d3 ) end do deallocate( has_data_2d, has_data_3d) end if #endif end do reg #ifdef with_ch4_emis if (has_ch4_emis) then deallocate( ch4_emis_2d ) deallocate( ch4_emis_3d ) end if #endif status = 0 end subroutine Emission_CH4_Done !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_CH4_DECLARE ! ! !DESCRIPTION: Opens, reads and evaluates input files (per month). ! Provides emissions on 2d/3d-arrays which are then added ! to mixing ratios in routine *apply. !\\ !\\ ! !INTERFACE: ! SUBROUTINE EMISSION_CH4_DECLARE( status ) ! ! !USES: ! use toolbox, only : coarsen_emission use dims, only : im, jm, lm, idate, sec_month, nlon360, nlat180, iglbsfc use chem_param, only : xmch4, ich4 use emission_data, only : msg_emis, LAR5BMB ! ---------------- CMIP6 - AR5 - EDGAR4 - ETC use emission_data, only : ch4_fixyear use emission_data, only : emis_input_year_ch4, emis_input_year_nat use emission_data, only : emis_input_dir_retro use emission_data, only : emis_input_dir_gfed use emission_data, only : emis_input_dir_ed4 #ifdef with_ch4_emis use emission_data, only : emis_input_dir_natch4 #endif use emission_read, only : read_cmip6_zch4 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_ed4_ReadSector use emission_read, only : emission_gfed_ReadSector use emission_read, only : emission_retro_ReadSector use emission_read, only : emission_lpj_ReadSector use emission_read, only : emission_hymn_ReadSector use emission_read, only : sectors_def, numb_sectors use emission_read, only : ar5_dim_3ddata use emission_read, only : ed42_ch4_sectors use Grid, only : FillGrid use meteodata, only : lli_z ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - adapted for AR5 ! 1 Dec 2011 - Narcisa Banda - added EDGAR 4.1 ! 28 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/emission_ch4_declare' integer :: region, hasData integer, parameter :: add_field=0 integer, parameter :: amonth=2 integer :: lsec, nyear, nmon integer :: target_year, target_month integer :: lmr, i1, i2, j1, j2 ! Emissions arrays real, dimension(:,:,:), allocatable :: field3d, field3d_p, field3d_n type(d3_data), dimension(nregions) :: emis3d, work, work3D type(emis_data) :: wrk2D(nregions) integer :: seccount2d, seccount3d ! --- begin ---------------------------------- #ifdef with_ch4_emis write(gol,'(" EMISS-INFO ------------- read CH4 emissions -------------")'); call goPr if(has_ch4_emis) then do region = 1, nregions do lsec=1,ch4_2dsec ch4_emis_2d(region,lsec)%surf = 0.0 end do do lsec=1,ch4_3dsec ch4_emis_3d(region,lsec)%d3 = 0.0 end do 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 end do ! global arrays for coarsening do region = 1, nregions if (isRoot)then allocate(work(region)%d3(im(region),jm(region),ar5_dim_3ddata)) 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_ch4 == sectors_def(lsec)%prov).eq.0) cycle if ((trim(sectors_def(lsec)%prov).eq.'ED42') .and. (count(ed42_ch4_sectors.eq.sectors_def(lsec)%name) .eq. 0)) cycle if (associated(sectors_def(lsec)%species)) then ! AR5 checks if (count('CH4'.eq.sectors_def(lsec)%species).eq.0) cycle if ((trim(sectors_def(lsec)%catname) .eq. 'biomassburning').and.(.not.LAR5BMB)) cycle endif ! Historical emissions of CH4 from aircraft are zero everywhere ! and it is assumed that this will also be the case for the future scenarios if ((trim(sectors_def(lsec)%prov).eq.'CMIP6') .and. (trim(sectors_def(lsec)%catname) .eq. 'aircraft')) cycle 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( 'CMIP6' ) call emission_cmip6_ReadSector( 'CH4', emis_input_year_ch4, idate(2), lsec, field3d, status ) IF_NOTOK_RETURN(status=1;deallocate(field3d)) case( 'CMIP6BMB' ) call emission_cmip6bmb_ReadSector( 'CH4', emis_input_year_ch4, idate(2), lsec, field3d, status ) IF_NOTOK_RETURN(status=1;deallocate(field3d)) case( 'AR5' ) ! Screen out solvent sector for CH4, ! because it is zero in the RCPs ! and not present in the historical files. if (trim(sectors_def(lsec)%name) .ne. 'emiss_slv') then call emission_ar5_ReadSector( 'CH4', emis_input_year_ch4, idate(2), lsec, field3d, status ) IF_NOTOK_RETURN(status=1) endif case( 'ED41' ) ! only transport sector (others provided by ED42) select case(trim(sectors_def(lsec)%name)) case ('1A3b_c_e','1A3d_SHIP','1A3d1') ! anthropogenic sources (only 2d) call emission_ed4_ReadSector( emis_input_dir_ed4, 'CH4', 'ch4', emis_input_year_ch4, idate(2), & lsec, trim(sectors_def(lsec)%prov) , 'kg / s', field3d, status ) IF_NOTOK_RETURN(status=1;deallocate(field3d)) end select case( 'ED42' ) ! biomass burning (GFED/RETRO/AR5BMB) and transport (ED41) are excluded through ED42_CH4_SECTORS definition call emission_ed4_ReadSector( emis_input_dir_ed4, 'CH4', 'ch4', emis_input_year_ch4, 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, 'ch4', emis_input_year_ch4, idate(2), & sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status ) IF_NOTOK_RETURN(status=1;deallocate(field3d)) case('RETRO') call emission_retro_ReadSector( emis_input_dir_retro, 'CH4', emis_input_year_ch4, idate(2), & sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status ) IF_NOTOK_RETURN(status=1;deallocate(field3d)) case( 'LPJ' ) ! this here is for natural sources (only 2d) call emission_lpj_ReadSector( trim(emis_input_dir_natch4)//'/LPJdata-jan2011', emis_input_year_nat, idate(2), & sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status ) IF_NOTOK_RETURN(status=1;deallocate(field3d)) case( 'HYMN' ) ! this here is for natural sources (only 2d) call emission_hymn_ReadSector( trim(emis_input_dir_natch4), & sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status ) IF_NOTOK_RETURN(status=1;deallocate(field3d)) case default write(gol,*) "Error in buidling list of providers USED_PROVIDERS_CH4"; 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 CH4 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 CH4 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 ! --------------------------------------- ! 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, 'CH4', xmch4, & sum(field3d) ) ! distribute to work arrays in regions call Coarsen_Emission( 'CH4 '//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 ) ch4_emis_3d(region,seccount3d)%d3 = ch4_emis_3d(region,seccount3d)%d3 + emis3d(region)%d3 end do else ! --------------------------- ! 2d data (Anthropogenic, Biomassburning, Natural) ! --------------------------- if (isRoot) then ! regrid ! write some numbers call msg_emis( amonth, trim(sectors_def(lsec)%prov), sectors_def(lsec)%name, 'CH4', xmch4, & sum(field3d(:,:,1)) ) call coarsen_emission( 'CH4 '//trim(sectors_def(lsec)%prov)//' '//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), ch4_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 ) deallocate( work3d(region)%d3 ) deallocate( emis3d(region)%d3 ) end do ! check sectors found if( seccount2d /= ch4_2dsec ) then write(gol,'(80("-"))') ; call goPr write(gol,'("ERROR: 2d sectors do not equal total number:",i4," /= ",i4," !")') seccount2d, ch4_2dsec ; call goErr write(gol,'(80("-"))') ; call goPr status=1; return end if if( seccount3d /= ch4_3dsec ) then write(gol,'(80("-"))') ; call goPr write(gol,'("ERROR: 3d sectors do not equal total number:",i4," /= ",i4," !")') seccount3d, ch4_3dsec ; call goErr write(gol,'(80("-"))') ; call goPr status=1; return end if end if #endif ! ! Read in either the zonal mean surface field for this year/month, or set to fixed value ! if (isRoot) then allocate( field3d(1,nlat180,1) ) allocate( field3d_p(1,nlat180,1) ) allocate( field3d_n(1,nlat180,1) ) field3d = 0.0 field3d_p = 0.0 field3d_n = 0.0 if (LCMIP6_CH4) then write(gol,'("Reading CMIP6 monthly zonal mean CH4 mixing ratios")'); call goPr ! previous month target_month = idate(2)-1 target_year = MIN(2500,MAX(emis_input_year_ch4,1850)) if (idate(2) .eq. 1) then target_month = 12 if (.not.ch4_fixyear) then target_year = MIN(2500,MAX(idate(1)-1,1850)) endif endif write (gol,*) ' Use year ', target_year,' and month ', target_month,' for previous month'; call goPr call read_cmip6_zch4(field3d_p(1,:,1),target_year,target_month, status) ! current month target_month = idate(2) target_year = MIN(2500,MAX(emis_input_year_ch4,1850)) write (gol,*) ' Use year ', target_year,' and month ', target_month,' for current month'; call goPr call read_cmip6_zch4(field3d(1,:,1),target_year,target_month, status) ! next month target_month = idate(2)+1 target_year = MIN(2500,MAX(emis_input_year_ch4,1850)) if (idate(2) .eq. 12) then target_month = 1 if (.not.ch4_fixyear) then target_year = MIN(2500,MAX(emis_input_year_ch4+1,1850)) endif endif write (gol,*) ' Use year ', target_year,' and month ', target_month,' for next month'; call goPr call read_cmip6_zch4(field3d_n(1,:,1),target_year,target_month, status) !write (gol,'("STOP AFTER TESTING READING OF CMIP6 CH4")') ; call goErr !status=1; TRACEBACK; return else if (emis_ch4_single) then ! fixed value ! TvN: bug fix ! unit should be in ppbv as in CMIP6 and NOAA/GMD files !field3d(:,:,:)=emis_ch4_fixed_ppb*1.0e-9 field3d(:,:,:)=emis_ch4_fixed_ppb else ! read NOAA/GMD monthly surface latitudinal distribution call read_bkglat_ch4(field3d(1,:,1),idate(1),idate(2), status) call prev_mon(idate(1),idate(2),nyear,nmon) call read_bkglat_ch4(field3d_p(1,:,1),nyear,nmon, status) call next_mon(idate(1),idate(2),nyear,nmon) call read_bkglat_ch4(field3d_n(1,:,1),nyear,nmon, status) endif ! coarsen or distribute to zoom regions: do region = 1, nregions call FillGrid( lli_z(region), 'n', wrk_c(region)%d23(:,:), & lli_z(iglbsfc), 'n', field3d(:,:,1), 'area-aver', status ) IF_NOTOK_RETURN(status=1) call FillGrid( lli_z(region), 'n', wrk_p(region)%d23(:,:), & lli_z(iglbsfc), 'n', field3d_p(:,:,1), 'area-aver', status ) IF_NOTOK_RETURN(status=1) call FillGrid( lli_z(region), 'n', wrk_n(region)%d23(:,:), & lli_z(iglbsfc), 'n', field3d_n(:,:,1), 'area-aver', status ) IF_NOTOK_RETURN(status=1) end do deallocate( field3d, field3d_p, field3d_n ) end if ! scatter along latitude direction, then broadcast to all longitudes DO region = 1, nregions call SCATTER_I_BAND( dgrid(region), zch4(region)%d2, wrk_c(region)%d23(1,:), status ) IF_NOTOK_RETURN(status=1) call SCATTER_I_BAND( dgrid(region), zch4_p(region)%d2, wrk_p(region)%d23(1,:), status ) IF_NOTOK_RETURN(status=1) call SCATTER_I_BAND( dgrid(region), zch4_n(region)%d2, wrk_n(region)%d23(1,:), status ) IF_NOTOK_RETURN(status=1) tmp2d(region)%d2 = zch4(region)%d2 call PAR_BROADCAST( tmp2d(region)%d2, status, row=.true. ) IF_NOTOK_RETURN(status=1) zch4(region)%d2 = tmp2d(region)%d2 tmp2d(region)%d2 = zch4_p(region)%d2 call PAR_BROADCAST( tmp2d(region)%d2, status, row=.true. ) IF_NOTOK_RETURN(status=1) zch4_p(region)%d2 = tmp2d(region)%d2 tmp2d(region)%d2 = zch4_n(region)%d2 call PAR_BROADCAST( tmp2d(region)%d2, status, row=.true. ) IF_NOTOK_RETURN(status=1) zch4_n(region)%d2 = tmp2d(region)%d2 END DO ! ok status = 0 END SUBROUTINE EMISSION_CH4_DECLARE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_CH4_APPLY ! ! !DESCRIPTION: Takes monthly emissions, and ! - split them vertically ! - apply time splitting factors ! - add them up (add_3d) !\\ !\\ ! !INTERFACE: ! SUBROUTINE EMISSION_CH4_APPLY( region, status) ! ! !USES: ! use binas, only : xmair use dims, only : itaur, nsrce, tref use dims, only : im, jm, lm, at, bt 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 emission_read, only : sectors_def, numb_sectors use emission_read, only : ed42_ch4_sectors use chem_param, only : ich4, xmch4 ! NB use chem_param, only : ch4_ps use meteodata, only : m_dat use global_data, only : region_dat, mass_dat use toolbox, only : lvlpress use partools, only : par_reduce_element #ifdef with_budgets use budget_global, only : budg_dat, nzon_vg use emission_data, only : sum_emission, budemi_dat #endif ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - adapted to AR5 emissions structures ! - added nudging if emissions ! 28 Jun 2012 - Ph. Le Sager - nudging in case of with_ch4_emis no more ! limited to 3x2 runs ! - optimized loop order for case of NO emissions ! - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! !EOP !--------------------------------------------------------------------------- !BOC character(len=*), parameter :: rname = mname//'/Emission_CH4_apply' integer, dimension(6) :: idater real :: dtime, fraction integer :: imr, jmr, lmr, lsec integer :: seccount2d, seccount3d type(d3_data) :: emis3d type(emis_data) :: emis2d integer :: i1, j1, i2, j2 real, dimension(:,:,:), allocatable :: field3d real, dimension(:,:,:,:), pointer :: rm real, dimension(:,:,:), pointer :: m real :: ch4ppb, toadd integer :: nzone_v real :: zvmr_obs, f_ratio, vmr_min integer :: i, j, l, lm_ch4, idateline real, allocatable :: vmr_mod(:,:), zvmrsum_local(:) real, allocatable :: zvmr_mod(:) real, parameter :: gnud_ch4 = 4.0e-6 ! ~3day time scale, see below ! --- begin --------------------------- call get_distgrid( dgrid(region), i_strt=i1, j_strt=j1, i_stop=i2, j_stop=j2) nullify(rm) nullify(m) rm => mass_dat(region)%rm m => m_dat(region)%data allocate(zvmr_mod(j1:j2)) !--------------------------------------------------------------- ! CH4 emissions !-------------------------------------------------------------- #ifdef with_ch4_emis if (has_ch4_emis) then call tau2date(itaur(region),idater) dtime=float(nsrce)/(2*tref(region)) !emissions are added in two steps...XYZECCEZYX. if(okdebug) then write(gol,*)'emission_ch4_apply in region ',region,' at date: ',idater, ' with time step:', dtime call goPr end if ! get a working structure lmr = lm(region) allocate( emis3d%d3 (i1:i2, j1:j2, lmr) ) ; emis3d%d3 = 0.0 allocate( emis2d%surf(i1:i2, j1:j2 ) ) ; emis2d%surf = 0.0 ! count 2d and 3d sectors seccount2d = 0 seccount3d = 0 ! cycle over sectors do lsec = 1, numb_sectors if (count(used_providers_ch4.eq.sectors_def(lsec)%prov).eq.0) cycle if ((trim(sectors_def(lsec)%prov).eq.'ED42') .and. (count(ed42_ch4_sectors.eq.sectors_def(lsec)%name) .eq. 0)) cycle if (associated(sectors_def(lsec)%species)) then ! AR5 checks if (count('CH4'.eq.sectors_def(lsec)%species).eq.0) cycle if ((trim(sectors_def(lsec)%catname) .eq. 'biomassburning').and.(.not.LAR5BMB)) cycle endif ! Historical emissions of CH4 from aircraft are zero everywhere ! and it is assumed that this will also be the case for the future scenarios if ((trim(sectors_def(lsec)%prov).eq.'CMIP6') .and. (trim(sectors_def(lsec)%catname) .eq. 'aircraft')) cycle ! 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 = ch4_emis_3d(region,seccount3d)%d3 else seccount2d = seccount2d + 1 if (.not.has_data_2d(seccount2d)) cycle emis2d%surf = ch4_emis_2d(region,seccount2d)%surf ! account for soil sink scale with respect to surface [CH4] over i,j if (trim(sectors_def(lsec)%name).eq.'soilconsumption') then do j=j1,j2 do i=i1,i2 if(emis2d%surf(i,j)>0.) & emis2d%surf(i,j) = ( - emis2d%surf(i,j)) * (rm(i,j,1,ich4)/m(i,j,1)/1.0e-6) * (xmair/xmch4) enddo enddo end if emis3d%d3 = 0.0 ! vertically distribute according to sector call emission_vdist_by_sector( sectors_def(lsec)%vdisttype, 'CH4', region, emis2d, 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, ich4, i1, j1, emis3d%d3, bb_cycle(region)%scalef, xmch4, xmch4, status,fraction ) IF_NOTOK_RETURN(status=1) else call do_add_3d( region, ich4, i1, j1, emis3d%d3, xmch4, xmch4, status, fraction ) IF_NOTOK_RETURN(status=1) end if end do deallocate( emis3d%d3 ) deallocate( emis2d%surf ) end if #endif /* with_ch4_emis */ !----------------------------------------------------- ! Fix methane mixing ratio !----------------------------------------------------- ! Set highest level for nudging !------------------------------ #ifdef with_ch4_emis !apply nuddging up to 550 hpa lm_ch4 = lvlpress(region,55000.,98400.) #else ! default: nudging in the lowest layer lm_ch4 = 1 if ( emis_ch4_single .and. emis_ch4_fix3d ) then lm_ch4 = lm(region) else ! TvN: in case of a shallow first layer (incl. L60, L91), nudge the lowest 2 layers if (at(2)+bt(2)*101325. > 101000.) lm_ch4 = 2 endif #endif ! In simulations without CH4 emissions, ! the CH4 mixing ratios are prescribed (not nudged). !--------------------------------------------------- #ifndef with_ch4_emis do l = 1, lm_ch4 do j = j1, j2 do i = i1, i2 ! budget for fixing: ch4ppb = 1e9 * (xmair/xmch4) * rm(i,j,l,ich4) / m(i,j,l) ! unit of zch4 is in ppb: toadd = (zch4(region)%d2(j) - ch4ppb) * 1e-9 * m(i,j,l) * 1e3 / xmair ! moles ch4 rm(i,j,l,ich4) = m(i,j,l)/xmair * zch4(region)%d2(j) *1e-9 * xmch4 #ifdef with_budgets nzone_v=nzon_vg(l) budemi_dat(region)%emi(i,j,nzone_v,ich4) = budemi_dat(region)%emi(i,j,nzone_v,ich4) + toadd if(ich4 ==1) sum_emission(region) = sum_emission(region) + toadd*xmch4/1e3 !in kg #endif end do end do end do #endif ! If CH4 emissions are used ! the CH4 mixing ratios are nudged, ! to either zonal means from CMIP6 or ! zonal background values from NOAA/GMD. !--------------------------------------- #ifdef with_ch4_emis if (LCMIP6_CH4) then ! Nudging for CMIP6 based on zonal means ! adapted from code for nudging to HALOE allocate(vmr_mod(i1:i2,j1:j2)) allocate(zvmrsum_local(j1:j2)) ! convert CH4 mass to volume mixing ratios (mol/mol) vmr_mod = (rm(i1:i2,j1:j2,1,ich4)/m(i1:i2,j1:j2,1))*(xmair/xmch4) zvmrsum_local=sum(vmr_mod,dim=1) ! zonal total for this subdomain CALL PAR_REDUCE_ELEMENT(zvmrsum_local, 'SUM', zvmr_mod, status, all=.true., row=.true.) IF_NOTOK_RETURN(status=1) ! zonal average zvmr_mod(:)=zvmr_mod(:)/im(1) deallocate(vmr_mod) deallocate(zvmrsum_local) else ! Nudging to zonal background values from NOAA/GMD ! nudging full 3D field using surface observations in background atmosphere ! the observations represent the minimum of rm in the zonal band ! the unit of zch4 in the netcdf files is in ppb ! zonal_mod is the zonal minimum volume mixing ratio (mol/mol) in the model ! use arbitary large initial value zvmr_mod(:)=3.0e-6 ! ! Use value at the dateline away from major CH4 sources ! idateline = 1 ! if current processor include dateline if (i1==1) then do j=j1,j2 vmr_min = (rm(idateline,j,1,ich4)/m(idateline,j,1))*(xmair/xmch4) zvmr_mod(j) = min(zvmr_mod(j),vmr_min) ! mimimum CH4 at latitude j end do endif ! broadcast (from proc with i1==1, ie from root of row_communicator, ie 0, ie default broadcaster) call par_broadcast(zvmr_mod, status, row=.true.) IF_NOTOK_RETURN(status=1) endif ! LCMIP6_CH4 ! nudging like in the stratosphere for ozone would be: ! rm = (rmold+dtime*gnud*zvmr_obs)/(1.+dtime*gnud) ! nudging used here is done without the approximation: ! 1 - exp(-dtime*gnud) = dtime*gnud ! we determine per zonal band the ratio between ! the CMIP6 mean/observations ! and the zonal mean/minimum surface model field: ! f_ratio = zvmr_obs / zvmr_mod ! we apply the surface ratio to the total 3d field ! with a slow nudging time scale to distribute the ! corrections at this latitude band over all heights ! we use gnud_ch4 = 4.e-7 (1/sec; for ~1 month nudging time scale ) ! The nudging equation is: ! vmr = vmrold * exp(-dtime*gnud) + zvmr_obs * (1 - exp(-dtime*gnud)) or: ! vmr = vmrold - (1 - exp(-dtime*gnud))(vmrold - vmrobs) or: ! vmr = vmrold + vmrold*(f_ratio-1)*(1 - exp(-dtime*gnud)) ! toadd = vmrold*(f_ratio-1)*(1 - exp(-dtime*gnud)) dtime=nsrce/(2*tref(region)) do j = j1, j2 ! interpolate in time and convert to volume mixing ratio (mol/mol) call interp_zch4(region,j,idater,itaur(region),zvmr_obs) ! trap bad cases (0. as initial conditions,...) if (zvmr_mod(j)==0.) then f_ratio=1. else f_ratio = zvmr_obs/zvmr_mod(j) ! f_ratio should be always very close to 1 ! (forcing f_ratio to 1 removes the nudging) end if ! expected behavior is nudging within 10% !PLS: commented since it appears too often and fills log file !if (f_ratio>1.1 .or. f_ratio<0.9) then !write(gol,*)"WARNING: CH4 nudging larger than 10% for region/ijl",region,i,j,l; call goPr !end if do l = 1, lm_ch4 do i = i1, i2 toadd = rm(i,j,l,ich4) * (f_ratio - 1.) * (1. - exp(- dtime*gnud_ch4)) *1e3 / xmair ! moles ch4 rm(i,j,l,ich4) = rm(i,j,l,ich4) + toadd * xmch4/1e3 #ifdef with_budgets nzone_v=nzon_vg(l) ! add the nudging adjustments to the emission budget (now full 3D field filled) ! the emissions are not added here budemi_dat(region)%emi(i,j,nzone_v,ich4) = budemi_dat(region)%emi(i,j,nzone_v,ich4) + toadd #endif end do end do end do #endif !----------------------------------------------------- ! Done !----------------------------------------------------- nullify(rm) nullify(m) deallocate(zvmr_mod) status = 0 END SUBROUTINE EMISSION_CH4_APPLY !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: READ_BKGLAT_CH4 ! ! !DESCRIPTION: Read ERA Interim background surface zonal CH4 for specified ! year/month !\\ !\\ ! !INTERFACE: ! SUBROUTINE READ_BKGLAT_CH4(field, year, month, status) ! ! !USES: ! USE MDF, ONLY : MDF_OPEN, MDF_NETCDF, MDF_READ, MDF_Get_Var, MDF_Close, MDF_Inq_VarID use emission_data, only : emis_zch4_fname use dims, only : nlat180 ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status real, dimension(nlat180), intent(out) :: field ! ! !INPUT PARAMETERS: ! integer, intent(in) :: year, month ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Read_BkgLat_CH4' character(len=256) :: filemon, emis_zch4_fullname integer :: fid, varid, year_valid integer, dimension(2), parameter :: era_interim_avail = (/1989, 2014/) ! valid year year_valid = MIN(era_interim_avail(2),MAX(year, era_interim_avail(1))) ! target file name with year e.g. surface_ch4_zm_1989.nc write (emis_zch4_fullname,'(a,"/surface_ch4_zm_",i4.4,".nc")') trim(emis_zch4_fname), year_valid CALL MDF_Open( TRIM(emis_zch4_fullname), MDF_NETCDF, MDF_READ, fid, status ) IF_NOTOK_RETURN(status=1) ! search for the record for requested month: if (month < 10) then write(filemon,'(a,i1)')'CH4_zm', month else write(filemon,'(a,i2)')'CH4_zm', month endif CALL MDF_Inq_VarID( fid, TRIM(filemon), varid, status ) IF_ERROR_RETURN(status=1) if ( varid < 0 ) then write (gol,'("WARNING - no surface CH4 concentrations `",a,"`")') trim(emis_zch4_fullname) ; call goPr status=1; TRACEBACK; return else CALL MDF_Get_Var( fid, varid, field, status ) IF_NOTOK_RETURN(status=1) endif CALL MDF_Close( fid, status ) IF_NOTOK_RETURN(status=1) status = 0 END SUBROUTINE READ_BKGLAT_CH4 !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: INTERP_ZCH4 ! ! !DESCRIPTION: ! Interpolates in time the latitudinal monthly background CH4 ! concentration. The data are considered representative of the ! 16th of each month, except February when it is 15th. !\\ !\\ ! !INTERFACE: ! SUBROUTINE INTERP_ZCH4(region,j,idater,itau, zch4_obs) ! ! !USES: ! use datetime, only : date2tau ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region,j integer, dimension(6), intent(in) :: idater integer(kind=8), intent(in) :: itau ! ! !OUTPUT PARAMETERS: ! real, intent(out) :: zch4_obs ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC integer, dimension(6) :: idate_c, idate_p, idate_n integer :: nmon, nyear integer(kind=8) :: itau_c, itau_p, itau_n idate_c(1) = idater(1) idate_c(2) = idater(2) if (idater(2) == 2) then idate_c(3) = 15 else idate_c(3) = 16 endif idate_c(4:6) = 0 call date2tau(idate_c,itau_c) call prev_mon(idater(1),idater(2),nyear,nmon) idate_p(1) = nyear idate_p(2) = nmon if (idate_p(2) == 2) then idate_p(3) = 15 else idate_p(3) = 16 endif idate_p(4:6) = 0 call date2tau(idate_p,itau_p) call next_mon(idater(1),idater(2),nyear,nmon) idate_n(1) = nyear idate_n(2) = nmon if (idate_n(2) == 2) then idate_n(3) = 15 else idate_n(3) = 16 endif idate_n(4:6) = 0 call date2tau(idate_n,itau_n) if (itau.lt.itau_c) then zch4_obs = zch4_p(region)%d2(j) * float(itau_c-itau)/float(itau_c-itau_p) + zch4(region)%d2(j) * float(itau-itau_p)/float(itau_c-itau_p) else zch4_obs = zch4(region)%d2(j) * float(itau_n-itau)/float(itau_n-itau_c) + zch4_n(region)%d2(j) * float(itau-itau_c)/float(itau_n-itau_c) endif zch4_obs = zch4_obs *1.e-9 END SUBROUTINE INTERP_ZCH4 !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: NEXT_MON ! ! !DESCRIPTION: Return month number and year of next month !\\ !\\ ! !INTERFACE: ! SUBROUTINE NEXT_MON(year,mon,nyear,nmon) ! ! !INPUT PARAMETERS: ! integer, intent(in) :: year, mon ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: nyear, nmon ! !EOP !------------------------------------------------------------------------ !BOC if (mon.lt.12) then nyear = year nmon = mon+1 else nyear = year+1 nmon = 1 endif END SUBROUTINE NEXT_MON !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: PREV_MON ! ! !DESCRIPTION: Return month number and year of previous month !\\ !\\ ! !INTERFACE: ! SUBROUTINE PREV_MON(year,mon,nyear,nmon) ! ! !INPUT PARAMETERS: ! integer, intent(in) :: year, mon ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: nyear, nmon ! !EOP !------------------------------------------------------------------------ !BOC if (mon.gt.1) then nyear = year nmon = mon-1 else nyear = year-1 nmon = 12 endif END SUBROUTINE PREV_MON !EOC END MODULE EMISSION_CH4