module cmaqmod
  use kinds
! ajl need to store some things here to pass to other places
  integer ncidcmaq
  real(r_kind), allocatable, dimension(:,:,:,:) :: cmaqpsol,cmaqpdash
  real(r_kind), allocatable, dimension(:,:) :: cmaqpbar,szcmaq
  real(r_kind), allocatable, dimension(:,:) :: cmaqpdashbar
  integer(i_kind), allocatable :: i2d(:),j2d(:),i2dg(:),j2dg(:)
  character *6 nameno2
  integer(i_kind),allocatable :: ncidcmaqd(:)
  real(r_kind), allocatable :: cmaqno2(:,:,:,:)
  real(r_kind) :: no2berrormult
  data nameno2/'NO2'/,no2berrormult/1.0/
  save nameno2,no2berrormult
end module cmaqmod
#define NODIAGCMAQ
    subroutine testnetcdfcmaq(fname,lnetcdf)
    use kinds
    use netcdf
    use cmaqmod
    use mpimod,only : mype,mpi_comm_world
    use mpeu_mpif, only : mpi_logical
    implicit none
    integer ierr
    logical lnetcdf
    character (*) fname
    character*20 cno2berrormult
    cno2berrormult=' '
    call getenv('NO2BERRORMULT',cno2berrormult)
    if(cno2berrormult.ne.' ')then
      read(cno2berrormult,*)no2berrormult
      if(mype.eq.0)then
        write(6,*)'no2berrormult ',no2berrormult
      endif
    endif
    if(mype.eq.0)then
      write(6,*)'fname',trim(fname)
      ierr=nf90_open(fname,nf90_write,ncidcmaq)
      if(ierr.eq.nf90_noerr)then
        write(6,*)'open cmaq file ',trim(fname)
        call flush(6)
        lnetcdf=.true.
      else
        lnetcdf=.false.
      endif
    endif
    call mpi_bcast(lnetcdf,1,mpi_logical,0,mpi_comm_world,ierr)
    return
    end subroutine testnetcdfcmaq
    subroutine readcmaqnetcdfdim(cmaq_time,nlon_cmaq,nlat_cmaq,nsig_cmaq)
    use kinds
    use netcdf
    use mpimod
    use cmaqmod
    implicit none
    integer cmaq_time(6)  ,rowdimid,coldimid,laydimid,nrowin,ncolin,nlayin
    integer nlon_cmaq,nlat_cmaq,nsig_cmaq,ierr,idimuse
    character*20 dimnames(3,2)
    data dimnames/'south-north','west-east','lev','ROW','COL','LAY'/
    idimuse=1 ! use wrf names for dims
!    idimuse=2 ! use CMAQ names for dims
!   get from cmaq netcdf file
    if(mype.eq.0)then
      ierr=nf90_inq_dimid(ncidcmaq,dimnames(1,idimuse),rowdimid)
      if(ierr.ne.nf90_noerr)then
        write(6,*)'error getting rowdimid cmaq file '
        write(6,*)trim(nf90_strerror(ierr))
      endif
      ierr=nf90_inquire_dimension(ncidcmaq,rowdimid,len=nrowin)
      if(ierr.ne.nf90_noerr)then
        write(6,*)'error getting nrow cmaq file '
        write(6,*)trim(nf90_strerror(ierr))
      endif
      ierr=nf90_inq_dimid(ncidcmaq,dimnames(2,idimuse),coldimid)
      if(ierr.ne.nf90_noerr)then
        write(6,*)'error getting rowdimid cmaq file '
        write(6,*)trim(nf90_strerror(ierr))
      endif
      ierr=nf90_inquire_dimension(ncidcmaq,coldimid,len=ncolin)
      if(ierr.ne.nf90_noerr)then
        write(6,*)'error getting nrow cmaq file '
        write(6,*)trim(nf90_strerror(ierr))
      endif
      ierr=nf90_inq_dimid(ncidcmaq,dimnames(3,idimuse),laydimid)
      if(ierr.ne.nf90_noerr)then
        write(6,*)'error getting rowdimid cmaq file '
        write(6,*)trim(nf90_strerror(ierr))
      endif
      ierr=nf90_inquire_dimension(ncidcmaq,laydimid,len=nlayin)
      if(ierr.ne.nf90_noerr)then
        write(6,*)'error getting nrow cmaq file '
        write(6,*)trim(nf90_strerror(ierr))
      endif
      nlat_cmaq=nrowin
      nlon_cmaq=ncolin
      nsig_cmaq=nlayin
      ierr=nf90_get_att(ncidcmaq,NF90_GLOBAL,'TIME',cmaq_time)
      write(6,*)'cmaq_time',cmaq_time
      call flush(6)
    endif
    call mpi_bcast(nlat_cmaq,1,mpi_integer4,0,mpi_comm_world,ierr)
    call mpi_bcast(nlon_cmaq,1,mpi_integer4,0,mpi_comm_world,ierr)
    call mpi_bcast(nsig_cmaq,1,mpi_integer4,0,mpi_comm_world,ierr)
    call mpi_bcast(cmaq_time,6,mpi_integer4,0,mpi_comm_world,ierr)
    end subroutine readcmaqnetcdfdim

    subroutine readcmaqnetcdfgrid(glat,dx_mc,glon,dy_mc,nlon_cmaq,nlat_cmaq)
    use kinds
    use netcdf
    use cmaqmod
    use mpimod
    implicit none
    integer(i_kind) :: nlon_cmaq,nlat_cmaq,ierr
    real(r_single),dimension(nlon_cmaq,nlat_cmaq) :: glat,dx_mc,glon,dy_mc
    integer idvar,nwds
    if(mype.eq.0)then
      ierr=nf90_inq_varid(ncidcmaq,'XLAT',idvar)
      if(ierr /= nf90_NOERR)then
        write(6,*)'XLAT error'
      endif
      ierr=nf90_get_var(ncidcmaq,idvar,glat)
      if(ierr /= nf90_NOERR)then
        write(6,*)'XLAT get error'
      endif
      ierr=nf90_inq_varid(ncidcmaq,'XLONG',idvar)
      if(ierr /= nf90_NOERR)then
        write(6,*)'GLONT error'
      endif
      ierr=nf90_get_var(ncidcmaq,idvar,glon)
      if(ierr /= nf90_NOERR)then
        write(6,*)'GLON GET error'
      endif
      ierr=nf90_inq_varid(ncidcmaq,'DXY_MC',idvar)
      if(ierr /= nf90_NOERR)then
        write(6,*)'DXY_MC error'
      endif
      ierr=nf90_get_var(ncidcmaq,idvar,dx_mc)
      if(ierr /= nf90_NOERR)then
        write(6,*)'dx_mc error'
      endif
    endif
    nwds=nlon_cmaq*nlat_cmaq
    call mpi_bcast(glat,nwds,mpi_real4,0,mpi_comm_world,ierr)
    call mpi_bcast(glon,nwds,mpi_real4,0,mpi_comm_world,ierr)
    call mpi_bcast(dx_mc,nwds,mpi_real4,0,mpi_comm_world,ierr)
    dy_mc=dx_mc
    return
    end subroutine readcmaqnetcdfgrid
subroutine readcmaq2d(data2d,name)
  use gridmod, only : nlon,nlat
  use mpimod, only : mype
  use kinds, only : r_kind,i_kind
  use netcdf
  use cmaqmod, only : ncidcmaq
  implicit none
  character *(*) name
  real(r_kind) data2d(nlat,nlon),data2din(nlon,nlat)
  integer(i_kind) ierr,ipnt,i,j
  if(mype.ne.0)return
!  write(6,*)'readcmaq2d',name
!  call flush(6)
  ierr=nf90_inq_varid(ncidcmaq,name,ipnt)
  if(ierr /= NF90_NOERR)then
    write(6,*)'readcmaq2d var not found',trim(name)
  endif
  ierr=nf90_get_var(ncidcmaq,ipnt,data2din)
  if(ierr /= NF90_NOERR)then
    write(6,*)'error reading ',trim(name)
  endif
  call flush(6)
  do i=1,nlon
    do j=1,nlat
      data2d(j,i)=data2din(i,j)
    end do
  end do
  return
