subroutine read_cmaq_files(mype)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    read_cmaq_files   same as read_files, but cmaq
!pagowski - based on read_wrf_mass_files.f90 but for cmaq
!   abstract: read cmaq input file, add analysis, and write out analysis
!   in cmaq intermediate format, based on other writers
!   all to be replaced when cmaq nemsio is ready
! along with changes to chemmod.f90
!
! abstract: figure out available time levels of background fields for 
!             later input. This is patterned after read_wrf_mass_files.
!
! program history log:
!   2010-09-07  pagowski
!   2015-09-17  Thomas  - add l4densvar to data selection procedure
!   
!   input argument list:
!     mype     - pe number
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block
  
  use kinds, only: i_kind,r_kind,r_single
  use mpimod, only: mpi_comm_world,ierror,mpi_rtype,npe
  use guess_grids, only: nfldsig,nfldsfc,ntguessig,ntguessfc,&
       ifilesig,ifilesfc,hrdifsig,hrdifsfc,create_gesfinfo
  use guess_grids, only: hrdifsig_all,hrdifsfc_all
  use gsi_4dvar, only: l4dvar, iwinbgn, winlen, nhr_assimilation,&
       l4densvar
  use gridmod, only: regional_time,regional_fhr
  use constants, only: zero,one,r60inv
  use obsmod, only: iadate,time_offset
  implicit none
  
! declare passed variables
  integer(i_kind),intent(in   ) :: mype
  
! declare local parameters
  real(r_kind),parameter:: r0_001=0.001_r_kind
  
! declare local variables
  logical(4) fexist
  character(6) filename
  integer(i_kind) :: i,j,iwan,npem1
  integer(i_kind) :: nhr_half
  integer(i_kind) :: nminanl,nmings,nming2,ndiff
  integer(i_kind),dimension(4):: idateg
  integer(i_kind),dimension(5):: idate5
  real(r_single) hourg4
  real(r_kind) hourg,temp,t4dv
  real(r_kind),dimension(202,2):: time_ges
  
! start read_cmaq_files here.
  nhr_half=nhr_assimilation/2
  if(nhr_half*2 < nhr_assimilation) nhr_half=nhr_half+1  
  npem1=npe-1

  do i=1,202
     time_ges(i,1) = 999_r_kind
     time_ges(i,2) = 999_r_kind
  end do
  

! let a single task query the guess files.
  if(mype==npem1) then
     
!    convert analysis time to minutes relative to fixed date
     call w3fs21(iadate,nminanl)
     write(6,*)'read_cmaq_files:  analysis date,minutes ',iadate,nminanl
     
!    check for consistency of times from sigma guess files.
     iwan=0
     do i=0,99
        write(filename,100)i
100     format('sigf',i2.2)
        inquire(file=filename,exist=fexist)
        if(fexist)then
           idateg(1)=regional_time(4)  !  hour
           idateg(2)=regional_time(2)  !  month
           idateg(3)=regional_time(3)  !  day
           idateg(4)=regional_time(1)  !  year
           hourg4= regional_fhr        !  fcst hour
           hourg = hourg4
           idate5(1)=idateg(4); idate5(2)=idateg(2)
           idate5(3)=idateg(3); idate5(4)=idateg(1); idate5(5)=regional_time(5)
           call w3fs21(idate5,nmings)
           nming2=nmings+60*hourg
           write(6,*)'read_cmaq_files:  sigma guess file, nming2 ',hourg,idateg,nming2
           t4dv=real((nming2-iwinbgn),r_kind)*r60inv
           if (l4dvar.or.l4densvar) then
              if (t4dv<zero .or. t4dv>winlen) cycle
           else
              ndiff=nming2-nminanl
              if(abs(ndiff) > 60*nhr_half ) cycle
           endif
           
           iwan=iwan+1
           time_ges(iwan,1) = real((nming2-iwinbgn),r_kind)*r60inv
           time_ges(iwan+100,1)=i+r0_001
        end if
     end do
     
     time_ges(201,1)=one
     time_ges(202,1)=one
     
     if(iwan > 1)then
        do i=1,iwan
           do j=i+1,iwan 
              if(time_ges(j,1) < time_ges(i,1))then
                 temp=time_ges(i+100,1)
                 time_ges(i+100,1)=time_ges(j+100,1)
                 time_ges(j+100,1)=temp
                 temp=time_ges(i,1)
                 time_ges(i,1)=time_ges(j,1)
                 time_ges(j,1)=temp
              end if
           end do
           if(abs(time_ges(i,1)-time_offset) < r0_001)time_ges(202,1) = i
        end do
     end if
     time_ges(201,1) = iwan+r0_001
     
