ecearth3/sources/tm5mp/base/budget_global.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"
#include "output.inc"
!
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:      BUDGET_GLOBAL
!
! !DESCRIPTION: This module holds most of the budget variables (see
!               wet_deposition for exceptions). 
! 
! Budgets are defined for a set of horizontal regions, split vertically, and
! available for each transported tracers. In the module:
! 
!    o  defines the budget output filename
!    o  defines the horizontal 'regions' for which a budgets are calculated
!            e.g. Europe, HNH, Africa, zonal bands
!    o  defines the vertical split of the budgets
!            e.g. boundary layer, free troposphere, stratosphere
!    o  routines to compute the transport budgets (advx, advy, advz, conv)
!    o  sets up the possibility to calculate transport fluxes in/out a
!        predefined XY, XZ, and YZ planes (2 of each => define a box) 
!    o  other budget terms (emission, deposition, etc.) can be added to the budget file
!
!\\
!\\
! !INTERFACE: 
!
MODULE BUDGET_GLOBAL
  !
  ! !USES:
  !
  use GO,           only : gol, goErr, goPr
  use dims,         only : im, jm, lm, idatei, idatee, nregions, region_name
  use dims,         only : itau, idate, dx, dy, xref, yref, xbeg, ybeg
  use chem_param,   only : ntracet, ntrace, names, ra
  use tm5_distgrid, only : dgrid, Get_DistGrid, reduce, gather

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: init_budget_global  ! initialize and allocate variables
  public :: done_budget_global  ! final computation, write budgets to file
  public :: budget_transportg   ! evaluate advection/convection budgets
  public :: budget_displayg     ! print out budget of tracer #1
  !
  ! !PUBLIC DATA MEMBERS:
  !
  character(len=256), public    :: budget_file_global  ! name of budget output file
  integer, public, parameter    :: nbudg   = 7         ! # horizontal budget regions (h-budgets)
  integer, public, parameter    :: nbud_vg = 4         ! # vertical   budget regions (v-budgets)
  logical, public               :: apply_budget_global      = .false. ! compute budgets     
  logical, public               :: apply_budget_global_flux = .false. ! compute extra fluxes
  logical, public               :: NHAB ! Non-Horizontally-Aggregated-Budgets (if true (default), writes [im,jm,nbud_vg] budgets
                                        !   which are available in the full chemistry (proj/cbm4 & proj/cb05)
  real, public, dimension(nbudg, nbud_vg, ntracet) :: budconvg ! convection budget

  real, public, dimension(nregions)  :: init_mass = 0.0     ! tracer #1 : initial mass (meaningful on root only)
  real, public, dimension(nregions)  :: sum_update = 0.0    ! tracer #1 : 
  real, public, dimension(nregions)  :: sum_advection = 0.0 ! tracer #1 : mass change thru advection (meaningful on root only)

  integer, public, dimension(lm(1))  :: nzon_vg  ! vertical budget code of each level

  integer, public, dimension(nregions)  :: iflux1, iflux2  ! I index of YZ plane where flux flowing in from west  
  integer, public, dimension(nregions)  :: jflux1, jflux2  ! J index of XZ plane where flux flowing in from south
  integer, public, dimension(nregions)  :: lflux1, lflux2  ! L index of XY plane where flux flowing in from below
  !
  ! !PUBLIC TYPES:
  !
  type, public :: budg_data
     ! holds horizontal budget code of each cell (1:SH 90-30, 2:S30-N30, 3:NH 30-90, 4: Europe...)
     integer, dimension(:,:), pointer    :: nzong
  end type budg_data

  type(budg_data),    public, dimension(nregions), target :: budg_dat    ! horizontal budget codes

  type, public :: budflux_data
     ! fluxes of all species through some predefined slices:
     real, dimension(:,:,:), pointer     :: flux_x1   ! (jm,lm,ntracet)
     real, dimension(:,:,:), pointer     :: flux_x2   ! (jm,lm,ntracet)
     real, dimension(:,:,:), pointer     :: flux_y1   ! (im,lm,ntracet)
     real, dimension(:,:,:), pointer     :: flux_y2   ! (im,lm,ntracet)
     real, dimension(:,:,:), pointer     :: flux_z1   ! (im,jm,ntracet)
     real, dimension(:,:,:), pointer     :: flux_z2   ! (im,jm,ntracet)
  end type budflux_data

  type(budflux_data), public, dimension(nregions), target :: budget_flux
  type(budflux_data),         dimension(nregions), target :: budget_flux_all
  !
  ! !PRIVATE DATA MEMBERS:
  ! 
  real, dimension(nbudg, nbud_vg, ntracet) :: rmoldg, rmnewg   ! work arrays (old/new rm)
  real, dimension(nbudg, nbud_vg, ntracet) :: budconcg         ! final   concentrations (end run)
  real, dimension(nbudg, nbud_vg, ntracet) :: budinig          ! initial concentrations (begin run)
  real, dimension(nbudg, nbud_vg, ntracet) :: budchngg = 0.0   ! change of concentration b/w end-begin
  real, dimension(nbudg, nbud_vg, ntracet) :: budadvxg         ! advection-x budget
  real, dimension(nbudg, nbud_vg, ntracet) :: budadvyg         ! advection-y budget
  real, dimension(nbudg, nbud_vg, ntracet) :: budadvzg         ! advection-z budget

  character(len=20), dimension(nbudg)      :: zone_nameg ! name of horizontal budget regions
  character(len=20), dimension(nbud_vg)    :: zone_namvg ! name of vertical   budget regions
                                            
  integer                                  :: lon_in, lon_out, lat_in, lat_out ! read box boundaries for extra fluxes
                                            
  character(len=*), parameter              :: mname = 'budget_global'  ! module name
  !
  ! !REVISION HISTORY: 
  !   28 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !  (1) All processors compute the same aggregated budgets in case of MPI, but 
  !       account only for the contribution of their own tile, if any. MPI
  !       reduction is used to sum up contributions.
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    INI_ZONEG
  !
  ! !DESCRIPTION: set structure of user-defined "budget regions"
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE INI_ZONEG ( status )
    !
    ! !USES:
    !
    use GO         , only : TrcFile, Init, Done, ReadRc
    use global_data, only : outdir,   rcfile
    use toolbox    , only : lvlpress
    use partools   , only : isRoot
#ifdef with_hdf4          
    use file_hdf,    only : THdfFile, TSds
    use file_hdf,    only : Init, Done, WriteAttribute, WriteData
#endif
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)   :: status
    !
    ! !REVISION HISTORY: 
    !   28 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