end subroutine readcmaq2d
subroutine readcmaq3d(data3d,nl,name)
  use gridmod, only : nlon,nlat
  use mpimod, only : mype
  use kinds, only : r_kind,i_kind,r_double,r_single
  use netcdf
  use cmaqmod, only : ncidcmaq
  implicit none
  character *(*) name
  real(r_kind) data3d(nlat,nlon,nl),data3din(nlon,nlat,nl)
  integer(i_kind) ierr,ipnt,nl,i,j,k
  if(mype.ne.0)return
  ierr=nf90_inq_varid(ncidcmaq,name,ipnt)
  if(ierr /= NF90_NOERR)then
    write(6,*)'readcmaq3d var not found',trim(name)
    call flush(6)
    call stop2(876)
 
  endif
  ierr=nf90_get_var(ncidcmaq,ipnt,data3din)
  if(ierr /= NF90_NOERR)then
    write(6,*)'error reading ',trim(name)
    call flush(6)
    call stop2(875)
  endif
!   write(6,*)'r_kind',r_kind,'r_double',r_double,'r_single',r_single
!   write(6,*)'shape',shape(data3din),'ipnt',ipnt
   write(6,*)name,'readcmaq3di data3din max',maxval(data3din),minval(data3din),'nlon',nlon,nlat,nl
   call flush(6)
! note 1 is bottom for global GSI and also cmaq data
!  write(6,*)'name',name

  do k=1,nl
!    write(6,*)'datain',minval(data3din(:,:,k)),maxval(data3din(:,:,k))
!    write(6,*)trim(name),k,data3din(nlon,nlat,k)
    do i=1,nlon
      do j=1,nlat
        data3d(j,i,k)=data3din(i,j,k)
      end do
    end do
  end do
!  write(6,*)'shape data3d',shape(data3d),'data3din',shape(data3din)
!  write(6,*)'data3d output',maxval(data3d),minval(data3d)
  return
end subroutine readcmaq3d
subroutine scat2dcmaq(datacmaq,datanode)
!  use cmaqmod, only : i2d,j2d,scattervsend,scattervdisp,itotscatter
  use cmaqmod, only : i2d,j2d
!  use cmaqmod, only : gathervrecv,gathervdisp,itotgatherv
  use kinds, only : i_kind,r_kind
  use mpimod, only : npe,mpi_comm_world,mpi_rtype,mype
  use gridmod
  implicit none
  integer(i_kind) ierr,nccmaq,nrcmaq,i,j,ii
  real(r_kind) send2d(itotsub)
! note datacmaq flipped already to be lat,lon
  real(r_kind) datanode(lat2,lon2),datacmaq(nlat,nlon)
!  write(300+mype,*)mype,'top scat2dcmaq lat2',lat2,lon2,'lat1',lat1,lat2
!  call flush(300+mype)
!  write(300+mype,*)'itotsub',itotsub,'itotscatterv',itotscatterv
!  call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
  if(mype.eq.0)then
    do ii=1,itotsub
      i=i2d(ii) ! lon
      j=j2d(ii) ! lat
      send2d(ii)=datacmaq(j,i)
!      if(ii<1000)then
!        write300+mype,*)'i',i,j,'ii',ii,'send2d',send2d(ii)
!      endif
    end do
  endif
!  write(300+mype,*)mype,'do scatterv ',shape(datanode),'shapei2d',shape(i2d),'scattervsend',scattervsend(mype+1)
!  write(300+mype,*)'size',size(datanode)
!  call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
 