!surface fields in cmaq are assigned to defaults
     
  end if
  
! broadcast guess file information to all tasks
  call mpi_bcast(time_ges,404,mpi_rtype,npem1,mpi_comm_world,ierror)
  nfldsig   = nint(time_ges(201,1))
  
  nfldsfc   = 1
  
! allocate space for guess information files
  call create_gesfinfo
  
  do i=1,nfldsig
     ifilesig(i) = -100
     hrdifsig(i) = zero
  end do
  
  do i=1,nfldsfc
     ifilesfc(i) = -100
     hrdifsfc(i) = zero
  end do
  
! load time information for sigma guess field sinfo into output arrays
  ntguessig = nint(time_ges(202,1))
  
  do i=1,nfldsig
     hrdifsig(i) = time_ges(i,1)
     ifilesig(i) = nint(time_ges(i+100,1))
     hrdifsig_all(i) = hrdifsig(i)
  end do
  
  if(mype == 0) write(6,*)'read_cmaq_files:  sigma fcst files used in analysis  :  ',&
       (ifilesig(i),i=1,nfldsig),(hrdifsig(i),i=1,nfldsig),ntguessig
  
  ntguessfc = ntguessig
  
  do i=1,nfldsfc
     hrdifsfc(i) = hrdifsig(ntguessig)
     ifilesfc(i) = ifilesig(ntguessig)
     hrdifsfc_all(i) = hrdifsfc(i)
  end do
  
  if(mype == 0) write(6,*)'read_cmaq_files:  surface fcst files used in analysis:  ',&
       (ifilesfc(i),i=1,nfldsfc),(hrdifsfc(i),i=1,nfldsfc),ntguessfc
  
  return
  
end subroutine read_cmaq_files

subroutine read_cmaq_guess(mype)
  
!$$$  subprogram documentation block
!                .      .    .                                       .
!read intermediate cmaq binary file that included aerosol species
!and meteorology, all input variables in cmaq binary are on A-grid
!   2010-08-02 pagowski - largely based on reading wrf intermediate binary file
!                         routine read_wrf_mass_netcdf_guess
!   2011-06-29  todling - no explict reference to internal bundle arrays
!   2013-10-19  todling - metguess now holds background
!
!   input argument list:
!     mype     - pe number
!
  
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  
  use kinds, only: r_kind,r_single,i_kind
  
  use mpimod, only: mpi_real4,mpi_comm_world,npe,ierror
  use guess_grids, only: nfldsig,ifilesig,ges_tsen
  use guess_grids, only: isli,fact10,sfct,dsfct,sno,veg_type,veg_frac,&
       soil_type,soil_temp,soil_moi,sfc_rough
  
  use gridmod, only: lat2,lon2,nlat_regional,nlon_regional,&
       nsig,ijn_s,displs_s,itotsub

  use constants, only: zero,one,fv
  use gsi_io, only: lendian_in
  
  use gsi_bundlemod, only : gsi_bundlegetpointer
  use gsi_metguess_mod, only: gsi_metguess_bundle
  use gsi_chemguess_mod, only: gsi_chemguess_bundle

  use constants, only : max_varname_length
  use chemmod, only : nmet2d_cmaq,nmet3d_cmaq,&
       naero_cmaq,pm2_5_guess,init_pm2_5_guess
    
  use mpeu_util, only: die
  implicit none
  
! declare passed variables
  integer(i_kind),intent(in):: mype
  integer(i_kind) :: status

  
!fields below are assigned defaults
!isli,fact10,sfct,dsfct,sno,veg_type,veg_frac,
!soil_type,soil_temp,soil_moi,sfc_rough

!for cmaq declare default values for constant/semi constant sfc fields
  character(len=*),parameter::myname='read_cmaq_guess'
  integer(i_kind),parameter :: isli_default=1
  real(r_kind),parameter :: fact10_default=one,&
       sno_default=zero,veg_type_default=2.0_r_kind,&
       veg_frac_default=0.1_r_kind,soil_type_default=4_r_kind,&
       dsfct_default=zero,soil_moi_default=0.3_r_kind,&
       rough_default=0.05_r_kind
