ecearth3/sources/tm5mp/proj/co2/emission_data__co2.F90

!
#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_DATA
!
! !DESCRIPTION: Provides general methods and fields needed for the emission
!               routines. 
!
!    emission_vdist_by_sector : distribute vertically emissions according to source sector
!    msg_emis                 : print formated emission total
!    do_add_2d                : add 2D emission to tracer mass array
!    do_add_3d                : add 3D emission to tracer mass array
!\\
!\\
! !INTERFACE: 
!
MODULE EMISSION_DATA
  !
  ! !USES:
  !
  use GO,           only : gol, goPr, goErr
  use Dims,         only : nregions, lm
  use global_types, only : emis_data, d2_data, d23_data, d3_data

#ifdef with_budgets
  use budget_global,only : budg_dat, nzon_vg
#endif

  IMPLICIT NONE
  !
  ! !PUBLIC DATA MEMBERS & TYPES:
  !
  public :: emis_data, d2_data, d23_data
  !
  type(emis_data), target    :: plandr(nregions)    ! Land fraction
  type(emis_data), target    :: emis2D(nregions)    ! 2D emiss prior vertical redistribution
  type(emis_data), target    :: emis_temp(nregions) ! temp array to store emission

  ! AR5
  logical            :: LAR5
  character(len=256) :: emis_input_dir
  character(len=256) :: emis_input_dir_ar5
  integer            :: emis_input_year

  ! biomass burning levels:
  integer, parameter ::  bb_nlev = 5
  integer, parameter ::  bb_lm   = lm(1)    ! to be reduced to strip stratosphere!
  integer, parameter ::  bl_lm   = lm(1)    ! to be reduced to the BL (around 8)
  real,    parameter ::  bb_hlow (0:bb_nlev) = (/ 0.,   100.,  500., 1000., 2000., 3000./)
  real,    parameter ::  bb_hhigh(0:bb_nlev) = (/ 100., 500., 1000., 2000., 3000., 6000./)

  ! biomass burning levels in tropics (-20 - 20)::
  integer, parameter ::  bb_nlev_trop = 3
  real,    parameter ::  bb_hlow_trop (0:bb_nlev_trop) = (/ 0.,   100.,  500., 1000./)
  real,    parameter ::  bb_hhigh_trop(0:bb_nlev_trop) = (/ 100., 500., 1000., 2000./)

  ! -----------------------------------------------------------------------------------
  ! Budgets
#ifdef with_budgets
  real, dimension(nregions)                    :: sum_emission

  type budemi_data
     real,dimension(:,:,:,:),pointer           :: emi    ! [i, j, nbud_vg, ntracet]
  end type budemi_data
  type(budemi_data),dimension(nregions),target :: budemi_dat
