subroutine getraqmsdate(filename)
  use raqmsmod, only : timeraq
  use kinds, only : r_kind,i_kind
  use mpimod, only : mpi_comm_world,mpi_rtype,mpi_integer4,mype
  implicit none
  include 'netcdf.inc'
  character *(*) filename,cdate*10
  integer(i_kind) :: ierr,ncidraq
  if(mype.eq.0)then
    write(6,*)'open filename ',trim(filename)
    call flush(6)
    ierr=nf_open(filename,NF_NOWRITE,ncidraq)
    if(ierr /= NF_NOERR)then
      write(6,*)'cant open file',trim(filename)
      call system('pwd ; ls -lt * ')
      call stop2(9999)
    endif
    ierr=nf_get_att_text(ncidraq,NF_GLOBAL,'CDATE',cdate)
    write(6,*)'cdate',cdate
    read(cdate,'(i4,i2,i2,i2)')timeraq(1),timeraq(2),timeraq(3),timeraq(4)
    timeraq(5)=0
    write(6,*)'timeraq',timeraq
    ierr=nf_close(ncidraq)
  endif
  call mpi_bcast(timeraq,5,mpi_integer4,0,mpi_comm_world,ierr)
  return
end subroutine getraqmsdate
subroutine raqmsset(mype)
  use kinds, only : i_kind,r_kind,r_single
  use guess_grids, only : nfldsig,ntguessfc,nfldsfc
  use obsmod, only: iadate,time_offset
  use constants, only: zero,one,r60inv
  use guess_grids, only : ifilesfc,ntguessig,ifilesig
  use gsi_4dvar, only: l4dvar, iwinbgn, winlen, nhr_assimilation
  use guess_grids, only : hrdifsig,hrdifsfc,hrdifsig_all,hrdifsfc_all,create_gesfinfo
  use mpimod, only : npe,ierror
  use raqmsmod, only : timeraq,wallstart
  use mpimod, only : mpi_comm_world,mpi_rtype,mpi_integer4
  use omimod,only : readapriori
  implicit none
  include 'netcdf.inc'
  integer mype
! declare local parameters
  real(r_kind),parameter:: r0_001=0.001_r_kind
! declare local variables
  integer(i_kind) iamana(3)
  integer(i_kind),parameter:: max_file = 100
  integer(i_kind),dimension(:),allocatable:: irec, fcst_hr_sig, &
     fcst_hr_sfc, fcst_hr_nst
  character(24)filename
  character(10)cdate
  logical(4) fexist
  integer(i_kind) :: i,j,iwan,npem1
  integer(i_kind) :: nhr_half,ierr
  integer(i_kind) :: nminanl,nmings,nming2,ndiff,ncidraq
  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
  real(r_kind),allocatable,dimension(:,:):: time_atm

! start read_raqms files here here.
  call getwall(wallstart)
  nhr_half=nhr_assimilation/2
  if(nhr_half*2 < nhr_assimilation) nhr_half=nhr_half+1
  npem1=npe-1
  if(mype.eq.0)print *,'raqmsset'
! find out number of uwh files
  nfldsig=0
  allocate( irec(max_file) )
  do i=0,min(99,max_file)
    write(filename,'("sigf",i2.2)')i
    inquire(file=filename,exist=fexist)
    if(fexist)then
      if(mype.eq.0)write(6,*)'ajl sigf ',filename,fexist
    endif
    if(fexist)then
      nfldsig=nfldsig+1 
      irec(nfldsig)=i
    endif
  end do
  if(mype.eq.0)write(6,*)'find ajl nfldsig',nfldsig
  allocate( fcst_hr_sig(nfldsig) )
  fcst_hr_sig(:) = irec(1:nfldsig)
  deallocate( irec )
  allocate(time_atm(nfldsig,2))
  if(mype.eq.0)then
!    write(6,*)'nfldsig ajl',nfldsig


!   convert analysis time to minutes relative to fixed date
    call w3fs21(iadate,nminanl)
!    write(6,*)'read_raqms_files:  analysis date,minutes ',iadate,nminanl
    iwan=0
    do i=1,nfldsig
      write(filename,'(''sigf'',i2.2)')fcst_hr_sig(i)
!     open uwh file and get date
      ierr=nf_open(filename,NF_NOWRITE,ncidraq)
      if(ierr /= NF_NOERR)then
         write(6,*)'cant open file',trim(filename)
         call stop2(9999)
      endif
      ierr=nf_get_att_text(ncidraq,NF_GLOBAL,'CDATE',cdate)
      ierr=nf_close(ncidraq)
!      write(6,*)'cdate',cdate
      read(cdate,'(i4,i2,i2,i2)')idateg(4),idateg(2),idateg(3),idateg(1)
      hourg4=zero
      hourg = hourg4
      idate5(1)=idateg(4); idate5(2)=idateg(2)
      idate5(3)=idateg(3); idate5(4)=idateg(1); idate5(5)=0
      if(mype.eq.0)write(6,*)'idate5 ajl',idate5
      call w3fs21(idate5,nmings)
      nming2=nmings+60*hourg
      write(6,*)'READ_RAQMS_FILES:  atm guess file, nming2 ',hourg,idateg,nming2
      t4dv=real((nming2-iwinbgn),r_kind)*r60inv
      if (l4dvar) then
        if (t4dv<zero .OR. t4dv>winlen) cycle
      else
        ndiff=nming2-nminanl
        if(mype.eq.0)write(6,*)'ndiff',ndiff,'nhr_half',nhr_half,'min',60*nhr_half
!       ajl want two have two time levels so change logic to test
         if(ndiff<0.0.and.abs(ndiff)>60*nhr_assimilation)cycle
         if(ndiff>0.0.and.abs(ndiff)>60*nhr_half)cycle
!        ajl below original
!        if(abs(ndiff) > 60*nhr_half ) cycle
       endif
       iwan=iwan+1
!       if(mype.eq.0)write(6,*)'nminanl',nminanl,'nming2',nming2,'iwan',iwan
!       if(mype.eq.0)write(6,*)'find ajl i',i,'ndiff',ndiff,'t4dv',t4dv
       if(nminanl==nming2) then
           iamana(1)=iwan
           timeraq=idate5
!           if(mype.eq.0)write(6,*)'set iamana ',iamana,'timeraq',timeraq
       endif
       time_atm(iwan,1) = t4dv
       time_atm(iwan,2) = fcst_hr_sig(i)+r0_001
!       write(6,*)'find ajl time_atm iwan ',iwan,time_atm(iwan,:)
    end do
    deallocate(fcst_hr_sig)
  endif
!  if(mype.eq.0)write(6,*)'timeraq',timeraq
!  call mpi_bcast(timeraq,5,mpi_integer4,0,mpi_comm_world,ierr)
!  write(6,*)mype,'timeraq',timeraq
  ntguessfc=1
  nfldsfc=1
! Broadcast guess file information to all tasks
  call mpi_bcast(time_atm,2*nfldsig,mpi_rtype,0,mpi_comm_world,ierror)
  call mpi_bcast(iamana,3,mpi_integer4,0,mpi_comm_world,ierror)
! Allocate space for guess information files

!  if(mype.eq.0)write(6,*)'create_gesfinfo'
  call create_gesfinfo
! Load time information for surface guess field info into output arrays
  ntguessfc = 1
  ntguessig = iamana(1)
!  write(6,*)mype,'ntguessig',ntguessig,'iamana',iamana,'nfldsig',nfldsig,'tim_atm',time_atm
!  call flush(6)
  do i=1,nfldsig
     hrdifsig(i) = time_atm(i,1)
     ifilesig(i) = nint(time_atm(i,2))
     hrdifsig_all(i) = hrdifsig(i)
!     write(100+mype,*)'find ajl raqmsset i',i,'hrdifsig',hrdifsig(i),'ifilesig',ifilesig(i)
!     if(mype.eq.0)write(6,*)'find ajl raqmsset i',i,'hrdifsig',hrdifsig(i),'ifilesig',ifilesig(i)
  end do
  if(mype.eq.0)write(6,*)'call readapriori'
  call readapriori

  return
