!!!!!!!!!!!!!!!!!!!!!?????????????????????????
!!!!!!!!!!!!!!!!!!! one line change here just  ????????
!!!!!!!!!!!!!!!!!for using subroutine outgrads1 ????
!!!!!!!!!!!!!!!!!!!!!?????????????????????????
subroutine antest_maps0(mype,theta0f,z0f)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    antest_maps0
!   prgmmr:
!
! abstract: this routine creates output maps of background error correlations
!
! program history log:
!   2009-09-18  lueken - added subprogram doc block
!   2010-03-30  zhu    - use nvars from control_vectors
!   2010-12-05  pondeca - change variable names: "tv" is now "t" , and "st" is "sf" 
!   2013-10-19  todling - metguess now holds background
!   2016-04-19  pondeca - updated to use bundle_work as argument to ansmoothrf
!
!   input argument list:
!    mype
!    theta0f
!    z0f
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind,r_single
  use anberror, only: kvar_start,kvar_end,var_names,pf2aP1,indices
  use gridmod, only: nsig,nsig1o,nlon,nlat,istart,jstart,lat2,lon2
  use constants, only: zero_single,zero,one,rd_over_cp,r100
  use mpimod, only: ierror,mpi_real4,mpi_real8,mpi_sum,mpi_comm_world
  use guess_grids, only: ntguessig,ges_prsl
  use fgrid2agrid_mod, only: fgrid2agrid
  use control_vectors, only: nvars
  use gsi_bundlemod, only: gsi_bundlecreate
  use gsi_bundlemod, only: gsi_bundle
  use gsi_bundlemod, only: gsi_bundlegetvar
  use gsi_bundlemod, only: gsi_bundleputvar
  use gsi_bundlemod, only: gsi_bundledestroy
  use gsi_bundlemod, only: gsi_grid
  use gsi_bundlemod, only: gsi_gridcreate
  use gsi_bundlemod, only: gsi_bundlegetpointer
  use gsi_metguess_mod, only: gsi_metguess_bundle
  use mpeu_util, only: die
  implicit none

  integer(i_kind),intent(in   ) :: mype
  type(gsi_bundle):: bundle_work
  type(gsi_grid) :: grid
  character(2) :: names2dwork(1),names3dwork(4)

  real(r_single) ,intent(in   ) :: theta0f(pf2aP1%nlatf,pf2aP1%nlonf,nsig1o)
  real(r_single) ,intent(in   ) :: z0f(pf2aP1%nlatf,pf2aP1%nlonf,nsig1o)

  character(len=*),parameter::myname='antest_maps0'
  real(r_kind),dimension(lat2,lon2,nsig):: twork,qwork,stwork,vpwork
  real(r_kind),dimension(lat2,lon2):: pwork


  real(r_kind) tempf(nlat,nlon),tempc(pf2aP1%nlatf,pf2aP1%nlonf)
  real(r_single) outwork(nlon,nlat),outwork0(nlon,nlat)

  character(80) ref_plotcor
  character(80) var_plotcor
  integer(i_kind) i_plotcor,j_plotcor,k_plotcor

  real(r_kind),dimension(:,:  ),pointer::ges_z_it =>NULL()
  real(r_kind),dimension(:,:,:),pointer::ges_tv_it=>NULL()

  real(r_kind)h00,h000
  integer(i_kind) lunin,i,j,k,ivar,iglob,jglob,ivar_plot
  integer(i_kind) it,mm1,ier,istatus
  integer(i_kind) iloc,jloc,kloc

  integer(i_kind):: iref

  integer(i_kind):: nsig_aux

  names2dwork(1)='ps'
  names3dwork(1)='sf'
  names3dwork(2)='vp'
  names3dwork(3)='t'
  names3dwork(4)='q'

!*********************************************************************
!          variable names expected for var_plotcor are
!
!    sf  -- stream function
!    vp  -- velocity potential
!    ps  -- surface pressure
!    t   -- virtual temperature
!    q   -- specific humidity
!    note: add more as desired
!*********************************************************************
! Make choices here!
  i_plotcor=nlat/2
  j_plotcor=nlon/2
  k_plotcor=min(nsig,25)
  iloc=i_plotcor-istart(mype+1)+2
  jloc=j_plotcor-jstart(mype+1)+2
  kloc=k_plotcor
  var_plotcor='sf'
