ecearth3/sources/tm5mp/proj/budget10/budget_global__detailed.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:      BUDGET_GLOBAL
!
! !DESCRIPTION: This module holds most of the budget variables (see
!               wet_deposition for exceptions). 
!
! ******************************************************************
! ***** THIS VERSION DIFFERS FROM THE ONE IN BASE IN THE        ****
! ***** DEFINITION OF THE HORIZONTAL AND VERTICAL BUDGET-REGION ****
! ******************************************************************
!
! 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
!            HERE: 10degree wide latitudinal bands
!    o  defines the vertical split of the budgets
!            HERE: each levels (instead of vertical 4 regions)
!    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   = 18        ! # horizontal budget regions (h-budgets)
  integer, public, parameter    :: nbud_vg = lm(1)     ! # 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)

  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
  !   16 Jul 2013 - P. Le Sager - version for "budget_detailed"
  !
  ! !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
    use file_hdf,    only : THdfFile, TSds
    use file_hdf,    only : Init, Done, WriteAttribute, WriteData
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)   :: status
    !
    ! !REVISION HISTORY: 
    !   28 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    type(THdfFile)       :: io
    type(TSds)           :: sds
    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 

    do i = 1, nbud_vg
       write( zone_namvg(i),'("Level_",i2.2)')i
    enddo

    do i=1,nbudg
       write( zone_nameg(i),'("i_Latitude_",i2.2)')i
    enddo
    
    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)

    do i=1,lm(1)
       nzon_vg(i) = i
    enddo
    

    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
          budg_dat(n)%nzong(:,j) = ( nint(lats)+90+10 ) / 10 
       end do


       ! --- 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)

       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

       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'

    if ( isRoot ) then
       !write(gol,'(A  )') '----------------------------------------------------'; call goPr
       write(gol,'(A,A)') 'Budget will be written to:', trim(budget_file_global); 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
    use file_hdf,     only : THdfFile, TSds
    use file_hdf,     only : Done, WriteAttribute, Init, WriteData, SetDim
    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
    type(THdfFile) :: io
    type(TSds)     :: sds

    ! 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)


    ! 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)
       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

    ! 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)

          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

          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)

          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

          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)

          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

          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)

          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

          call GATHER( dgrid(region), budget_flux(region)%flux_z1, budget_flux_all(region)%flux_z1, 0, status )
          IF_NOTOK_RETURN(status=1)

          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

          call GATHER( dgrid(region), budget_flux(region)%flux_z2, budget_flux_all(region)%flux_z1, 0, status )
          IF_NOTOK_RETURN(status=1)

          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

          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)
          deallocate( budg_dat(region)%nzong)

       end do
    endif   ! apply_budget_flux

    if (isRoot) then 
       call Done(io, status)
       IF_NOTOK_RETURN(status=1)
    end if

    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. 
  !\\
  !\\
  ! !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
    use file_hdf,    only : THdfFile
    use file_hdf,    only : Init, Done, WriteAttribute, WriteData
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)      ::  status
    !
    ! !REVISION HISTORY: 
    !   28 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    ! --- 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 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(region); call goPr
          write (gol,frmt) 'advection      :            ',sum_advection(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(region) +sum_update(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

    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, status)
       IF_NOTOK_RETURN(status=1)

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

       call Done(io,status)
       IF_NOTOK_RETURN(status=1)
    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