#endif

  ! ----------------------------------------------------
  ! Vertical splitting : used now for all sectors
  ! (tables in Dentener, 2006, ACP)
  ! 
  ! ATTENTION: changes here have impact on the routine emission_vdist_by_sector (contained)
  ! Note: vdist_* variable could be hold private
  !
  integer, parameter                      ::  vd_class_name_len = 12
  ! biomassburning
  integer, parameter                      ::  vdist_bb_nlev  = 6
  real,dimension(vdist_bb_nlev),parameter ::  vdist_bb_hlow  = (/ 0.,   100.,  500., 1000., 2000., 3000./)
  real,dimension(vdist_bb_nlev),parameter ::  vdist_bb_hhigh = (/ 100., 500., 1000., 2000., 3000., 6000./)
  ! other sources (related to surface)
  integer, parameter                      ::  vdist_nn_nlev  = 6
  real,dimension(vdist_nn_nlev),parameter ::  vdist_nn_hlow  = (/ 0.,  30. , 100., 300.,  600., 1000. /)
  real,dimension(vdist_nn_nlev),parameter ::  vdist_nn_hhigh = (/ 30., 100., 300., 600., 1000., 2000. /)
  !
  ! !PRIVATE TYPES:
  !
  character(len=*), parameter, private  ::  mname = 'emission_data'
  !
  ! !REVISION HISTORY: 
  !    1 Oct 2010 - Achim Strunk - adapted for AR5
  !   28 Jun 2012 - Ph. Le Sager - made everything public by default
  !                              - adapted for lon-lat MPI domain decomposition
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_VDIST_BY_SECTOR
  !
  ! !DESCRIPTION: Return vertically distributed emissions depending on given
  !               sector and constituent.
  !               Distribution is done as 'compromise' between 
  !               1)  Dentener et al., 2006, ACP, slightly modified for
  !                   emissions supposed to enter the surface layer. 
  !               2)  De Meij et al., ACP, 2006
  !               3)  EMEP model (docu, model code, publications)
  !               4)  Bieser et al., GMDD, 2010 
  !               5)  Small updates to biomassburning following 
  !                   Huijnen et al., GMDD, 2010
  !\\
  !\\
  ! !INTERFACE:
  !  
  SUBROUTINE EMISSION_VDIST_BY_SECTOR( splitname, constituent, region, emis2D, emis3D, status )
    !
    ! !USES:
    !
    use toolbox,         only : distribute_emis2D
    !
    ! !INPUT PARAMETERS:
    !
    character(len=vd_class_name_len), intent(in)          :: splitname
    character(len=*),                 intent(in)          :: constituent
    integer,                          intent(in)          :: region
    type(emis_data),                  intent(in)          :: emis2D
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    type(d3_data),                    intent(inout)       :: emis3D
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,                          intent(out)         :: status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - v0
    !
    ! !REMARKS: Splitting depending on constituent is still to code. 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/emission_vdist_by_sector'

    integer                         :: ilev
    character(len=2)                :: splittype
    integer, parameter              :: maxvdist_nlev = max(vdist_bb_nlev,vdist_nn_nlev)
    real, dimension(maxvdist_nlev)  :: sfglobal, sfintrop, sfintemp

    ! --- begin --------------------------------------
    ! control implicit settings here in this routine
    if( vdist_bb_nlev /= 6  .or. vdist_nn_nlev /= 6 ) then 
       write(gol, '("ERROR: wrong number of layers: bb_nlev /= 6  or  nn_lev /= 6")'); call goErr
       status=1; return
    end if

    ! initialise
    sfglobal = 0.0 
    sfintrop = 0.0 
    sfintemp = 0.0 

    ! -------------------------------------------------
    ! This here is selection by emission source sectors 
    ! -------------------------------------------------
    select case( trim(splitname) )

    case( 'combenergy' ) 
       ! Combustion from energy production and transformation, power-plants
       ! no zonal differences, no dependency on species
       ! --
       ! assumed to be related to EMEP SNAP category [1]
       ! AEROCOM equivalent: power-plants
       ! --
       splittype = 'nn' 
       sfglobal  = (/ 0.0, 0.1, 0.7, 0.2, 0.0, 0.0 /) 

    case( 'combrescom' )
       ! Residential and commercial combustion
       ! no zonal differences, no dependency on species
       ! --
       ! assumed to be related to EMEP SNAP category [2]
       ! AEROCOM equivalent: domestic/industry
       ! --
       splittype = 'nn' 
       sfglobal  = (/ 0.4, 0.4, 0.2, 0.0, 0.0, 0.0 /) 

    case( 'industry' )
       ! Industrial processes and combustion
       ! no zonal differences, no dependency on species
       ! --
       ! assumed to be an aggregate of EMEP SNAP categories [3,4,5]
       ! AEROCOM equivalent: industry
       ! --
       splittype = 'nn' 
       sfglobal  = (/ 0.1, 0.2, 0.6, 0.1, 0.0, 0.0 /)

    case( 'waste' )
       ! Waste treatment and disposal
       ! no zonal differences, no dependency on species
       ! --
       ! assumed to be related to EMEP SNAP category [9]
       ! AEROCOM equivalent: - 
       ! --
       splittype = 'nn' 
       sfglobal  = (/ 0.1, 0.2, 0.4, 0.3, 0.0, 0.0 /)

    case( 'nearsurface' )
       ! Near surface emissions 
       ! no zonal differences, no dependency on species
       ! --
       ! AR5: solvent production and use, agricultural waste burning, &
       !      ships, grassland fire, mineral dust (AEROCOM or Tegen-Vignati)
       ! EMEP equivalents: SNAP categories [6,8]
       ! AEROCOM equivalents: ships, agricultural waste
       ! --
       splittype = 'nn' 
       sfglobal  = (/ 0.8, 0.2, 0.0, 0.0, 0.0, 0.0 /)

    case( 'surface' )
       ! Surface emissions
       ! no zonal differences, no dependency on species
       ! --
       ! AR5: agriculture, transport, soil, oceanic, biogenic
       ! EMEP equivalents: SNAP categories [7,8,10]
       ! AEROCOM equivalents: surface emissions from road/off-road, ships, agricultural waste
       ! --
       splittype = 'nn' 
       sfglobal  = (/ 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)

    case( 'forestfire' )
       ! forest fires: 6 layer model, different splitting for different regions
       !               dont distinguish between boreal Europe and boreal Canada, use Europe globally
       ! Dentener et al., modified in tropics (up to 2 km) with respect to Huijnen et al., GMD 2010, \
       !  who are citing Labonne et al., 2007 ofr new information based on satellite data.
       splittype = 'bb' 
       sfglobal  = (/ 0.1, 0.1, 0.2, 0.2, 0.4, 0.0 /)
       sfintemp  = (/ 0.2, 0.2, 0.2, 0.4, 0.0, 0.0 /)
       sfintrop  = (/ 0.2, 0.2, 0.2, 0.4, 0.0, 0.0 /)

    case default
       write(gol, '("ERROR: wrong specification of emission sector ",a,"in vdist_by_sector called for constituent ",a)') &
            trim(splitname), trim(constituent); call goErr
       status=1; return

    end select

    ! --------------------
    ! do the splitting 
    ! --------------------
    select case( splittype )
    case( 'nn' )
       do ilev = 1, vdist_nn_nlev
          call distribute_emis2D( region, emis2D, emis3D, vdist_nn_hlow(ilev), vdist_nn_hhigh(ilev), status, sfglobal(ilev) )
          IF_NOTOK_RETURN(status=1)
       end do
    case( 'bb' )
       do ilev = 1, vdist_bb_nlev
          ! boreal 
          call distribute_emis2D( region, emis2D, emis3D, vdist_bb_hlow(ilev), vdist_bb_hhigh(ilev), status, sfglobal(ilev),  -90.,-61.)
          IF_NOTOK_RETURN(status=1)
          call distribute_emis2D( region, emis2D, emis3D, vdist_bb_hlow(ilev), vdist_bb_hhigh(ilev), status, sfglobal(ilev),   60., 90.)
          IF_NOTOK_RETURN(status=1)
          ! temperate 
          call distribute_emis2D( region, emis2D, emis3D, vdist_bb_hlow(ilev), vdist_bb_hhigh(ilev), status, sfintemp(ilev),  -60.,-31.)
          IF_NOTOK_RETURN(status=1)
          call distribute_emis2D( region, emis2D, emis3D, vdist_bb_hlow(ilev), vdist_bb_hhigh(ilev), status, sfintemp(ilev),   30., 59.)
          IF_NOTOK_RETURN(status=1)
          ! tropics
          call distribute_emis2D( region, emis2D, emis3D, vdist_bb_hlow(ilev), vdist_bb_hhigh(ilev), status, sfintrop(ilev),  -30., 29.)
          IF_NOTOK_RETURN(status=1)
       end do
    case default
       write(gol, '("ERROR: wrong specification of splitting type in vdist_by_sector.")'); call goErr
       status=1; return
    end select

    ! OK 
    status=0

  END SUBROUTINE EMISSION_VDIST_BY_SECTOR
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    MSG_EMIS
  !
  ! !DESCRIPTION: Provide output to verify the amount of emissions 
  !               entering the calculations
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE MSG_EMIS(year_month, msg, msg_comp, msg_mass, summed_emissions)
    !
    ! !USES:
    !
    use dims, only: idate
    !
    ! !INPUT PARAMETERS:
    !
    integer,          intent(in) :: year_month      ! 1: year, 2: monthly amount
    character(len=*), intent(in) :: msg, msg_comp
    real,             intent(in) :: msg_mass, summed_emissions
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - protex; increase space for output
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=7),dimension(2), parameter:: ym=(/'Yearly ','Monthly'/)