!  call mpi_scatterv(send2d,ijn_s,displs_s,mpi_rtype, &
  call mpi_scatterv(send2d,scattervsend,scattervdisp,mpi_rtype, &
                     datanode ,scattervsend(mype+1),mpi_rtype,0,mpi_comm_world,ierr)
!  write(300+mype,*)mype,'ddid scatterv ',shape(datanode),'shapei2d',shape(i2d)
!  call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
  return
end subroutine scat2dcmaq
subroutine scat3dcmaq(datacmaq,datanode)
!  use cmaqmod, only : i2d,j2d,scattervsend,scattervdisp
  use cmaqmod, only : i2d,j2d
  use kinds, only : i_kind,r_kind
  use mpimod, only : npe,mpi_comm_world,mpi_rtype,mype
  use gridmod
  implicit none
!  integer(i_kind) i2d(itotsub),j2d(itotsub),ierr,nccmaq,nrcmaq,i,j,ii,k
  integer(i_kind) ierr,nccmaq,nrcmaq,i,j,ii,k
  integer(i_kind) isend3d(npe),disp3d(npe)
  real(r_kind) send3d(nsig,itotsub)
! note datacmaq flipped already to be lat,lon
  real(r_kind),dimension(grd_a%lat2,grd_a%lon2,grd_a%nsig) :: datanode
  real(r_kind) datacmaq(nlat,nlon,grd_a%nsig)
!  real(r_kind), dimension(grd%lat2,grd%lon2,grd%nsig), intent(out) :: g_u,g_v, g_q,g_oz,g_tv
  real(r_kind) datarecv(grd_a%nsig,grd_a%lat2,grd_a%lon2)
!  isend3d=grd_a%nsig*ijn_s
!  disp3d=grd_a%nsig*displs_s

  isend3d=grd_a%nsig*scattervsend
  disp3d=grd_a%nsig*scattervdisp
!  write(300+mype,*)mype,'top scat3dcmaq nsig',grd_a%nsig
!  call flush(300+mype)
!  write(300+mype,*)mype,nsig*lon2*lat2,'datarecv',shape(datarecv),'isend3d',isend3d(mype+1)
!  call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
  if(mype.eq.0)then
    do ii=1,itotsub
      i=i2d(ii) ! lon
      j=j2d(ii) ! lat
      send3d(:,ii)=datacmaq(j,i,:)
    end do
  endif
!  write(300+mype,*)mype,'did fill',isend3d(mype+1),disp3d(mype+1),'shape',shape(datarecv),'send',shape(send3d)
!  write(300+mype,*)'isend3d',isend3d(mype+1),'size',size(datarecv)
!  call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
 
  call mpi_scatterv(send3d,isend3d,disp3d,mpi_rtype, &
                     datarecv ,isend3d(mype+1),mpi_rtype,0,mpi_comm_world,ierr)
!  write(300+mype,*)mype,'did scatter'
!  call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
  do k=1,nsig
    do i=1,lon2
      do j=1,lat2
        datanode(j,i,k)=datarecv(k,j,i)
      end do
    end do
  end do
!  write(300+mype,*)'did node fill'
!  call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
  return
end subroutine scat3dcmaq
subroutine scat3dcmaqpr(datacmaq,datanode,cvar)
  use cmaqmod, only : i2d,j2d
!  scattervsend,scattervdisp
  use kinds, only : i_kind,r_kind
  use mpimod, only : npe,mpi_comm_world,mpi_rtype,mype
  use gridmod
  implicit none
  character *(*) cvar
!  integer(i_kind) i2d(itotsub),j2d(itotsub),ierr,nccmaq,nrcmaq,i,j,ii,k
  integer(i_kind) ierr,nccmaq,nrcmaq,i,j,ii,k
  integer(i_kind) isend3d(npe),disp3d(npe)
  real(r_kind) send3d(nsig,itotsub)
! note datacmaq flipped already to be lat,lon
  real(r_kind),dimension(grd_a%lat2,grd_a%lon2,grd_a%nsig) :: datanode
  real(r_kind) datacmaq(nlat,nlon,grd_a%nsig)
!  real(r_kind), dimension(grd%lat2,grd%lon2,grd%nsig), intent(out) :: g_u,g_v, g_q,g_oz,g_tv
  real(r_kind) datarecv(grd_a%nsig,grd_a%lat2,grd_a%lon2)
!  isend3d=grd_a%nsig*ijn_s
!  disp3d=grd_a%nsig*displs_s

  isend3d=grd_a%nsig*scattervsend
  disp3d=grd_a%nsig*scattervdisp
!  write(300+mype,*)mype,'top scat3dcmaq nsig',grd_a%nsig,cvar
!  call flush(300+mype)
!  !write(300+mype,*)mype,nsig*lon2*lat2,'datarecv',shape(datarecv),'isend3d',isend3d(mype+1),cvar
!  call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
  if(mype.eq.0)then
    do ii=1,itotsub
      i=i2d(ii) ! lon
      j=j2d(ii) ! lat
      send3d(:,ii)=datacmaq(j,i,:)
    end do
  endif
!  write(300+mype,*)mype,'did fill',isend3d(mype+1),disp3d(mype+1),'shape',shape(datarecv),'send',shape(send3d)
!  write(300+mype,*)'isend3d',isend3d(mype+1),'size',size(datarecv),cvar
!  call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
 
  call mpi_scatterv(send3d,isend3d,disp3d,mpi_rtype, &
                     datarecv ,isend3d(mype+1),mpi_rtype,0,mpi_comm_world,ierr)
!  write(300+mype,*)mype,'did scatter',cvar
!  call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
  do k=1,nsig
    do i=1,lon2
      do j=1,lat2
        datanode(j,i,k)=datarecv(k,j,i)
      end do
    end do
  end do
!  write(300+mype,*)'bottom scat3dmaqpr did node fill',cvar
!  call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
  return
end subroutine scat3dcmaqpr
subroutine scat3dpcmaq(datacmaq,datanode)
  use cmaqmod, only : i2d,j2d
! ,scattervsend,scattervdisp
  use kinds, only : i_kind,r_kind
  use mpimod, only : npe,mpi_comm_world,mpi_rtype,mype
  use gridmod
  implicit none
!  integer(i_kind) i2d(itotsub),j2d(itotsub),ierr,nccmaq,nrcmaq,i,j,ii,k
  integer(i_kind) ierr,nccmaq,nrcmaq,i,j,ii,k
  integer(i_kind) isend3d(npe),disp3d(npe)
  real(r_kind) send3d(nsig+1,itotsub)
! note datacmaq flipped already to be lat,lon
  real(r_kind),dimension(grd_a%lat2,grd_a%lon2,grd_a%nsig+1) :: datanode
  real(r_kind) datacmaq(nlat,nlon,grd_a%nsig+1)
!  real(r_kind), dimension(grd%lat2,grd%lon2,grd%nsig), intent(out) :: g_u,g_v, g_q,g_oz,g_tv
  real(r_kind) datarecv(grd_a%nsig+1,grd_a%lat2,grd_a%lon2)
!  isend3d=(grd_a%nsig+1)*ijn_s
  isend3d=(grd_a%nsig+1)*scattervsend
!  disp3d=(grd_a%nsig+1)*displs_s
  disp3d=(grd_a%nsig+1)*scattervdisp
!  write(300+mype,*)'top scat3dpcamq ',isend3d(mype+1)
!  call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
  if(mype.eq.0)then
    do ii=1,itotsub
      i=i2d(ii) ! lon
      j=j2d(ii) ! lat
      send3d(:,ii)=datacmaq(j,i,:)
    end do
  endif
!  write(300+mype,*)'fill send3d scat3dpcamq ',isend3d(mype+1),disp3d(mype+1)
!  !call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
 
  call mpi_scatterv(send3d,isend3d,disp3d,mpi_rtype, &
                     datarecv ,isend3d(mype+1),mpi_rtype,0,mpi_comm_world,ierr)
!  write(300+mype,*)'did mpi_scatterv',isend3d(mype+1),disp3d(mype+1)
!  write(300+mype,*)'size',size(datarecv)
!  call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
  do k=1,nsig+1
    do i=1,lon2
      do j=1,lat2
        datanode(j,i,k)=datarecv(k,j,i)
      end do
    end do
  end do
!  write(300+mype,*)'bottom scat3dcmaq'
!  call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
  return
end subroutine scat3dpcmaq
subroutine read_cmaq_guess_netcdf(mype)
  use kinds, only : i_kind,r_kind
  use gridmod, only : nlon,nlat,grd_a
  use guess_grids, only : ges_z,ges_ps,ges_u,ges_v,ges_tv,ges_q,ges_oz,nfldsig,ifilesig,ntguessig
  use guess_grids, only : ges_tsen
!  use guess_grids, only : ges_no2
  use guess_grids, only : ifilesfc
  use guess_grids, only : hrdifsig,hrdifsfc,hrdifsig_all,hrdifsfc_all
  use constants, only : zero
  use cmaqmod, only : cmaqpsol,ncidcmaqd
  use chemmod, only : laeroana_gocart
  implicit none
  integer(i_kind)mype,iret,i,k
  character(24)filename
  integer(i_kind) :: it,ncidcmaqdum
!  ntguessig=1
!  nfldsig=1
!  write(6,*)'readcmaq set nfldsig',nfldsig
!  call flush(6)
!  call create_gesfinfo
  hrdifsfc=zero
!  hrdifsig(1)=0.
!  hrdifsig_all(1)=hrdifsig(1)
  hrdifsfc_all(1)=hrdifsfc(1)
!  ifilesig(1)=0
!  write(6,*)mype,'nfldsig',nfldsig,'sz',shape(ges_z),'sp',shape(ges_ps)
!  call flush(6)
!  ifilesig=0
  ifilesfc=0
  ncidcmaqdum=0
!  if(mype==0) write(6,*)'find ajl read_cmaq ',(hrdifsig(i),i=1,nfldsig)
  
  do it=1,nfldsig
    write(filename,'("sigf",i2.2)')ifilesig(it)
!    if(mype.eq.0)write(6,*)'ifilesig',it,ifilesig(it)
    call get_cmaq_netcdf(grd_a,filename,mype,ges_z(1,1,it),ges_ps(1,1,it),ges_u(1,1,1,it), &
    ges_v(1,1,1,it),ges_tsen(1,1,1,it),ges_tv(1,1,1,it),ges_q(1,1,1,it),it,iret)
!    ges_v(1,1,1,it),ges_tv(1,1,1,it),ges_q(1,1,1,it),ges_oz(1,1,1,it),it,iret)
!    ges_v(1,1,1,it),ges_tv(1,1,1,it),ges_q(1,1,1,it),ges_oz(1,1,1,it),ges_no2(1,1,1,it),it,iret)
  end do
  return
end subroutine read_cmaq_guess_netcdf
subroutine get_cmaq_netcdf(grd,filename,mype,g_z,g_ps,g_u,g_v,g_tsen,g_tv,g_q,it,iret_read)
! subroutine to read in cmaq met and ozone now no2 also
! define dimension names
  use kinds, only : r_kind,i_kind
  use netcdf
  use mpimod, only : npe,mpi_comm_world,mpi_rtype,mpi_integer4
  use gridmod, only : nlon,nlat,nsig,itotsub,iglobal,jstart,istart,jlon1,ilat1,lat1,lon1
  use gridmod, only : ijn,displs_g,displs_s,ltosi,ltosj,strip
  use gridmod, only : ird_s,ijn_s,irc_s,isd_g,isc_g,lat2,lon2
  use guess_grids, only : nfldsig
  use constants, only : r1000,rd_over_cp
  use general_sub2grid_mod, only : sub2grid_info
  use constants, only : one,fv,qmin,zero
!  use cmaqmod, only :i2d,j2d,cmaqpdash,cmaqpbar,i2dg,j2dg,szraq,cmaqoz,cmaqpsol
  use cmaqmod, only :i2d,j2d,cmaqpdash,cmaqpbar,i2dg,j2dg,szcmaq,cmaqpsol
  use cmaqmod, only : cmaqpdashbar
!  use cmaqmod, only : ncidcmaqd,nameoz
  use cmaqmod, only : ncidcmaqd,ncidcmaq
!  use cmaqmod, only : nameno2,cmaqno2,itotscatter
  use cmaqmod, only : nameno2,cmaqno2
  use gridmod, only :itotscatterv
  use gridmod, only : scattervsend,scattervdisp
  use gridmod, only : gathervrecv,gathervdisp,itotgatherv
!  use cmaqmod, only : aodobsgrd,wtaodobsgrd
  use ozinfo, only : ihave_oz
  use no2info, only : ihave_no2
  use coinfo, only : ihave_co
  implicit none
! Declare passed variables
  type(sub2grid_info)                          ,intent(in) :: grd
  character(*), intent(in) :: filename
  integer(i_kind),intent(in) :: mype
! these are local grid on nodes
  real(r_kind), dimension(grd%lat2,grd%lon2) , intent(out) :: g_z,g_ps
  real(r_kind), dimension(grd%lat2,grd%lon2,grd%nsig), intent(out) :: g_u,g_v, g_q,g_tv,g_tsen
!  real(r_kind), dimension(grd%lat2,grd%lon2,grd%nsig), intent(out) :: g_no2
  integer(i_kind) iret_read