!Note: Must call this subroutine from anprewgt_reg.f90
!Make sure statement has been uncommented!
! End of choice section
!*********************************************************************

  call gsi_gridcreate(grid,lat2,lon2,nsig)
  ref_plotcor='theta'
  ier=0
  it=ntguessig
  call gsi_bundlegetpointer (gsi_metguess_bundle(it),'z' ,ges_z_it,   istatus)
  ier=ier+istatus
  call gsi_bundlegetpointer (gsi_metguess_bundle(it),'tv',ges_tv_it,  istatus)
  ier=ier+istatus
  if(ier/=0) call die(myname,'missing fields, ier= ', ier)

  do 200 iref=1,5

     if (iref==1) var_plotcor='sf'
     if (iref==2) var_plotcor='vp'
     if (iref==3) var_plotcor='ps'
     if (iref==4) var_plotcor='t'
     if (iref==5) var_plotcor='q'

     ref_plotcor=var_plotcor
     lunin=1
     if(mype==0) then
        open(lunin,file="cormaps_"//trim(var_plotcor),form='unformatted')
        rewind lunin
     end if
     ivar_plot=0
     do ivar=1,nvars
        if(trim(var_names(ivar))==trim(var_plotcor)) then
           ivar_plot=ivar
           exit
        end if
     end do
     if(ivar_plot==0) then
        write(0,*)' in antest_maps0, variable ',trim(var_plotcor),'  not found.  program stops'
        call mpi_finalize(ierror)
        stop
     end if

     twork=zero
     pwork=zero
     qwork=zero
     stwork=zero
     vpwork=zero

     if(ivar_plot==1) stwork(iloc,jloc,kloc)=one
     if(ivar_plot==2) vpwork(iloc,jloc,kloc)=one
     if(ivar_plot==3) pwork(iloc,jloc)=one
     if(ivar_plot==4) twork(iloc,jloc,kloc)=one
     if(ivar_plot==5) qwork(iloc,jloc,kloc)=one

     call gsi_bundlecreate(bundle_work,grid,'bundle work',ier, &
              names2d=names2dwork,names3d=names3dwork,bundle_kind=r_kind)

     call gsi_bundleputvar(bundle_work,'ps',pwork,ier)
     call gsi_bundleputvar(bundle_work,'sf',stwork,ier)
     call gsi_bundleputvar(bundle_work,'vp',vpwork,ier)
     call gsi_bundleputvar(bundle_work,'t' ,twork,ier)
     call gsi_bundleputvar(bundle_work,'q' ,qwork,ier)

     call ansmoothrf(bundle_work)

     call gsi_bundlegetvar(bundle_work,'ps',pwork,ier)
     call gsi_bundlegetvar(bundle_work,'sf',stwork,ier)
     call gsi_bundlegetvar(bundle_work,'vp',vpwork,ier)
     call gsi_bundlegetvar(bundle_work,'t' ,twork,ier)
     call gsi_bundlegetvar(bundle_work,'q' ,qwork,ier)

     call gsi_bundledestroy(bundle_work,ier)

     if(mype==0) write(0,*) ' refvar= ',trim(ref_plotcor),' corvar= ',trim(var_plotcor), &
            '  i,j,k_plotcor =', j_plotcor,i_plotcor,k_plotcor, ' nlon,nlat,nsig=', &
                  nlon,nlat,kvar_end(ivar_plot)-kvar_start(ivar_plot)+1

     if(ivar_plot==1) h00=stwork(iloc,jloc,kloc)
     if(ivar_plot==2) h00=vpwork(iloc,jloc,kloc)
     if(ivar_plot==3) h00=pwork(iloc,jloc)
     if(ivar_plot==4) h00=twork(iloc,jloc,kloc)
     if(ivar_plot==5) h00=qwork(iloc,jloc,kloc)

     call mpi_allreduce(h00,h000,1,mpi_real8,mpi_sum,mpi_comm_world,ierror)

     stwork=stwork/h000
     vpwork=vpwork/h000
     pwork=pwork/h000
     twork=twork/h000
     qwork=qwork/h000

!     output original pot temp  (slow way to get full 2d field)  -- this is reference field

     mm1=mype+1
     do k=1,nsig
        outwork=zero_single
        do j=2,lon2-1
           jglob=jstart(mm1)-2+j
           do i=2,lat2-1
              iglob=istart(mm1)-2+i
              outwork(jglob,iglob)=ges_tv_it(i,j,k)/(ges_prsl(i,j,k ,it)/r100)**rd_over_cp
           end do
        end do
        call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
        if(mype==0) write(lunin) outwork0
     end do

!             output "smoothed pot temp"

     do k=kvar_start(ivar_plot),kvar_end(ivar_plot)
        outwork=zero_single
        if(k>=indices%kps.and.k<=indices%kpe) then
           do j=1,pf2aP1%nlonf
              do i=1,pf2aP1%nlatf
                 tempc(i,j)=theta0f(i,j,k-indices%kps+1)
              end do
           end do
           call fgrid2agrid(pf2aP1,tempc,tempf)
           do j=1,nlon
              do i=1,nlat
                 outwork(j,i)=tempf(i,j)
              end do
           end do
        end if
        call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
        if(mype==0) write(lunin) outwork0
     end do

     if(ivar_plot==3) then 
        nsig_aux=1
      else
        nsig_aux=nsig
     endif

     mm1=mype+1
     do k=1,nsig_aux
        outwork=zero_single
        do j=2,lon2-1
           jglob=jstart(mm1)-2+j
           do i=2,lat2-1
              iglob=istart(mm1)-2+i
              if(ivar_plot==1) outwork(jglob,iglob)=stwork(i,j,k)
              if(ivar_plot==2) outwork(jglob,iglob)=vpwork(i,j,k)
              if(ivar_plot==3) outwork(jglob,iglob)=pwork(i,j)
              if(ivar_plot==4) outwork(jglob,iglob)=twork(i,j,k)
              if(ivar_plot==5) outwork(jglob,iglob)=qwork(i,j,k)
           end do
        end do
!             very slow way to move field from local processor to processor 0

        call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
        if(mype==0) write(lunin) outwork0
     end do

!     output original terrain (slow way to get full 2d field)  -- this is reference field

     mm1=mype+1
     outwork=zero_single
     do j=2,lon2-1
        jglob=jstart(mm1)-2+j
        do i=2,lat2-1
           iglob=istart(mm1)-2+i
           outwork(jglob,iglob)=ges_z_it(i,j)
        end do
     end do
     call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
     if(mype==0) write(lunin) outwork0

!             output "smoothed terrain"

     if(mype==0)PRINT*,'IN ANPREWGT_REG,KPS,KPE=',indices%KPS,indices%KPE
     do k=1, 1  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!kvar_start(ivar_plot),kvar_end(ivar_plot)
        outwork=zero_single
        if(k>=indices%kps.and.k<=indices%kpe) then
           do j=1,pf2aP1%nlonf
              do i=1,pf2aP1%nlatf
                 tempc(i,j)=z0f(i,j,k)!  theta0f(i,j,k-indices%kps+1)
              end do
           end do
           call fgrid2agrid(pf2aP1,tempc,tempf)
           do j=1,nlon
              do i=1,nlat
                 outwork(j,i)=tempf(i,j)
              end do
           end do
        end if
        call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
        if(mype==0) write(lunin) outwork0
     end do

     close(lunin)
!    if(mype>-1000) then
!       call mpi_finalize(i)
!       stop
!    end if

200 continue

end subroutine antest_maps0
!-------------------------------------------------------------------------------------

subroutine antest_maps0_subdomain_option(mype,theta0f,z0f)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    antest_maps0_subdomain_option
!   prgmmr:
!
! abstract: this routine creates output maps of background error correlations 
!
! program history log:
!   2009-09-18  lueken - added subprogram doc block
!   2013-01-15  parrish - ansmoothrf_reg_subdomain_option now uses bundle input.
!                             Add changes to update from old to new version.
!   2013-10-19  todling - metguess now holds background
!
!   input argument list:
!    mype
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind,r_single
  use constants, only: max_varname_length
  use anberror, only: kvar_start,kvar_end,var_names,levs_jdvar,indices,pf2aP1
  use gridmod, only: nsig,nlon,nlat,istart,jstart,lat2,lon2
  use constants, only: zero_single,zero,one,rd_over_cp,r100
  use mpimod, only: ierror,mpi_real4,mpi_real8,mpi_sum,mpi_comm_world
  use guess_grids, only: ntguessig,ges_prsl
  use control_vectors, only: nvars
  use gsi_bundlemod, only: gsi_bundlecreate
  use gsi_bundlemod, only: gsi_bundle
  use gsi_bundlemod, only: gsi_bundlegetvar
  use gsi_bundlemod, only: gsi_bundleputvar
  use gsi_bundlemod, only: gsi_bundledestroy
  use gsi_bundlemod, only: gsi_grid
  use gsi_bundlemod, only: gsi_gridcreate
  use gsi_bundlemod, only: gsi_bundlegetpointer
  use gsi_metguess_mod, only: gsi_metguess_bundle
  use mpeu_util, only: die

  implicit none

  integer(i_kind),intent(in   ) :: mype
  real(r_single) ,intent(in   ) :: theta0f(lat2,lon2,nsig)
  real(r_single) ,intent(in   ) :: z0f(lat2,lon2,nsig)

  character(len=*),parameter::myname='antest_maps0_subdomain_option'
  real(r_kind),dimension(lat2,lon2,nsig):: twork,qwork,stwork,vpwork
  real(r_kind),dimension(lat2,lon2):: pwork
  real(r_single) outwork(nlon,nlat),outwork0(nlon,nlat)
  character(80) ref_plotcor
  character(80) var_plotcor
  character(80) plotname
  integer(i_kind) i_plotcor,j_plotcor,k_plotcor
  integer(i_kind) iloc,jloc,kloc,ier,istatus

  real(r_kind),dimension(:,:  ),pointer::ges_z_it =>NULL()
  real(r_kind),dimension(:,:,:),pointer::ges_tv_it=>NULL()

  real(r_kind) h00,h000
  integer(i_kind) lunin,i,j,k,ivar,iglob,jglob,ivar_plot
  integer(i_kind) it,mm1
  integer(i_kind) lvar
  integer(i_kind):: ips,ipe,jps,jpe,kps,kpe
  integer(i_kind):: nlatf,nlonf
  type(gsi_bundle):: bundle_work
  type(gsi_grid) :: grid
  character(max_varname_length) :: names2dwork(1),names3dwork(4)

  names2dwork(1)='ps'
  names3dwork(1)='sf'
  names3dwork(2)='vp'
  names3dwork(3)='t'
  names3dwork(4)='q'
  ips=indices%ips; ipe=indices%ipe
  jps=indices%jps; jpe=indices%jpe
  kps=indices%kps; kpe=indices%kpe
  nlatf=pf2aP1%nlatf
  nlonf=pf2aP1%nlonf

!*********************************************************************
!          variable names expected for var_plotcor are
!
!    sf  -- stream function
!    vp  -- velocity potential
!    ps  -- surface pressure
!    t   -- virtual temperature
!    q   -- specific humidity
!    oz  -- ozone
!    sst -- sea surface temperature
!    stl -- skin temp over land
!    sti -- skin temp over ice
!    cw  -- cloud water
!*********************************************************************
! Make choices here!
  i_plotcor=460 !440!  430
  j_plotcor=265!  250
  k_plotcor=1
  iloc=i_plotcor-istart(mype+1)+2
  jloc=j_plotcor-jstart(mype+1)+2
!  var_plotcor='sf'
!Note: Must call this subroutine from anprewgt_reg.f90
!Make sure statement has been uncommented!
! End of choice section
!*********************************************************************

     call gsi_gridcreate(grid,lat2,lon2,nsig)
  ref_plotcor='theta'
  ier=0
  it=ntguessig
  call gsi_bundlegetpointer (gsi_metguess_bundle(it),'z' ,ges_z_it,   istatus)
  ier=ier+istatus
  call gsi_bundlegetpointer (gsi_metguess_bundle(it),'tv',ges_tv_it,  istatus)
  ier=ier+istatus
  if(ier/=0) call die(myname,'missing fields, ier= ', ier)
  
  do 200 lvar=1,5
     if (lvar==1)  var_plotcor='sf'
     if (lvar==2)  var_plotcor='vp'
     if (lvar==3)  var_plotcor='ps'
     if (lvar==4)  var_plotcor='t'
     if (lvar==5)  var_plotcor='q'
     lunin=1
     if(mype==0) then
!       open(lunin,file='cormaps_'//trim(var_plotcor),form='unformatted')
!       rewind lunin
     end if
     ivar_plot=0
     do ivar=1,nvars
        if(trim(var_names(ivar))==trim(var_plotcor)) then
           ivar_plot=ivar
           exit
        end if
     end do
     if(ivar_plot==0) then
        write(6,*)' in antest_maps0, variable ',trim(var_plotcor),'  not found.  program stops'
        call mpi_finalize(ierror)
        stop
     end if
     kloc=k_plotcor
     twork=zero
     pwork=zero
     qwork=zero
     stwork=zero
     vpwork=zero
     if(i_plotcor>=ips.and.i_plotcor<=ipe.and. &
        j_plotcor>=jps.and.j_plotcor<=jpe) then
        if(ivar_plot==1)     stwork(iloc,jloc,kloc)=one
        if(ivar_plot==2) vpwork(iloc,jloc,kloc)=one
        if(ivar_plot==3) pwork(iloc,jloc)=one
        if(ivar_plot==4) twork(iloc,jloc,kloc)=one
        if(ivar_plot==5) qwork(iloc,jloc,kloc)=one
     end if

     call gsi_bundlecreate(bundle_work,grid,'bundle work',ier, &
              names2d=names2dwork,names3d=names3dwork,bundle_kind=r_kind)
     call gsi_bundleputvar(bundle_work,'ps',pwork,ier)
     call gsi_bundleputvar(bundle_work,'sf',stwork,ier)
     call gsi_bundleputvar(bundle_work,'vp',vpwork,ier)
     call gsi_bundleputvar(bundle_work,'t' ,twork,ier)
     call gsi_bundleputvar(bundle_work,'q' ,qwork,ier)

     call ansmoothrf_reg_subdomain_option(bundle_work)
     call gsi_bundlegetvar(bundle_work,'ps',pwork,ier)
     call gsi_bundlegetvar(bundle_work,'sf',stwork,ier)
     call gsi_bundlegetvar(bundle_work,'vp',vpwork,ier)
     call gsi_bundlegetvar(bundle_work,'t' ,twork,ier)
     call gsi_bundlegetvar(bundle_work,'q' ,qwork,ier)
     call gsi_bundledestroy(bundle_work,ier)
 !   if(mype==0) write(lunin) ref_plotcor,var_plotcor,j_plotcor,i_plotcor,k_plotcor, &
 !                nlon,nlat,kvar_end(ivar_plot)-kvar_start(ivar_plot)+1
     if(mype==0) write(6,*) ' refvar= ',trim(ref_plotcor),' corvar= ',trim(var_plotcor), &
            '  i,j,k_plotcor =', j_plotcor,i_plotcor,k_plotcor, ' nlon,nlat,nsig=', &
                  nlon,nlat,kvar_end(ivar_plot)-kvar_start(ivar_plot)+1

    !---------------------in case we haven't normalized, divide by value of correlation point
     h00=zero
     if(i_plotcor>=ips.and.i_plotcor<=ipe.and. &
        j_plotcor>=jps.and.j_plotcor<=jpe) then
        if(ivar_plot==1)     h00=stwork(iloc,jloc,kloc)
        if(ivar_plot==2) h00=vpwork(iloc,jloc,kloc)
        if(ivar_plot==3) h00=pwork(iloc,jloc)
        if(ivar_plot==4) h00=twork(iloc,jloc,kloc)
        if(ivar_plot==5) h00=qwork(iloc,jloc,kloc)
     end if
     call mpi_allreduce(h00,h000,1,mpi_real8,mpi_sum,mpi_comm_world,ierror)
     stwork=stwork/h000
     vpwork=vpwork/h000
     pwork=pwork/h000
     twork=twork/h000
     qwork=qwork/h000


!     output original pot temp  (slow way to get full 2d field)  -- this is reference field

     mm1=mype+1
     do k=1,nsig
        outwork=zero_single
        do j=2,lon2-1
           jglob=jstart(mm1)-2+j
           do i=2,lat2-1
              iglob=istart(mm1)-2+i
              outwork(jglob,iglob)=ges_tv_it(i,j,k)/(ges_prsl(i,j,k ,it)/r100)**rd_over_cp
           end do
        end do
        call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
        !if(mype==0) write(lunin) outwork0
        if(mype==0) call outgrads1(outwork0,nlon,nlat,'theta')
     end do

!             output "smoothed pot temp"

     do k=1,nsig
        outwork=zero_single
        do j=2,lon2-1
           jglob=jstart(mm1)-2+j
           do i=2,lat2-1
              iglob=istart(mm1)-2+i
              outwork(jglob,iglob)=theta0f(i,j,k)
           end do
        end do
        call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
        !if(mype==0) write(lunin) outwork0
        if(mype==0) call outgrads1(outwork0,nlon,nlat,'sm_theta')
     end do

     do k=kvar_start(ivar_plot),kvar_end(ivar_plot)
        kloc=levs_jdvar(k)
        outwork=zero_single
        do j=2,lon2-1
           jglob=jstart(mm1)-2+j
           do i=2,lat2-1
              iglob=istart(mm1)-2+i
              if(ivar_plot==1)     outwork(jglob,iglob)=stwork(i,j,kloc)
              if(ivar_plot==2) outwork(jglob,iglob)=vpwork(i,j,kloc)
              if(ivar_plot==3) outwork(jglob,iglob)=pwork(i,j)
              if(ivar_plot==4) outwork(jglob,iglob)=twork(i,j,kloc)
              if(ivar_plot==5) outwork(jglob,iglob)=qwork(i,j,kloc)
           end do
        end do
!             very slow way to move field from local processor to processor 0

        call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
        !if(mype==0) write(lunin) outwork0
        if(mype==0) then
           write(plotname,'("sub_",a)')trim(var_plotcor)
           call outgrads1(outwork0,nlon,nlat,plotname)
        end if
     end do

!     output original terrain (slow way to get full 2d field)  -- this is reference field

     mm1=mype+1
     outwork=zero_single
     do j=2,lon2-1
        jglob=jstart(mm1)-2+j
        do i=2,lat2-1
           iglob=istart(mm1)-2+i
           outwork(jglob,iglob)=ges_z_it(i,j)
        end do
     end do
     call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
     !if(mype==0) write(lunin) outwork0
     if(mype==0) call outgrads1(outwork0,nlon,nlat,'z')

!             output "smoothed terrain"

     if(mype == 0)PRINT*,'IN ANPREWGT_REG_subdomain_option,KPS,KPE=',KPS,KPE

     do k=1,1
        outwork=zero_single
        do j=2,lon2-1
           jglob=jstart(mm1)-2+j
           do i=2,lat2-1
              iglob=istart(mm1)-2+i
              outwork(jglob,iglob)=z0f(i,j,k)
           end do
        end do
        call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
        !if(mype==0) write(lunin) outwork0
        if(mype==0) call outgrads1(outwork0,nlon,nlat,'sm_z')
     end do

     !close(lunin)
200 continue

end subroutine antest_maps0_subdomain_option

!-------------------------------------------------------------------------------------

subroutine outgrads1(f,nx,ny,label)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    outgrads1
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-09-18  lueken - added subprogram doc block
!   2012-12-11  parrish - assign np a value.
!
!   input argument list:
!    label
!    nx,ny
!    f
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block
  use kinds, only: i_kind,r_single
  implicit none

  character(*)   ,intent(in   ) :: label
  integer(i_kind),intent(in   ) :: nx,ny
  real(r_single) ,intent(in   ) :: f(nx,ny)

  integer(i_kind) i,l,next,last,np,ntime,ioutdat,ioutcor,koutmax
  real(r_single) rlonmap0,undef,dlonmap,pinc,startp,rlatmap0,dlatmap
  character(80) dsdes,dsdat
  character(80) datdes(1000)
  character(1) blank
  data blank/' '/
  data undef/-9.99e33_r_single/

  ioutcor=10
  ioutdat=11
  np=1

  write(dsdes,'(a,".des")')trim(label)
  write(dsdat,'(a,".dat")')trim(label)
  open(unit=ioutcor,file=dsdes,form='formatted')
  open(unit=ioutdat,file=dsdat,form='unformatted')
  ntime=1
  rlonmap0=1._r_single
  dlonmap=1._r_single
  rlatmap0=1._r_single
  dlatmap=1._r_single
  startp=1._r_single
  pinc=1._r_single
  koutmax=1
  do i=1,1000
     write(datdes(i),'(80a1)')(blank,l=1,80)
  end do
  write(datdes(1),'("DSET ",a)')trim(dsdat)
  write(datdes(2),'("options big_endian sequential")')
  write(datdes(3),'("TITLE ",a)')trim(label)
  write(datdes(4),'("UNDEF ",e11.2)')undef
  write(datdes(5),'("XDEF ",i5," LINEAR ",f7.2,f7.2)')nx,rlonmap0,dlonmap
  write(datdes(6),'("YDEF ",i5," LINEAR ",f7.2,f7.2)')ny,rlatmap0,dlatmap
  next=7
  write(datdes(next),'("ZDEF ",i5," LINEAR ",f7.2,f7.2)')np,startp,pinc
  next=next+1
  write(datdes(next),'("TDEF ",i5," LINEAR 0Z23may1992 24hr")')koutmax
  next=next+1
  write(datdes(next),'("VARS 1")')
  next=next+1
  write(datdes(next),'("f   ",i5," 99 f   ")')np
  next=next+1
  write(datdes(next),'("ENDVARS")')
  last=next
  write(ioutcor,'(a80)')(datdes(i),i=1,last)
  close(ioutcor)

  write(ioutdat) f
  close(ioutdat)

end subroutine outgrads1