1111 format('EMISS-INFO -   ',a7,' Emission category ',a34,1x,a6,' [kg], mw:',f5.1,', month=',i2.2,', Sum',1x,1pe10.4)

    write (gol,1111) ym(year_month), msg, msg_comp, msg_mass, idate(2), summed_emissions; call goPr

  END SUBROUTINE MSG_EMIS
  !EOC
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DO_ADD_2D
  !
  ! !DESCRIPTION: General routine to add a 2D emission field (given in kg
  !                /month) into tracer mass array. The mass increase is done
  !                at level L_LEVEL, and for tracer I_TRACER.
  !                
  !               Update accordingly the budget.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DO_ADD_2D( region,   i_tracer,   l_level, is, js, emis_field, &
                        mol_mass, mol_mass_e, status,  xfrac               )
    !
    ! !USES:
    !
    use dims,          only : CheckShape, tref, nsrce
    use dims,          only : sec_month, im, jm, okdebug
    use global_data,   only : mass_dat, region_dat
    use chem_param,    only : ntracet, names
    use tm5_distgrid,  only : get_distgrid, dgrid
#ifdef with_budgets
    use budget_global, only : budg_dat, nzon_vg
#endif

    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)          :: region  
    integer, intent(in)          :: i_tracer 
    integer, intent(in)          :: l_level
    integer, intent(in)          :: is, js
    real, intent(in)             :: emis_field(is:,js:)
    real, intent(in)             :: mol_mass    ! mol mass of tracer field
    real, intent(in)             :: mol_mass_e  ! mol mass of emission field (e.g. NOx emission ar in kgN/month)
    real, intent(in), optional   :: xfrac       ! partial addition 
    !
    ! !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//'/do_add_2d'

    integer :: i, j, nzone, nzone_v, i1, i2, j1, j2
    real    :: x, efrac, dtime, month2dt
    real    :: mbef, maft, sume  ! mass BEFore, AFTer, SUM Emiss (debug)
    
    real, dimension(:,:,:,:), pointer :: rm, rzm   

    ! --- begin --------------------------------------

    call GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )          
    
    call CheckShape( (/i2-i1+1,j2-j1+1/), shape(emis_field), status )
    IF_NOTOK_RETURN(status=1)
    
    dtime=float(nsrce)/(2*tref(region))  ! timestep emissions
    month2dt=dtime/sec_month             ! conversion to emission per timestep
   
    rm  => mass_dat(region)%rm
    rzm => mass_dat(region)%rzm

    if(present(xfrac)) then
       efrac = xfrac
    else
       efrac = 1.0
    endif

    if (okdebug) then 
       write (gol, '("--------------------------------------------------------------")')            ; call goPr
       write (gol, '(" DO ADD 2D debug in region:",i3) '    ) region                                ; call goPr
       write (gol, '(" itracer                  :",i3,a8)'  ) i_tracer, ' :'//names(i_tracer)       ; call goPr
       write (gol, '(" xfrac                    :",f6.1)'   ) efrac                                 ; call goPr
       write (gol, '(" mol_mass                 :",f6.1 ) ' ) mol_mass                              ; call goPr
       write (gol, '(" mol_mass_e               :",f6.1 )'  ) mol_mass_e                            ; call goPr
       write (gol, '(" emisfield                :",2e12.4 )') minval(emis_field),maxval(emis_field) ; call goPr
       sume =  sum(emis_field)*month2dt*(mol_mass/mol_mass_e)*efrac
       write (gol, '(" sum emmis                :",e12.4 )' ) sume                                  ; call goPr
       mbef = sum(rm(i1:i2,j1:j2,:,i_tracer))
    endif

    do j=j1,j2
       do i=i1,i2
          x = emis_field(i,j)*month2dt*(mol_mass/mol_mass_e)*efrac
          rm(i,j,l_level,i_tracer) = rm(i,j,l_level,i_tracer) + x