! declare local variables

  
  character(len=max_varname_length) :: cvar

! declare local variables
  
  real(r_kind),parameter:: r0_001=0.001_r_kind
  
! other internal variables
  real(r_single) tempa(itotsub)
  real(r_single),allocatable::temp1(:,:)
  real(r_single),allocatable::all_loc(:,:,:)
  real(r_kind),pointer,dimension(:,:,:)::ptr3d=>NULL()
  character(6) filename 
  integer(i_kind) irc_s_reg(npe),ird_s_reg(npe)
  integer(i_kind) ifld,im,jm,lm,num_cmaq_fields
  integer(i_kind) num_all_fields,num_loc_groups,num_all_pad
  integer(i_kind) i,icount,icount_prev,it,j,k
  integer(i_kind) i_0,i_psfc,i_fis,i_t,i_q,i_u,i_v
  
  integer(i_kind) :: kfis,kpsfc,kt,kq,ku,kv,ier,istatus
  
  integer(i_kind) :: nskip
  
  real(r_kind),dimension(:,:  ),pointer::ges_ps_it=>NULL()
  real(r_kind),dimension(:,:  ),pointer::ges_z_it =>NULL()
  real(r_kind),dimension(:,:,:),pointer::ges_u_it =>NULL()
  real(r_kind),dimension(:,:,:),pointer::ges_v_it =>NULL()
  real(r_kind),dimension(:,:,:),pointer::ges_tv_it=>NULL()
  real(r_kind),dimension(:,:,:),pointer::ges_q_it =>NULL()

  if(mype==0) write(6,*)' at 0 in read_cmaq_guess'
  
! big section of operations done only on first outer iteration
  
  im=nlon_regional
  jm=nlat_regional
  lm=nsig
  
  i_fis=1
  i_psfc=i_fis+1
  i_t=i_psfc+1
  i_q=i_t+nsig
  i_u=i_q+nsig
  i_v=i_u+nsig
  
!  CMAQ input grid dimensions in module gridmod
!      These are the following:
!          im -- number of x-points on C-grid
!          jm -- number of y-points on C-grid
!          lm -- number of vertical levels ( = nsig for now)
  
  
  num_cmaq_fields=nmet2d_cmaq+(nmet3d_cmaq+naero_cmaq)*nsig
  
  num_all_fields=num_cmaq_fields*nfldsig
  num_loc_groups=num_all_fields/npe
  
  if(mype==0) then 
     write(6,'(" at 1 in read_cmaq_guess, lm            =",i6)')lm
     write(6,'(" at 1 in read_cmaq_guess, num_cmaq_fields=",i6)')num_cmaq_fields
     write(6,'(" at 1 in read_cmaq_guess, nfldsig       =",i6)')nfldsig
     write(6,'(" at 1 in read_cmaq_guess, num_all_fields=",i6)')num_all_fields
     write(6,'(" at 1 in read_cmaq_guess, npe           =",i6)')npe
     write(6,'(" at 1 in read_cmaq_guess, num_loc_groups=",i6)')num_loc_groups
  endif

  do 
     num_all_pad=num_loc_groups*npe
     if(num_all_pad >= num_all_fields) exit
     num_loc_groups=num_loc_groups+1
  end do
  
  if(mype==0) then
     write(6,'(" at 2 in read_cmaq_guess, num_all_pad   =",i6)')num_all_pad
     write(6,'(" at 2 in read_cmaq_guess, num_loc_groups=",i6)')num_loc_groups
  endif

  allocate(all_loc(lat2,lon2,num_all_pad))
  
  allocate(temp1(im,jm))
  
  do i=1,npe
     irc_s_reg(i)=ijn_s(mype+1)
  end do
  ird_s_reg(1)=0
  do i=1,npe
     if(i /= 1) ird_s_reg(i)=ird_s_reg(i-1)+irc_s_reg(i-1)
  end do
  
!read cmaq intermediate binary
  