! declare local variables
  integer(i_kind) ndim,isum1,isum2,ibnd,jbnd,ip,ic,ii,kk,it
  parameter (ndim=4)
  character*12 dimname(ndim)
!  data dimname/'lon','lat','lev','time'/
  data dimname/'west-east','south-north','lev','time'/
!  character *6 namesp,namesz,nameu,namev,nameq,nameptrop
  character *6 namesp,namesz,nameu,namev,nameq
  data namesp/'PSFC'/,namesz/'SZ'/,nameu/'U'/,namev/'V'/,nameq/'QV'/
!  data nameptrop/'ptrop'/
  integer(i_kind) dimid(ndim),ierr,dimlen(ndim),nccmaq,nrcmaq,nlcmaq,i,j,k
  real(r_kind), allocatable :: pdash(:,:,:),u(:,:,:),v(:,:,:),q(:,:,:),oz(:,:,:)
  real(r_kind), allocatable :: temp(:,:,:)
!  real(r_kind), allocatable :: no2(:,:,:)
!  real(r_kind), allocatable :: co(:,:,:)
  real(r_kind), allocatable :: pdashin(:,:,:),psgather(:)
!  real(r_kind), allocatable :: delp(:,:,:),psol(:,:,:)
  real(r_kind), allocatable :: psol(:,:,:)
!  real(r_kind), allocatable :: sp(:,:),sz(:,:),ptrop(:,:)
  real(r_kind), allocatable :: sp(:,:),sz(:,:)
!  integer(i_kind) iu,iv,ioz,ipdash,idim,isz,isp,iq,idelp
  integer(i_kind) iu,iv,ioz,ipdash,idim,isz,isp,iq
  real(r_kind), allocatable :: tempps(:,:),tempsz(:,:),tempu(:,:,:),tempv(:,:,:),tempoz(:,:,:),t(:,:,:)
!  real(r_kind), allocatable :: tempno2(:,:,:)
!  real(r_kind), allocatable :: tempco(:,:,:)
  real(r_kind), allocatable :: tempq(:,:,:),tv(:,:,:)
  integer(i_kind), allocatable :: iremap(:,:,:),jremap(:,:,:)
  real(r_kind), allocatable :: send2d(:),send3d(:,:)
  real(r_kind), allocatable :: recv2d(:),recv3d(:,:,:),ps11(:,:)
  real(r_kind) :: div
  integer(i_kind) :: idateid,numbergridb,numbergridg,isumpoints
  integer *4 idate,idate2,ipnt
  character *10 cdate
!  integer(i_kind), allocatable :: i2d(:),j2d(:)
!  real(r_kind), dimension(grd%lat2,grd%lon2) , intent(out) :: g_z,g_ps
!  write(6,*)'mype',mype,'shape g_x',shape(g_z),'ps',shape(g_ps),'gridlat2',grd%lat2,grd%lon2
!  call flush(6)
  go to 888
!  if(.not.allocated(scattervsend))then
!    allocate (scattervsend(0:npe-1),scattervdisp(0:npe-1))
!    allocate (gathervrecv(0:npe-1),gathervdisp(0:npe-1))
    numbergridb=lat2*lon2
    numbergridg=lat1*lon1
!    write(300+mype,*)'lat2',lat2,lon2,'prod',lat2*lon2,'lat1',lat1,lat2,lat1*lat2
!    write(300+mype,*)'ilat1',ilat1(mype+1),jlon1(mype+1),'NUMBER',numbergridb
    call mpi_gather(numbergridb,1,mpi_integer4,scattervsend,1,mpi_integer4,0,mpi_comm_world,ierr)
    call mpi_gather(numbergridg,1,mpi_integer4,gathervrecv,1,mpi_integer4,0,mpi_comm_world,ierr)
    if(mype.eq.0)then
      itotscatterv=0
      itotgatherv=0
      scattervdisp(0)=0
      gathervdisp(0)=0
      itotscatterv=scattervsend(npe-1)
      itotgatherv=gathervrecv(npe-1)
      do i=0,npe-2
!        write(6,*)'scattervsend',i,scattervsend(i),'disp',scattervdisp(i)
!        write(6,*)'gathervrecv',i,gathervrecv(i),'disp',gathervdisp(i)
        scattervdisp(i+1)=scattervdisp(i)+scattervsend(i)
        gathervdisp(i+1)=gathervdisp(i)+gathervrecv(i)
        itotscatterv=itotscatterv+scattervsend(i)
        itotgatherv=itotgatherv+gathervrecv(i)
      end do
      write(6,*)'scattervsend',npe-1,scattervsend(npe-1),'disp',scattervdisp(npe-1)
      write(6,*)'gathervrecv',npe-1,gathervrecv(npe-1),'disp',gathervdisp(npe-1)
      write(6,*)'itotsub',itotsub,'itotscatterv',itotscatterv,'itotgatherv',itotgatherv
    endif
    call mpi_bcast(scattervdisp,npe,mpi_integer4,0,mpi_comm_world,ierr)
    call mpi_bcast(scattervsend,npe,mpi_integer4,0,mpi_comm_world,ierr)
    call mpi_bcast(itotscatterv,1,mpi_integer4,0,mpi_comm_world,ierr)
    call mpi_bcast(gathervdisp,npe,mpi_integer4,0,mpi_comm_world,ierr)
    call mpi_bcast(gathervrecv,npe,mpi_integer4,0,mpi_comm_world,ierr)
    call mpi_bcast(itotgatherv,1,mpi_integer4,0,mpi_comm_world,ierr)
!  endif
888 continue
!  call mpi_barrier(mpi_comm_world,ierr)
  nccmaq=nlon
  nrcmaq=nlat
  nlcmaq=nsig
!  write(6,*)mype,'lat1',lat1,lon1,'lat2',lat2,lon2,'istart',istart(mype+1),jstart(mype+1)
!  call flush(6)
  if(mype.eq.0)then
!    write(6,*)'find ajl getcmaqnetcdf filename',filename,'it',it
!    call flush(6)
!    do i=1,npe
!      write(6,*)i,'istart',istart(i),jstart(i),'ilat1',ilat1(i),jlon1(i)
!    end do
    if(.not.allocated(ncidcmaqd))then
      allocate(ncidcmaqd(nfldsig))
      ncidcmaqd=0
    endif
    if(it.eq.nfldsig)then
      ierr=nf90_open(filename,NF90_WRITE,ncidcmaqd(it))
!      write(6,*)'find ajl ncidcmaq it 2 ',it,ncidcmaqd(it),'file',trim(filename)
    else
      ierr=nf90_open(filename,NF90_NOWRITE,ncidcmaqd(it))
!      write(6,*)'find ajl ncidcmaq it 1 ',it,ncidcmaqd(it),'file',trim(filename)
    endif
    ncidcmaq=ncidcmaqd(it)
    do idim=1,ndim
      ierr=nf90_inq_dimid(ncidcmaq,dimname(idim),dimid(idim))
      if(ierr /= NF90_NOERR)then
        write(6,*)'can t get id for dim ',idim,dimname(idim)
        call stop2(9999)
        
      endif
      ierr=nf90_inquire_dimension(ncidcmaq,dimid(idim),len=dimlen(idim))
      if(ierr /= NF90_NOERR)then
        write(6,*)'can t get dimlen  ',idim,dimname(idim)
        call stop2(9999)
      endif
    end do
!   check if match gsi dimensions
    if(dimlen(1).ne.nlon.or.dimlen(2).ne.nlat.or.dimlen(3).ne.nsig)then
      write(6,*)'dimension mismatch dimlen',dimlen
      write(6,*)'nlon',nlon,'nlat',nlat,'nsig',nsig
      call stop2(999)
    endif
!    ierr=nf90_inq_varid(ncidcmaq,'IDATE',idateid)
!    if(ierr /= nf90_noerr)then
!      write(6,*)'can t get idateid  ',idateid,dimname(idim)
!      call stop2(9999)
!    endif
    idate2=2015123112
!    ierr=nf90_get_var(ncidcmaq,idateid,idate)
!    ierr=nf90_get_att(ncidcmaq,NF90_GLOBAL,'CDATE',cdate)
!    write(6,*)'cdate',cdate
    