#ifdef with_hdf4          
    type(THdfFile)       :: io
    type(TSds)           :: sds
#endif
    integer              :: i1, i2, j1, j2, i, j, n, nix, niy, region
    integer              :: ilvl1, ilvl2, ilvl3
    real                 :: psurf, plev, rlat, rlon
    real                 :: dyr, dxr, lats  ! made real instead of integer to handle hi-res regions
    integer, allocatable :: gbl_data(:,:)

    type(TrcFile)       ::  rcF
    character(len=256)  ::  fname

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

    ! ------- Compute budgets and extra fluxes ?
#ifdef with_budgets
    apply_budget_global = .true.
#endif

    if (.not.apply_budget_global) return

    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)

    call ReadRc(rcF, 'apply.budget.global.flux', apply_budget_global_flux, status, default=.false.)
    IF_ERROR_RETURN(status=1)

    if (apply_budget_global_flux) then
       
       ! Example Europe: lat_in = 34, lat_out = 62, lon_in = -12, lon_out  = 36

       call ReadRc(rcF, 'budget.global.flux.lon_in', lon_in, status )
       IF_NOTOK_RETURN(status=1)
       call ReadRc(rcF, 'budget.global.flux.lon_out', lon_out, status )
       IF_NOTOK_RETURN(status=1)
       call ReadRc(rcF, 'budget.global.flux.lat_in', lat_in, status )
       IF_NOTOK_RETURN(status=1)
       call ReadRc(rcF, 'budget.global.flux.lat_out', lat_out, status )
       IF_NOTOK_RETURN(status=1)
    endif

    call ReadRc(rcF, 'write.column.budgets', NHAB, status, default=.true.)
    IF_ERROR_RETURN(status=1)
    
    call Done( rcF, status )
    IF_NOTOK_RETURN(status=1)

    
    ! -------- Fill budgets name, code and vertical indices
    region = 1 

    zone_namvg(1)='BL'
    zone_namvg(2)='850-500hPa'                         
    zone_namvg(3)='500-100hPa'                        
    zone_namvg(4)='100-0hPa'

    zone_nameg(1)='SH 90-30'
    zone_nameg(2)='30S-30N '
    zone_nameg(3)='NH 90-30'
    zone_nameg(4)='Europe'
    zone_nameg(5)='N-america'
    zone_nameg(6)='Asia'
    zone_nameg(7)='Africa'

    if ( isRoot ) then
       write(gol,'(A)') '-----------------------------------------------'; call goPr
       write(gol,'(A)') 'Horizontal definition of budget regions        '; call goPr
       write(gol,'(A)') '-----------------------------------------------'; call goPr
       do i=1,nbudg
          write(gol,'(A,i4,4x,A)') 'Budget region:',i,zone_nameg(i); call goPr
       enddo
       !write(gol,'(A)') '-----------------------------------------------'; call goPr
    endif

    ! Associate a v-budget code to each level
    psurf = 98400.0
    plev  = 85000.0 ; ilvl1 = lvlpress(region,plev, psurf)
    plev  = 50000.0 ; ilvl2 = lvlpress(region,plev, psurf)
    plev  = 10000.0 ; ilvl3 = lvlpress(region,plev, psurf)
    
    nzon_vg(1:ilvl1)      = 1 !boundary layer 
    nzon_vg(ilvl1+1:ilvl2)= 2 !850-->500 hPa
    nzon_vg(ilvl2+1:ilvl3)= 3 !500-->100 hPa
    nzon_vg(ilvl3+1:lm(1))= 4 !100 hPa-->TOA

    if ( isRoot ) then
       write(gol,'(A)') '-----------------------------------------------'; call goPr
       write(gol,'(A)') 'Distribution models levels over budget levels  '; call goPr
       write(gol,'(A)') '-----------------------------------------------'; call goPr
       write(gol,'(A,A20,i4,A3,i4)') 'Budget level 1:', zone_namvg(1), 1,      '...',ilvl1 ; call goPr
       write(gol,'(A,A20,i4,A3,i4)') 'Budget level 2:', zone_namvg(2), ilvl1+1,'...',ilvl2 ; call goPr
       write(gol,'(A,A20,i4,A3,i4)') 'Budget level 3:', zone_namvg(3), ilvl2+1,'...',ilvl3 ; call goPr
       write(gol,'(A,A20,i4,A3,i4)') 'Budget level 4:', zone_namvg(4), ilvl3+1,'...',lm(1) ; call goPr
       write(gol,'(A)') '-----------------------------------------------'; call goPr
    endif

    if (apply_budget_global_flux) then ! Limited to free troposphere

       lflux1(:) = ilvl1+1    ! where to calculate the vertical flux (1 = BL->FT)
       lflux2(:) = ilvl2+1    ! where to calculate the vertical flux (2 = FT->UFT)

    endif

    !---------- Allocate/fill arrays w/r/t horizontal definitions
    REG: do n=1,nregions

       call Get_DistGrid( dgrid(n), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
       
       ! allocate horizontal budget code :
       allocate( budg_dat(n)%nzong(i1:i2,j1:j2))

       ! allocate flux data:
       if (apply_budget_global_flux) then
          allocate( budget_flux(n)%flux_x1(j1:j2,lm(n),ntracet))
          allocate( budget_flux(n)%flux_y1(i1:i2,lm(n),ntracet))
          allocate( budget_flux(n)%flux_z1(i1:i2,j1:j2,ntracet))
          allocate( budget_flux(n)%flux_x2(j1:j2,lm(n),ntracet))
          allocate( budget_flux(n)%flux_y2(i1:i2,lm(n),ntracet))
          allocate( budget_flux(n)%flux_z2(i1:i2,j1:j2,ntracet))

          if (isRoot) then
             allocate( budget_flux_all(n)%flux_x1(jm(n),lm(n),ntracet))
             allocate( budget_flux_all(n)%flux_y1(im(n),lm(n),ntracet))
             allocate( budget_flux_all(n)%flux_z1(im(n),jm(n),ntracet))
          else
             allocate( budget_flux_all(n)%flux_x1(1,1,1))
             allocate( budget_flux_all(n)%flux_y1(1,1,1))
             allocate( budget_flux_all(n)%flux_z1(1,1,1))
          end if
          
          budget_flux(n)%flux_x1(:,:,:) = 0.0
          budget_flux(n)%flux_y1(:,:,:) = 0.0
          budget_flux(n)%flux_z1(:,:,:) = 0.0
          budget_flux(n)%flux_x2(:,:,:) = 0.0
          budget_flux(n)%flux_y2(:,:,:) = 0.0
          budget_flux(n)%flux_z2(:,:,:) = 0.0
       endif
       
       
       dyr = dy/yref(n)
       dxr = dx/xref(n)

       ! Fill horizontal budget codes for zonal bands
       do j=j1,j2          
          lats = ybeg(n) + (j-1)*dyr
       
! PREV: not correct with resolutions less than 1 degree
!           select case(lats)
!           case(-90:-31)
!              budg_dat(n)%nzong(:,j) = 3  !SH 90-30
!           case(-30:29)
!              budg_dat(n)%nzong(:,j) = 2  !-30->30
!           case(30:90)
!              budg_dat(n)%nzong(:,j) = 1  !NH 30-90
!           case default
!              write (gol,'(A)') 'ini_zoneg: Wrong latitude in budget routine' ; call goErr
!              IF_NOTOK_RETURN(status=1)
!           end select

          if ( lats >= -90. .and. lats < -30.) then
             budg_dat(n)%nzong(:,j) = 3  !SH 90-30
          elseif ( lats >= -30. .and. lats < 30.) then
             budg_dat(n)%nzong(:,j) = 2  !-30->30
          elseif ( lats >= 30. .and. lats < 90.) then
             budg_dat(n)%nzong(:,j) = 1  !NH 30-90
          else
             write (gol,'(A)') 'ini_zoneg: Wrong latitude in budget routine' ; call goErr
             IF_NOTOK_RETURN(status=1)
          end if

       end do

       
       ! Fill horizontal budget codes for continents : Europe, NAM, Asia, Africa
       do j=j1,j2
          rlat = ybeg(n) + (float(j)-0.5)*dyr
          do i=i1,i2
             rlon = xbeg(n) + (float(i)-0.5)*dxr
             if(rlon > lon_in  .and. rlon < lon_out .and. rlat > lat_in .and. rlat < lat_out)  then
                budg_dat(n)%nzong(i,j) = 4   ! EUROPE
             elseif(rlon > -126. .and. rlon < -66. .and. rlat > 22. .and. rlat < 62.)  then
                budg_dat(n)%nzong(i,j) = 5   ! NAM
             elseif(rlon > 72. .and. rlon < 150. .and. rlat > -10. .and. rlat < 50.)  then
                budg_dat(n)%nzong(i,j) = 6   ! Asia
             elseif(rlon > -18. .and. rlon < 48. .and. rlat > -34. .and. rlat < 34.)  then
                budg_dat(n)%nzong(i,j) = 7   ! Africa
             endif
          enddo
       enddo

       ! --- Calculate in each region where to calculate the flux:
       !
       ! Get the index of the flux planes. The plane is always exactly between cells.
       ! The cell i-index east of the YZ-plane will be put in iflux(n) (for each region).
       ! iflux1 is the western plan (lon_in) and iflux2 is the eastern plane (lon_out).
       ! The same way jflux1 and jflux2 will become the j-index of the cells just north
       ! of the XZ-planes, etc.
       ! 
       if (apply_budget_global_flux) then
          jflux1(n)=0 ; jflux2(n)=0
          iflux1(n)=0 ; iflux2(n)=0
          do j=1,jm(n)
             ! rlat = ybeg(n) + (float(j)-0.5)*dyr
             ! Example Europe: lat_in = 34, lat_out = 62, lon_in = -12, lon_out  = 36
             if(  nint(ybeg(n) + float(j-1)*dyr) == lat_in) jflux1(n) = j
             if(  nint(ybeg(n) + float(j-1)*dyr) == lat_out) jflux2(n) = j
             !if(  nint(ybeg(n) + float(j)*dyr) == 62)   jflux2(n) = j !CMK corrected
          enddo
          do i=1,im(n)
             ! rlon = xbeg(n) + (float(i)-0.5)*dxr
             if(  nint(xbeg(n) + float(i-1)*dxr) == lon_in) iflux1(n) = i
             if(  nint(xbeg(n) + float(i-1)*dxr) == lon_out)  iflux2(n) = i   ! corrected CMK
             !cmk if(  nint(xbeg(n) + float(i)*dxr) == 36)  iflux2(n) = i   ! corrected CMK
          enddo
       endif

       
       ! --- Write the budgets regional definition to HDF file:
       if ( isRoot ) then
          allocate( gbl_data(im(n),jm(n)))
       else
          allocate( gbl_data(1,1))
       end if
       call gather( dgrid(n), budg_dat(n)%nzong, gbl_data, 0, status)

#ifdef with_hdf4
       if ( isRoot ) then
          
          write (fname,'(a,"/regionsg_",a,".hdf")') trim(outdir), trim(region_name(n))
          call Init(io,fname,'create',status)
          IF_NOTOK_RETURN(status=1)

          call WriteAttribute(io, 'im', im(n), status)
          call WriteAttribute(io, 'jm', jm(n), status)
          call WriteAttribute(io, 'lm', lm(n), status)

          call Init(Sds, io, 'region_global', (/im(n),jm(n)/), 'integer(4)', status)
          IF_NOTOK_RETURN(status=1)

          call WriteData(Sds, gbl_data, status)
          IF_NOTOK_RETURN(status=1)

          call Done(Sds,status)
          IF_NOTOK_RETURN(status=1)

          call Done(io,status)
          IF_NOTOK_RETURN(status=1)
       end if
#endif
       
       deallocate(gbl_data)

    end do REG

    
    ! Filename for budget data (include time range: budget_2001010100_2001020100_global.hdf)
    write (budget_file_global,'(a,"/budget_",i4.4,3i2.2,"_",i4.4,3i2.2,"_",a,".hdf")') &
         trim(outdir), idatei(1:4), idatee(1:4), 'global'

#ifdef with_hdf4          
    if ( isRoot ) then
       !write(gol,'(A  )') '----------------------------------------------------'; call goPr
       write(gol,'(A,A)') 'Budget output file:', trim(budget_file_global); call goPr
       write(gol,'(A  )') '----------------------------------------------------'; call goPr
    endif
#else
    write(gol,'(A)') 'Budget will **NOT** be written - Requires HDF4'; call goPr
    write(gol,'(A)') '----------------------------------------------'; call goPr
#endif    
    ! OK
    status = 0

  END SUBROUTINE INI_ZONEG
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    INIT_BUDGET_GLOBAL
  !
  ! !DESCRIPTION: initialize budgets
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE INIT_BUDGET_GLOBAL ( status )
    !
    ! !USES:
    !
    use ParTools,    only : isRoot
    use global_data, only : region_dat, mass_dat
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)  :: status
    !
    ! !REVISION HISTORY: 
    !   28 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real,dimension(:,:,:,:),pointer    :: rm
    integer,dimension(:,:) ,pointer    :: zoomed, edge
    integer                            :: nzone,  nzone_v, region
    integer                            :: i,j,l,n, i1,i2,j1,j2
    
    character(len=*), parameter ::  rname = mname//'/init_budget_global'

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

    !  first set up the budget regions, and extra info from rcfile
    call ini_zoneg( status )
    IF_NOTOK_RETURN(status=1)

    ! reset budget variables
    budconcg = 0.
    budinig  = 0.
    budconvg = 0.0
    budadvxg = 0.0
    budadvyg = 0.0
    budadvzg = 0.0

    REG: do region = 1, nregions

       rm     => mass_dat(region)%rm
       zoomed => region_dat(region)%zoomed
       edge   => region_dat(region)%edge

       call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )


       ! Tracer #1 init budgets
       call REDUCE( dgrid(region), rm(:,:,:,1), mass_dat(region)%halo, 'SUM', init_mass(region), status)
       IF_NOTOK_RETURN(status=1)
       
       sum_update(region)       = 0.0
       sum_advection(region)    = 0.0
       !
       ! initialise budinig
       !
       do n=1,ntracet
       do l=1,lm(region)
             
          nzone_v=nzon_vg(l)  ! vertical bud code
          
          do j=j1,j2
          do i=i1,i2

             !skip zoomed part of the region
             if ( zoomed(i,j) /= region .or. edge(i,j) >= 1 ) cycle
                
             ! include the edges....parent takes care!
             nzone = budg_dat(region)%nzong(i,j)
                
             budinig(nzone, nzone_v, n) = budinig(nzone, nzone_v, n) + &
                                          rm(i,j,l,n)/ra(n)*1000.
          end do !i
          end do !j
       end do  !l
       end do  !n 

       nullify(rm)
       nullify(zoomed)
       nullify(edge)
       
    enddo REG

    ! ok
    status = 0

  END SUBROUTINE INIT_BUDGET_GLOBAL
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DONE_BUDGET_GLOBAL
  !
  ! !DESCRIPTION: compute concentration change b/w beginning and end of the
  !               run; write budgets to output file.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DONE_BUDGET_GLOBAL ( status )
    !
    ! !USES:
    !
    use ParTools,     only : isRoot, par_reduce_element, par_barrier
    use global_data,  only : region_dat, mass_dat
