! #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). ! ! ****************************************************************** ! ***** 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 #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 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) #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 will be written to:', 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. !\\ !\\ ! !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