!    write(6,*)'idate',idate,idate2
!   get surface variables
    allocate (sp(nrcmaq,nccmaq),sz(nrcmaq,nccmaq))
!    allocate (ptrop(nrcmaq,nccmaq))
    call readcmaq2d(sp,namesp)
!   sp is in pascals convert to mb
    call readcmaq2d(sz,namesz)
!    call readcmaq2d(ptrop,nameptrop)
!    write(6,*)'did sz'
!    call flush(6)
!    do j=1,10
!      write(6,*)'spin j',j,(sp(j,i),i=1,5)
!    end do
!   calculate temperature
!   need to read in theta and convert to temperature
!   first need to read in pressure on model layers
!    allocate(theta(nrcmaq,nccmaq,nlcmaq),t(nrcmaq,nccmaq,nlcmaq),tv(nrcmaq,nccmaq,nlcmaq),q(nrcmaq,nccmaq,nlcmaq))
    allocate(t(nrcmaq,nccmaq,nlcmaq),tv(nrcmaq,nccmaq,nlcmaq),q(nrcmaq,nccmaq,nlcmaq))
!    allocate(delp(nrcmaq,nccmaq,nlcmaq),psol(nrcmaq,nccmaq,nlcmaq+1),pdash(nrcmaq,nccmaq,nlcmaq))
    allocate(psol(nrcmaq,nccmaq,nlcmaq+1),pdash(nrcmaq,nccmaq,nlcmaq))
!   #define DIAGOZONE
#ifdef DIAGOZONE
    if(.not.allocated(cmaqoz))then
      allocate(cmaqoz(nccmaq,nrcmaq,nlcmaq,nfldsig))
    endif
!    if(.not.allocated(cmaqno2))then
!      allocate(cmaqno2(nccmaq,nrcmaq,nlcmaq,nfldsig))
!    endif
#endif
!    write(6,*)'call readcmaq3d'
!    call flush(6)
!   note readcmaq3d puts 1 at bottom like gsi wants
    call readcmaq3d(pdash,nlcmaq,'PDASH')
    call readcmaq3d(psol(1,1,2),nlcmaq,'PSOL')
!    call readcmaq3d(delp,nlcmaq,'delp')
!    psol=0.0
    do i=1,nccmaq
      do j=1,nrcmaq
        psol(j,i,1)=sp(j,i)
        if(psol(j,i,1)<psol(j,i,2))then
           write(6,*)mype,'ps error ',j,i,psol(j,i,2),sp(j,i)
        endif
      end do
    end do
    sp=sp*.01
!   add to read cmaq temperature 6/7/2016
!   test if have temperature field
    ierr=nf90_inq_varid(ncidcmaq,'TA',ipnt)
!   have cmaq temperature use it instead of using theta to calculate temperature
!    if(mype.eq.0)then
!      write(6,*)'have cmaq temperature use it'
!      call flush(6)
!    endif
    allocate(temp(nrcmaq,nccmaq,nlcmaq))
    call readcmaq3d(temp,nlcmaq,'TA')
    do k=1,nlcmaq
      do i=1,nccmaq
        do j=1,nrcmaq
          t(j,i,k)=temp(j,i,k)
!          psol(j,i,k+1)=psol(j,i,k)-delp(j,i,k)
        end do
      end do
    end do
    deallocate (temp)
    if(ihave_oz)then
      allocate (oz(nrcmaq,nccmaq,nlcmaq))
    endif
    allocate (u(nrcmaq,nccmaq,nlcmaq),v(nrcmaq,nccmaq,nlcmaq))
    call readcmaq3d(q,nlcmaq,nameq)
!    call readcmaq3d(ncidcmaq,u,nlcmaq,nameu)
    !call readcmaq3d(ncidcmaq,v,nlcmaq,namev)
     u=0.
     v=0.
!    if(ihave_no2)then
!      allocate (no2(nrcmaq,nccmaq,nlcmaq))
!      call readcmaq3d(ncidcmaq,no2,nlcmaq,'NO2')
!    endif
    write(6,*)'haveoz',ihave_oz
    call flush(6)
!    if(ihave_oz)then
!      call readcmaq3d(ncidcmaq,oz,nlcmaq,nameoz)
!    endif
    if(it.ne.nfldsig)then
      ierr=nf90_close(ncidcmaq)
    endif
!    print *,'maxoz cmaq',maxval(oz),minval(oz)
!    if(ihave_no2.and.ihave_oz)then
!      do k=1,nlcmaq
!        do i=1,nccmaq
!          do j=1,nrcmaq
!            q(j,i,k)=max(qmin,q(j,i,k))
!            tv(j,i,k)=t(j,i,k)*(one+fv*q(j,i,k))
!            !oz(j,i,k)=oz(j,i,k)*(47.98/28.97) ! convert vmr to mixing ratio !  close to 1/constoz time 1.e6
!            no2(j,i,k)=no2(j,i,k)*(46.0055/28.97) ! convert vmr to no2 mixing ratio !  
!#ifdef DIAGOZONE
!            cmaqno2(i,j,k,it)=no2(j,i,k)
!            cmaqoz(i,j,k,it)=oz(j,i,k)
!#endif
!          end do
!        end do
!      end do
!    elseif(ihave_no2)then
!      do k=1,nlcmaq
!        do i=1,nccmaq
!          do j=1,nrcmaq
!            q(j,i,k)=max(qmin,q(j,i,k))
!            tv(j,i,k)=t(j,i,k)*(one+fv*q(j,i,k))
!            oz(j,i,k)=oz(j,i,k)*(47.98/28.97) ! convert vmr to mixing ratio !  close to 1/constoz time 1.e6
!            no2(j,i,k)=no2(j,i,k)*(46.0055/28.97) ! convert vmr to no2 mixing ratio !  
!#ifdef DIAGOZONE
!            cmaqno2(i,j,k,it)=no2(j,i,k)
!#endif
!          end do
!        end do
!      end do
!    else
      do k=1,nlcmaq
        do i=1,nccmaq
          do j=1,nrcmaq
            q(j,i,k)=max(qmin,q(j,i,k))
            tv(j,i,k)=t(j,i,k)*(one+fv*q(j,i,k))
            if(ihave_oz)then
             oz(j,i,k)=oz(j,i,k)*(47.98/28.97) ! convert vmr to mixing ratio !  close to 1/constoz time 1.e6
            endif
#ifdef DIAGOZONE
            cmaqoz(i,j,k,it)=oz(j,i,k)
#endif
          end do
        end do
      end do
!    endif
!    write(6,*)'q',maxval(q),'u',maxval(u),'oz',maxval(oz),'tv',maxval(tv),'t',maxval(t),'v',maxval(v)
    write(6,*)'tv',maxval(tv),minval(tv)
    write(6,*)'q',maxval(q),minval(q)
!    call flush(6)
  else
    allocate (psol(1,1,1),pdash(1,1,1))
    allocate (sp(1,1),sz(1,1),q(1,1,1),tv(1,1,1))
    if(ihave_oz)then
      allocate (oz(1,1,1))
    endif
    allocate (u(1,1,1),v(1,1,1))
!    allocate (no2(1,1,1))
!    allocate (ptrop(1,1))
  endif
! scatter values to other processors
! set up big array on masterproc and scatterv them
  isum1=0
  isum2=0 
  if(iglobal.ne.itotgatherv)then
    write(6,*)mype,'iglobal',iglobal,itotgatherv
    call flush(6)
  endif
  if(mype.eq.0)then
    if(itotscatterv>itotsub)then
      allocate (send2d(itotscatterv),send3d(nlcmaq,itotscatterv))
      if(.not.allocated(i2d))then
        allocate (i2d(itotscatterv),j2d(itotscatterv))
      endif
    else
      allocate (send2d(itotsub),send3d(nlcmaq,itotsub))
      if(.not.allocated(i2d))then
        allocate (i2d(itotsub),j2d(itotsub))
      endif
    endif
 
!    allocate (i2dg(iglobal),j2dg(iglobal))
    allocate (i2dg(itotgatherv),j2dg(itotgatherv))
    i2dg=-999
    j2dg=-999
  else
    allocate (send2d(1),send3d(1,1))
    if(.not.allocated(i2d))then
      allocate (i2d(1),j2d(1))
    endif
    allocate (i2dg(1),j2dg(1))
  endif
  if(.not.allocated(cmaqpdash))then
    allocate (cmaqpdash(grd%lat2,grd%lon2,grd%nsig,nfldsig),cmaqpsol(grd%lat2,grd%lon2,grd%nsig+1,nfldsig))
    allocate (cmaqpbar(grd%nsig+1,nfldsig),cmaqpdashbar(grd%nsig,nfldsig))
  endif
  if(.not.allocated(szcmaq))allocate (szcmaq(nlon,nlat))