end subroutine raqmsset
subroutine read_raqms(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_no2
  use guess_grids, only : ifilesfc
  use guess_grids, only : hrdifsig,hrdifsfc,hrdifsig_all,hrdifsfc_all
  use constants, only : zero
  use raqmsmod, only : raqpsol
  implicit none
  integer(i_kind)mype,iret,i
  character(24)filename
  integer(i_kind) :: it
!  ntguessig=1
!  nfldsig=1
!  write(6,*)'readraqms 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
!  if(mype==0) write(6,*)'find ajl read_raqms ',(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_raqms_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_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_raqms
subroutine read_raqmssfc(filename,mype,fact10,sfct,sno,isli,sfc_rough,terrain)
  use kinds, only : i_kind,r_kind,r_single
  use gridmod, only : nlat_sfc,nlon_sfc
  use raqmsmod, only : szraq
  use mpimod, only : mpi_comm_world,mpi_integer4,mpi_real4,mpi_real8
  implicit none
  character *(*), intent(in) :: filename
  integer(i_kind), intent(in) :: mype
  integer(i_kind), dimension(nlat_sfc,nlon_sfc),intent(out) :: isli
  real(r_kind), dimension(nlat_sfc,nlon_sfc),intent(out) :: fact10,sfct,sno,sfc_rough,terrain
! local r_single variable to read and bcast
  real(r_kind), dimension(nlon_sfc,nlat_sfc,4) :: data2din
  integer(i_kind), dimension(nlon_sfc,nlat_sfc) :: idata2din
  integer(i_kind) ncidsfc,ierr,ndim,idim,i,j
  include 'netcdf.inc'
  integer(i_kind) :: nsfc
  parameter (nsfc=5)
  character *10 csfc(nsfc)
  integer(i_kind) :: ivar(nsfc)
  data csfc/'F10M','TSEA','SHELEG','ZORL','SLMSK'/
  parameter (ndim=2)
  character *10 dimname(ndim)
  data dimname/'lon','lat'/
  integer(i_kind) :: dimlen(ndim),dimid(ndim),iv
! read
  if(mype.eq.0)then
!    write(6,*)'sfc filename',trim(filename)
    ierr=nf_open(filename,NF_NOWRITE,ncidsfc)
    if(ierr /= NF_NOERR)then
       write(6,*)'cant open file',trim(filename)
       call stop2(9999)
    endif
    do idim=1,ndim
      ierr=nf_inq_dimid(ncidsfc,dimname(idim),dimid(idim))
      if(ierr /= NF_NOERR)then
        write(6,*)'can t get id for dim ',idim,dimname(idim)
        call stop2(9999)
        
      endif
      ierr=nf_inq_dimlen(ncidsfc,dimid(idim),dimlen(idim))
      if(ierr /= NF_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_sfc.or.dimlen(2).ne.nlat_sfc)then
      write(6,*)'dimension mismatch dimlen',dimlen
      write(6,*)'nlon',nlon_sfc,'nlat_sfc',nlat_sfc
      call stop2(999)
    endif
    do iv=1,nsfc
      ierr=nf_inq_varid(ncidsfc,csfc(iv),ivar(iv))
    end do
  endif
  do iv=1,nsfc
    if(mype.eq.0)then
      if(csfc(iv).eq.'SLMSK')then
        ierr=nf_get_var_int(ncidsfc,ivar(iv),idata2din)
!         write(6,*)'mask',idata2din(1,1)

      else
        ierr=nf_get_var_double(ncidsfc,ivar(iv),data2din(1,1,iv))
!        write(6,*)csfc(iv),maxval(data2din(:,:,iv)),minval(data2din(:,:,iv))
      endif
    endif
    if(csfc(iv).eq.'SLMSK')then
      call mpi_bcast(idata2din,nlon_sfc*nlat_sfc,mpi_integer4,0,mpi_comm_world,ierr)
    else
      call mpi_bcast(data2din(1,1,iv),nlon_sfc*nlat_sfc,mpi_real8,0,mpi_comm_world,ierr)
    endif
  end do
  ierr=nf_close(ncidsfc)
!  data csfc/'F10M','TSEA','SHELEG','ZORL','SLMSK'/

  do i=1,nlon_sfc
    do j=1,nlat_sfc
      fact10(j,i)=data2din(i,j,1)
      sfct(j,i)=data2din(i,j,2)
      sno(j,i)=data2din(i,j,3)
      sfc_rough(j,i)=data2din(i,j,4)
      terrain(j,i)=szraq(i,j)
      isli(j,i)=idata2din(i,j)
    end do
  end do
  return
end subroutine read_raqmssfc


subroutine get_raqms_netcdf(grd,filename,mype,g_z,g_ps,g_u,g_v,g_tv,g_q,g_oz,g_no2,it,iret_read)
! subroutine to read in raqms met and ozone now no2 also
! define dimension names
  use kinds, only : r_kind,i_kind
  use mpimod, only : npe,mpi_comm_world,mpi_rtype
  use gridmod
  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 raqmsmod, only :i2d,j2d,raqpdash,raqpbar,i2dg,j2dg,szraq,raqmsoz,raqpsol
  use raqmsmod, only : ncidraqd,nameoz
  use raqmsmod, only : nameno2,raqmsno2,raqmsptrop
  use ozinfo, only : ihave_oz
  use no2info, only : ihave_no2
  implicit none
  include 'netcdf.inc'
! 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_oz,g_tv
  real(r_kind), dimension(grd%lat2,grd%lon2,grd%nsig), intent(out) :: g_no2
  integer(i_kind) iret_read,ncidraq
! declare local variables
  integer(i_kind) ndim,isum1,isum2,ibnd,jbnd,ip,ic,ii,kk,it
  parameter (ndim=4)
  character*4 dimname(ndim)
  data dimname/'lon','lat','lev','time'/
  character *6 namesp,namesz,nameu,namev,nameq,nameptrop
  data namesp/'psfc'/,namesz/'zsfc'/,nameu/'uu'/,namev/'vv'/,nameq/'iq'/
  data nameptrop/'ptrop'/
  integer(i_kind) dimid(ndim),ierr,dimlen(ndim),ncraq,nrraq,nlraq,i,j,k
  real(r_kind), allocatable :: pdash(:,:,:),theta(:,:,:),u(:,:,:),v(:,:,:),q(:,:,:),oz(:,:,:)
  real(r_kind), allocatable :: no2(:,:,:)
  real(r_kind), allocatable :: pdashin(:,:,:),psgather(:)
  real(r_kind), allocatable :: delp(:,:,:),psol(:,:,:)
  real(r_kind), allocatable :: sp(:,:),sz(:,:),ptrop(:,:)
  integer(i_kind) iu,iv,ioz,ipdash,itheta,idim,isz,isp,iq,idelp
  real(r_kind), allocatable :: tempps(:,:),tempsz(:,:),tempu(:,:,:),tempv(:,:,:),tempoz(:,:,:),t(:,:,:)
  real(r_kind), allocatable :: tempno2(:,:,:)
  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
  integer *4 idate,idate2
  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)
  if(mype.eq.0)then
    write(6,*)'find ajl getraqmsnetcdf filename',filename,'it',it
    call flush(6)
    if(.not.allocated(ncidraqd))then
      allocate(ncidraqd(nfldsig))
      ncidraqd=0
    endif
    if(it.eq.nfldsig)then
      ierr=nf_open(filename,NF_WRITE,ncidraqd(it))
      write(6,*)'find ajl ncidraq it 2 ',it,ncidraqd(it),'file',trim(filename)
    else
      ierr=nf_open(filename,NF_NOWRITE,ncidraqd(it))
      write(6,*)'find ajl ncidraq it 1 ',it,ncidraqd(it),'file',trim(filename)
    endif
    ncidraq=ncidraqd(it)
    do idim=1,ndim
      ierr=nf_inq_dimid(ncidraq,dimname(idim),dimid(idim))
      if(ierr /= NF_NOERR)then
        write(6,*)'can t get id for dim ',idim,dimname(idim)
        call stop2(9999)
        
      endif
      ierr=nf_inq_dimlen(ncidraq,dimid(idim),dimlen(idim))
      if(ierr /= NF_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
    ncraq=nlon
    nrraq=nlat
    nlraq=nsig
    ierr=nf_inq_varid(ncidraq,'IDATE',idateid)
    if(ierr /= NF_NOERR)then
      write(6,*)'can t get idateid  ',idateid,dimname(idim)
      call stop2(9999)
    endif
    idate2=2015123112
    ierr=nf_get_var_int(ncidraq,idateid,idate)
    ierr=nf_get_att_text(ncidraq,NF_GLOBAL,'CDATE',cdate)
!    write(6,*)'cdate',cdate
    
!    write(6,*)'idate',idate,idate2
!   get surface variables
    allocate (sp(nrraq,ncraq),sz(nrraq,ncraq))
    allocate (ptrop(nrraq,ncraq))
    call readraqms2d(ncidraq,sp,namesp)
    call readraqms2d(ncidraq,sz,namesz)
    call readraqms2d(ncidraq,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(nrraq,ncraq,nlraq),t(nrraq,ncraq,nlraq),tv(nrraq,ncraq,nlraq),q(nrraq,ncraq,nlraq))
    allocate(delp(nrraq,ncraq,nlraq),psol(nrraq,ncraq,nlraq+1),pdash(nrraq,ncraq,nlraq))
    if(.not.allocated(raqmsoz))then
      allocate(raqmsoz(ncraq,nrraq,nlraq,nfldsig))
    endif
    if(.not.allocated(raqmsno2))then
      allocate(raqmsno2(ncraq,nrraq,nlraq,nfldsig))
    endif
!    write(6,*)'call readraqms3d'
!    call flush(6)
!   note readraqms3d puts 1 at bottom like gsi wants
    call readraqms3d(ncidraq,pdash,nlraq,'pdash')
    call readraqms3d(ncidraq,delp,nlraq,'delp')
    call readraqms3d(ncidraq,theta,nlraq,'ttheta')
!    write(6,*)'ncraq',ncraq,nrraq,nlraq
    psol=0.0
    do i=1,ncraq
      do j=1,nrraq
        psol(j,i,1)=sp(j,i)
      end do
    end do
!    do k=1,nlraq
!      write(6,*)'pdash sp',pdash(1,1,k),'delp',delp(1,1,k)
!    end do
    do k=1,nlraq
!    do k=nlraq,1,-1
      do i=1,ncraq
        do j=1,nrraq
          t(j,i,k)=theta(j,i,k)*(pdash(j,i,k)/r1000)**rd_over_cp
          psol(j,i,k+1)=psol(j,i,k)-delp(j,i,k)
!          if(j.eq.1.and.i.eq.1)print *,'psol',k,psol(j,i,k:k+1),'delp',k,delp(j,i,k)
        end do
      end do
!      write(6,*)'temp',k,maxval(t(:,:,k)),minval(t(:,:,k)),'psol',maxval(psol(:,:,k+1))
      !call flush(6)
    end do
    deallocate (theta)
    allocate (oz(nrraq,ncraq,nlraq),u(nrraq,ncraq,nlraq),v(nrraq,ncraq,nlraq))
    allocate (no2(nrraq,ncraq,nlraq))
    call readraqms3d(ncidraq,q,nlraq,nameq)
    call readraqms3d(ncidraq,u,nlraq,nameu)
    call readraqms3d(ncidraq,v,nlraq,namev)
    if(ihave_no2)then
      call readraqms3d(ncidraq,no2,nlraq,'ino2')
    endif
    if(ihave_oz)then
      call readraqms3d(ncidraq,oz,nlraq,nameoz)
    endif
    if(it.ne.nfldsig)then
      ierr=nf_close(ncidraq)
    endif
!    print *,'maxoz raqms',maxval(oz),minval(oz)
    if(ihave_no2.and.ihave_oz)then
      do k=1,nlraq
        do i=1,ncraq
          do j=1,nrraq
            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 !  
            raqmsno2(i,j,k,it)=no2(j,i,k)
          end do
        end do
      end do
    elseif(ihave_no2)then
      do k=1,nlraq
        do i=1,ncraq
          do j=1,nrraq
            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 !  
            raqmsno2(i,j,k,it)=no2(j,i,k)
          end do
        end do
      end do
    else
      do k=1,nlraq
        do i=1,ncraq
          do j=1,nrraq
            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
            raqmsoz(i,j,k,it)=oz(j,i,k)
          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)
!    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))
    allocate (oz(1,1,1),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(mype.eq.0)then
    allocate (send2d(itotsub),send3d(nlraq,itotsub))
    allocate (i2d(itotsub),j2d(itotsub))
    allocate (i2dg(iglobal),j2dg(iglobal))
  else
    allocate (send2d(1),send3d(1,1))
    allocate (i2d(1),j2d(1))
    allocate (i2dg(1),j2dg(1))
  endif
  if(.not.allocated(raqpdash))then
    allocate (raqpdash(grd%lat2,grd%lon2,grd%nsig,nfldsig),raqpsol(grd%lat2,grd%lon2,grd%nsig+1,nfldsig))
    allocate (raqpbar(grd%nsig+1,nfldsig))
  endif
  if(.not.allocated(szraq))allocate (szraq(nlon,nlat))
  if(.not.allocated(raqmsptrop))allocate (raqmsptrop(grd%lat2,grd%lon2))
!   write(500+mype,*)'shape raqmsptrop',shape(raqmsptrop)
!  write(500+mype,*)mype
!  do i=1,npe
!    write(500+mype,*)i,ijn_s(i),displs_s(i)
!  end do
!  call flush(500+mype)

  if(mype.eq.0)then
#ifdef DIAGRAQ
    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 DIAGRAQ
      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+ncraq
        if(i>ncraq)ibnd=i-ncraq
        do j=istart(ip)-1,istart(ip)+ilat1(ip)
          jbnd=max(1,min(nrraq,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 DIAGRAQ
      write(6,*)'ip',ip,'ic',ic,'ncraq',ncraq,nrraq
#endif
    end do
#ifdef DIAGRAQ
    do ic=1,100
      write(6,*)'ic',ic,'i2d',i2d(ic),'j2d',j2d(ic)
    end do
    call flush(6)
#endif
    ic=0
    do ip=1,npe
#ifdef DIAGRAQ
      write(6,*)'j limlat',istart(ip)-1,istart(ip)+ilat1(ip)
      call flush(6)
#endif
      do i=jstart(ip),jstart(ip)+jlon1(ip)-1
        ibnd=i
        if(i<1)ibnd=i+ncraq
        if(i>ncraq)ibnd=i-ncraq
        do j=istart(ip),istart(ip)+ilat1(ip)-1
          jbnd=max(1,min(nrraq,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
          i2dg(ic)=ibnd ! lon
          j2dg(ic)=jbnd ! lat
        end do
      end do
#ifdef DIAGRAQ
      write(6,*)'ip',ip,'icg',ic,'ncraq',ncraq,nrraq
#endif
    end do
#ifdef DIAGRAQ
    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>ncraq.or.j<1.or.j>nrraq)then
         write(6,*)'ij out of bounds',i,j,ncraq,nrraq
         call flush(6)
     endif
     end do
     sp=sp*.1 ! cb
   
     pdash=pdash*.1 ! cb
     psol=psol*.1 ! cb
     raqpbar(:,it)=0.0
!     write(6,*)mype,'do pbar'
!     call flush(6)
     do k=1,nsig+1
       do i=1,nlon
         do j=1,nlat
           raqpbar(k,it)=raqpbar(k,it)+psol(j,i,k)
!           if(i.eq.1.and.j.eq.1)print *,'raqpbar',k,raqpbar(k,it),'psol',psol(j,i,k)
         end do
       end do
     end do
 !    write(6,*)mype,'did raqpbar'
 !    call flush(6)
     div=one/float(nlon*nlat)
     raqpbar(:,it)=raqpbar(:,it)*div
!     do k=1,nsig+1
!       write(6,*)'raqpbar',k,raqpbar(k)
!     end do
     
  endif
  call mpi_bcast(raqpbar(1,it),nsig+1,mpi_rtype,0,mpi_comm_world,ierr)
!  if(mype.eq.1)then
!    do k=1,nsig+1
!       write(6,*)'raqpbar',k,raqpbar(k)
!       call flush(6)
!    end do
!  endif
  call scat2d(sp,g_ps)
  call scat2d(sz,g_z)
  call scat2d(ptrop,raqmsptrop)
  call scat3d(tv,g_tv)
  call scat3d(q,g_q)
  call scat3d(u,g_u)
  call scat3d(v,g_v)
  if(ihave_oz)then
    call scat3d(oz,g_oz)
  endif
  if(ihave_no2)then
    call scat3d(no2,g_no2)
  endif
  call scat3d(pdash,raqpdash(1,1,1,it))
  call scat3dp(psol,raqpsol(1,1,1,it))
  if(mype.eq.0)then
   szraq=sz
  endif
 ! write(6,*)'do bast',mype
 ! call flush(6)
  call mpi_bcast(szraq,nlon*nlat,mpi_rtype,0,mpi_comm_world,ierr)
#ifdef DIAGRAQ
  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 DIAGRAQ
    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
! 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)
  if(mype.eq.0)then
    do ii=1,iglobal
!      i=ltosi(ii)
!      j=ltosj(ii)
      i=i2dg(ii)
      j=j2dg(ii)
      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)

  endif
  deallocate (psgather,sp,oz,u,v,ps11)
  deallocate (sz,tv,q,psol,pdash)
  deallocate(send2d,send3d,i2d,j2d,i2dg,j2dg)
  if(mype.eq.0)then
    write(6,*)'end get_raqms_netcdf'
    call flush(6)
  endif
  return
end subroutine get_raqms_netcdf
subroutine readraqms2d(ncidraq,data2d,name)
  use gridmod, only : nlon,nlat
  use mpimod, only : mype
  use kinds, only : r_kind,i_kind
  implicit none
  include 'netcdf.inc'
  character *(*) name
  real(r_kind) data2d(nlat,nlon),data2din(nlon,nlat)
  integer(i_kind) ierr,ipnt,ncidraq,i,j
  if(mype.ne.0)return
!  write(6,*)'readraqms2d',name
!  call flush(6)
  ierr=nf_inq_varid(ncidraq,name,ipnt)
  if(ierr /= NF_NOERR)then
    write(6,*)'readraqms2d var not found',trim(name)
  endif
  ierr=nf_get_var_double(ncidraq,ipnt,data2din)
  if(ierr /= NF_NOERR)then
    write(6,*)'error reading ',trim(name)
  endif
#ifdef DIAGRAQ
  write(6,*)name,'shape',shape(data2din),shape(data2d)
#endif
  call flush(6)
  do i=1,nlon
    do j=1,nlat
      data2d(j,i)=data2din(i,j)
    end do
  end do
#ifdef DIAGRAQ
  write(6,*)'did readraqms2d ',trim(name)
  call flush(6)
#endif
  return
end subroutine readraqms2d
subroutine readraqms3d(ncidraq,data3d,nl,name)
  use gridmod, only : nlon,nlat
  use mpimod, only : mype
  use kinds, only : r_kind,i_kind
  implicit none
  include 'netcdf.inc'
  character *(*) name
  real(r_kind) data3d(nlat,nlon,nl),data3din(nlon,nlat,nl)
  integer(i_kind) ierr,ipnt,ncidraq,nl,i,j,k
  if(mype.ne.0)return
#ifdef DIAGRAQ
  write(6,*)'readraqms3d',name
  call flush(6)
#endif
  ierr=nf_inq_varid(ncidraq,name,ipnt)
  if(ierr /= NF_NOERR)then
    write(6,*)'readraqms2d var not found',trim(name)
  endif
  ierr=nf_get_var_double(ncidraq,ipnt,data3din)
  if(ierr /= NF_NOERR)then
    write(6,*)'error reading ',trim(name)
  endif
! note 1 is bottom for global GSI
  do k=1,nl
!    write(6,*)trim(name),k,data3din(nlon,nlat,k)
    do i=1,nlon
      do j=1,nlat
        data3d(j,i,k)=data3din(i,j,nl-k+1)
      end do
    end do
  end do
  return
end subroutine readraqms3d
!subroutine scat2d(i2d,j2d,dataraq,datanode)
subroutine scat2d(dataraq,datanode)
  use raqmsmod, 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,ncraq,nrraq,i,j,ii
  integer(i_kind) ierr,ncraq,nrraq,i,j,ii
  real(r_kind) send2d(itotsub)
! note dataraq flipped already to be lat,lon
  real(r_kind) datanode(lat2,lon2),dataraq(nlat,nlon)
#ifdef DIAGRAQ
  write(6,*)mype,'top scat2d'
  call flush(6)
#endif
  if(mype.eq.0)then
#ifdef DIAGRAQ
    write(6,*)'nlon',nlon,nlat,'dataraq',shape(dataraq)
#endif
    do ii=1,itotsub
      i=i2d(ii) ! lon
      j=j2d(ii) ! lat
      send2d(ii)=dataraq(j,i)
!      if(ii<1000)then
!        write(6,*)'i',i,j,'ii',ii,'send2d',send2d(ii)
!      endif
    end do
  endif
 
  call mpi_scatterv(send2d,ijn_s,displs_s,mpi_rtype, &
                     datanode ,ijn_s(mype+1),mpi_rtype,0,mpi_comm_world,ierr)
  return
end subroutine scat2d
subroutine scat3d(dataraq,datanode)
  use raqmsmod, 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,ncraq,nrraq,i,j,ii,k
  integer(i_kind) ierr,ncraq,nrraq,i,j,ii,k
  integer(i_kind) isend3d(npe),disp3d(npe)
  real(r_kind) send3d(nsig,itotsub)
! note dataraq flipped already to be lat,lon
  real(r_kind),dimension(grd_a%lat2,grd_a%lon2,grd_a%nsig) :: datanode
  real(r_kind) dataraq(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
  if(mype.eq.0)then
    do ii=1,itotsub
      i=i2d(ii) ! lon
      j=j2d(ii) ! lat
      send3d(:,ii)=dataraq(j,i,:)
    end do
  endif
 
  call mpi_scatterv(send3d,isend3d,disp3d,mpi_rtype, &
                     datarecv ,isend3d(mype+1),mpi_rtype,0,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
  return
end subroutine scat3d
subroutine scat3dp(dataraq,datanode)
  use raqmsmod, 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,ncraq,nrraq,i,j,ii,k
  integer(i_kind) ierr,ncraq,nrraq,i,j,ii,k
  integer(i_kind) isend3d(npe),disp3d(npe)
  real(r_kind) send3d(nsig+1,itotsub)
! note dataraq flipped already to be lat,lon
  real(r_kind),dimension(grd_a%lat2,grd_a%lon2,grd_a%nsig+1) :: datanode
  real(r_kind) dataraq(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
  disp3d=(grd_a%nsig+1)*displs_s
  if(mype.eq.0)then
    do ii=1,itotsub
      i=i2d(ii) ! lon
      j=j2d(ii) ! lat
      send3d(:,ii)=dataraq(j,i,:)
    end do
  endif
 
  call mpi_scatterv(send3d,isend3d,disp3d,mpi_rtype, &
                     datarecv ,isend3d(mype+1),mpi_rtype,0,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
  return
end subroutine scat3dp
subroutine write_raqms(mype)
use raqmsmod, only : raqmsoz,ncidraqd,nameoz,wallstart,wallend
use raqmsmod, only : raqmsno2,nameno2
use guess_grids, only : ges_oz,nfldsig
use guess_grids, only : ges_no2
use gridmod,only : iglobal,itotsub,lat1,lon1,nsig,ijn,displs_g,nlon,nlat,strip
use kinds, only: r_kind,i_kind
use mpimod, only : mpi_comm_world,mpi_integer4,mpi_real4,mpi_rtype,mpi_real8
use ozinfo, only : ihave_oz
use no2info, only : ihave_no2
implicit none
integer mype,mm1
  include 'netcdf.inc'
  real(r_kind), dimension(lat1*lon1,nsig) :: ozsm,no2sm
  real(r_kind), dimension(max(iglobal,itotsub)) :: work,workno2
  real(r_kind), dimension(nlon,nlat) :: gridoz,gridno2
  real(r_kind), dimension(nlon,nlat,nsig) :: newoz,diffoz
  real(r_kind), dimension(nlon,nlat,nsig) :: newno2,diffno2
  integer(i_kind) k,ierror,i,j,idoz,ierr,kk,ncidraq,idno2
! need to gather ges_oz and difference with raqmsoz
!  write(100+mype,*)'ges_oz ajl',maxval(ges_oz),minval(ges_oz)
!  if(maxval(ges_oz(:,:,:,nfldsig))>0.0)write(6,*)mype,'maxges_oz',maxval(ges_oz)
  call strip(ges_oz(1,1,1,nfldsig),ozsm,nsig)
  call strip(ges_no2(1,1,1,nfldsig),no2sm,nsig)
!  write(6,*)mype,'ozsmmax',maxval(ozsm),minval(ozsm)
  mm1=mype+1
!          oz(j,i,k)=oz(j,i,k)*(47.98/28.97) ! convert vmr to mixing ratio
  do k=1,nsig
    kk=nsig-k+1
    work=0.0
    workno2=0.0
    if(ihave_oz)then
    call mpi_gatherv(ozsm(1,k),ijn(mm1),mpi_rtype, &
      work,ijn,displs_g,mpi_rtype, &
      0,mpi_comm_world,ierror)
    endif
    if(ihave_no2)then
      call mpi_gatherv(no2sm(1,k),ijn(mm1),mpi_rtype, &
      workno2,ijn,displs_g,mpi_rtype, &
      0,mpi_comm_world,ierror)
    endif
    if(mype.eq.0)then
!      write(6,*)'work ajl',k,maxval(work),minval(work)
!      write(6,*)'shape',shape(work),'gridoz',shape(gridoz)
      if(ihave_oz)then
        call reorderraqms(work,gridoz)
      endif
      if(ihave_no2)then
        call reorderraqms(workno2,gridno2)
      endif
!      write(6,*)'ajl gridoz',maxval(gridoz),minval(gridoz)
      if(ihave_no2)then
        do j=1,nlat
          do i=1,nlon
            newno2(i,j,kk)=gridno2(i,j)
            diffno2(i,j,k)=newoz(i,j,kk)-raqmsno2(i,j,k,nfldsig)
            newno2(i,j,kk)=newno2(i,j,kk)*28.97/46.0055
          end do
        end do
      endif
      if(ihave_oz)then
        do j=1,nlat
          do i=1,nlon
            newoz(i,j,kk)=gridoz(i,j)
            diffoz(i,j,k)=newoz(i,j,kk)-raqmsoz(i,j,k,nfldsig)
            newoz(i,j,kk)=newoz(i,j,kk)*28.97/47.98
          end do
        end do
      endif
!          if(newoz(i,j,k).ne.0.and.diffoz(i,j,k).ne.0)then
!             write(6,*)'diffoz',i,j,k,diffoz(i,j,k),'new',newoz(i,j,k),'raqmsoz',raqmsoz(i,j,k)
!          endif
!      write(6,*)'diffoz',k,maxval(diffoz(:,:,k)),minval(diffoz(:,:,k))
    endif
  end do
  if(mype.eq.0)then
    ncidraq=ncidraqd(nfldsig)
!    write(6,*)'ajl find nldsig',nfldsig,'ncidraq',ncidraq,'ncidraqd',ncidraqd
    if(ihave_oz)then
      ierr=nf_inq_varid(ncidraq,nameoz,idoz)
!      write(6,*)'idoz',idoz
      if(ierr /= NF_NOERR)then
        write(6,*)'error varid ',nameoz,nf_strerror(ierr)
      endif
      write(6,*)'newoz',maxval(newoz),minval(newoz)
      ierr=nf_put_var_double(ncidraq,idoz,newoz)
      if(ierr /= NF_NOERR)then
        write(6,*)'error write ',nameoz,nf_strerror(ierr)
      endif
    endif
    if(ihave_no2)then
      ierr=nf_inq_varid(ncidraq,nameno2,idno2)
!    write(6,*)'idoz',idoz
      if(ierr /= NF_NOERR)then
        write(6,*)'error varid ',nameno2,nf_strerror(ierr)
      endif
!      write(6,*)'newoz',maxval(newoz),minval(newoz)
      ierr=nf_put_var_double(ncidraq,idno2,newno2)
      if(ierr /= NF_NOERR)then
        write(6,*)'error write ',nameno2,nf_strerror(ierr)
      endif
      write(6,*)'newno2',maxval(newno2),minval(newno2)
      call flush(6)
     endif
    ierr=nf_close(ncidraq)
    if(ierr /= NF_NOERR)then
      write(6,*)'error close',nf_strerror(ierr)
    endif
!    write(6,*)'close netcdf file'
  endif
  call getwall(wallend)
  if(mype.eq.0)then
       write(6,*)'diff wall ',wallend-wallstart, 'min ',(wallend-wallstart)/60.
 endif
  
return
end subroutine write_raqms
subroutine reorderraqms(grid_in,grid_out)
use kinds, only: r_kind,i_kind
use gridmod, only : itotsub,ltosi,ltosj,nlat,nlon,iglobal
use mpimod, only : mype
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
   
   do k=1,iglobal
      i=ltosi(k)
      j=ltosj(k)
!      if(k<100)write(6,*)'k',k,'i',i,'j',j
      grid_out(j,i) = grid_in(k)
!      if(grid_in(k).ne.0)then
!        write(6,*)'grid_out',i,j,grid_out(j,i)
!      endif
   end do
return
end subroutine reorderraqms
      subroutine getwall(wall)
      use mpimod, only : mype
      implicit none 
      include 'mpif.h'
!      real wall,mpi_wtime
      real*8 wall 
      wall=mpi_wtime()
      return
      end  

subroutine intrp3ozdp1(f,g,dx,dy,obstime,dpoz,dp,sfcp,mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    intrp3oz    space-time linear interpolation for ozone
!   prgmmr: parrish          org: np22                date: 2013-01-26
! added dp and sfc press and dpoz(k) for omi apriori adjustment Lenzen 2015-05-10
!
! abstract: same as intrp3oz but for special case n=1, with n argument removed.
!            This has been created to solve problem of type mismatch debug
!            compile
!            error on WCOSS.
!
! program history log:
!   2013-01-26  parrish
!
!   input argument list:
!     f        - input interpolator (gridded guess ozone fields)
!     dx,dy - input x,y-coords of interpolation points (grid units)
!     obstime  - observation times
!     mype     - mpi task id
!
!   output argument list:
!     g        - output interpolatees (guess ozone at observation location)
!     dp       - delta pressure at observation location
!     dpoz     - delta pressure * layer ozone at observation location
!     sfcp     - surface pressure at observation location
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
!--------
  use kinds, only: r_kind,i_kind
  use guess_grids, only: nfldsig,hrdifsig,ges_prsi
  use gridmod, only: lat2,lon2,nlat,nlon,nsig,lon1,istart,jstart
  use constants, only: zero, one, rozcon
  implicit none

! Declare passed variables
  integer(i_kind)                               ,intent(in   ) :: mype
  real(r_kind)                                  ,intent(in   ) :: dx,dy,obstime
  real(r_kind),dimension(lat2,lon2,nsig,nfldsig),intent(in   ) :: f
  real(r_kind)                                  ,intent(  out) :: g
  real(r_kind),dimension(nsig)                 ,intent(out)    :: dp,dpoz
  real(r_kind)                                 ,intent(out)    :: sfcp
  

! Declare local variables
  integer(i_kind) j,k,ix,ix1,iy,iy1,kk,itsig,itsigp,iz1,iz2
  integer(i_kind) ixp,iyp,mm1
  real(r_kind) w00,w01,w10,w11,delx,dely,delx1,dely1
  real(r_kind) dtsig,dtsigp,pob
  real(r_kind) delp1,delp2,delp3,delp4,delp5,delp6,delp7,delp8


!*************************************************************************
! Initialize variables
  g=zero
  dp=zero
  dpoz=zero
  sfcp=zero
  mm1=mype+1


! Loop over number of observations.

!    Get horizontal interpolation information.  This information includes
!    the (i,j) indices of the grid points surrounding the observation,
!    plus the corresponding interpolation weights between these points
!    and the observation.

     ix1=dx; iy1=dy
     ix1=max(1,min(ix1,nlat))
     delx=dx-ix1; dely=dy-iy1; delx=max(zero,min(delx,one))
     ix=ix1-istart(mm1)+2; iy=iy1-jstart(mm1)+2
     if(iy<1) then
        iy1=iy1+nlon
        iy=iy1-jstart(mm1)+2
     end if
     if(iy>lon1+1) then
        iy1=iy1-nlon
        iy=iy1-jstart(mm1)+2
     end if
     ixp=ix+1; iyp=iy+1
     if(ix1==nlat) then
        ixp=ix
     end if
     delx1=one-delx; dely1=one-dely
     w00=delx1*dely1; w10=delx*dely1; w01=delx1*dely; w11=delx*dely
        
!    Set the weights for linear time iterpolation from the guess to 
!    the observation time.
!     if(mype.eq.0)then
!          write(6,*)'obstime',obstime,'hrdiffsig',hrdifsig,'nfldsig',nfldsig
!           call flush(6)
!     endif

     if(obstime > hrdifsig(1) .and. obstime < hrdifsig(nfldsig))then
        do j=1,nfldsig-1
           if(obstime > hrdifsig(j) .and. obstime <= hrdifsig(j+1))then
              itsig=j
              itsigp=j+1
              dtsig=((hrdifsig(j+1)-obstime)/(hrdifsig(j+1)-hrdifsig(j)))
           end if
        end do
     else if(obstime <=hrdifsig(1))then
        itsig=1
        itsigp=1
        dtsig=one
     else
        itsig=nfldsig
        itsigp=nfldsig
        dtsig=one
     end if
     dtsigp=one-dtsig
     
!
!    Perform spatial and temporal interpolation for the total column 
!    ozone observation
!     write(100+mype,*)'ix',ix,iy,itsig,itsigp,'prsi',ges_prsi(ix,iy,1,itsig),ges_prsi(ix,iy,nsig,itsig)
     do kk=1,nsig
        delp1=ges_prsi(ix ,iy ,kk,itsig )-ges_prsi(ix ,iy ,kk+1,itsig )
        delp2=ges_prsi(ixp,iy ,kk,itsig )-ges_prsi(ixp,iy ,kk+1,itsig )
        delp3=ges_prsi(ix ,iyp,kk,itsig )-ges_prsi(ix ,iyp,kk+1,itsig )
        delp4=ges_prsi(ixp,iyp,kk,itsig )-ges_prsi(ixp,iyp,kk+1,itsig )
        delp5=ges_prsi(ix ,iy ,kk,itsigp)-ges_prsi(ix ,iy ,kk+1,itsigp)
        delp6=ges_prsi(ixp,iy ,kk,itsigp)-ges_prsi(ixp,iy ,kk+1,itsigp)
        delp7=ges_prsi(ix ,iyp,kk,itsigp)-ges_prsi(ix ,iyp,kk+1,itsigp)
        delp8=ges_prsi(ixp,iyp,kk,itsigp)-ges_prsi(ixp,iyp,kk+1,itsigp)
        dp(kk)=dp(kk)+(w00*delp1+w10*delp2+w01*delp3+w11*delp4)*dtsig + &
                      (w00*delp5+w10*delp6+w01*delp7+w11*delp8)*dtsigp
        dpoz(kk)=dpoz(kk)+ &
              (f(ix ,iy ,kk,itsig )*w00*delp1 &
             + f(ixp,iy ,kk,itsig )*w10*delp2 &
             + f(ix ,iyp,kk,itsig )*w01*delp3 &
             + f(ixp,iyp,kk,itsig )*w11*delp4)*dtsig + &
              (f(ix ,iy ,kk,itsigp)*w00*delp5 &
             + f(ixp,iy ,kk,itsigp)*w10*delp6 &
             + f(ix ,iyp,kk,itsigp)*w01*delp7 &
             + f(ixp,iyp,kk,itsigp)*w11*delp8)*dtsigp
        dpoz(kk)=dpoz(kk)*rozcon
!        write(100+mype,*)kk,'dopoz',dpoz(kk),'rozcon',rozcon
        g=g + dpoz(kk)
!        g=g + &
!              (f(ix ,iy ,kk,itsig )*w00*delp1 &
!             + f(ixp,iy ,kk,itsig )*w10*delp2 &
!             + f(ix ,iyp,kk,itsig )*w01*delp3 &
!             + f(ixp,iyp,kk,itsig )*w11*delp4)*dtsig + &
!              (f(ix ,iy ,kk,itsigp)*w00*delp5 &
!             + f(ixp,iy ,kk,itsigp)*w10*delp6 &
!             + f(ix ,iyp,kk,itsigp)*w01*delp7 &
!             + f(ixp,iyp,kk,itsigp)*w11*delp8)*dtsigp
     enddo
     sfcp=(ges_prsi(ix,iy,1,itsig)*w00+ges_prsi(ixp,iy,1,itsig)*w10 + &
           ges_prsi(ix,iyp,1,itsig)*w01+ges_prsi(ixp,iyp,1,itsig)*w11)*dtsig + &
          (ges_prsi(ix,iy,1,itsigp)*w00+ges_prsi(ixp,iy,1,itsigp)*w10 + &
           ges_prsi(ix,iyp,1,itsigp)*w01+ges_prsi(ixp,iyp,1,itsigp)*w11)*dtsigp 
!     write(100+mype,*)'sfcp',sfcp,dtsig,dtsigp
!     call flush(100+mype)

! End of routine
  return
end subroutine intrp3ozdp1
subroutine omiadjust(xdata1,rrlat,dpoz,dp,sfcp)
  use kinds, only: r_kind,i_kind
  use omimod, only : ncol,nprs,apreff,aprprb,aprprt,aprprs,aprlat,aprcol,aprmon,nmon
  use omimod, only : nlat,apriori,pintapr,aprdp,mapriori
  use obsmod, only: iadate,ianldate
  use gridmod, only : nsig
  use mpimod, only : mype
  implicit none
  real(r_kind),dimension(nsig) :: dpoz,dp,eff,apr,psav,ydata2
  real(r_kind),dimension(nsig+1) :: presi
  real(r_kind) :: apri0,apri1,xdata1,rrlat,sfcp,frac2,frac,frac1,xdata2
  real(r_kind) :: mindex0,mindex1,xscale,yscale,cscale,spob,apr0,apr1,mscale,pscale
  real(r_kind) :: yapriori(ncol,nprs)
  integer(i_kind) :: iday,ifailed,yindex0,yindex1,cindex0,cindex1,ii,jj,ll,llev,iaprlev
  integer(i_kind) :: imonth,ihr,iyear,idat,k,kbot
  integer(i_kind),dimension(nsig) :: filled
  real(r_kind) :: intapr(0:nprs+1),praqtop,praqbot,afbot,aprintbot,aftop,aprinttop
  logical :: gotbot
! construct TOMS apriori based on RAQMS total column
!  if(mype.eq.0)write(6,*)'top omiajust',xdata1
  presi(1)=sfcp*10. ! make mb
  do k=1,nsig
    dp(k)=dp(k)*10. ! makd mb
    presi(k+1)=presi(k)-dp(k)
    psav(k)=.5*(presi(k)+presi(k+1))
!    write(100+mype,*)k,dp(k),'psav',psav(k)
  end do
  if(rrlat.lt.aprlat(1)) then
    yindex0=1
    yindex1=1
    yscale=0
  endif
  do jj=2,nlat
    if((rrlat.ge.aprlat(jj-1)).and.(rrlat.lt.aprlat(jj))) then
      yindex0=jj-1
      yindex1=jj
      yscale=(rrlat-aprlat(jj-1))/(aprlat(jj)-aprlat(jj-1))
    endif
  enddo
  if(rrlat.ge.aprlat(nlat)) then
    yindex0=nlat
    yindex1=nlat
    yscale=0
  endif
!  write(100+mype,*)'rrlat',rrlat,'yindex0',yindex0,yindex1,'yscale',yscale
!  write(100+mype,*)'mapriori',maxval(mapriori),minval(mapriori),'rrlat',rrlat,'yindex0',yindex0,yindex1,yscale
! construct interpolated monthly-latitudinal apriori
  do ii=1,ncol
    do ll=1,nprs
      apr0=mapriori(yindex0,ii,ll)
      apr1=mapriori(yindex1,ii,ll)
      if((apr0.lt.999.).and.(apr1.lt.999.)) yapriori(ii,ll)=apr0+yscale*(apr1-apr0)
      if((apr0.eq.999.).and.(apr1.lt.999.)) yapriori(ii,ll)=apr1
      if((apr0.lt.999.).and.(apr1.eq.999.)) yapriori(ii,ll)=apr0
      if((apr0.eq.999.).and.(apr1.eq.999.)) then
!        write(100+mype,*)'latitudinal apriori faile at col ',aprcol(ii)
!       print,'latitudinal apriori interp failed at (col) ',aprcol(ii)
        yapriori(ii,ll)=999.
      endif
    enddo
  enddo
  ifailed=0
  if(xdata1.lt.aprcol(1)) then
    cindex0=1
    cindex1=1
    cscale=0
  endif
! construct interpolated monthly-latitudinal-column apriori
  do ii=2,ncol
    if((xdata1.ge.aprcol(ii-1)).and.(xdata1.lt.aprcol(ii))) then
     cindex0=ii-1
     cindex1=ii
     cscale=(xdata1-aprcol(ii-1))/(aprcol(ii)-aprcol(ii-1))
    endif
  enddo
  if(xdata1.ge.aprcol(ncol)) then
    cindex0=ncol
    cindex1=ncol
    cscale=0
  endif
  do ll=1,nprs
    apr0=yapriori(cindex0,ll)
    apr1=yapriori(cindex1,ll)
    if((apr0.lt.999.).and.(apr1.lt.999.)) apr(ll)=apr0+cscale*(apr1-apr0)
    if((apr0.eq.999.).and.(apr1.lt.999.)) apr(ll)=apr1
    if((apr0.lt.999.).and.(apr1.eq.999.)) apr(ll)=apr0
    if((apr0.eq.999.).and.(apr1.eq.999.)) then
     apr(ll)=999.
! print *,'apriori interp failed at (i,j) ',ilon,ilat
!     write(100+mype,*)'failed level ',ll
!     if(mype.eq.0)write(6,*)'failed',ll,'apr',apr0,apr1,'cscale',cscale
     ifailed=1
    endif
  enddo
!
  filled=0
  ydata2=0.0

  intapr(nprs+1)=0.0
  do ll=nprs,1,-1
    intapr(ll)=intapr(ll+1)+apr(ll)
  end do
  intapr(0)=intapr(1)+apr(1)
!
!  llevold=-1
  gotbot=.false.
RAQL: do llev=1,nsig ! bottom up  here for  GSI
    if(psav(llev).lt.aprprs(nprs))eff(llev)=apreff(nprs)
    do ll=1,nprs-1
      if((aprprs(ll).ge.psav(llev)).and.(aprprs(ll+1).le.psav(llev))) then
        pscale=(psav(llev)-aprprs(ll+1))/(aprprs(ll)-aprprs(ll+1))
        eff(llev)=apreff(ll+1)+pscale*(apreff(ll)-apreff(ll+1))
      endif
    enddo
    if(psav(llev).gt.aprprs(1))eff(llev)=apreff(1)
!    write(100+mype,*)'eff',eff(llev),'llev',llev,'presi',presi(llev),'nprs',nprs
!    call flush(6)
    if(eff(llev).ne.1.0)then
      praqtop=presi(llev+1)
      praqbot=presi(llev)
!      dp=dpsav(llev)
!      write(100+mype,*)'praqbot',llev,praqbot,praqtop,'pintaprp',pintapr(nprs+1)
      if(praqbot<=pintapr(nprs+1))then
        ydata2(llev)=0.0
        cycle RAQL
      endif
      kbot=0
APRL: do iaprlev=kbot,nprs
        if(.not.gotbot)then
          if(pintapr(iaprlev+1)>=praqbot)then
            cycle APRL
          endif
!         in iaprlev box
          if(praqbot<pintapr(iaprlev+1).or.praqtop>pintapr(iaprlev))then
            write(100+mype,*)'error praqbot llev',llev,'iaprlev',iaprlev
            call flush(100+mype)
            write(6,*)'error praqbot llev',llev,'iaprlev',iaprlev
          call flush(6)
            call stop2(9999)
          endif
          afbot=(praqbot-pintapr(iaprlev+1))/aprdp(iaprlev)
          if(afbot<0.0.or.afbot>1.0)then
            write(100+mype,*)'afbot error',afbot,llev,'iaprlev',iaprlev, &
             'praq',praqtop,praqbot,'apr',aprprb(iaprlev),aprprt(iaprlev), &
            'pintapr',iaprlev+1,pintapr(iaprlev+1),'aprdp',aprdp(iaprlev)
            call flush(100+mype)
            write(6,*)'afbot error',afbot,llev,'iaprlev',iaprlev, &
             'praq',praqtop,praqbot,'apr',aprprb(iaprlev),aprprt(iaprlev), &
            'pintapr',iaprlev+1,pintapr(iaprlev+1),'aprdp',aprdp(iaprlev)
          call flush(6)
            call stop2(9999)
          endif
          aprintbot=intapr(iaprlev+1)*(1.-afbot)+intapr(iaprlev)*afbot
!          if(aprintbot<0.0)then
!            write(100+mype,*)'aprintbot',aprintbot,'afbot',afbot,'iaprlev',iaprlev, & 
!              'intpartop',intapr(iaprlev+1),intapr(iaprlev)
!            call flush(100+mype)
!          endif
!          if(iaprlev.eq.0)then
!            write(6,*)'below 1000 mb aprintbot ',aprintbot,' afbot ',afbot
!            call flush(6)
!          endif

          gotbot=.true.
        endif
        if(gotbot)then
          if(pintapr(iaprlev+1)>=praqtop)then
            cycle APRL
          endif
          if(praqtop<pintapr(iaprlev+1).or.praqtop>pintapr(iaprlev))then
            write(100+mype,*)'error praqtop llev',llev,'iaprlev',iaprlev,'praqtop',praqtop,praqbot, & 
               'pintapr',pintapr(iaprlev:iaprlev+1)
            call flush(100+mype)
            write(6,*)'error praqtop llev',llev,'iaprlev',iaprlev,'praqtop',praqtop,praqbot, & 
               'pintapr',pintapr(iaprlev:iaprlev+1)
          call flush(6)
            call stop2(9999)
          endif
          aftop=(praqtop-pintapr(iaprlev+1))/aprdp(iaprlev)
          if(aftop<0.0.or.aftop>1.0)then
            write(100+mype,*)'aftop error',aftop,llev,'iaprlev',iaprlev,pintapr(iaprlev), &
       'praq',praqtop,praqbot,'apr',aprprb(iaprlev),aprprt(iaprlev),'iaprlevp',iaprlev+1,pintapr(iaprlev+1), & 
        'aprdp',aprdp(iaprlev)
            call flush(100+mype)
            write(6,*)'aftop error',aftop,llev,'iaprlev',iaprlev,pintapr(iaprlev), &
       'praq',praqtop,praqbot,'apr',aprprb(iaprlev),aprprt(iaprlev),'iaprlevp',iaprlev+1,pintapr(iaprlev+1), & 
        'aprdp',aprdp(iaprlev)
          call flush(6)
            call stop2(9999)
          endif
          aprinttop=intapr(iaprlev+1)*(1.-aftop)+intapr(iaprlev)*aftop
          if(aprinttop<=0.0)then
            write(100+mype,*)mype,'aprinttop',aprinttop,'aftop',aftop,'iaprlev',iaprlev, & 
                  'intapr',intapr(iaprlev:iaprlev+1),'llev',llev,'rrlat',rrlat
            call flush(100+mype)
            write(6,*)mype,'aprinttop',aprinttop,'aftop',aftop,'iaprlev',iaprlev, & 
                  'intapr',intapr(iaprlev:iaprlev+1),'llev',llev,'rrlat',rrlat
          call flush(6)
            call stop2(9999)
          endif
        endif

        ydata2(llev)=aprintbot-aprinttop
        aprintbot=aprinttop
        !llevold=llev-1
        if(ydata2(llev)<0.0)Then
          write(100+mype,*)'aprintbot',aprintbot,aprinttop
            call flush(100+mype)
          write(6,*)mype,'aprintbot',aprintbot,aprinttop
          call flush(6)
          call stop2(9999)
        endif
        kbot=iaprlev
        exit APRL
      end do APRL
    endif
  enddo RAQL ! llev
!
! incorporate TOMS apriori and averaging kernel
!
  xdata2=0.0 ! ajl added 9/6/2015
  do llev=1,nsig
! tks  1-17-07
    if(eff(llev).ne.1.)then
!  ajl ydata2 only needed if eff ne 1
      xdata2=xdata2+ydata2(llev)*(1.-eff(llev))+dpoz(llev)*eff(llev)
    else
      xdata2=xdata2+dpoz(llev)
    endif
  enddo
!  write(100+mype,*)'omiadj',xdata2-xdata1,'orig',xdata1,'new',xdata2
! incorporate TOMS apriori and averaging kernel
  if(ifailed.eq.1)then
!    xdata2=xdata1
    return
  endif
  xdata1=xdata2
  return

end subroutine omiadjust

subroutine intrp3no2dp1(f,dpno2,dx,dy,obstime,mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    intrp3oz    space-time linear interpolation for ozone
!   prgmmr: parrish          org: np22                date: 2013-01-26
!  dpno2(k) for omi apriori adjustment Lenzen 2015-05-10
!
! abstract: same as intrp3oz but for special case n=1, with n argument removed.
!            This has been created to solve problem of type mismatch debug
!            compile
!            error on WCOSS.
!
! program history log:
!   2013-01-26  parrish
!
!   input argument list:
!     f        - input interpolator (gridded guess ozone fields)
!     dx,dy - input x,y-coords of interpolation points (grid units)
!     obstime  - observation times
!     mype     - mpi task id
!
!   output argument list:
!     g        - output interpolatees (guess ozone at observation location)
!     dpno2     - delta pressure * layer ozone at observation location
!     sfcp     - surface pressure at observation location
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
!--------
  use kinds, only: r_kind,i_kind
  use guess_grids, only: nfldsig,hrdifsig,ges_prsi
  use gridmod, only: lat2,lon2,nlat,nlon,nsig,lon1,istart,jstart
  use constants, only: zero, one, rno2con
  implicit none

! Declare passed variables
  integer(i_kind)                               ,intent(in   ) :: mype
  real(r_kind)                                  ,intent(in   ) :: dx,dy,obstime
  real(r_kind),dimension(lat2,lon2,nsig,nfldsig),intent(in   ) :: f
!  real(r_kind)                                  ,intent(  out) :: g
  real(r_kind),dimension(nsig)                 ,intent(out)    :: dpno2
  

! Declare local variables
  integer(i_kind) j,k,ix,ix1,iy,iy1,kk,itsig,itsigp,iz1,iz2
  integer(i_kind) ixp,iyp,mm1
  real(r_kind) w00,w01,w10,w11,delx,dely,delx1,dely1
  real(r_kind) dtsig,dtsigp,pob
  real(r_kind) delp1,delp2,delp3,delp4,delp5,delp6,delp7,delp8


!*************************************************************************
! Initialize variables
!  g=zero
!  dp=zero
  dpno2=zero
!  sfcp=zero
  mm1=mype+1


! Loop over number of observations.

!    Get horizontal interpolation information.  This information includes
!    the (i,j) indices of the grid points surrounding the observation,
!    plus the corresponding interpolation weights between these points
!    and the observation.

     ix1=dx; iy1=dy
     ix1=max(1,min(ix1,nlat))
     delx=dx-ix1; dely=dy-iy1; delx=max(zero,min(delx,one))
     ix=ix1-istart(mm1)+2; iy=iy1-jstart(mm1)+2
     if(iy<1) then
        iy1=iy1+nlon
        iy=iy1-jstart(mm1)+2
     end if
     if(iy>lon1+1) then
        iy1=iy1-nlon
        iy=iy1-jstart(mm1)+2
     end if
     ixp=ix+1; iyp=iy+1
     if(ix1==nlat) then
        ixp=ix
     end if
     delx1=one-delx; dely1=one-dely
     w00=delx1*dely1; w10=delx*dely1; w01=delx1*dely; w11=delx*dely
        
!    Set the weights for linear time iterpolation from the guess to 
!    the observation time.
!     if(mype.eq.0)then
!          write(6,*)'obstime',obstime,'hrdiffsig',hrdifsig,'nfldsig',nfldsig
!           call flush(6)
!     endif

     if(obstime > hrdifsig(1) .and. obstime < hrdifsig(nfldsig))then
        do j=1,nfldsig-1
           if(obstime > hrdifsig(j) .and. obstime <= hrdifsig(j+1))then
              itsig=j
              itsigp=j+1
              dtsig=((hrdifsig(j+1)-obstime)/(hrdifsig(j+1)-hrdifsig(j)))
           end if
        end do
     else if(obstime <=hrdifsig(1))then
        itsig=1
        itsigp=1
        dtsig=one
     else
        itsig=nfldsig
        itsigp=nfldsig
        dtsig=one
     end if
     dtsigp=one-dtsig
     
!
!    Perform spatial and temporal interpolation for the total column 
!    ozone observation
!     write(100+mype,*)'ix',ix,iy,itsig,itsigp,'prsi',ges_prsi(ix,iy,1,itsig),ges_prsi(ix,iy,nsig,itsig)
     do kk=1,nsig
        delp1=ges_prsi(ix ,iy ,kk,itsig )-ges_prsi(ix ,iy ,kk+1,itsig )
        delp2=ges_prsi(ixp,iy ,kk,itsig )-ges_prsi(ixp,iy ,kk+1,itsig )
        delp3=ges_prsi(ix ,iyp,kk,itsig )-ges_prsi(ix ,iyp,kk+1,itsig )
        delp4=ges_prsi(ixp,iyp,kk,itsig )-ges_prsi(ixp,iyp,kk+1,itsig )
        delp5=ges_prsi(ix ,iy ,kk,itsigp)-ges_prsi(ix ,iy ,kk+1,itsigp)
        delp6=ges_prsi(ixp,iy ,kk,itsigp)-ges_prsi(ixp,iy ,kk+1,itsigp)
        delp7=ges_prsi(ix ,iyp,kk,itsigp)-ges_prsi(ix ,iyp,kk+1,itsigp)
        delp8=ges_prsi(ixp,iyp,kk,itsigp)-ges_prsi(ixp,iyp,kk+1,itsigp)
!        dp(kk)=dp(kk)+(w00*delp1+w10*delp2+w01*delp3+w11*delp4)*dtsig + &
!                      (w00*delp5+w10*delp6+w01*delp7+w11*delp8)*dtsigp
        dpno2(kk)=dpno2(kk)+ &
              (f(ix ,iy ,kk,itsig )*w00*delp1 &
             + f(ixp,iy ,kk,itsig )*w10*delp2 &
             + f(ix ,iyp,kk,itsig )*w01*delp3 &
             + f(ixp,iyp,kk,itsig )*w11*delp4)*dtsig + &
              (f(ix ,iy ,kk,itsigp)*w00*delp5 &
             + f(ixp,iy ,kk,itsigp)*w10*delp6 &
             + f(ix ,iyp,kk,itsigp)*w01*delp7 &
             + f(ixp,iyp,kk,itsigp)*w11*delp8)*dtsigp
        dpno2(kk)=dpno2(kk)*rno2con
     enddo
!     write(500+mype,*)'bottom intrp3no2dp1'
!     call flush(500+mype)

! End of routine
  return
end subroutine intrp3no2dp1
subroutine raqmsno2col(dpno2,no2ges,mype,ptrop,prsitmp,scatwt,amf,nsig,ntrop,scatwtmodel)
  use kinds, only: r_kind,i_kind
  use omimod, only : nlevscatwt,presscatwt,pintscatwt,pintscatwtcb
  implicit none
  integer(i_kind)                               ,intent(in   ) :: mype
  integer(i_kind)                               ,intent(in   ) :: nsig
  real(r_kind), dimension(nlevscatwt) , intent(in) :: scatwt
  real(r_kind)                    ,intent(in)  :: amf
! dpno2 in millidobson
  real(r_kind),dimension(nsig)                 ,intent(in)    :: dpno2
  real(r_kind)                                 ,intent(out)   :: no2ges
! ptrop in mb
  real(r_kind)                                 ,intent(in)    :: ptrop
! prsitmp in cb
  real(r_kind),dimension(0:nsig)               ,intent(in)    :: prsitmp
  real(r_kind),dimension(0:nlevscatwt)            :: aintscatwt
  real(r_kind), dimension(0:nsig)      :: aintscatwtmodel
  real(r_kind), dimension(nsig)    ,intent(out)    :: scatwtmodel
  real(r_kind)                         :: ptropcb,af,prskk
  integer(i_kind)                :: ntop,ntrop,k,kk,kbot
  integer entry
  data entry/0/
  save entry
  entry=entry+1
  
! need to integrate scatterweight
  aintscatwt(0)=0.
! need to find aintscatwt at model interfaces
! need to find which layer ptrop is in
  ptropcb=ptrop/10.
! need highest model layer containing ptrop
  ntrop=nsig
!  write(540+mype,*)'prsitmp sfc',prsitmp(0)
  do k=nsig,1,-1
!    write(540+mype,*)k,'ptropcb',ptropcb,'prsitmp',k-1,prsitmp(k-1),'ptrop',ptrop,'prsitmpk',prsitmp(k)
    if(ptropcb<prsitmp(k-1))then
!    if(ptropcb>=prsitmp(k).and.ptropcb<prsitmp(k-1))then
!     so prsitmp(ntrop) is highest model interface need pintscalwtcb to bracket
!     it so can interp to it
      ntrop=k
      exit ! since not bracket any more
    end if
  end do
  ntop=nlevscatwt
  do k=nlevscatwt,1,-1
!    if(prsitmp(ntrop)<=pintscatwt(k))then
    if(prsitmp(ntrop)<pintscatwt(k-1))then
      ntop=k
!      ntop=k+1
      exit
    endif
  end do
  do k=1,ntop
    aintscatwt(k)=aintscatwt(k-1)+scatwt(k)*(pintscatwtcb(k-1)-pintscatwtcb(k))
!    if(entry<10)then
!      write(540+mype,*)k,'aintscatwt',aintscatwt(k),'scatwt',scatwt(k),'pint',pintscatwtcb(k)
!    endif
  end do
!    if(entry<10)then
      !write(540+mype,*)'ntop',ntop,'ntrop',ntrop
    !endif
! in between pintscatwt(ntop) and pintscatwt(ntop+1)
! find index for aintscatwt on pintscatwt
  kbot=1
  do kk=0,ntrop
LP: do k=kbot,ntop
!      if(kk.eq.ntrop)then
!        !write(540+mype,*)'kk',kk,'prsitmp',prsitmp(kk)
!        !write(540+mype,*)'k ',k-1,k,pintscatwtcb(k-1:k),'kbot',kbot,'ntop',ntop
!      endif
      if(kk.eq.ntrop)then
        prskk=ptropcb
      else
        prskk=prsitmp(kk)
      endif
!      if(prsitmp(kk)<pintscatwtcb(k-1).and.prsitmp(kk)>pintscatwtcb(k))then
      if(prskk<pintscatwtcb(k-1).and.prskk>pintscatwtcb(k))then
!        kint(kk)=k
!        if(kk.eq.ntrop)then
!          af=(max(ptropcb,prsitmp(kk))-pintscatwtcb(k))/ &
          !(pintscatwtcb(k-1)-pintscatwtcb(k))
!          write(540+mype,*)k,'aftrop',af,'ptropcb',ptropcb,pintscatwtcb(k-1:k)
!          write(540+mype,*)'ntop',ntop,'ntrop',ntrop,'prsitmp',prsitmp(kk-1:kk+1)
!          write(540+mype,*)'aintscatwt',k-1,k,aintscatwt(k-1:k)
!        else
          af=(prskk-pintscatwtcb(k))/ &
          (pintscatwtcb(k-1)-pintscatwtcb(k))
!        endif
        if(af<0.0.or.af>1.000)then
          write(540+mype,*)'af out of bounds kk',kk,'ntrop',ntrop,'k',k, &
          'ptropcb',ptropcb,'pintscatwtcb',k,pintscatwtcb(k-1:k),'prsitmp',kk,prsitmp(kk)
        endif
        aintscatwtmodel(kk)=aintscatwt(k-1)*af+aintscatwt(k)*(1.-af)
        kbot=k
        exit LP
      endif 
    end do LP
  end do
!  write(540+mype,*)'amf',amf
  do k=1,ntrop
    scatwtmodel(k)=(aintscatwtmodel(k)-aintscatwtmodel(k-1))/ &
!      (prsitmp(k-1)-prsitmp(k))/amf
      (prsitmp(k-1)-prsitmp(k))
!     this adjusts for fact use smaller deltap for top layer add /amf so get a
!     normalized wt
!    write(540+mype,*)'scatwtmodel',k,scatwtmodel(k), 'scaled',scatwtmodel(k)/amf, &
!    .5*(prsitmp(k)+prsitmp(k-1))*10.,'aintscatwtmodel',aintscatwtmodel(k-1:k)
    scatwtmodel(k)=scatwtmodel(k)/amf
  end do
  no2ges=0.0
  do k=1,ntrop
    if(k.eq.ntrop)then
!     only part of layer ?
!      if(entry<10)then
!        write(580+mype,*)'ptropcb',ptropcb,prsitmp(k-1:k)
!      endif
      
      no2ges=no2ges+scatwtmodel(k)*dpno2(k)
    else
      no2ges=no2ges+scatwtmodel(k)*dpno2(k)
    endif
!    if(entry<10)then
!      write(580+mype,*)k,'scatwtmodel',scatwtmodel(k),'dpno2',dpno2(k)
!    endif
  end do
!  no2ges=no2ges/amf
  return
end subroutine raqmsno2col