!    this is done by reading in parallel from every pe, and redistributing
!    to local domains once for every npe fields read in, using 
!    mpi_all_to_allv
  
  icount=0
  icount_prev=1
  
  do it=1,nfldsig
     write(filename,'("sigf",i2.2)')ifilesig(it)
     open(lendian_in,file=filename,form='unformatted') 
     rewind(lendian_in)
     
     if (mype==0) &
          write(6,*)'read_cmaq_guess:  open cmaq input=',&
          lendian_in,' to file=',filename
     
     read(lendian_in)nskip
     
     do ifld=1,nskip !need to skip grid variables
        read(lendian_in)
     enddo
     
     do ifld=1,num_cmaq_fields
        icount=icount+1
        if(mype==mod(icount-1,npe)) then
           read(lendian_in)((temp1(i,j),i=1,im),j=1,jm)
           call fill_mass_grid2t(temp1,im,jm,tempa,1)
        else
           read(lendian_in)
        end if
!          distribute to local domains everytime we have npe fields
        if(mod(icount,npe) == 0.or.icount==num_all_fields) then
           call mpi_alltoallv(tempa,ijn_s,displs_s,mpi_real4, &
                all_loc(1,1,icount_prev),irc_s_reg,ird_s_reg,mpi_real4,mpi_comm_world,ierror)
           icount_prev=icount+1
        end if
     end do
     close(lendian_in)
  end do
  
  do it=1,nfldsig
     i_0=(it-1)*num_cmaq_fields
     kfis=i_0+i_fis
     kpsfc=i_0+i_psfc
     kt=i_0+i_t-1 !since +1 in sigma loop
     kq=i_0+i_q-1 
     ku=i_0+i_u-1
     kv=i_0+i_v-1

     ier=0
     call gsi_bundlegetpointer (gsi_metguess_bundle(it),'ps',ges_ps_it,  istatus)
     ier=ier+istatus
     call gsi_bundlegetpointer (gsi_metguess_bundle(it),'z' ,ges_z_it ,  istatus)
     ier=ier+istatus
     call gsi_bundlegetpointer (gsi_metguess_bundle(it),'u' ,ges_u_it ,  istatus)
     ier=ier+istatus
     call gsi_bundlegetpointer (gsi_metguess_bundle(it),'v' ,ges_v_it ,  istatus)
     ier=ier+istatus
     call gsi_bundlegetpointer (gsi_metguess_bundle(it),'tv',ges_tv_it,  istatus)
     ier=ier+istatus
     call gsi_bundlegetpointer (gsi_metguess_bundle(it),'q' ,ges_q_it ,  istatus)
     ier=ier+istatus
     if(ier/=0) call die(myname,'missing fields, ier= ', ier)

     do k=1,nsig
        kt=kt+1
        kq=kq+1
        ku=ku+1
        kv=kv+1
        do i=1,lon2
           do j=1,lat2
              ges_u_it(j,i,k) = real(all_loc(j,i,ku),r_kind)
              ges_v_it(j,i,k) = real(all_loc(j,i,kv),r_kind)
              ges_tsen(j,i,k,it) = real(all_loc(j,i,kt),r_kind)
              ges_q_it(j,i,k) = real(all_loc(j,i,kq),r_kind)
!                convert guess mixing ratio to specific humidity
              ges_q_it(j,i,k) = ges_q_it(j,i,k)/(one+ges_q_it(j,i,k))
              ges_tv_it(j,i,k) = ges_tsen(j,i,k,it) * &
                   (one+fv*ges_q_it(j,i,k))
           end do
        end do

     enddo

     do i=1,lon2
        do j=1,lat2
           ges_z_it(j,i) = real(all_loc(j,i,kfis),r_kind)
           ges_ps_it(j,i)=r0_001*real(all_loc(j,i,kpsfc),r_kind)! convert from Pa to cb
        end do
     end do

     do i=1,lon2
        do j=1,lat2
           isli(j,i,it)=isli_default
           fact10(j,i,it)=fact10_default
           sfct(j,i,it)=ges_tsen(j,i,1,it)
           dsfct(j,i,it)=dsfct_default
           sno(j,i,it)=sno_default
           veg_type(j,i,it)=veg_type_default
           veg_frac(j,i,it)=veg_frac_default
           soil_type(j,i,it)=soil_type_default
           soil_temp(j,i,it)=ges_tsen(j,i,1,it)
           soil_moi(j,i,it)=soil_moi_default
           sfc_rough(j,i,it)=rough_default
        end do
     end do