#ifdef with_hdf4          
    use file_hdf,     only : THdfFile, TSds
    use file_hdf,     only : Done, WriteAttribute, Init, WriteData, SetDim
#endif
    use tm5_distgrid, only : gather_i_band, gather_j_band
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !   28 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real,    dimension(:,:,:,:), pointer :: rm
    integer, dimension(:,:) ,    pointer :: zoomed, edge
    
    integer        :: nzone, nzone_v, region
    integer        :: i,j,l,n, i1,i2,j1,j2
#ifdef with_hdf4          
    type(THdfFile) :: io
    type(TSds)     :: sds
#endif
    ! Need these arrays to deal with buggy MPI implementations: if send and receive buffers
    ! are the same, you can get this error on some system:
    ! Fatal error in PMPI_Reduce: Invalid buffer pointer, error stack:
    ! PMPI_Reduce(1701): MPI_Reduce(sbuf=0x1bae7a0 rbuf=0x1bae7a0, count=3060, MPI_DOUBLE_PRECISION, MPI_SUM, root=0, MPI_COMM_WORLD) failed
    ! PMPI_Reduce(1511): Buffers must not be aliased
    ! For MPICH2, this is a bug that was fixed in 2009 (http://lists.mcs.anl.gov/pipermail/mpich-discuss/2009-October/005804.html)
    ! but this is not the case for Intel MPI (the version used at KNMI) 
    real, dimension(nbudg,nbud_vg,ntracet) :: budconcg_all
    real, dimension(nbudg,nbud_vg,ntracet) :: budchngg_all
    real, dimension(nbudg,nbud_vg,ntracet) :: budinig_all
    real, dimension(nbudg,nbud_vg,ntracet) :: budadvxg_all 
    real, dimension(nbudg,nbud_vg,ntracet) :: budadvyg_all 
    real, dimension(nbudg,nbud_vg,ntracet) :: budadvzg_all 
    real, dimension(nbudg,nbud_vg,ntracet) :: budconvg_all 

    ! --- const -------------------------------
    
    character(len=*), parameter ::  rname = mname//'/done_budget_global'

    integer,          parameter :: nproces=4   ! x,y,z, conv

    character(len=10), dimension(nproces),parameter :: proces_name=(/'names     ','names     ','names     ','names     '/)
    integer,           dimension(nproces),parameter :: proces_nr = (/ ntracet, ntracet, ntracet, ntracet /)


    ! --- begin --------------------------------
    if (.not.apply_budget_global) return


    ! final concentrations
    budconcg = 0.0

    do region = 1, nregions

       call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

       zoomed => region_dat(region)%zoomed
       edge   => region_dat(region)%edge
       rm     =>   mass_dat(region)%rm

       do n=1,ntracet
       do l=1,lm(region)

          nzone_v = nzon_vg(l)
             
          do j=j1,j2
          do i=i1,i2

             if (zoomed(i,j)/=region .or. edge(i,j) >= 1) cycle  
             nzone = budg_dat(region)%nzong(i,j)
             
             budconcg(nzone, nzone_v, n) = budconcg(nzone, nzone_v, n) + &
                                           rm(i,j,l,n)/ra(n)*1000.
          end do !i
          end do !j
       end do !l
       end do !n

       nullify(rm)
       nullify(zoomed)
       nullify(edge)

    enddo


    ! Calculate concentrations change
    budchngg = budconcg - budinig

    
    ! Sum up contribution from each proc into root array    
    call PAR_REDUCE_ELEMENT( budconcg, 'SUM', budconcg_all, status )
    IF_NOTOK_RETURN(status=1)

    call PAR_REDUCE_ELEMENT( budchngg, 'SUM', budchngg_all, status )
    IF_NOTOK_RETURN(status=1)

    call PAR_REDUCE_ELEMENT( budinig,  'SUM', budinig_all, status )
    IF_NOTOK_RETURN(status=1)

    call PAR_REDUCE_ELEMENT( budadvxg, 'SUM', budadvxg_all, status )
    IF_NOTOK_RETURN(status=1)                 

    call PAR_REDUCE_ELEMENT( budadvyg, 'SUM', budadvyg_all, status )
    IF_NOTOK_RETURN(status=1)                 

    call PAR_REDUCE_ELEMENT( budadvzg, 'SUM', budadvzg_all, status )
    IF_NOTOK_RETURN(status=1)

    call PAR_REDUCE_ELEMENT( budconvg, 'SUM', budconvg_all, status )
    IF_NOTOK_RETURN(status=1)

#ifdef with_hdf4          
    
    ! Write out budget file
    if ( isRoot ) then

       call Init(io, budget_file_global, 'create', status)
       IF_NOTOK_RETURN(status=1)
       
       call WriteAttribute(io, 'nregions',    nregions,    status)
       ! lib cannot handle kind=8
       ! call WriteAttribute(io, 'itau',        itau,        status)
       call WriteAttribute(io, 'nbud',        nbudg,       status)
       call WriteAttribute(io, 'nbud_v',      nbud_vg,     status)
       call WriteAttribute(io, 'ntrace',      ntrace,      status)
       call WriteAttribute(io, 'ntracet',     ntracet,     status)
       call WriteAttribute(io, 'idate',       idate,       status)
       call WriteAttribute(io, 'idatei',      idatei,      status)
       call WriteAttribute(io, 'idatee',      idatee,      status)
       call WriteAttribute(io, 'proces_nr',   proces_nr,   status)
       call WriteAttribute(io, 'proces_name', proces_name, status)
       call WriteAttribute(io, 'zone_name',   zone_nameg,  status)
       call WriteAttribute(io, 'zone_namv',   zone_namvg,  status)
       call WriteAttribute(io, 'names',       names,       status)
       if (apply_budget_global_flux) then
          call WriteAttribute(io, 'iflux1' , iflux1, status)
          call WriteAttribute(io, 'iflux2' , iflux2, status)
          call WriteAttribute(io, 'jflux1' , jflux1, status)
          call WriteAttribute(io, 'jflux2' , jflux2, status)
          call WriteAttribute(io, 'lflux1' , lflux1, status)
          call WriteAttribute(io, 'lflux2' , lflux2, status)
       endif
       IF_NOTOK_RETURN(status=1)
       
       call Init(Sds, io, 'budconc', (/nbudg, nbud_vg, ntracet/), 'real(8)',      status)
       call SetDim(Sds, 0, 'nbudg', 'horizontal budget', (/(j, j=1,nbudg)/),   status)
       call SetDim(Sds, 1, 'nbud_vg', 'vertical budget', (/(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, budconcg_all, status)       
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)

       call Init(Sds, io, 'budchng', (/nbudg, nbud_vg, ntracet/), 'real(8)', status)
       call SetDim(Sds, 0, 'nbudg', 'horizontal budget', (/(j, j=1,nbudg)/), status)
       call SetDim(Sds, 1, 'nbud_vg', 'vertical budget', (/(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, budchngg_all, status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)

       call Init(Sds, io, 'budini', (/nbudg, nbud_vg, ntracet/), 'real(8)', status)
       call SetDim(Sds, 0, 'nbudg', 'horizontal budget', (/(j, j=1,nbudg)/), status)
       call SetDim(Sds, 1, 'nbud_vg', 'vertical budget', (/(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, budinig_all, status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)

       call Init(Sds, io, 'budadvx', (/nbudg, nbud_vg, ntracet/), 'real(8)', status)
       call SetDim(Sds, 0, 'nbudg', 'horizontal budget', (/(j, j=1,nbudg)/), status)
       call SetDim(Sds, 1, 'nbud_vg', 'vertical budget', (/(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, budadvxg_all, status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)

       call Init(Sds, io, 'budadvy', (/nbudg, nbud_vg, ntracet/), 'real(8)', status)
       call SetDim(Sds, 0, 'nbudg', 'horizontal budget', (/(j, j=1,nbudg)/), status)
       call SetDim(Sds, 1, 'nbud_vg', 'vertical budget', (/(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, budadvyg_all, status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)

       call Init(Sds, io, 'budadvz', (/nbudg, nbud_vg, ntracet/), 'real(8)', status)
       call SetDim(Sds, 0, 'nbudg', 'horizontal budget', (/(j, j=1,nbudg)/), status)
       call SetDim(Sds, 1, 'nbud_vg', 'vertical budget', (/(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, budadvzg_all, status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)

       call Init(Sds, io, 'budconv_nocp', (/nbudg, nbud_vg, ntracet/), 'real(8)', status)
       call SetDim(Sds, 0, 'nbudg', 'horizontal budget', (/(j, j=1,nbudg)/), status)
       call SetDim(Sds, 1, 'nbud_vg', 'vertical budget', (/(j, j=1,nbud_vg)/), status)
       call SetDim(Sds, 2, 'ntracet', 'tracer number', (/(j, j=1,ntracet)/), status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(Sds, 'comment', 'Convection budget not corrected for cp removal', status)
       IF_NOTOK_RETURN(status=1)
       call WriteData(Sds, budconvg_all, status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)

    end if ! root
#endif
    ! PLS: old unclear comment:
       ! Jottum. We have the chemistry. We have it before dry deposition.
       ! Let us not write it down. That is an idea
       ! Chemistry is done by budget_global
       ! Dry deposition is done by chemistry
       ! That is logical. We do not write chemistry here, because 'chemistry', the one
       ! about dry deposition can access it. And we need dry deposition (that is not accessible).
       ! because chemistry has been infested by dry deposition.

    ! Gather and write the extra flux through 2D slices
    if (apply_budget_global_flux) then

       do region=1,nregions

          call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

          
          call GATHER_I_BAND( dgrid(region), budget_flux(region)%flux_x1, budget_flux_all(region)%flux_x1, status, iflux1(region) )
          IF_NOTOK_RETURN(status=1)

#ifdef with_hdf4          
          if (isRoot) then
             call Init(Sds, io, 'budget_flux_x1_'//region_name(region), (/jm(region),lm(region),ntracet/), 'real(4)', status)
             call SetDim(Sds, 0, 'jm_'//region_name(region), 'latitude index', (/(j, j=1,jm(region))/), status)
             call SetDim(Sds, 1, 'lm_'//region_name(region), 'layer index   ', (/(j, j=1,lm(region))/), status)
             call SetDim(Sds, 2, 'ntracet', 'tracer number', (/(j, j=1,ntracet)/), status)
             IF_NOTOK_RETURN(status=1)
             call WriteData(Sds, budget_flux_all(region)%flux_x1, status)
             IF_NOTOK_RETURN(status=1)
             call Done(Sds, status)
             IF_NOTOK_RETURN(status=1)
          end if
#endif
          call GATHER_I_BAND( dgrid(region), budget_flux(region)%flux_x2, budget_flux_all(region)%flux_x1, status, iflux2(region) )
          IF_NOTOK_RETURN(status=1)

#ifdef with_hdf4          
          if (isRoot) then
             call Init(Sds, io, 'budget_flux_x2_'//region_name(region), (/jm(region),lm(region),ntracet/), 'real(4)', status)
             call SetDim(Sds, 0, 'jm_'//region_name(region), 'latitude index', (/(j, j=1,jm(region))/), status)
             call SetDim(Sds, 1, 'lm_'//region_name(region), 'layer index   ', (/(j, j=1,lm(region))/), status)
             call SetDim(Sds, 2, 'ntracet', 'tracer number', (/(j, j=1,ntracet)/), status)
             IF_NOTOK_RETURN(status=1)
             call WriteData(Sds, budget_flux_all(region)%flux_x1, status)
             IF_NOTOK_RETURN(status=1)
             call Done(Sds, status)
             IF_NOTOK_RETURN(status=1)
          end if
#endif
          call GATHER_J_BAND( dgrid(region), budget_flux(region)%flux_y1, budget_flux_all(region)%flux_y1, status, jflux1(region) )
          IF_NOTOK_RETURN(status=1)

#ifdef with_hdf4          
          if (isRoot) then
             call Init(Sds, io, 'budget_flux_y1_'//region_name(region), (/im(region),lm(region),ntracet/), 'real(4)', status)
             call SetDim(Sds, 0, 'im_'//region_name(region), 'longitude index', (/(j, j=1,im(region))/), status)
             call SetDim(Sds, 1, 'lm_'//region_name(region), 'layer index    ', (/(j, j=1,lm(region))/), status)
             call SetDim(Sds, 2, 'ntracet', 'tracer number', (/(j, j=1,ntracet)/), status)
             IF_NOTOK_RETURN(status=1)
             call WriteData(Sds, budget_flux_all(region)%flux_y1, status)
             IF_NOTOK_RETURN(status=1)
             call Done(Sds, status)
             IF_NOTOK_RETURN(status=1)
          end if
#endif
          call GATHER_J_BAND( dgrid(region), budget_flux(region)%flux_y2, budget_flux_all(region)%flux_y1, status, jflux2(region) )
          IF_NOTOK_RETURN(status=1)

#ifdef with_hdf4          
          if (isRoot) then
             call Init(Sds, io, 'budget_flux_y2_'//region_name(region), (/im(region),lm(region),ntracet/), 'real(4)', status)
             call SetDim(Sds, 0, 'im_'//region_name(region), 'longitude index', (/(j, j=1,im(region))/), status)
             call SetDim(Sds, 1, 'lm_'//region_name(region), 'layer index    ', (/(j, j=1,lm(region))/), status)
             call SetDim(Sds, 2, 'ntracet', 'tracer number', (/(j, j=1,ntracet)/), status)
             IF_NOTOK_RETURN(status=1)
             call WriteData(Sds, budget_flux_all(region)%flux_y1, status)
             IF_NOTOK_RETURN(status=1)
             call Done(Sds, status)
             IF_NOTOK_RETURN(status=1)
          end if
#endif
          call GATHER( dgrid(region), budget_flux(region)%flux_z1, budget_flux_all(region)%flux_z1, 0, status )
          IF_NOTOK_RETURN(status=1)

#ifdef with_hdf4          
          if (isRoot) then
             call Init(Sds, io, 'budget_flux_z1_'//region_name(region), (/im(region),jm(region),ntracet/), 'real(4)', status)
             call SetDim(Sds, 0, 'im_'//region_name(region), 'longitude index', (/(j, j=1,im(region))/), status)
             call SetDim(Sds, 1, 'jm_'//region_name(region), 'latitude index ', (/(j, j=1,jm(region))/), status)
             call SetDim(Sds, 2, 'ntracet', 'tracer number', (/(j, j=1,ntracet)/), status)
             IF_NOTOK_RETURN(status=1)
             call WriteData(Sds, budget_flux_all(region)%flux_z1, status)
             IF_NOTOK_RETURN(status=1)
             call Done(Sds, status)
             IF_NOTOK_RETURN(status=1)
          end if
#endif
          call GATHER( dgrid(region), budget_flux(region)%flux_z2, budget_flux_all(region)%flux_z1, 0, status )
          IF_NOTOK_RETURN(status=1)

#ifdef with_hdf4          
          if (isRoot) then
             call Init(Sds, io, 'budget_flux_z2_'//region_name(region), (/im(region),jm(region),ntracet/), 'real(4)', status)
             call SetDim(Sds, 0, 'im_'//region_name(region), 'longitude index', (/(j, j=1,im(region))/), status)
             call SetDim(Sds, 1, 'jm_'//region_name(region), 'latitude index ', (/(j, j=1,jm(region))/), status)
             call SetDim(Sds, 2, 'ntracet', 'tracer number', (/(j, j=1,ntracet)/), status)
             IF_NOTOK_RETURN(status=1)
             call WriteData(Sds, budget_flux_all(region)%flux_z1, status)
             IF_NOTOK_RETURN(status=1)
             call Done(Sds, status)
             IF_NOTOK_RETURN(status=1)
          end if
#endif
          deallocate( budget_flux(region)%flux_x1)
          deallocate( budget_flux(region)%flux_y1)
          deallocate( budget_flux(region)%flux_z1)
          deallocate( budget_flux(region)%flux_x2)
          deallocate( budget_flux(region)%flux_y2)
          deallocate( budget_flux(region)%flux_z2)
          deallocate( budget_flux_all(region)%flux_x1)
          deallocate( budget_flux_all(region)%flux_y1)
          deallocate( budget_flux_all(region)%flux_z1)
          !KLUDGE (variable used in chemistr_done!) deallocate( budg_dat(region)%nzong)
          
       end do
    endif   ! apply_budget_flux

#ifdef with_hdf4          
    if (isRoot) then 
       call Done(io, status)
       IF_NOTOK_RETURN(status=1)
    end if
#endif
    call budget_displayg ( status )
    IF_NOTOK_RETURN(status=1)

    ! OK
    status = 0

  END SUBROUTINE DONE_BUDGET_GLOBAL
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    BUDGET_TRANSPORTG
  !
  ! !DESCRIPTION: evaluate advection/convection budgets
  !
  ! Routine is called before and after a transport step. The difference is used to
  ! calculate the budget. Problems occur with zoome regions. Here advection is
  ! applied by writing x-adges and y-edges (in put_xedges, put_yedges).
  ! Knowledge of the previous step is required to calculate the budget....
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE BUDGET_TRANSPORTG(region, stat, proces, prev_step)
    !
    ! !USES:
    !
    use dims,        only : okdebug
    use global_data, only : region_dat, mass_dat
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)   :: stat      ! 0 before process, 1 after
    integer, intent(in)   :: region
    character(len=*)      :: proces    ! current process (xyzvcs)
    character(len=*)      :: prev_step ! previous process.
    !
    ! !REVISION HISTORY: 
    !   28 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real,    dimension(:,:,:,:), pointer :: rm
    integer, dimension(:,:),     pointer :: zoomed, edge
    logical, dimension(:,:), allocatable :: skip_flag

    integer :: i,j,l,n, nzone, nzone_v, lmr, i1,j1,i2,j2

    ! ------ start
    
    if (.not.apply_budget_global) return

    rm     =>   mass_dat(region)%rm    
    zoomed => region_dat(region)%zoomed
    edge   => region_dat(region)%edge

    lmr = lm(region)

    if ( okdebug ) then
       write(gol,*) 'budget_transportg: region, stat, proces, prev_step = ',&
            region, stat, proces,' ', prev_step ; call goPr
    end if

    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )


    !---------- flag

    allocate(skip_flag(i1:i2,j1:j2))

    skip_flag = .false.

    where ( zoomed /= region .or. edge >= 1 )
       skip_flag = .true.
    end where

    if ( proces(1:4) == 'advx' .and. prev_step /= 'y' ) then
       ! include y-edges if XYZ sequence
       where(edge >= 2) 
          skip_flag = .false.    !only remaining problem: reduced grid!   
       end where
    else if ( proces(1:4) == 'advy' .and. prev_step /= 'x' ) then
       where ( edge == 1 .or. edge == 3 ) 
          skip_flag = .false.
       end where
    else if ( proces(1:4) == 'advz' .and. prev_step /= 'y' ) then
       ! zyx sequence: include interface
       where ( edge >= 1 ) 
          skip_flag = .false.
       end where
    else if ( proces(1:4) == 'conv' ) then   ! always interface
       where ( edge >= 1 ) 
          skip_flag = .false.
       end where
    end if

    !----------- Calculate budget
    
    if ( stat == 0 ) then ! save current tracer concentrations

       rmoldg = 0.0
       do n=1, ntracet
       do l=1,lmr

          nzone_v = nzon_vg(l)
          do j=j1, j2
          do i=i1, i2
             if (skip_flag(i,j)) cycle
             nzone = budg_dat(region)%nzong(i,j)
             rmoldg(nzone,nzone_v,n) = &
                  rmoldg(nzone,nzone_v,n) + rm(i,j,l,n)
          end do
          end do

       end do
       end do

    else ! new concentration and fill budget

       rmnewg = 0.0
       do n=1, ntracet
       do l=1,lmr

          nzone_v = nzon_vg(l)
          do j=j1,j2
          do i=i1,i2
             if (skip_flag(i,j)) cycle
             nzone=budg_dat(region)%nzong(i,j)
             rmnewg(nzone,nzone_v,n) = rmnewg(nzone,nzone_v,n) + rm(i,j,l,n)
          end do
          end do

       end do
       end do

       select CASE(proces(1:4))
       case('advx')
          do n=1, ntracet
             budadvxg(:,:,n) = budadvxg(:,:,n) + &
                  (rmnewg(:,:,n)-rmoldg(:,:,n))/ra(n)*1e3
          end do
       case('advy')
          do n=1, ntracet
             budadvyg(:,:,n) = budadvyg(:,:,n) + &
                  (rmnewg(:,:,n)-rmoldg(:,:,n))/ra(n)*1e3
          end do
       case('advz')
          do n=1, ntracet
             budadvzg(:,:,n) = budadvzg(:,:,n) + &
                  (rmnewg(:,:,n)-rmoldg(:,:,n))/ra(n)*1e3
          end do
       case('conv')
          do n=1, ntracet
             budconvg(:,:,n) = budconvg(:,:,n) + &
                  (rmnewg(:,:,n)-rmoldg(:,:,n))/ra(n)*1e3
          end do
       case default
          write(gol,*)'budget_transportg: no valid value for transport budget,', &
               ' budget not updated!'; call goPr
       end select
    end if

    deallocate(skip_flag)
    nullify(rm)
    nullify(zoomed)
    nullify(edge)

  END SUBROUTINE BUDGET_TRANSPORTG
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    BUDGET_DISPLAYG
  !
  ! !DESCRIPTION: print out budget for tracer #1
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE BUDGET_DISPLAYG( status )
    !
    ! !USES:
    !
    use dims,        only : nregions
    use chem_param,  only : names
    use global_data, only : mass_dat
    use partools,    only : isRoot, par_reduce
#ifdef with_hdf4          
    use file_hdf,    only : THdfFile
    use file_hdf,    only : Init, Done, WriteAttribute, WriteData
#endif
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)      ::  status
    !
    ! !REVISION HISTORY: 
    !   28 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

#ifdef with_hdf4          
    type(THdfFile)              :: io
#endif
    character(len=40)           :: frmt
    integer                     :: region, imr, jmr, lmr
    real                        :: mass(nregions)
    character(len=*), parameter :: rname = mname//'/budget_displayg'

    real, dimension(nregions) :: sum_update_all, sum_advection_all
    
    ! --- begin --------------------------------

    if (.not.apply_budget_global) return

    if ( isRoot ) then 
       write (gol,'("--------------------------------------------------------")'); call goPr
       write (gol,'(" Budget of tracer ",a," (kg) :     ")') trim(names(1))      ; call goPr
       write (gol,'("--------------------------------------------------------")'); call goPr
    end if

    do region = 1, nregions 

       call REDUCE( dgrid(region), mass_dat(region)%rm(:,:,:,1), mass_dat(region)%halo, 'SUM', mass(region), status)
       IF_NOTOK_RETURN(status=1)

       CALL PAR_REDUCE(sum_update(region), 'SUM', sum_update_all(region), status)
       IF_NOTOK_RETURN(status=1)
       
       CALL PAR_REDUCE(sum_advection(region), 'SUM', sum_advection_all(region), status)
       IF_NOTOK_RETURN(status=1)

       
       call get_format(mass(region),frmt)

       if ( isRoot ) then 
          !write (gol,'("Region             ",i3)') region; call goPr
          write (gol,frmt) 'Initial mass   :            ',init_mass(region); call goPr
          !          write (gol,frmt) 'mass emitted   :            ',sum_emission(region); call goPr
          !          write (gol,frmt) 'chemistry      :            ',sum_chemistry_all; call goPr
          !          write (gol,frmt) 'wet chemistry  :            ',-sum_wetchem_all; call goPr
          write (gol,frmt) 'update_p       :            ',sum_update_all(region); call goPr
          write (gol,frmt) 'advection      :            ',sum_advection_all(region); call goPr
          !          write (gol,frmt) 'deposition     :            ',sum_deposition(region); call goPr
          !          write (gol,frmt) 'stratosphere   :            ',sum_stratosphere(region); call goPr
          write (gol,frmt) 'Final mass     :            ',mass(region); call goPr
          write (gol,frmt) 'Budget (error) :            ', &
               init_mass(region) +sum_advection_all(region) +sum_update_all(region)-mass(region) ; call goPr
          !          write (gol,frmt) 'Budget (error) :            ', &
          !               init_mass(region)  + sum_emission(region) &
          !               +sum_advection(region) + sum_chemistry_all &
          !               +sum_update(region)  &
          !               +sum_deposition(region) + sum_stratosphere(region) - sum_wetchem_all - mass; call goPr
          write (gol,'("--------------------------------------------------------")'); call goPr
       end if

    end do

#ifdef with_hdf4          
    if ( isRoot ) then 

       call Init(io, budget_file_global, 'write', status)
       IF_NOTOK_RETURN(status=1)
       
       call WriteAttribute(io, 'initial mass',      init_mass, status)
       IF_NOTOK_RETURN(status=1)

       call WriteAttribute(io, 'final mass',             mass, status)
       IF_NOTOK_RETURN(status=1)

       call WriteAttribute(io, 'sum_advection', sum_advection_all, status)
       IF_NOTOK_RETURN(status=1)

       call WriteAttribute(io, 'sum_update',       sum_update_all, status)
       IF_NOTOK_RETURN(status=1)

       call Done(io,status)
       IF_NOTOK_RETURN(status=1)
    endif
#endif
    ! ok
    status = 0

  END SUBROUTINE BUDGET_DISPLAYG
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    GET_FORMAT
  !
  ! !DESCRIPTION: returns convenient format to write a real 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE GET_FORMAT( mm, frmt)
    !
    ! !INPUT PARAMETERS:
    !
    real, intent(in)              :: mm
    !
    ! !OUTPUT PARAMETERS:
    !
    character(len=*), intent(out) :: frmt
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer :: is

    if ( mm > 0.0 ) then
       is = int(alog10(mm))
    else
       is = 0
    end if

    select case(is)
    case(1:10)
       write (frmt,'(a9,i1,a1)') '(a30,f14.', 10-is, ')'
    case default
       frmt =  '(a30,1pe12.5)'
    end select

  END SUBROUTINE GET_FORMAT
  !EOC

END MODULE BUDGET_GLOBAL