!  if(.not.allocated(aodobsgrd))then
!    allocate (aodobsgrd(grd%lat2,grd%lon2),wtaodobsgrd(grd%lat2,grd%lon2))
!    aodobsgrd=-999.
!    wtaodobsgrd=0.0
!  endif
!  if(.not.allocated(cmaqptrop))allocate (cmaqptrop(grd%lat2,grd%lon2))
!   write(500+mype,*)'shape cmaqptrop',shape(cmaqptrop)
!  write(500+mype,*)mype
!  do i=1,npe
!    write(500+mype,*)i,ijn_s(i),displs_s(i)
!  end do
!  call flush(500+mype)
!   write(6,*)'before mype',mype
!   call flush(6)
!   write(300+mype,*)'before mype',mype
!   call flush(300+mype)
!   call mpi_barrier(mpi_comm_world,ierr)

  if(mype.eq.0)then
#ifdef DIAGCMAQ
    write(6,*)'itotsub',itotsub,'npe',npe
    call flush(6)
    write(6,*)'ilat1',shape(ilat1),shape(jlon1)
    write(6,*)'ird_s',shape(ird_s)
    call flush(6)
    write(6,*)'ijn_s',shape(ijn_s)
    call flush(6)
    write(6,*)'ijn',shape(ijn)
    call flush(6)
    do i=1,npe
      write(6,*)'i,ilat1',ilat1(i),'jlon1',jlon1(i),'irc_s',irc_s(i)
      write(6,*)'proid',ilat1(i)*jlon1(i)
      call flush(67)
      write(6,*)i,'ird_s',ird_s(i),'isd_g',isd_g(i)
      write(6,*)i,'ijn_s',ijn_s(i),'ijn',ijn(i)
      isum1=isum1+ilat1(i)*jlon1(i)
      isum2=isum2+ijn_s(i)
      if(ilat1(i)*jlon1(i).ne.ijn(i))then
       write(6,*)'ne i',i
      endif
      write(6,*)'displs_s',displs_s(i),'g',displs_g(i)
      write(6,*)'isc_g',isc_g(i),'ltosi',ltosi(i),'ltosj',ltosj(i)
    end do
    write(6,*)'isum1',isum1,isum2,'itotsub',itotsub
#endif
    ic=0
    do ip=1,npe
#ifdef DIAGCMAQ
      write(6,*)'j limlat',istart(ip)-1,istart(ip)+ilat1(ip)
      call flush(6)
#endif
      do i=jstart(ip)-1,jstart(ip)+jlon1(ip)
        ibnd=i
!        if(i<1)ibnd=i+nccmaq
        if(i<1)ibnd=1
!        if(i>nccmaq)ibnd=i-nccmaq
        if(i>nccmaq)ibnd=nccmaq
        do j=istart(ip)-1,istart(ip)+ilat1(ip)
          jbnd=max(1,min(nrcmaq,j))
!         write(6,*)'j',j,'jbnd',jbnd,'limlon',jstart(ip)-1,jstart(ip)+jlon1(ip)
!         call flush(6)
          ic=ic+1
          if(ic>itotsub)then
            write(6,*)'ic out of bounds'
            call flush(6)
          endif
          i2d(ic)=ibnd ! lon
          j2d(ic)=jbnd ! lat
        end do
      end do
#ifdef DIAGCMAQ
      write(6,*)'ip',ip,'ic',ic,'nccmaq',nccmaq,nrcmaq
#endif
    end do
#ifdef DIAGCMAQ
      write(6,*)'ip',ip,'ic',ic,'nccmaq',nccmaq,nrcmaq
    do ic=1,100
      write(6,*)'ic',ic,'i2d',i2d(ic),'j2d',j2d(ic)
    end do
    call flush(6)
#endif
    ic=0
    isumpoints=0
    do ip=1,npe
      isumpoints=isumpoints+ilat1(ip)*jlon1(ip)
    end do
!    write(6,*)'isumpoints',isumpoints,'itotgatherv=',itotgatherv
!    call flush(6)
    do ip=1,npe
#ifdef DIAGCMAQ
      write(6,*)'ip',ip,'ic',ic,'nccmaq',nccmaq,'nrcmaq',nrcmaq
      write(6,*)'j limlat',istart(ip)-1,istart(ip)+ilat1(ip),'ilat1',ilat1(ip)
      write(6,*)'i limlon',jstart(ip)-1,jstart(ip)+jlon1(ip),'jlon1',jlon1(ip)
      write(6,*)'nlocal',ilat1(ip)*jlon1(ip)

      call flush(6)
#endif
      do i=jstart(ip),jstart(ip)+jlon1(ip)-1
        ibnd=i
!        if(i<1)ibnd=i+nccmaq
        if(i<1)ibnd=1
!        if(i>nccmaq)ibnd=i-nccmaq
        if(i>nccmaq)ibnd=nccmaq
        do j=istart(ip),istart(ip)+ilat1(ip)-1
          jbnd=max(1,min(nrcmaq,j))
!         write(6,*)'j',j,'jbnd',jbnd,'limlon',jstart(ip)-1,jstart(ip)+jlon1(ip)
!         call flush(6)
          ic=ic+1
!          if(ic>itotsub)then
!          if(ic>itotgatherv)then
!            write(6,*)ip,'ic out of bounds',ic,itotgatherv
!            call flush(6)
!          endif
          i2dg(ic)=ibnd ! lon
          j2dg(ic)=jbnd ! lat
!          if(ic<100)then
!            write(6,*)'ic',ic,'i',i,j,'ip',ip,'ibnd',ibnd,jbnd
!            call flush(6)
!          endif
        end do
      end do
!      write(6,*)'ic final',ic
#ifdef DIAGCMAQ
      write(6,*)'ip',ip,'icg',ic,'nccmaq',nccmaq,nrcmaq
#endif
    end do
#ifdef DIAGCMAQ
    write(6,*)'ic',ic,'iglobal',iglobal
    call flush(6)
#endif
        
          
!   ijn_s number to send displs_s displace for scatter for ilat2*jlon2
!   ijnm number to gather displs_g displace for gather for ilat1*jlon1
    send2d=-10.e20
    do ii=1,itotsub
      i=i2d(ii)
      j=j2d(ii)
      if(i<1.or.i>nccmaq.or.j<1.or.j>nrcmaq)then
         write(6,*)'ij out of bounds',i,j,nccmaq,nrcmaq
         call flush(6)
     endif
     end do
!    cmaq pressure in pascals
     sp=sp*.001 ! cb
   
     pdash=pdash*.001 ! cb
     psol=psol*.001 ! cb
     cmaqpbar(:,it)=0.0
     cmaqpdashbar(:,it)=0.0
     write(6,*)mype,'do pbar'
     call flush(6)
     do k=1,nsig+1
       do i=1,nlon
         do j=1,nlat
           cmaqpbar(k,it)=cmaqpbar(k,it)+psol(j,i,k)
!           if(i.eq.1.and.j.eq.1)print *,'cmaqpbar',k,cmaqpbar(k,it),'psol',psol(j,i,k)
         end do
       end do
     end do
     do k=1,nsig
       do i=1,nlon
         do j=1,nlat
           cmaqpdashbar(k,it)=cmaqpdashbar(k,it)+pdash(j,i,k)
         end do
       end do
     end do
!    write(6,*)mype,'did cmaqpbar'
!   call flush(6)
     div=one/float(nlon*nlat)
     cmaqpbar(:,it)=cmaqpbar(:,it)*div
     cmaqpdashbar(:,it)=cmaqpdashbar(:,it)*div
!    cb to mb so need to make mb
!    if(mype.eq.0)then
!       write(6,*)'cmaqpdashbar',cmaqpdashbar,' pa?'
!    endif
     cmaqpbar=cmaqpbar*10.
     cmaqpdashbar=cmaqpdashbar*10.