!do the bundle assignment
     
!the last variable written before aero is "v"
!assign icount for aerosols
     
     call init_pm2_5_guess

     icount=kv

     do ifld=1,naero_cmaq
        
        do k=1,nsig
           
!aerosols come after v - just do the counts
           icount=icount+1
     
           do i=1,lon2
              do j=1,lat2
                 pm2_5_guess(j,i,k)=pm2_5_guess(j,i,k)+real(all_loc(j,i,icount),r_kind)
              enddo
           enddo
        enddo

     enddo

     cvar='pm2_5'

     call gsi_bundlegetpointer (gsi_chemguess_bundle(it), &
          cvar, ptr3d, status)
     
     if(status/=0) then 
        write(6,*)'problem with bundle pointers in cmaq_routines'
        call stop2(425)
     endif
     
     do k=1,nsig
        do i=1,lon2
           do j=1,lat2
              ptr3d(j,i,k)=pm2_5_guess(j,i,k)
           enddo
        enddo
     enddo
     
  enddo

  deallocate(all_loc)
  deallocate(temp1)

  return
  
end subroutine read_cmaq_guess

subroutine make_sigf
  
!since cmaq file is read without conversions just link cmaq inout to sigf
!temporary routine until cmaq input is ready in nemsio 
  
!   2010-09-15  pagowski - link cmaq_input to sigf
  
  use gsi_4dvar, only: nhr_assimilation
  use chemmod, only : maxstr,in_fname
  use kinds, only: i_kind
  
  implicit none
  
  character(len=maxstr) :: cmaq_infile_name

  character(len=6) :: filename
  character(len=maxstr) :: command
  
  logical :: fexist
  
  integer(i_kind) :: ihr

  cmaq_infile_name=in_fname
  
! get regional constants
  
  ihr=-999
  
  inquire(file=trim(cmaq_infile_name),exist=fexist)
  if (fexist) then
     ihr=nhr_assimilation
     write(filename,100) nhr_assimilation
100  format('sigf',i2.2)
     command='if [ ! -r '//trim(filename)//' ] ;then ln -sf '//trim(cmaq_infile_name)//' '//trim(filename)//' ;fi'
     call system(command)
  endif
  
end subroutine make_sigf