#ifdef slopes
          ! injection of emissions at surface, thus vertical slope more negative
          ! (keep concentration at top of layer the same as before emission)
          if ( i_tracer <= ntracet ) then
             rzm(i,j,l_level,i_tracer) = rzm(i,j,l_level,i_tracer) - x
          end if
#endif

#ifdef with_budgets
          nzone   = budg_dat(region)%nzong(i,j)    !global budget
          nzone_v = nzon_vg(l_level) 
          budemi_dat(region)%emi(i,j,nzone_v,i_tracer) = &
               budemi_dat(region)%emi(i,j,nzone_v,i_tracer) + x/mol_mass*1e3
          if(i_tracer == 1) sum_emission(region) = sum_emission(region) + x  !in kg
#endif
       enddo
    enddo

    if (okdebug) then
       maft = sum(rm(i1:i2,j1:j2,:,i_tracer))
       write (gol, '(" Added amount             :",e12.4 )' ) maft-mbef                ; call goPr
       write (gol, '(" Added amount can be different when inner zoom is present!" )' ) ; call goPr
       if (maft-mbef-sume > 1e-8*sume) then
          write (gol, '(" Warning: error in emission ")' )  ; call goErr
       endif
    endif
    nullify(rm)
    nullify(rzm)

    status=0

  END SUBROUTINE DO_ADD_2D
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DO_ADD_3D
  !
  ! !DESCRIPTION: General routine to add a 3D emission field (given in kg
  !                /month) into tracer mass array. The mass increase is done
  !                for tracer I_TRACER.
  !                
  !               Update accordingly the budget.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DO_ADD_3D( region, i_tracer, is, js, emis_field, mol_mass, mol_mass_e, status, xfrac)
    !
    ! !USES:
    !
    use dims,          only : CheckShape, nsrce, tref
    use dims,          only : sec_month, im, jm, lm, okdebug
    use global_data,   only : mass_dat, region_dat
    use chem_param,    only : ntracet, names
    use tm5_distgrid,  only : get_distgrid, dgrid