!     if(mype.eq.0)then
!         write(6,*)'cmaqpadhbar',cmaqpdashbar,' mb'
!      endif
!     do k=1,nsig+1
!       write(6,*)'cmaqpbar',k,cmaqpbar(k)
!     end do
!      do k=1,nsig
!        write(6,*)'cmaqpbar ',k,cmaqpdashbar(k,1)
!      end do
 
     
  endif
!  write(6,*)mype,'mpibcastpbsr'
!  call flush(6)
!  write(300+mype,*)mype,'mpibcastpbar'
!  call flush(300+mype)
  call mpi_barrier(mpi_comm_world,ierr)
  call mpi_bcast(cmaqpbar(1,it),nsig+1,mpi_rtype,0,mpi_comm_world,ierr)
  call mpi_bcast(cmaqpdashbar(1,it),nsig,mpi_rtype,0,mpi_comm_world,ierr)
     call mpi_barrier(mpi_comm_world,ierr)
!  if(mype.eq.1)then
!    do k=1,nsig+1
!       write(6,*)'cmaqpbar',k,cmaqpbar(k)
!       call flush(6)
!    end do
!  endif
!  write(6,*)'sp',shape(sp)
  !call flush(6)
!  write(6,*)'g_ps',shape(g_ps)
!  call flush(6)
!  write(300+mype,*)mype,'call scat2dsp'
!  call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
  call scat2dcmaq(sp,g_ps)
!  write(300+mype,*)mype,'did sp'
!  call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
 
!  write(300+mype,*)'sz',shape(sz)
!  call flush(300+mype)
!  write(300+mype,*)mype,'g_z',shape(g_z)
!  call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
  call scat2dcmaq(sz,g_z)
!  write(300+mype,*)mype,'did sz'
!  call flush(300+mype)
!  !call mpi_barrier(mpi_comm_world,ierr)
!  call scat2d(ptrop,cmaqptrop)
  call scat3dcmaq(tv,g_tv)
  call scat3dcmaq(t,g_tsen)
  call scat3dcmaq(q,g_q)
!  write(6,*)mype,'did q',maxval(g_q),minval(g_q)
!  write(6,*)mype,'did tv',maxval(g_tv),minval(g_tv)

!  call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
!  call scat3d(u,g_u)
!  call scat3d(v,g_v)
   g_u=0.0
   g_v=0.0
!   g_q=0.0
!  if(ihave_no2)then
!    call scat3d(no2,g_no2)
!  endif
!  write(300+mype,*)mype,'do pdash'
!  call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
! ajl pdash is in pascals convert to millibars
!  if(mype.eq.0)then
!    do k=1,25
!      print *,'befpdash',k,maxval(pdash(:,:,k)),'psol',maxval(psol(:,:,k))
!    end do
!  endif
  pdash=pdash
  call scat3dcmaq(pdash,cmaqpdash(1,1,1,it))
!  write(300+mype,*)mype,'did pdash'
  !call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
! ajl psol is in pascals convert to millibars
  psol=psol
  call scat3dpcmaq(psol,cmaqpsol(1,1,1,it))
!  if(mype.eq.0)then
!    do k=1,25
!      print *,'pdash',k,maxval(pdash(:,:,k)),'psol',maxval(psol(:,:,k))
!    end do
!  endif
  if(mype.eq.0)then
   szcmaq=sz
  endif
 ! write(6,*)'do bast',mype
 ! call flush(6)
  call mpi_bcast(szcmaq,nlon*nlat,mpi_rtype,0,mpi_comm_world,ierr)
#ifdef DIAGCMAQ
  write(6,*)mype,'maxps',maxval(g_ps),minval(g_ps)
  write(6,*)mype,'maxsz',maxval(g_z),minval(g_z)
  write(6,*)mype,'tv',maxval(g_tv),minval(g_tv)
  call flush(6)
#endif
! now gather to test transfer
  allocate (ps11(lat1,lon1))
  call strip(g_ps,ps11,1)
  if(mype.eq.0)then
    allocate(psgather(nlat*nlon))
#ifdef DIAGCMAQ
    do j=1,lat2
      write(6,*)'g_ps',j,(g_ps(j,i),i=1,5)
    end do
    do j=1,lat1
      write(6,*)'ps11',j,(ps11(j,i),i=1,5)
    end do
#endif
  else
    allocate(psgather(1))
  endif
  if(size(ps11).ne.gathervrecv(mype+1))then
     write(6,*)mype,'sizeps11',size(ps11),'gathervrecv',gathervrecv(mype+1),'disp',gathervdisp(mype+1)
  endif
! write(6,*)'gather',mype
! call flush(6)
!  call mpi_gatherv(ps11,ijn(mype+1),mpi_rtype,psgather,ijn,displs_g,mpi_rtype,0,mpi_comm_world,ierr)
  call mpi_gatherv(ps11,gathervrecv(mype+1),mpi_rtype,psgather,gathervrecv,gathervdisp,mpi_rtype,0,mpi_comm_world,ierr)
  if(mype.eq.0)then
!    write(6,*)'at 44'
!    call flush(6)
!    do ii=1,iglobal
    do ii=1,itotgatherv
!      i=ltosi(ii)
!      j=ltosj(ii)
      i=i2dg(ii)
      j=j2dg(ii)
      if(i<1.or.j<1)then
        call stop2(888)
      endif
      if(abs(sp(j,i)-psgather(ii))>1)then
        write(6,*)'i',i,j,'sp',sp(j,i),'psgather',psgather(ii),'ii',ii,'ltosi',ltosi(ii),ltosj(ii)
      endif
    end do
!    deallocate (delp,t)
    if(allocated(t))then
    deallocate (t)
    endif

  endif
  if(ihave_oz)then
    deallocate (oz)
  endif
  deallocate (psgather,sp,u,v,ps11)
  deallocate (sz,tv,q,psol,pdash)
!  deallocate(send2d,send3d,i2d,j2d,i2dg,j2dg)
!  deallocate(send2d,send3d,i2dg,j2dg)
  deallocate(send2d,send3d)
!  if(ihave_co)then
!    call get_cmaq_co(ncidcmaq)
!  endif
  if(ihave_no2)then
    call get_cmaq_no2
  endif
!  g_tv=0.0
  return
end subroutine get_cmaq_netcdf
subroutine get_cmaq_no2
  use kinds, only : i_kind,r_kind,r_single
  use gsi_chemguess_mod, only : GSI_ChemGuess_Bundle, gsi_chemguess_get
  use gsi_bundlemod, only: GSI_BundleGetPointer
  use mpimod,only : mype,mpi_comm_world
  use gridmod, only : nlon,nlat,nsig
!  use cmaqmod, only : ncidcmaqd,ncidcmaq
  use mpimod, only : mype
  implicit none
  integer(i_kind) :: ino2,ier,it,k,ierr,istatus
  integer(i_kind) :: ncidcmaq
  character(len=5), allocatable :: cvar(:)
  real(r_kind),pointer :: no2array(:,:,:)
  real(r_kind),allocatable :: no2tmp(:,:,:)
!  real(r_kind) :: mwco=28.0101,scale
  real(r_kind) :: mwno2=46.0055
  real(r_kind) :: nwair=28.97
! units of co in GSI is ppv
!  scale=mwco/28.97
!  scale=1.e6 ! ppv to ppmv
  it=1
!  write(6,*)mype,'call chemguess'
!  call flush(6)
  call gsi_chemguess_get ( 'var::no2',ino2,ier)
!  write(6,*)mype,'ino2',ino2,'nlat',nlat,nlon,nsig
!  call flush(6)
  allocate (cvar(ino2))
  if(mype.eq.0)then
!    write(6,*)'find nlat',nlat,nlon,nsig
    allocate(no2tmp(nlat,nlon,nsig))
  else
    allocate(no2tmp(1,1,1))
  endif
  if(mype.eq.0)then
!    write(6,*)'call readcmaq3d no2'
!    call flush(6)
    call readcmaq3d(no2tmp,nsig,'NO2')
!    force to be non negative
    no2tmp=max(no2tmp,1.e-18)
!    write(6,*)'did readcmaq3d no2',maxval(no2tmp),minval(no2tmp)
!    call flush(6)
    no2tmp=no2tmp*mwno2/nwair*1.e-6 ! ppmv to ppv to mxratio
!    cotmp=cotmp*scale
!    do k=1,nsig
!      write(6,*)'co in',k,maxval(cotmp(:,:,k)),minval(cotmp(:,:,K))
!    end do
  endif