subroutine write_cmaq(mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    write_cmaq - write out cmaq intermediate binary file
!   prgmmr: pagowski   date: 2010-09-22
!
! program history log:
!   2010-09-23  pagowski 
!   2013-10-24  todling - revisit strip interface
!
!   input argument list:
!     mype     - pe number
!
!   output argument list:
!     no output arguments
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,r_single,i_kind
  use guess_grids, only: ntguessig,ifilesig
  use mpimod, only: mpi_comm_world,ierror,mpi_real4
  use gridmod, only: lat2,iglobal,itotsub,&
       lon2,nsig,lon1,lat1,nlon_regional,nlat_regional,ijn,displs_g,&
       strip
  use constants, only: zero_single,&
       tiny_single,max_varname_length
  use gsi_io, only: lendian_in, lendian_out
  use chemmod, only : ngrid2d_cmaq,nmet2d_cmaq,nmet3d_cmaq,&
       naero_cmaq,aeronames_cmaq,pm2_5_guess,diag_incr,maxstr,out_fname,&
       incr_fname
  use gsi_chemguess_mod, only: gsi_chemguess_bundle
  use gsi_bundlemod, only : gsi_bundlegetpointer

  implicit none
  
! declare passed variables
  integer(i_kind),intent(in   ) :: mype
  
! declare local parameters

!wish this could be combined to a module rather then repeat

! declare local variables
  integer(i_kind) im,jm,lm
  real(r_single),allocatable::temp1(:),tempa(:),tempb(:)
  real(r_single),allocatable::all_loc(:,:,:),ratio(:,:,:)
  real(r_single),allocatable::incr(:,:,:)
  real(r_single),allocatable::strp(:)
  real(r_kind),pointer,dimension(:,:,:)::ptr3d=>NULL()
  character(6) filename
  integer(i_kind) :: i,j,k,it
  integer(i_kind) :: ncmaqin,nskip
  integer(i_kind) :: regional_time0(6),nlon_regional0,nlat_regional0,nsig0
  real(r_single) aeta10(nsig),eta10(nsig+1),aeta20(nsig),eta20(nsig+1),pt0,pdt0
  integer(i_kind) :: ifld
  integer :: status
  character(len=max_varname_length) :: cvar
  character(len=maxstr) :: cmaq_outfile_name,cmaq_incrementfile_name

  cmaq_outfile_name=out_fname
  cmaq_incrementfile_name=incr_fname

  ncmaqin=ngrid2d_cmaq+nmet2d_cmaq+nmet3d_cmaq*nsig

  im=nlon_regional
  jm=nlat_regional
  lm=nsig
  
  allocate(all_loc(lat2,lon2,nsig),ratio(lat2,lon2,nsig))
  allocate(incr(lat2,lon2,nsig))
  allocate(strp(lat1*lon1))
  
  allocate(temp1(im*jm))
  allocate(tempa(itotsub),tempb(itotsub))

  if(mype == 0) then
     write(filename,'("sigf",i2.2)')ifilesig(ntguessig)
     open (lendian_in,file=filename,form='unformatted')
     open (lendian_out,file=trim(cmaq_outfile_name),form='unformatted')
     rewind lendian_in ; rewind lendian_out
  end if
  
! convert analysis variables to cmaq variables
  it=ntguessig
  
  if(mype == 0) then
     read(lendian_in) nskip
     read(lendian_in) regional_time0,nlon_regional0,nlat_regional0,nsig0,pt0,pdt0
     write(lendian_out) regional_time0,nlon_regional0,nlat_regional0,&
          nsig0

     read(lendian_in) aeta10,aeta20
     read(lendian_in) eta10,eta20
     
     do i=1,ncmaqin
        read(lendian_in)temp1
     enddo

  end if

  cvar='pm2_5'

  call gsi_bundlegetpointer (gsi_chemguess_bundle(it), &
       cvar, ptr3d, status)

  if(status/=0) then
     write(6,*)'problem with bundle pointers in cmaq_routines'
     call stop2(426)
  endif

  do k=1,nsig
     do i=1,lon2
        do j=1,lat2
           incr(j,i,k)=ptr3d(j,i,k)-pm2_5_guess(j,i,k)
           ratio(j,i,k)=max(min(1.0_r_single+incr(j,i,k)/pm2_5_guess(j,i,k),10.0_r_single),&
                tiny_single)
           all_loc(j,i,k)=ratio(j,i,k)
        enddo
     enddo
  enddo

  do ifld=1,naero_cmaq
     
     if (mype==0) then 
        write(lendian_out)aeronames_cmaq(ifld)
     endif

     cvar=trim(aeronames_cmaq(ifld))
     
     do k=1,nsig
        
        call strip(all_loc(:,:,k),strp)
        
        call mpi_gatherv(strp,ijn(mype+1),mpi_real4,&
             tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
        
        if(mype == 0) then

           read(lendian_in)temp1

           call fill_mass_grid2t(temp1,im,jm,tempb,2)
        
           do i=1,iglobal
              tempa(i)=tempb(i)*tempa(i)
           end do

           do i=1,im*jm
              temp1(i)=zero_single
           enddo

           call unfill_mass_grid2t(tempa,im,jm,temp1)

           write(lendian_out)temp1

        end if

     end do
     
  enddo
  
  if (mype==0) then
     close(lendian_in)
     close(lendian_out)
  endif

  if (diag_incr) then

     open(unit=lendian_out,file=trim(incr_fname),form='unformatted')

     if (mype==0)  write(lendian_out)int(im),int(jm),int(nsig)
     
     do k=1,nsig

        call strip(incr(:,:,k),strp)

        call mpi_gatherv(strp,ijn(mype+1),mpi_real4,&
             tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)

        if(mype == 0) then
           do i=1,im*jm
              temp1(i)=zero_single
           enddo
           call unfill_mass_grid2t(tempa,im,jm,temp1)
           write(lendian_out)temp1
        endif

     enddo

     if (mype == 0) close(lendian_out)
     
  endif
  
  call mpi_barrier(mpi_comm_world,ierror)
  
  deallocate(temp1)  
  deallocate(all_loc)
  deallocate(strp)
  deallocate(tempa)
  deallocate(tempb)
  deallocate(incr)

end subroutine write_cmaq