! #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" ! !----------------------------------------------------------------------------- ! TM5 ! !----------------------------------------------------------------------------- !BOP ! ! !MODULE: EMISSION_NOX ! ! !DESCRIPTION: hold data and methods for NOx emissions. ! ----------------- ! AR5 emissions ! ----------------- ! For each month, arrays emis_nox_2d/3d have to be filled. ! It follows the following settins: ! - take emiss_ene/dom/ind/wst/agr/awb/slv/tra/shp/ ! /air/forestfire/grassfire from AR5 data sets (NO GFED!!) ! - use natural emissions from MACC data sets ! (emiss_nat/soil/bio/oc) ! - vertical distribution is done via emission_vdist_by_sector ! (emission_data.F90) ! - lightning is done online (eminox_lightning) !\\ !\\ ! !INTERFACE: ! MODULE EMISSION_NOX ! ! !USES: ! use GO, only : gol, goErr, goPr, goBug use tm5_distgrid, only : dgrid, get_distgrid, scatter, gather use dims, only : nregions, idate, dy, 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_NOx_Init public :: Emission_NOx_Done public :: Emission_NOx_Declare public :: Emission_NOx_bb_daily_cycle #ifndef without_convection public :: lightningNOX #endif public :: nox_emis_3d, nox_emis_3d_bb_app public :: eminox_lightning ! ! !DATA MEMBERS: ! character(len=*), parameter :: mname = 'emission_nox' type(d3_data), dimension(nregions), target :: nox_emis_3d, nox_emis_3d_bb, nox_emis_3d_bb_app type(d3_data), dimension(nregions), target :: eminox_lightning integer :: nox_2dsec, nox_3dsec real :: fscalelig ! scaling used in lightning NOX production to get 5.98 Tg for 2006 ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - overhaul for AR5 ! 1 Dec 2011 - Narcisa Banda - added EDGAR 4 ! 27 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! NOx emissions are added directly in chemistry, instead of apart from it. ! !EOP !------------------------------------------------------------------------ CONTAINS !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_NOX_INIT ! ! !DESCRIPTION: allocate memory !\\ !\\ ! !INTERFACE: ! subroutine Emission_NOx_Init( status ) ! ! !USES: ! use dims, only : lm use emission_read, only : providers_def, numb_providers, ed42_nsect_nox 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 #ifndef without_convection use meteodata, only : set, gph_dat, temper_dat, cp_dat use emission_data, only : use_tiedkte #endif ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - adapted for AR5 ! 27 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition ! 10 Jul 2013 - Ph. Le Sager - init lightning when no inventory is selected ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Emission_NOx_Init' integer :: region, i1, j1, i2, j2 integer :: imr, jmr, lmr, lsec, lprov ! --- begin -------------------------------------- status = 0 #ifndef without_convection ! Meteo used for LightningNOx do region=1,nregions call Set( temper_dat(region), status, used=.true. ) call Set( gph_dat(region), status, used=.true. ) call Set( cp_dat(region), status, used=.true. ) enddo ! ! Scaling parameter for LiNOx ! ! Set to get 5.98 Tg N with 2006 EI met fields. This is resolution ! and met fields dependent. The factor has been estimated for: ! - @3x2 and 34 levels, Tiedkte : fscalelig=13.715 ! - @1x1 and 34 levels, Tiedkte : fscalelig=17.051 ! - @3x2 and 34 levels, EI conv : fscalelig=13.715*0.786 ! - @1x1 and 34 levels, EI conv : fscalelig=17.051*0.649 ! if (use_tiedkte) then ! convective fluxes computed from T/rh/wind (Tiedkte) fscalelig=13.715 ! 3x2-34L, Tiedkte scheme if (dy == 1) fscalelig=17.051 ! 1x1-34L, Tiedkte scheme else fscalelig=10.78 ! 3x2-34L, EI convective fluxes if (dy == 1) fscalelig=11.066 ! 1x1-34L, EI convective fluxes endif #endif ! 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) allocate(eminox_lightning(region)%d3(i1:i2,j1:j2,lmr) ) eminox_lightning(region)%d3=0. enddo ! Check if any inventory is used if(.not. has_emis) return ! nb of sectors nox_2dsec = 0 nox_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 nox_2dsec = nox_2dsec + n_ar5_ant_sec*count('NO'.eq.ar5_cat_ant) + & n_ar5_shp_sec*count('NO'.eq.ar5_cat_shp) if (LAR5BMB) nox_2dsec = nox_2dsec + n_ar5_bmb_sec*count('NO'.eq.ar5_cat_bmb) nox_3dsec = nox_3dsec + n_ar5_air_sec*count('NO'.eq.ar5_cat_air) ! nox_2dsec = nox_2dsec + providers_def(lprov)%nsect2d ! nox_3dsec = nox_3dsec + count('NO'.eq.ar5_cat_air) elseif (trim(providers_def(lprov)%name) .eq. 'ED42') then nox_2dsec = nox_2dsec + ed42_nsect_nox ! no 3d sectors in EDGAR 4.2 else nox_2dsec = nox_2dsec + providers_def(lprov)%nsect2d nox_3dsec = nox_3dsec + providers_def(lprov)%nsect3d endif endif enddo ! 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) allocate( nox_emis_3d (region)%d3(i1:i2,j1:j2,lmr) ) allocate( nox_emis_3d_bb (region)%d3(i1:i2,j1:j2,lmr) ) allocate( nox_emis_3d_bb_app(region)%d3(i1:i2,j1:j2,lmr) ) enddo status = 0 END SUBROUTINE EMISSION_NOX_INIT !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_NOX_DONE ! ! !DESCRIPTION: Free memory !\\ !\\ ! !INTERFACE: ! SUBROUTINE EMISSION_NOX_DONE( status ) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - adapted for AR5 ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Emission_NOx_Done' integer :: region, lsec ! --- begin --------------------------------- status = 0 if(.not. has_emis) return do region = 1, nregions deallocate( nox_emis_3d (region)%d3 ) deallocate( nox_emis_3d_bb (region)%d3 ) deallocate( nox_emis_3d_bb_app(region)%d3 ) deallocate( eminox_lightning (region)%d3 ) end do status = 0 END SUBROUTINE EMISSION_NOX_DONE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_NOX_DECLARE ! ! !DESCRIPTION: Opens, reads and evaluates input files (per month). ! Provides emissions on 2d/3d-arrays which are then added ! in the chemistry routine (no *apply !). ! Vertically distribute the 2D dataset according to sector. !\\ !\\ ! !INTERFACE: ! SUBROUTINE EMISSION_NOX_DECLARE( status ) ! ! !USES: ! use toolbox, only : coarsen_emission use partools, only : isRoot, par_broadcast use dims, only : im, jm, lm, nlon360, nlat180, iglbsfc use dims, only : newsrun, idate, sec_month use chem_param, only : xmn, xmno2, xmno use emission_data, only : msg_emis, emission_vdist_by_sector, LAR5BMB ! ---------------- AR5 - EDGAR 4 - ETC. -------------------- use emission_data, only : emis_input_year use emission_data, only : emis_input_dir_mac, emis_input_dir_ed4 use emission_data, only : emis_input_dir_retro, emis_input_dir_gfed use emission_read, only : emission_ar5_regrid_aircraft use emission_read, only : emission_ar5_ReadSector, emission_macc_ReadSector use emission_read, only : emission_ed4_ReadSector, 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_nox_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 ! 27 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition ! 25 Feb 2014 - Jason Williams - separate array for BMB so that burning daily cycle can be applied ! ! !REMARKS: ! (1) Because we do not use an apply method, the vertical distribution ! is done here. However this is a bug, since this is time dependent. ! Possible solution: do vert dist in chemistry like the BMB cycle, ! or in the more general BMB cycle routine. ! ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/emission_nox_declare' integer :: region, hasData integer,parameter :: add_field=0 integer,parameter :: amonth=2 integer :: imr, jmr, lmr, lsec, i1, i2, j1, j2 ! AR5 temporary arrays real, dimension(:,:,:), allocatable :: field3d !, field3d2 type(d3_data), dimension(nregions), target :: emis3d, work(nregions), work3d(nregions) type(emis_data), dimension(nregions), target :: emis2d, wrk2D(nregions) ! defensive integer :: seccount2d, seccount3d ! --- begin ----------------------------------------- status = 0 if(.not. has_emis) return write(gol,'(" EMISS-INFO ------------- read NOx emissions -------------")'); call goPr ! reset emissions, allocate work array do region = 1, nregions nox_emis_3d(region)%d3 = 0.0 ; nox_emis_3d_bb(region)%d3 = 0.0 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( emis2d(region)%surf(i1:i2,j1:j2) ) ; emis2d(region)%surf = 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.eq.sectors_def(lsec)%prov).eq.0) cycle if ((trim(sectors_def(lsec)%prov).eq.'ED42') .and. (count(ed42_nox_sectors.eq.sectors_def(lsec)%name) .eq. 0)) cycle if (associated(sectors_def(lsec)%species)) then ! AR5 check if (count('NO'.eq.sectors_def(lsec)%species).eq.0) cycle if ((trim(sectors_def(lsec)%catname) .eq. 'biomassburning').and.(.not.LAR5BMB)) cycle endif if( sectors_def(lsec)%f3d ) then seccount3d = seccount3d + 1 else seccount2d = seccount2d + 1 end if field3d = 0.0 #ifdef with_online_nox ! skip natural nox in case it is calculated online if( trim(sectors_def(lsec)%namecat == 'natural') ) then write (gol,'(80("-"))'); call goPr write (gol,'("INFO: skipping sector `",a,"` due to `with_online_nox` ")') trim(sectors_def(lsec)%name); call goPr cycle end if #endif root: if (isRoot) then ! READ select case( trim(sectors_def(lsec)%prov) ) case( 'AR5' ) ! AR5 emissions included NO and NO2 aircraft emissions, but they are duplicates. So ! only take into account one of the sets. (in TM5: NO, skip NO2). ! Screen out solvent sector for NO, ! 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( 'NO', emis_input_year, idate(2), lsec, field3d, status ) IF_NOTOK_RETURN(status=1) ! convert from (kg NO)/s to (kg N)/s field3d = field3d * xmn / xmno endif case( 'MACC' ) ! screen out sectors w/o NOx (bio, oc, nat) if ( (trim(sectors_def(lsec)%name) .eq. 'emiss_soil' ) .or. & (trim(sectors_def(lsec)%name) .eq. 'emiss_anthro') .or. & (trim(sectors_def(lsec)%name) .eq. 'emiss_air' ) ) then call emission_macc_ReadSector( emis_input_dir_mac, 'NO', emis_input_year, idate(2), & '0.5x0.5_kg.nc', sectors_def(lsec)%name, 'kg NO / s', field3d, status ) IF_NOTOK_RETURN(status=1) ! convert from (kg NO)/s to (kg N)/s field3d = field3d * xmn / xmno endif case( 'ED41' ) select case(trim(sectors_def(lsec)%name)) case ('1A3a','1A3b_c_e','1A3d_SHIP','1A3d1') ! anthropogenic sources call emission_ed4_ReadSector( emis_input_dir_ed4, 'NOx','nox', emis_input_year, 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_NOX_SECTORS definition call emission_ed4_ReadSector( emis_input_dir_ed4, 'NOx', 'nox', 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, 'nox', emis_input_year, idate(2), & sectors_def(lsec)%name, 'kg NO2 / s', field3d(:,:,1), status ) IF_NOTOK_RETURN(status=1) ! convert from (kg NO2)/s to (kg N)/s field3d = field3d * xmn / xmno2 case('RETRO') call emission_retro_ReadSector( emis_input_dir_retro, 'NOX', emis_input_year, idate(2), & sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status ) IF_NOTOK_RETURN(status=1) ! in the file kg(species)/m2/s - what does this mean?? by the numbers I assume kg NO2 ! convert from (kg NO2)/s to (kg N)/s field3d = field3d * xmn / xmno2 case('MEGAN') ! ! use soil emissions from MACC ! case('DUMMY') case default write(gol,*) "Error in building list of providers USED_PROVIDERS"; call goErr status=1; TRACEBACK; return end select ! verbose if(sum(field3d) < 100.*TINY(1.0) ) then if (okdebug) then write(gol,'("EMISS-INFO - no NOx emissions found for ",a," ",a," for month ",i2 )') & trim(sectors_def(lsec)%prov), sectors_def(lsec)%name, idate(2) ; call goPr endif hasData=0 else if (okdebug) then write(gol,'("EMISS-INFO - found NOx emissions for ",a," ",a," for month ",i2 )') & trim(sectors_def(lsec)%prov), sectors_def(lsec)%name, idate(2) ; call goPr endif ! scale from kg/s to kg/month field3d = field3d * sec_month ! kg / month hasData=1 endif end if root call Par_broadcast(hasData, status) IF_NOTOK_RETURN(status=1) if (hasData == 0) cycle sec ! reset temporary arrays do region = 1, nregions emis3d(region)%d3 = 0.0 work3d(region)%d3 = 0.0 emis2d(region)%surf = 0.0 end do ! distinguish 2d/3d sectors if( sectors_def(lsec)%f3d ) then ! --------------------------- ! 3d data (AIRCRAFT) ! --------------------------- if (isRoot) then ! write some numbers call msg_emis( amonth, trim(sectors_def(lsec)%prov), sectors_def(lsec)%name, 'NOx', xmn, sum(field3d) ) ! distribute to work arrays in regions call Coarsen_Emission( 'NOX '//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 ) nox_emis_3d(region)%d3 = nox_emis_3d(region)%d3 + emis3d(region)%d3 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, 'NOx', xmn, sum(field3d(:,:,1)) ) call coarsen_emission( 'NOx '//sectors_def(lsec)%name, nlon360, nlat180, field3d(:,:,1), & wrk2D, add_field, status ) IF_NOTOK_RETURN(status=1) end if ! scatter, distribute vertically according to sector, and sum up on target array do region = 1, nregions call scatter(dgrid(region), emis2d(region)%surf, work(region)%d3(:,:,1), 0, status) IF_NOTOK_RETURN(status=1) call emission_vdist_by_sector( sectors_def(lsec)%vdisttype, 'NOx', region, emis2d(region), emis3d(region), status ) IF_NOTOK_RETURN(status=1) if ( trim(sectors_def(lsec)%catname) .eq. 'biomassburning') then nox_emis_3d_bb(region)%d3 = nox_emis_3d_bb(region)%d3 + emis3d(region)%d3 else nox_emis_3d(region)%d3 = nox_emis_3d(region)%d3 + emis3d(region)%d3 endif end do end if ! sectors_def end do sec ! sectors deallocate( field3d ) do region = 1, nregions if (associated(wrk2D(region)%surf)) nullify(wrk2D(region)%surf) deallocate( emis3d(region)%d3, emis2d(region)%surf ) deallocate( work(region)%d3 ) deallocate( work3d(region)%d3 ) end do ! check sectors found if( seccount2d /= nox_2dsec ) then write(gol,'(80("-"))') ; call goPr write(gol,'("ERROR: 2d sectors do not equal total number:",i4," /= ",i4," !")') seccount2d, nox_2dsec ; call goErr write(gol,'(80("-"))') ; call goPr status=1; return end if if( seccount3d /= nox_3dsec ) then write(gol,'(80("-"))') ; call goPr write(gol,'("ERROR: 3d sectors do not equal total number:",i4," /= ",i4," !")') seccount3d, nox_3dsec ; call goErr write(gol,'(80("-"))') ; call goPr status=1; return end if ! ok status = 0 END SUBROUTINE EMISSION_NOX_DECLARE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_NOX_BB_DAILY_CYCLE ! ! !DESCRIPTION: Impose daily burning cycle to BMB NOx emissions for current ! time step. !\\ !\\ ! !INTERFACE: ! SUBROUTINE EMISSION_NOX_BB_DAILY_CYCLE( status ) ! ! !USES: ! use dims, only : itaur, nsrce, tref, lm use dims, only : dx, xref, xbeg, yref, ybeg, ndyn_max use partools, only : myid use emission_data, only : emis_bb_trop_cycle, bb_cycle use datetime, only : tau2date ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 23 Jan 2014 - Jason Williams - V0 ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/emission_nox_bb_daily_cycle' integer :: i,j,l,region, lmr, itim, ntim integer :: i1, i2, j1, j2, ipos, sec_in_day integer, dimension(6) :: idater real :: dtime, dtime2, xlon, xlat ! REG: do region = 1, nregions ! ! Re-initialize the bb NOx array ! nox_emis_3d_bb_app(region)%d3 = 0.0 call tau2date(itaur(region),idater) dtime = float(nsrce)/(2*tref(region)) ! emissions are added in two steps...XYZECCEZYX. dtime2 = float(ndyn_max)/(2*tref(region)) ntim = 86400/nint(dtime2) ! number of timesteps in 24 hours. sec_in_day = idater(4)*3600 + idater(5)*60 + idater(6) ! elapsed seconds this day itim = sec_in_day/nint(dtime2)+1 ! time interval if(okdebug) then write(gol,*)'emission_nox_bb_daily_cycle in region ',region,' at date: ',idater, ' with time step:', dtime,' on ',myid call goPr end if CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) lmr = lm(region) if( emis_bb_trop_cycle) then do l=1,lmr do j=j1,j2 do i=i1,i2 xlat = ybeg(region) + (j-0.5)*dy/yref(region) if (xlat .gt. -20 .and. xlat .lt. 20) then ! apply emission cycle in tropics only ! itim = 1 and lon = -180 --->position 1 ! itim = ntim ant lon = -180 --->position ntim, etc. ! itim = 1 and lon = 0 ---->position ntim/2 xlon = xbeg(region) + (i-0.5)*dx/xref(region) ipos = 1 + mod(int((xlon+180.)*ntim/360.) + (itim-1),ntim) !position in array depending on time and lon. nox_emis_3d_bb_app(region)%d3(i,j,l) = nox_emis_3d_bb(region)%d3(i,j,l)*bb_cycle(region)%scalef(ipos) else nox_emis_3d_bb_app(region)%d3(i,j,l) = nox_emis_3d_bb(region)%d3(i,j,l) endif enddo enddo enddo else nox_emis_3d_bb_app(region)%d3 = nox_emis_3d_bb(region)%d3 endif end do REG if(okdebug) then write(gol,*) 'end of emission_nox_bb_daily_cycle'; call goPr end if status=0 END SUBROUTINE EMISSION_NOX_BB_DAILY_CYCLE #ifndef without_convection !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: lightningNOx ! ! !DESCRIPTION: Calculates NOx emissions from lightning as input for ! photochemistry module. NOx lightning is calculated using a linear ! relationship between lightning flashes and convective precipitation ! ! * total annual production is approximately 5 Tg(N)/yr ! * marine lightning is ten times less active ! * fraction of cloud-to-ground over total flashes is determined by ! 4th order polynomial fit of the cold cloud thickness (Price and ! Rind, GRL 1993). ! * NOx production per IC and CG flash is according to Price et al, ! JGR, 1997. ! * vertical NOx profile is an approximation of the 'outflow' profile ! adopted from Pickering et al., JGR 1998. ! ! Calculate distribution of lightning using cloudtop heights ! of deep convection, cloud cover and convective precipitation ! ! Reference: E. Meijer, KNMI. ! Physics and Chemistry of the earth, Manuscript ST6.03-4 !\\ !\\ ! !INTERFACE: ! SUBROUTINE LIGHTNINGNOX(region, I1, J1, emilig, status) ! ! !USES: ! USE dims, only : im,jm,lm,ybeg,yref use Dims, only : CheckShape USE Binas, only : Avog use chem_param, only : xmn use partools, only : isRoot, par_reduce USE toolbox, only : ltropo, lvlpress USE meteodata, only : m_dat, phlb_dat use meteodata, only : temper_dat, gph_dat, cp_dat USE global_data, only : region_dat, conv_dat USE emission_data, only : plandr ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region, i1, j1 ! ! !OUTPUT PARAMETERS: ! real, intent(out) :: emilig(i1:,j1:,:) ! lighting emissions (kg N/s) integer, intent(out) :: status ! ! !REVISION HISTORY: ! ? ??? 2001 - Ernst Meijer - Set up ! ? ??? 2002 - Olivie, van Weele - Revisions ! ? Jul 2002 - Frank Dentener - adapted for TM5 ! ? Jan 2003 - Maarten krol - adapted for NEW TM5 ! 1 Oct 2010 - Achim Strunk - protex ! 27 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! 21 Aug 2013 - Ph. Le Sager - update fscalelig for 1x1 ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/lightningNOx' real, dimension(:,:,:), pointer :: phlb ! pressure (Pa) (1:lm+1) real, dimension(:,:,:), pointer :: m ! air mass (kg) real, dimension(:,:,:), pointer :: gph ! height (incl. oro) integer,dimension(:,:), pointer :: lutop ! cloud_top integer,dimension(:,:), pointer :: lubottom ! cloud base real,dimension(:), pointer :: dxyp ! area (m) real,dimension(:,:,:), pointer :: t ! temperature (K) real,dimension(:,:,:), pointer :: cp ! convective precipitation (mm/day) real,dimension(:,:), pointer :: plandregion ! landfraction (0-1) ! real,dimension(:),allocatable :: gphx,tx real :: top ! cloudtop (km) real :: cldd ! cold cloud extension (km) (0 deg - tep) integer :: ibase ! base of cloud layer integer :: itropo ! tropopuase layer integer :: itop ! top of cloud layer integer :: itmin15 ! layer with the t=-15 isotherm integer :: itmin24 ! layer with the t=-24 isotherm integer :: it0 ! layer with the t= 0 isotherm ! real :: cpc ! convective precipitation => (m s^-1) ! real :: fl ! flash rate (1/s) real :: cg ! cloud-to-ground flashe rate (1/s) real :: ic ! intra-cloud flash rate (1/s) real :: pnocg,pnoic ! molecules NO produced per CG, IC flash ! DIAGNOSTIC variables real :: flashglob ! global flash frequency (1/s) real :: flashtrop ! tropical flash frequency (1/s) real :: ic_flashglob ! global flash frequency (1/s) real :: ic_flashtrop ! tropical flash frequency (1/s) real :: noxltrop ! tropical NOxL (kg/s) ! logical :: cld ! flag for initialisation and cloud ! integer :: i,j,l,ll,nlay,i2,j2 real :: dsum real :: airmass integer :: level_pblh real :: xlat, dlat, tot_emilig logical :: lightning_output = .false. ! switch for diagnostic output logical :: lightning_output_2 = .false. ! for buggy MPI (see comment in budget_global.F90 for details) real :: flashglob_all, ic_flashglob_all, flashtrop_all, ic_flashtrop_all, noxltrop_all, tot_emilig_all ! --- begin -------------------------------- CALL GET_DISTGRID( dgrid(region), I_STOP=i2, J_STOP=j2 ) call CheckShape( (/i2-i1+1,j2-j1+1,lm(region)/), shape(emilig), status ) IF_NOTOK_RETURN(status=1) m => m_dat (region)%data phlb => phlb_dat (region)%data cp => cp_dat (region)%data plandregion => plandr (region)%surf dxyp => region_dat (region)%dxyp lubottom => conv_dat (region)%cloud_base lutop => conv_dat (region)%cloud_top t => temper_dat (region)%data gph => gph_dat (region)%data ! ! allocate(gphx(0:lm(region))) ! note now from 0-->lm allocate(tx(lm(region))) ! ! FD region coordinates are needed for determining statistics in tropics ! (-30 N,30N) and for excluding polar lightning (75N, 75 S) the emissions ! in parent regions containing a zoom are only set to zero after budget ! calculations using zoomed ! ! Initialising statistics ! flashglob = 0. flashtrop = 0. ic_flashglob = 0. ic_flashtrop = 0. noxltrop = 0. ! ! initialising lightning emission ! emilig(:,:,:) = 0. ! (im,jm,lm) dlat = dy/yref(region) ! do j=j1,j2 xlat = float(ybeg(region)) + j*dlat ! southern edge of gridbox if(xlat < -75.0 .or. (xlat+dlat) > 75.0) cycle ! exclude poles.... do i=i1,i2 ! fl = 0. cldd = 0. cg = 0. ic = 0. ! ibase = 0 itop = 0 itmin24 = 0 itmin15 = 0 it0 = 0 ! cld = .false. ! cpc = cp(i,j,1) !old data cpc = cp(i,j,1) / 1000./86400. ! mm/day => m/s ! if (cpc.gt.0.) then tx(:)=t(i,j,:) gphx(0:lm(region))=gph(i,j,1:lm(region)+1) !note the bounds ibase = lubottom(i,j) itop = lutop(i,j) itropo = ltropo(region,tx,gphx,lm(region)) if (ibase.ne.0.and.itop.ne.0) then do l = itop, ibase, -1 if (tx(l).le.249.15) itmin24=l if (tx(l).le.258.15) itmin15=l if (tx(l).le.273.15) it0 = l enddo !l top = (gphx(itop)+gphx(itop-1))/2000. ! cloud top (km) if (itmin24.ne.0.and.top.gt.5) then cld = .true. ! IF CLOUD REGIONS, IT IS A DEEP CLOUD if (itop.gt.itropo) itop = itropo if (it0 .gt.itropo) it0 = itropo endif ! itmin24 endif !ibase ne 0 ! if (cld) then !fd top = (gphx(itop)+gphx(itop-1))/2000. ! cloud top (km) cldd= top - (gphx(it0)+gphx(it0-1))/2000. ! cold top (km) fl = fscalelig *4.e6 * cpc * dxyp(j)*1.e-12 fl = (0.9*plandregion(i,j) + 0.1) * fl if (cldd.ge.5.5) then cg = fl / (.021*cldd**4-.648*cldd**3+7.493*cldd**2-36.54*cldd+64.09) ic = fl - cg else ! changed from [0.;fl] to [fl; fl-cg] (TvN, PLS, 22-04-2013) ! this increases the LiNOx by ~8.2 TgN/yr cg = fl ic = fl-cg endif !cldd ! changed pnocg factor from 6.7e9 to 6.7e8 (TvN, PLS, 22-04-2013) ! this is from more recent estimate from observations [Ott et al., 2010] pnocg = 1e17*6.7e8*cg *xmn*1e-3/Avog pnoic = 1e17*6.7e8*ic *xmn*1e-3/Avog ! DISTRIBUTION of LNOx over the COLUMN ! ! assume all IC-LNOx and 70% of CG-LNOx betweem t=-15 and cloudtop; ! assume 10% of CG-LNOx between EARTH SURFACE and t=-15 ! assume 20% of CG-LNOx in BOUNDARY LAYER ! ! To avoid dependency on the vertical resolution : ! - surface emission in lowest layers : boundary layer height; ! - LNOx within one of these three regions is distributed proportional to the mass of each layer. ! distributing all IC LNOx and 70% of CG LNOx BETWEEN t=-15 and CLOUDTOP dsum = 0.7*pnocg+pnoic ! determining nlay itop = min(itop,itropo-1) nlay = itop - itmin15 if (nlay.le.0) then itmin15=itop-1 nlay = 1 if (lightning_output_2) write(6,*) 'WARNING noxlight_cvp: itmin15>=itropo: ',i,j,itropo,itmin15 endif !nlay le 0 ! distributing according to airmass airmass = 0. do l=itmin15+1,itop airmass = airmass + m(i,j,l) enddo do l=itmin15+1,itop emilig(i,j,l) = dsum*m(i,j,l)/airmass enddo ! distributing 10% of CG LNOx between EARTH SURFACE and t=-15 dsum = 0.1 * pnocg ! distributing according to air mass airmass = 0. do l=1,itmin15 airmass = airmass + m(i,j,l) enddo do l=1,itmin15 emilig(i,j,l) = emilig(i,j,l) + dsum*m(i,j,l)/airmass enddo ! distributing 20% of CG LNOx between ground pressure and 0.8*ground pressure dsum = 0.2 * pnocg level_pblh = lvlpress(region,0.8*phlb(i,j,1),phlb(i,j,1)) ! distributing according to airmass airmass = 0. do l=1,level_pblh airmass = airmass + m(i,j,l) enddo do l=1,level_pblh emilig(i,j,l) = emilig(i,j,l) + dsum*m(i,j,l)/airmass enddo if (lightning_output) then ! CALCULATE GLOBAL/TROPICAL flash, NOxL rates flashglob = flashglob + fl ic_flashglob = ic_flashglob + ic select case(nint(xlat)) ! xlat is the southern edge of box j.... case(-30:29) flashtrop = flashtrop + fl ic_flashtrop = ic_flashtrop + ic do l = 1, lm(region) noxltrop = noxltrop + emilig(i,j,l) enddo case default end select endif ! lightning_output endif !cld = .true. endif !cpc.gt.0. enddo !i enddo !j if (lightning_output) then call par_reduce( flashglob , 'sum', flashglob_all , status ) call par_reduce( ic_flashglob , 'sum', ic_flashglob_all , status ) call par_reduce( flashtrop , 'sum', flashtrop_all , status ) call par_reduce( ic_flashtrop , 'sum', ic_flashtrop_all , status ) call par_reduce( noxltrop , 'sum', noxltrop_all , status ) tot_emilig = sum(emilig) call par_reduce( tot_emilig, 'sum', tot_emilig_all, status ) if (isRoot) then write(gol,*) 'EMISS-INFO - global lightning frequency = ',flashglob_all,' s-1' ; call goPr write(gol,*) 'EMISS-INFO - ic global lightning frequency = ',ic_flashglob_all,' s-1' ; call goPr write(gol,*) 'EMISS-INFO - tropical lightning frequency = ',flashtrop_all,' s-1' ; call goPr write(gol,*) 'EMISS-INFO - ic tropical lightning frequency= ',ic_flashtrop_all,' s-1' ; call goPr write(gol,*) 'EMISS-INFO - global lightning emission : ',tot_emilig_all*1.e-9*365.*86400.,' Tg[N]/a' ; call goPr write(gol,*) 'EMISS-INFO - tropical lightning emission : ',noxltrop_all*1.e-9*365.*86400.,' Tg[N]/a' ; call goPr end if endif nullify(m) nullify(phlb) nullify(gph) nullify(cp) nullify(plandregion) nullify(lubottom) nullify(lutop) nullify(dxyp) nullify(t) deallocate(gphx) deallocate(tx) status=0 END SUBROUTINE LIGHTNINGNOX #endif END MODULE EMISSION_NOX