#ifdef with_budgets
    use budget_global, only : budg_dat, nzon_vg
#endif
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)           :: region   ! region #
    integer, intent(in)           :: i_tracer ! id of tracer to increase
    integer, intent(in)           :: is, js     
    real,    intent(in)           :: emis_field(is:,js:,:)
    real,    intent(in)           :: mol_mass, mol_mass_e
    real,    intent(in), optional :: xfrac                 !partial addition 
    !
    ! !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//'/do_add_3d'

    integer                             :: i, j, l, nzone, nzone_v, i1,j1,i2,j2
    integer, dimension(3)               :: ubound_e
    real                                :: x, efrac, dtime, month2dt
    real, dimension(:,:,:,:), pointer   :: rm
    
    real :: mbef, maft, sume ! debug

    ! --- begin --------------------------------------

    call GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    status = 0
    
    ! check arguments
    ubound_e = ubound(emis_field)
    if(ubound_e(1) /= i2) then
       write(gol,'(A)') 'do_add_3d: im emission not consistent' ; call goPr
       status = 1
    endif
    if(ubound_e(2) /= j2) then
       write(gol,'(A)') 'do_add_3d: jm emission not consistent' ; call goPr
       status = 1
    endif
    if(ubound_e(3) > lm(region)) then
       write(gol,'(A)') 'do_add_3d: lm emission not consistent' ; call goPr
       status = 1
    endif
    IF_NOTOK_RETURN(status=1)
    
    rm => mass_dat(region)%rm

    dtime=float(nsrce)/(2*tref(region))  ! timestep emissions
    month2dt=dtime/sec_month             ! conversion to emission per timestep

    if(present(xfrac)) then
       efrac = xfrac
    else
       efrac = 1.0
    endif

!     if (okdebug) then 
!        write (gol, '("--------------------------------------------------------------")')             ; call goPr
!        write (gol, '(" DO ADD 3D debug in region:",i3) '    ) region                                 ; call goPr
!        write (gol, '(" itracer                  :",i3,a8)'  ) i_tracer, ' :'//names(i_tracer)        ; call goPr
!        write (gol, '(" xfrac                    :",f6.1)'   ) efrac                                  ; call goPr
!        write (gol, '(" mol_mass                 :",f6.1 ) ' ) mol_mass                               ; call goPr
!        write (gol, '(" mol_mass_e               :",f6.1 )'  ) mol_mass_e                             ; call goPr
!        write (gol, '(" emisfield                :",2e12.4 )') minval(emis_field), maxval(emis_field) ; call goPr
!        write (gol, '(" ubound(3)                :",i3 )'    ) ubound_e(3)                            ; call goPr
!        sume = sum(emis_field)*month2dt*(mol_mass/mol_mass_e)*efrac
!        write (gol, '(" sume                     :", e12.4 )') sume                                   ; call goPr
!        mbef = sum(rm(i1:i2,j1:j2,:,i_tracer))
!     endif

    do l=1,ubound_e(3)
       do j=j1,j2
          do i=i1,i2
             x = emis_field(i,j,l)*month2dt*(mol_mass/mol_mass_e)*efrac
             rm(i,j,l,i_tracer) = rm(i,j,l,i_tracer) + x

#ifdef with_budgets
             ! budget___
             nzone=budg_dat(region)%nzong(i,j)    !global budget
             nzone_v=nzon_vg(l) 
             budemi_dat(region)%emi(i,j,nzone_v,i_tracer) = &
                  budemi_dat(region)%emi(i,j,nzone_v,i_tracer) + x/mol_mass*1e3
             if(i_tracer ==1) sum_emission(region) = sum_emission(region) + x  !in kg
#endif
          enddo
       enddo
    enddo  !levels

!     if (okdebug) then
!        maft = sum(rm(i1:i2,j1:j2,:,i_tracer))
!        write (gol, '(" after-before             :", e12.4 )' ) maft-mbef ; call goPr
!        write (gol, '(" END debug output          ")' )                   ; call goPr
!     endif

    nullify(rm)

  END SUBROUTINE DO_ADD_3D
  !EOC

END MODULE EMISSION_DATA