!  write(6,*)mype,'did no2tmp'
!  call flush(6)
!  call mpi_barrier(mpi_comm_world,ierr)
  call gsi_chemguess_get('gsinames',cvar,ier)
!  write(6,*)mype,'cvar',cvar
!  write(300+mype,*)mype,'cvar',cvar
!  call flush(300+mype)
!  call flush(6)
!  call mpi_barrier(mpi_comm_world,ierr)
  call GSI_BundleGetPointer(GSI_ChemGuess_Bundle(it),cvar(1),no2array,istatus)
!  write(6,*)mype,'call scat3dpr no2tmp',size(no2array),'shape',shape(no2array)
!  write(300+mype,*)mype,'call scat3dpr no2tmp',size(no2array),'shape',shape(no2array)
!  call flush(300+mype)
!  call flush(6)
!  call mpi_barrier(mpi_comm_world,ierr)
  call scat3dcmaqpr(no2tmp,no2array,'NO2') 
!  write(300+mype,*)mype,'did call scat3d no2tmp'
!  call flush(300+mype)
!  write(6,*)mype,'did call scat3d no2tmp'
!  call flush(6) 
! call mpi_barrier(mpi_comm_world,ierr)
  deallocate (cvar,no2tmp)
!  write(6,*)mype,'bottom get_cmaq-no2'
!  write(300+mype,*)mype,'bottom get_cmaq-no2'
!  call flush(300+mype)
!  call mpi_barrier(mpi_comm_world,ierr)
  return  
end subroutine get_cmaq_no2
subroutine write_cmaq_netcdf
  use cmaqmod, only : cmaqno2,nameno2,ncidcmaqd,ncidcmaq
  use no2info, only : ihave_no2
  use netcdf
  use mpimod, only : mype
  use guess_grids, only : nfldsig
  implicit none
  integer ierr
!  if(mype.eq.0)then
!    ncidcmaq=ncidcmaqd(nfldsig)
!  endif
  if(mype.eq.0)then
     write(6,*)'write_cmaq_no2'
     call flush(6)
  endif

  if(ihave_no2)then
    call write_cmaQ_no2
!    write(6,*)mype,'did write_cmaq_no2'
!    call flush(6)
  endif
  if(mype.eq.0)then
    ierr=nf90_close(ncidcmaq)
    if(ierr /= NF90_NOERR)then
      write(6,*)'error close',nf90_strerror(ierr)
    endif
!    write(6,*)'did close'
!    call flush(6)
  endif
  return
end subroutine write_cmaq_netcdf
subroutine write_cmaq_no2
  use kinds, only : i_kind,r_kind,r_single
  use gsi_chemguess_mod, only : GSI_ChemGuess_Bundle, gsi_chemguess_get
  use gsi_bundlemod, only: GSI_BundleGetPointer
  use gridmod, only : nlon,nlat,nsig,gathervrecv,gathervdisp,lat1,lon1,iglobal,itotsub,strip
  use mpimod,only : mype
  use mpimod, only : mpi_comm_world,mpi_integer4,mpi_real4,mpi_rtype,mpi_real8
!  use gridmod, only : nlon,nlat,nsig
  use cmaqmod, only : ncidcmaqd,ncidcmaq
  use mpimod, only : mype
  use netcdf
  implicit none
  integer(i_kind) :: ino2,ier,it
  character(len=5), allocatable :: cvar(:)
  real(r_kind),pointer :: no2array(:,:,:)
  real(r_kind),allocatable :: no2tmp(:,:,:)
  real(r_kind), dimension(lat1*lon1,nsig) :: no2sm
  real(r_kind), dimension(max(iglobal,itotsub)) :: work
  real(r_kind), dimension(nlon,nlat) :: gridno2
  real(r_kind), dimension(nlon,nlat,nsig) :: newno2
!  real(r_kind) :: mwco=28.0101,scale
  real(r_kind) :: mwno2=46.0055
  real(r_kind) :: mwair=28.97
  integer(i_kind) :: k,kk,istatus,mm1,ierror,i,j,idno2,ierr
! units of no2 in GSI is ppmv
!  scale=mwco/28.97
!  scale=1.e-6 ! ppmv to ppv
  it=1
  mm1=mype+1
  call gsi_chemguess_get ( 'var::no2',ino2,ier)
  allocate (cvar(ino2))
  if(mype.eq.0)then
!    write(6,*)'find nlat',nlat,nlon,nsig
    allocate(no2tmp(nlat,nlon,nsig))
  else
    allocate(no2tmp(1,1,1))
  endif
  call gsi_chemguess_get('gsinames',cvar,ier)
  call GSI_BundleGetPointer(GSI_ChemGuess_Bundle(it),cvar(1),no2array,istatus)
!  write(6,*)mype,'call strip'
!  call flush(6)
  call strip(no2array,no2sm,nsig)
  do k=1,nsig
!    kk=nsig-k+1
    work=0.0
!    call mpi_gatherv(no2sm(1,k),ijn(mm1),mpi_rtype, &
!      work,ijn,displs_g,mpi_rtype, &
    call mpi_gatherv(no2sm(1,k),gathervrecv(mm1),mpi_rtype, &
      work,gathervrecv,gathervdisp,mpi_rtype, &
      0,mpi_comm_world,ierror)
    if(mype.eq.0)then
!      write(6,*)'call reordercmaq'
!      call flush(6)
      call reordercmaq(work,gridno2)
!      write(6,*)'did call reordercmaq'
!      call flush(6)
      do j=1,nlat
        do i=1,nlon
          newno2(i,j,k)=gridno2(i,j)*mwair/mwno2*1.e6 ! put back top ppmv

!             leave ppv
!          newco(i,j,kk)=newco(i,j,kk)*scale
        end do
      end do
    endif
  end do ! k
  if(mype.eq.0)then
   ierr=nf90_inq_varid(ncidcmaq,'NO2',idno2)
   if(ierr /= NF90_NOERR)then
     write(6,*)'error varid no2 ',nf90_strerror(ierr)
   endif
!   write(6,*)'did newno2','no2',maxval(newno2),minval(newno2)
!   call flush(6)
   ierr=nf90_put_var(ncidcmaq,idno2,newno2)
   if(ierr /= NF90_NOERR)then
     write(6,*)'error write no2',nf90_strerror(ierr)
   endif
 endif
  !call scat3d(cotmp,coarray)
  deallocate (cvar,no2tmp)
return
end subroutine write_cmaq_no2
subroutine reordercmaq(grid_in,grid_out)
use kinds, only: r_kind,i_kind
use gridmod, only : itotsub,ltosi,ltosj,nlat,nlon,iglobal,itotgatherv
use mpimod, only : mype
use cmaqmod, only : i2dg,j2dg
implicit none
real(r_kind), dimension(itotsub) , intent(in) :: grid_in ! input grid
real(r_kind), dimension(nlon,nlat),intent(out) :: grid_out ! output grid
integer(i_kind) i,j,k
!  Transfer input 1d array to output 2d array
!   write(6,*)'maxltosi',maxval(ltosi),'ltosj',maxval(ltosj)
!  note ltosi in north south and ltosj is east west
!  write(6,*)'reorderraqms work',mype,maxval(grid_in),minval(grid_in),' shape
!  ',shape(grid_in),'itotsub',itotsub,'iglobal',iglobal
!    write(6,*)mype,'reodercmaq',iglobal,'itotgatherv',itotgatherv
!    call flush(6)
!   write(6,*)'grid_in',shape(grid_in),'size',size(grid_in)
!   write(6,*)'grid_out',shape(grid_out),'size',size(grid_out)
!   write(6,*)'iglobal',iglobal
!   call flush(6)
   do k=1,iglobal
!      write(6,*)'k ',k,'i2dg',i2dg(k),j2dg(k)
!      call flush(6)
      i=i2dg(k)
      j=j2dg(k)
!      i=ltosi(k)
!      j=ltosj(k)
!      if(k<100)write(6,*)'k',k,'i',i,'j',j
      grid_out(i,j) = grid_in(k)
!      if(grid_in(k).ne.0)then
!        write(6,*)'grid_out',i,j,grid_out(j,i)
!      endif
   end do
!   write(6,*)'did reordercmaq'
!   call flush(6)
return
end subroutine reordercmaq
