  implicit none
  include 'netcdf.inc'
  integer nvar,nv,nvarp
  parameter (nvar=35,nvarp=6)
  integer ivar(nvar),nvarin,ivarp(nvarp)
  character *8 cvar(nvar),cvarp(nvarp),filesurface*256,filesp*256,filepress*256,user*20
  character *256 pathgrid
  integer iret,ncid,nt,timdimid,latdimid,londimid,nc,nr,latsdimid,lonsdimid
  integer ivdiml(2),istartl(2),icountl(2),idate,iwrf
  integer ivdim2(3),icount2(3),istart2(3),i,j,isizdat,isizdat8,icountu2(3),icountv2(3)
  integer istart2i(2),icount2i(2),ivdims(2),ivdimss(3)
  integer icount3(4),istart3(4),levdimid,ivdim3(4)
  integer layid,jdatein11,jtimein11,jdatein12,jtimein12
  integer itop,timeid,time,it,ii,ncin,nrin,ntime,k,itsave(6),ntimep
  integer z8id,xlatid,xlonid,szid,hr,ihr,nskip,psfcid,ioff,hour
  integer xlatuid,xlonuid,xlatvid,xlonvid
  integer dxy_mcid,pdashid
  integer iymd,iymdp,juldate,juldatep,julhrs,julhrsp,skip,gjuldate,sigmaid
!  integer mapfacmid(3)
  character *60 exp,cymdh*20,chour*10,cymdhres*20,hourrestart*2
  character *256 metdir,cmaqrundir,ctime*20
  integer ymdh,id,dday,yr,md,mon,day,daysave,monsave,itime
  integer days(12),itflag(2),ierr,lent
  data days/31,28,31,30,31,30,31,31,30,31,30,31/
  real rx,rn,rs
  logical fexist,save,first
  save first
  data first/.true./
  integer icol,jdatein,jtimein,nvarsin,nlayin,nlay,nvarsinp,hrrestart
  character *256 filen,date*19,cin*4,name*16,fileno2*256
  ! JDE:
  character *8 metfiledate
!  parameter (nc=151,nr=142,nt=3*24)
!  parameter (nc=279,nr=279,nt=49)
! JDE why are nrin and ncol hardcoded
! JDE parameter (ncin=442,nrin=265,nc=ncin,nr=nrin,nlayin=35,nlay=nlayin)
  parameter (ncin=187,nrin=187,nc=ncin,nr=nrin,nlayin=44,nlay=nlayin)
  real*4, allocatable,dimension(:,:,:) :: alat,alon,landmask,psfc,alatu,alatv,alonu,alonv
  real*4, allocatable,dimension(:,:) :: msfx2
  real*4, allocatable,dimension(:,:) :: mapfacm,sz
  real*4, allocatable,dimension(:,:,:,:) :: p,no2,zh,zf,dens,densj,pdash,temp,qv
  real*4, allocatable,dimension(:,:,:) :: psolidm
  real*4, allocatable, dimension(:,:) :: dxy_mc
!  real, allocatable,dimension(:,:) :: sz8
  real*4, allocatable, dimension(:,:) :: sigma
!  real top(nc,nr,nt),wrf(nc,nr,nt),range(2),topmin
!  real col(nc,nr,nt)
!  real*4 input(nc,nr)
  real*4 dxy,pdashnew,dx,dy,no2max,converttoppmv
  data dxy/108000./
  real pdashbar,delptop,delpbot,delplay,grav
  real rmx,rmn,rms
  data cvarp/'P','DENS','ZF','ZH','TA','QV'/
  data cvar/'NO2','O3','CO','PAN','SO2',&
            'ASO4J','ASO4I','ANH4J','ANH4I','ANO3J',&
            'ANO3I','AORGAJ','AORGAI','AORGPAJ','AORGPAI',&
            'AORGBJ','AORGBI','AECJ','AECI','A25J',&
            'A25I','ACORS','ASOIL','AH2OJ','AH2OI',&
            'ANAJ','ANAI','ACLJ','ACLI','ANAK',&
            'ACLK','ASO4K','ANO3K','AH2OK','NH3'/
!            'ACLK','ASO4K','ANO3K','AH2OK','NH3',&
!            'EXT_Mie','EXT_Recon'/
!  character *8 yymd,yymdp,hrs*2,cpath*256,yyyyjjj*7
  character *8 yymdp,hrs*2,cpath*256,cpathout*256
  character *10 yymd,yyyyjjj*10,yymdres
  grav=9.80616

!
  write(*,*) 'IN, BEFORE GETTING DATES'
  call flush(6)
!

! read in yyyjjj and yymd from the environment
!
! they are set by script as julian_date=> yyyyjjj and normal_date=> yymd

! JDE think this isnt needed, was used for met file date
! JDE call getenv( 'HOUR', chour )
! JDE read(chour,*)hour
! JDE write(6,*)'hour',hour
! JDE call flush(6)

call getenv( 'METFILEDATE', metfiledate )
write(6,*)'metfiledate ', metfiledate
call flush(6)

!  write(*,*) '++++++++++++++++++++++++++++++++++++++++++++++++ '
!  write(*,*) 'DATES '
!  write(*,*) 'hour =', hour
!  write(*,*) '++++++++++++++++++++++++++++++++++++++++++++++++ '

!  exp=' '
  write(6,*)'get arg1'
  call flush(6)
  call getarg(1,cymdh)
!  call getarg(2,cymdhres)
!  yymdres=cymdhres(1:8)
!  hourrestart=cymdhres(9:10)
  write(6,*)'cymdh ',trim(cymdh) 
  call flush(6)
  read(cymdh,'(i10)')ymdh
  write(6,*)'ymdh',ymdh
  call flush(6)
! now doing every 3 hours need to figure out ioff

!  ioff=25 ! use last hour
    allocate (p(nc,nr,0:nlay,1),sigma(nlay,1),psfc(nc,nr,1))
    allocate (zf(nc,nr,0:nlay,1),zh(nc,nr,nlay,1))
    allocate (pdash(nc,nr,nlay,1),psolidm(nc,nr,0:nlay))
    allocate (dxy_mc(nc,nr))
  
  dday=ymdh
  yr=dday/1000000
  md=dday-1000000*yr
  mon=md/10000
  day=md-10000*mon
  hr=mod(day,100)
  day=day/100
  write(6,*)'yr',yr,'mon',mon,'day',day
  call flush(6)
  if(mon>12)stop
  
20 if(day > days(mon) )then
     day=day-days(mon)
     mon=mon+1
     go to 20
   endif
  iymd=(yr*100+mon)*100+day 
  write(yymd,'(i8.8)')iymd
  juldate=gjuldate(iymd)
!  user='lenzen'
!  cpath='/odyssey/proxy/mrogal/cmaqerr/'//yymd(1:8)//'00'
!  cpath='/odyssey/proxy/'//trim(user)//'/cmaqerr/'//trim(yymd)//'00'
!  cpathout='/iliad/proxy/'//trim(user)//'/cmaqgsi/'//trim(yymd)//'00'
  cpathout=' '
  call getenv('GSIOUT',cpathout)
  write(6,*)'cpathout ',trim(cpathout)
  call flush(6)
!  call system('mkdir /odyssey/proxy/mrogal/cmaqerr')
  call system('mkdir -p '//trim(cpathout))
!  call system('mkdir /odyssey/proxy/mrogal/GSI_cmaq/cmaqerr')
!  call system('mkdir /odyssey/proxy/mrogal/cmaqerr/'//trim(yymd)//'00')
!  call system('mkdir -p /odyssey/proxy/'//trim(user)//'/cmaqerr/'//trim(yymd)//'00')
!  call system('mkdir /odyssey/proxy/mrogal/GSI_cmaq/cmaqerr/'//yymd(1:8)//'00')


!  yyyyjjj='2011180'
  write(yyyyjjj,'(i4.4,i3.3)')yr,juldate
!  yymd='20110619'
  write(yymd,'(i4.4,i2.2,i2.2)')yr,mon,day
  metdir=' '
  call getenv('METDIR',metdir)
  write(6,*)'metdir ',trim(metdir)
  ctime=' '
  call getenv('ITIME',ctime)
  if(ctime.ne.'')then
    write(6,*)'ctime ',trim(ctime)
    read(ctime,*)itime
    call flush(6)
  else
    write(6,*)'no ITIME'
    stop

  endif
  call getenv('PATHGRID',pathgrid)
  write(6,*)'pathgrid ',pathgrid
!  cmaqrundir=' '
!  call getenv('CMAQRUNDIR',cmaqrundir)
!   write(6,*)'cmaqrundir',trim(cmaqrundir)
  if(metdir.eq.'') then
    write(*,*)'no METDIR env var, STOPPING'
    stop
  endif
!  filesurface=trim(metdir)//'/GRIDCRO2D.'//trim(yyyyjjj)//'.12US2.Z25.CB5.gmt'
   !filesp=trim(metdir)//'/METCRO2D.108NHEMI2.44L.'//trim(metfiledate)!//'.nc' !JDE  e.g. METCRO2D_20140709.nc
   filesp=trim(metdir)//'/METCRO2D_20'//trim(metfiledate)//'.nc4' !JDE  e.g. METCRO2D_20140709.nc
   !filesp=trim(metdir)//'/METCRO2D_'//trim(metfiledate)//'.nc' !JDE  e.g. METCRO2D_20140709.nc
   filepress=trim(metdir)//'/METCRO3D_20'//trim(metfiledate)//'.nc4' !JDE  e.g. METCRO3D_20140709.nc
   !filepress=trim(metdir)//'/METCRO3D_'//trim(metfiledate)//'.nc' !JDE  e.g. METCRO3D_20140709.nc
!JDE   filesp=trim(metdir)//'aqm.t'//trim(cHOUR)//'z.metcro2d.ncf'
!JDE   filepress=trim(metdir)//'aqm.t'//trim(cHOUR)//'z.metcro3d.ncf'
!  filesp=trim(metdir)//'/METCRO2D.'//trim(yyyyjjj)//'.12US2.Z25.CB5.gmt'
!  filepress=trim(metdir)//'/METCRO3D.'//trim(yyyyjjj)//'.12US2.Z25.CB5.gmt'
!  fileno2=trim(cmaqrundir)//'/CCTM_July2011_LADCO.CGRID.'//trim(yymdres)//'_'//hourrestart//'.nc'
  fileno2=' '
! JDE NO2 file is file from opGSI.F.
  call getenv('CMAQRESTART',fileno2)
  if(fileno2.eq.'') then
    write(*,*)'no CMAQRESTART env var, STOPPING'
    stop
  endif
  write(6,*)'fileno2',trim(fileno2)
!JDE  filesurface=trim(pathgrid)//'aqm.t'//trim(chour)//'z.grdcro2d.ncf'
  filesurface=trim(pathgrid)//'/GRIDCRO2D_20'//trim(metfiledate)//'.nc4'! !JDE eg GRIDCRO2D_2010709
!  filesurface=trim(pathgrid)//'/GRIDCRO2D.108NHEMI2.44L.'//trim(metfiledate)!//'.nc' !JDE eg GRIDCRO2D_2010709
!  filesurface='/odyssey/proxy/mharkey/CMAQv5/July2011_LADCO/mcip/GRIDCRO2D.'//trim(yyyyjjj)//'.12US2.Z25.CB5.gmt'
!  filesp='/odyssey/proxy/mharkey/CMAQv5/July2011_LADCO/mcip/METCRO2D.'//trim(yyyyjjj)//'.12US2.Z25.CB5.gmt'
!  filepress='/odyssey/proxy/mharkey/CMAQv5/July2011_LADCO/mcip/METCRO3D.'//trim(yyyyjjj)//'.12US2.Z25.CB5.gmt'
!  fileno2='/odyssey/proxy/mharkey/CMAQv5/July2011_LADCO/cctm/reduced_allNOx/CCTM_July2011_LADCO.CONC.'//trim(yymd)//'_00.nc'
  print *,'iymd',iymd

!
! Let's check what gets passed
!
  write(*,*) '                                                 '
  write(*,*) '++++++++++++++++++++++++++++++++++++++++++++++++ '
  write(*,*) 'DATES '
  write(*,*) 'filepress = ',filepress
  write(*,*) 'fileno2 = ', fileno2
  write(*,*) '++++++++++++++++++++++++++++++++++++++++++++++++ '
   call flush(6)


!  exp='raqms'
! read in cmaq lat lon
  if(first)then
    if(.not.allocated(alat))then
      allocate (msfx2(ncin,nrin))
      write(6,*)'find ncin',ncin,nrin

      allocate (alat(ncin,nrin,1),alon(ncin,nrin,1),sz(ncin,nrin))
      allocate (mapfacm(ncin,nrin))
    endif
    call readsurface(filesurface,ncin,nrin,alat,alon,sz,msfx2,dx,dy)
    write(6,*)'JDE dx: ',dx
    write(6,*)'JDE dy: ',dy
    zf(:,:,0,1)=sz(:,:)
    mapfacm=sqrt(msfx2)
    do j=1,nrin
      do i=1,ncin
        dxy_mc(i,j)=dx/mapfacm(i,j)
      end do
    end do
    first=.false.
  endif
!  stop
  write(6,*)'ncin ',ncin,nrin,nlayin
  call flush(6)
!  range(1)=0.
!  range(2)=500.
  istart2(1:2)=1
  istart2i=1
  icount2i(1)=nc
  icount2i(2)=nr
  icount2(1)=nc
  icount2(2)=nr
  icount2(3)=1
  icountu2(1)=nc+1
  icountu2(2)=nr
  icountu2(3)=1
  icountv2(1)=nc
  icountv2(2)=nr+1
  icountv2(3)=1
  istart2(3)=1
  istart3(1:3)=1
  icount3(1)=nc
  icount3(2)=nr
  icount3(3)=nlay
  icount3(4)=1
  istart3(4)=1
    write(yymd,'(i8.8)')iymd
  if(.not.allocated(no2))then
    allocate (no2(nc,nr,nlay,1),temp(nc,nr,nlay,1),qv(nc,nr,nlay,1))
  endif
  allocate (dens(nc,nr,nlay,1),densj(nc,nr,nlay,1))
  ihr=hr+1 !JDE +1 for nf90_get_var
  write(6,*)'hour ',ihr
  call readsp(filesp,nc,nr,psfc,ihr)
  print *,'psfc',maxval(psfc)
  p(:,:,0,1)=psfc(:,:,1)
! looks like p is pdash and not psolid
  call readpressure(filepress,nc,nr,nlay,pdash(1,1,1,1),dens,densj,zh,zf(1,1,1,1),temp,qv,ihr)
  print *,'p',maxval(p),minval(p)
  call readcmaqno2(fileno2,nc,nr,nlay,no2,itflag)
  no2max=maxval(no2)
  print *,'no2 max ',maxval(no2)
  if(no2max>4.)then
!   must be ppbV input need to convert to ppmv
    write(6,*)'convert no2 in ppbV to ppmV for GSI'
    no2=.001*no2
    print *,'no2 in ppmv',maxval(no2)
  elseif(no2max>.004)then
!   must be ppmv ok
    write(6,*)'no2 from NMMB-CMAQ is in ppmv good for GSI'
  elseif(no2max<10e-6)then
!   must be ppv input need to convert to ppmv
    write(6,*)'convert no2 in ppv to ppmv for GSI'
    no2=no2*1.e6
    print *,'no2 in ppmv',maxval(no2)
  endif
! calculate pdash two times and average
!  write(6,*)'ioff',ioff
!  do k=1,nlay
!    write(6,*)'pdash',k,maxval(pdash(:,:,k,1)),minval(pdash(:,:,k,1))
!  end do
  write(6,*)'psfc',maxval(psfc),minval(psfc)
  do j=1,nr
    do i=1,nc
      p(i,j,nlay,1)=pdash(i,j,nlay,1)-dens(i,j,nlay,1)*grav* (zf(i,j,nlay,1)-zh(i,j,nlay,1))
!      if(i.eq.nc/2.and.j.eq.nr/2)then
!       write(6,*)'ptop',p(i,j,nlay,1),'pdashtop',pdash(i,j,nlay,1)
!       write(6,*)'dens',dens(i,j,nlay,1),'grav',grav,'dz',zf(i,j,nlay,1)-zh(i,j,nlay,1)
!       write(6,*)'delp',dens(i,j,nlay,1)*grav* (zf(i,j,nlay,1)-zh(i,j,nlay,1))
!      endif
    end do
  end do
  do k=1,nlay-1
    do j=1,nr
      do i=1,nc
!        pdash(i,j,k,1)=.5*(p(i,j,k,1)+p(i,j,k-1,1))+dens(i,j,k,1)*grav*(.5*(zf(i,j,k,1)+zf(i,j,k-1,1))-zh(i,j,k,1))
        p(i,j,k,1)=.5* (pdash(i,j,k,1)+pdash(i,j,k+1,1)+grav* & 
        (dens(i,j,k+1,1)* (zh(i,j,k+1,1)-zf(i,j,k,1)) -dens(i,j,k,1)* (zf(i,j,k,1)-zh(i,j,k,1)) ) )
      end do
    end do
!    print *,'pdash',k,maxval(pdash(:,:,k,1)),minval(pdash(:,:,k,1))
!    print *,'psol',k,maxval(p(:,:,k,1)),minval(p(:,:,k,1))
  end do
  do j=1,nr
    do i=1,nc
      if(p(i,j,1,1)>psfc(i,j,1))then
       write(6,*)'surface problem',i,j,p(i,j,1,1),'psfc',psfc(i,j,1)
      endif
    end do
  end do
  print *,'psol sfc',maxval(p(:,:,0,1))
!  do k=1,nlay
!    print *,'pdash',k,maxval(pdash(:,:,k,1))
!    print *,'psol ',k,maxval(p(:,:,k,1))
!  end do
   
  
  

  time=0
  daysave=day
  monsave=mon
  ntime=1
!  itsave=0
! mrogal changed hr to hours to reflect the hours we are doing and enable
! passing of the hour as a variable
!  hour=hr
!  hr=hour
  write(6,*)'ntime ',ntime
  call flush(6)
  do it=1,ntime
    save=.true.
    iymd=(yr*100+mon)*100+day 
    itsave(1)=yr
    itsave(2)=mon
    itsave(3)=day
    itsave(4)=hr
    itsave(5)=0
    itsave(6)=0
!    itsave(3)=day-1
!    itsave(4)=22
    write(6,*)'after hr=',hr,'day',day,'mon',mon
    write(yymd,'(i8.8)')iymd
    write(hrs,'(i2.2)')hr
    print *,'yymd',yymd,' hr ',hrs
    write(6,*)'it',it
!   read over p stuff and get psfc
  
!   netcdf stuff here for each file
!    if(save)then
!     start save
      date=yymd(1:4)//'-'//yymd(5:6)//'-'//yymd(7:8)//'_'//hrs//':00:00'
      lent=len_trim(cpathout)
      if(cpathout(lent:lent).eq.'/')then
        filen=trim(cpathout)//'/cmaq_d01_'//yymd(1:4)//'-'//yymd(5:6)// &
        '-'//yymd(7:8)//'_'//hrs//':00:00'
      else
        filen=trim(cpathout)//'/cmaq_d01_'//yymd(1:4)//'-'//yymd(5:6)// &
        '-'//yymd(7:8)//'_'//hrs//':00:00'
      endif
      write(6,*)'fileout cmaqtogsifile filen',trim(filen)
      call flush(6)
!     allow to write over old one
      isizdat=NF_FLOAT
      isizdat8=NF_DOUBLE
      call system('mkdir '//trim(cpathout))
      write(6,*)'pathout',trim(cpathout)
      write(6,*)'create ',trim(filen)
      iret=nf_create(filen,NF_CLOBBER,ncid)
      if(iret.ne.nf_noerr)then
        write(6,*)'error open filen'
        write(6,*)trim(nf_strerror(ierr))
        stop
      endif
      write(6,*)'find nc',nc,nr
      call flush(6)
      iret=nf_def_dim(ncid,'south-north',nr,latdimid)
      iret=nf_def_dim(ncid,'south-north_stag',nr+1,latsdimid)
      iret=nf_def_dim(ncid,'west-east',nc,londimid)
      iret=nf_def_dim(ncid,'west-east_stag',nc+1,lonsdimid)
      iret=nf_def_dim(ncid,'lev',nlay,levdimid)
      iret=nf_def_dim(ncid,'time',NF_UNLIMITED,timdimid)
      iret=nf_def_dim(ncid,'DateStrLen',19,idate)
      ivdims(1)=levdimid
      ivdims(2)=timdimid
      ivdiml(1)=londimid
      ivdiml(2)=latdimid
      ivdim2(1)=londimid
      ivdim2(2)=latdimid
      ivdim2(3)=timdimid
      !iret=nf_def_var(ncid,'XLAT',isizdat,2,ivdiml,xlatid)
!      iret=nf_def_var(ncid,'XLONG',isizdat,2,ivdiml,xlonid)
      iret=nf_def_var(ncid,'XLAT',isizdat,3,ivdim2,xlatid)
      iret=nf_def_var(ncid,'XLONG',isizdat,3,ivdim2,xlonid)
      iret=nf_def_var(ncid,'SZ',isizdat,3,ivdim2,szid)
      iret=nf_def_var(ncid,'PSFC',isizdat,3,ivdim2,psfcid)
      iret=nf_def_var(ncid,'ZNU',isizdat,2,ivdims,sigmaid)
!      iret=nf_def_var(ncid,'MAPFAC_M',isizdat,3,ivdim2,mapfacmid(1))
      ivdimss(1)=londimid
      ivdimss(2)=latsdimid
      ivdimss(3)=timdimid
      iret=nf_def_var(ncid,'DXY_MC',isizdat,3,ivdimss,dxy_mcid)
  
      iret=nf_put_att_text(ncid,szid,'coordinates',10,'XLONG XLAT')
      iret=nf_put_att_real(ncid,NF_GLOBAL,'DX',NF_FLOAT,1,dxy)
      iret=nf_put_att_real(ncid,NF_GLOBAL,'DY',NF_FLOAT,1,dxy)
      iret=nf_put_att_int(ncid,NF_GLOBAL,'TIME',NF_INT,6,itsave)
      iret=nf_put_att_text(ncid,psfcid,'coordinates',10,'XLONG XLAT')
      iret=nf_put_att_text(ncid,xlatid,'coordinates',10,'XLONG XLAT')
      iret=nf_put_att_text(ncid,xlonid,'coordinates',10,'XLONG XLAT')
      iret=nf_put_att_text(ncid,xlonid,'description',27,'LONGITUDE, WEST IS NEGATIVE')
      iret=nf_put_att_text(ncid,xlonid,'units',11,'degree_east')
      iret=nf_put_att_text(ncid,xlonid,'stagger',0,'')
      iret=nf_put_att_text(ncid,xlatid,'description',27,'LATITUDE, SOUTH IS NEGATIVE')
      iret=nf_put_att_text(ncid,xlatid,'units',12,'degree_north')
      iret=nf_put_att_text(ncid,xlatid,'stagger',0,'')
!      iret=nf_def_var(ncid,'lay',isizdat,1,levdimid,layid)
!      write(6,*)'laynid',layid
!      call flush(6)
      iret=nf_def_var(ncid,'time',NF_INT,1,timdimid,timeid)
      write(6,*)'latdimid',latdimid,londimid,timdimid
      call flush(6)
!     isizdat=NF_DOUBLE
      write(6,*)'NF_FLOAT',NF_FLOAT
      write(6,*)'isizdat',isizdat
      call flush(6)
      ivdim3(1)=londimid
      ivdim3(2)=latdimid
      ivdim3(3)=levdimid
      ivdim3(4)=timdimid
      write(6,*)'filled ivdim'
      call flush(6)
!     do nv=1,nvarp
      nv=1
!      iret=nf_def_var(ncid,cvarp(nv),isizdat,4,ivdim3,ivarp(nv))
      iret=nf_def_var(ncid,'PSOL',isizdat,4,ivdim3,ivarp(nv))
     
      iret=nf_put_att_text(ncid,ivarp(nv),'coordinates',10,'XLONG XLAT')
      iret=nf_def_var(ncid,'PDASH',isizdat,4,ivdim3,pdashid)
      iret=nf_put_att_text(ncid,pdashid,'coordinates',10,'XLONG XLAT')

      write(6,*)'cvarp',nv,cvarp(nv),'ivarp',ivarp(nv)
      call flush(6)
!     end do
      do nv=5,6
        iret=nf_def_var(ncid,cvarp(nv),isizdat,4,ivdim3,ivarp(nv))
        iret=nf_put_att_text(ncid,ivarp(nv),'coordinates',10,'XLONG XLAT')
        write(6,*)'cvar',nv,cvar(nv),ivar(nv)
      end do
!     do nv=1,nvar
      nv=1 ! NO2
      write(6,*)'cvar nv',nv,cvar(nv),'ivar',ivar(nv)
      write(6,*)'isizdat',isizdat,'ivdim3',ivdim3
      call flush(6)
      iret=nf_def_var(ncid,cvar(nv),isizdat,4,ivdim3,ivar(nv))
      iret=nf_put_att_text(ncid,ivar(nv),'coordinates',10,'XLONG XLAT')
      write(6,*)'cvar',nv,cvar(nv),ivar(nv)
      call flush(6)
!     end do
      write(6,*)xlatid,xlonid,timeid
      call flush(6)
      write(6,*)'end definitions ncendf',ncid
      call ncendf(ncid,iret)
      write(6,*)'iret',iret,'nlay',nlay
      istartl=1
      icountl(1)=nlay
      icountl(2)=1
      istart2(3)=1
      write(6,*)'output xlatid',xlatid
      call flush(6)
      call ncvpt(ncid,xlatid,istart2,icount2,alat,ierr)
      print *,'istart2',istart2,'icount2',icount2
      !ierr=nf_put_var_real(ncid,xlatid,alat,istart2,icount2)
      ierr=nf_put_var_real(ncid,xlatid,alat)
      if(ierr.ne.nf_noerr)then
        write(6,*)'error put alat'
        write(6,*)trim(nf_strerror(ierr))
        stop
      endif
      ierr=nf_put_var_real(ncid,xlonid,alon)
      if(ierr.ne.nf_noerr)then
        write(6,*)'error put alon'
        write(6,*)trim(nf_strerror(ierr))
        stop
      endif
      call ncvpt(ncid,dxy_mcid,istart2,icount2,dxy_mc,iret)
     
       call ncvpt(ncid,szid,istart2,icount2,sz,iret)
      write(6,*)'did sz'
      call flush(6)
!     end netcdf stuff here for each file
      istart3(4)=1
      write(6,*)'itsave',itsave,'timeid',timeid,'ncid',ncid
      call flush(6)
!      iret=nf_put_var1_int(ncid,timeid,6,itsave)
      write(6,*)'did it'
      call flush(6)
!    endif ! end save
!   output pressure
!     do k=1,nlay
!        write(6,*)'out psol ',k,maxval(p(:,:,k,1))
!     end do
        call ncvpt(ncid,ivarp(1),istart3,icount3,p(1,1,1,1),iret)
        call ncvpt(ncid,pdashid,istart3,icount3,pdash,iret)
!       calculate pdash bar for sigma use 1013mb
!        write(6,*)'nc',nc,nr,nlay
!        call flush(6)
        do k=1,nlay
!          write(6,*)'k',k
!          call flush(6)
          pdashbar=0.0
          do j=1,nr
            do i=1,nc
              pdashbar=pdashbar+pdash(i,j,k,1)
            end do
          end do
          pdashbar=pdashbar/float(nc*nr)
          sigma(k,1)=pdashbar/1013.e2
!          print *,'sigma',k,sigma(k,1)
        end do
        call ncvpt(ncid,sigmaid,istartl,icountl,sigma,iret)
!        write(6,*)'did sigma'
!        call flush(6)
!      endif
!    end do
    if(save)then
!      write(6,*)'psfcid',psfcid,'istart2',istart2,icount2
      call ncvpt(ncid,psfcid,istart2,icount2,psfc,iret)
!      write(6,*)'did psfc',iret
!      call flush(6)
    endif
    
!      call mxmn(o3_diff(1,1,nlay),nc,nr*nlay,'DIFF',it,1,rx,rn,rs)
!     endif
     write(6,*)'no2',maxval(no2),minval(no2)
!    stop

!    do nv=1,nvarin
     nv=1 ! NO2
        call ncvpt(ncid,ivar(nv),istart3,icount3,no2,iret)
!    end do
    nv=5
        call ncvpt(ncid,ivarp(nv),istart3,icount3,temp,iret)
    nv=6
        call ncvpt(ncid,ivarp(nv),istart3,icount3,qv,iret)
    if(save)call ncclos(ncid,iret)
   end do
   stop
  end
      integer function gjuldate(mbdate   )

       implicit none
! ------------------------ input/output variables -----------------
! input
      integer mbdate       ! base date (yymmdd format) of run (e.g., 000901)
      integer mbsec        ! base seconds of base date (e.g., 0)
 
! output
      integer icalday          ! current julian day including fraction
! -----------------------------------------------------------------

! ---------------------- local variables --------------------------
      integer nsecs        ! run time in seconds
      integer ndays        ! ndays + nyears = number of day changes since start
      integer nyears       ! ndays + nyears = number of day changes since start
      integer mcyear       ! current year of current date (0-99)
      integer mcmnth       ! current month of current date (1-12)
      integer mcday        ! current day of current date (1-31)
      integer jday         ! Julian day (1-365)
      integer mbmnth       ! base month (1-12) validity check
      integer mbday        ! base day of base date (1-31) validity check
      integer ndm(12)      ! number of days per month
      integer jdcon(12)    ! convert month index to julian day
      save ndm,jdcon
      data ndm/31,28,31,30,31,30,31,31,30,31,30,31/
      data jdcon/0,31,59,90,120,151,181,212,243,273,304,334/
! -----------------------------------------------------------------

! check validity of input data 
      print *,'mbdate',mbdate
 
      mbmnth = mod(mbdate,10000)/100
      print *,'mbmnth',mbmnth
      mbday = mod(mbdate,100)
      print *,'mbday',mbday
   
      if (mbmnth.lt.1 .or. mbmnth.gt.12) then
         write(6,*)'CALENDR: invalid base month input:',mbmnth
        stop
      end if
      if (mbday.lt.1 .or. mbday.gt.ndm(mbmnth)) then
         write(6,*)'CALENDR: invalid base day of base date input:',mbday
         stop
      end if
      nsecs = 0
 
! first: current day, seconds
 
 
 
 
 
! initialize current year, month, day
 
      mcyear = mbdate/10000 
      mcmnth = mod(mbdate,10000)/100
      mcday = mod(mbdate,100) 
      print *,'mcyear',mcyear,'mcmnth',mcmnth,'mcday',mcday
 
! now loop through months, converting yy, mm, and ddd to yymmdd.
! e.g.: 791235 becomes 800104. 190001370 becomes 19010105.
 
   10 if (mcday.gt.ndm(mcmnth)) then
         mcday = mcday - ndm(mcmnth)
         mcmnth = mcmnth + 1
         if (mcmnth.eq.13) then            ! add a year
            mcyear = mcyear + 1
            mcmnth = 1
         end if
         go to 10
      end if
!     mcdate = mcyear*10000 + mcmnth*100 + mcday
 
! convert current month, day, seconds to Julian day + fraction 
     print *,'mcnmth',mcmnth,'mcday',mcday 
      jday = jdcon(mcmnth) + mcday
      print *,'jday=',jday
      gjuldate = jday 
      return
      end 
  subroutine readsp(filesur,nc,nr,psurf,ihr)
  use netcdf
  implicit none
  integer nc,nr,ncid,idvar,ierr,ioff,ihr
  real(kind=4),dimension(nc,nr,1,1) :: psurf
  integer start(4),count(4),sdate,stime,tstep
  character *(*) filesur
  start=1
  count(1)=nc
  count(2)=nr
  count(3)=1
  count(4)=1
  ierr=nf90_open(filesur,0,ncid)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error reading sur file ',trim(filesur)
    write(6,*)trim(nf90_strerror(ierr))
  endif
! get time attribute to calculate ioff
  ierr=nf90_get_att(ncid,NF90_GLOBAL,'SDATE',sdate)
  ierr=nf90_get_att(ncid,NF90_GLOBAL,'STIME',stime)
  ierr=nf90_get_att(ncid,NF90_GLOBAL,'TSTEP',tstep)
  print *,'sdate',sdate,stime,tstep
!  ioff=ihr-(stime/10000)+1
  ioff=ihr
  print *,'ioff',ioff
  start(4)=ioff
! now getvariables
  ierr=nf90_inq_varid(ncid,'PRSFC',idvar)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting prsfc id file '
    write(6,*)trim(nf90_strerror(ierr))
  endif
  print *,'start',start,'count',count,'id',idvar
  print *,'shape',shape(psurf)
  ierr=nf90_get_var(ncid,idvar,psurf,start,count)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting psurf surf file '
    write(6,*)trim(nf90_strerror(ierr))
  endif
  write(6,*)'filesur',trim(filesur)
  write(6,*)'nc',nc,nr,'ioff',ioff,'psurf366,11',psurf(366,11,1,1)
  ierr=nf90_close(ncid)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error close'
    call flush(6)
    stop
  endif
  return
  end subroutine readsp
  subroutine readpressure(filepres,nc,nr,nl,pres,dens,densj,zh,zf,ta,qv,ihr)
  use netcdf
  implicit none
  integer nc,nr,nl,ncid,idvar,ierr,ioff
  real(kind=4) :: pres(nc,nr,nl,1),zh(nc,nr,nl,1),zf(nc,nr,nl,1)
  real(kind=4) :: dens(nc,nr,nl,1),densj(nc,nr,nl,1),ta(nc,nr,nl,1),qv(nc,nr,nl,1)
  integer start(4),count(4),ihr
  character *(*) filepres
  integer sdate,stime,tstep
  start=1
  count(1)=nc
  count(2)=nr
  count(3)=nl
  count(4)=1
  ierr=nf90_open(filepres,0,ncid)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error reading pres file ',trim(filepres)
    write(6,*)trim(nf90_strerror(ierr))
  endif
! get time attribute to calculate ioff
  ierr=nf90_get_att(ncid,NF90_GLOBAL,'SDATE',sdate)
  ierr=nf90_get_att(ncid,NF90_GLOBAL,'STIME',stime)
  ierr=nf90_get_att(ncid,NF90_GLOBAL,'TSTEP',tstep)
  print *,'sdate',sdate,stime,tstep
!  ioff=ihr-(stime/10000)+1
  ioff=ihr
  print *,'ioff',ioff
  start(4)=ioff
! now getvariables
  ierr=nf90_inq_varid(ncid,'PRES',idvar)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting pres id file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_get_var(ncid,idvar,pres,start,count)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting press surf file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_inq_varid(ncid,'DENS',idvar)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting dens id file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_get_var(ncid,idvar,dens,start,count)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting dens surf file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_inq_varid(ncid,'DENSA_J',idvar)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting densj id file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_get_var(ncid,idvar,densj,start,count)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting densj surf file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_inq_varid(ncid,'ZH',idvar)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting zh id file '
!    write(6,*)trim(nf90_strerror(ierr))
    write(6,*)nf90_strerror(ierr)
  endif
  ierr=nf90_get_var(ncid,idvar,zh,start,count)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting zh file '
    write(6,*)trim(nf90_strerror(ierr))
  endif
  ierr=nf90_inq_varid(ncid,'ZF',idvar)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting zf id file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_get_var(ncid,idvar,zf,start,count)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting zf file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_inq_varid(ncid,'TA',idvar)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting TA id file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_get_var(ncid,idvar,ta,start,count)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting ta file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_inq_varid(ncid,'QV',idvar)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting QV id file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_get_var(ncid,idvar,qv,start,count)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting qv file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_close(ncid)
  end subroutine readpressure
  subroutine readsurface(filesur,nc,nr,lat,lon,sz,msfx2,dx,dy)
  use netcdf
  implicit none
  integer nc,nr,ncid,idvar,rowdimid,coldimid,ierr
  integer start(3),count(3)
  integer nrowin,ncolin
  real(kind=4) :: lat(nc,nr,1),lon(nc,nr,1),sz(nc,nr,1),msfx2(nc,nr,1),dx,dy
  character *(*) filesur
  logical exist
  start=1
  count(1)=nc
  count(2)=nr
  count(3)=1
  inquire(file=filesur,exist=exist)
  if(.not.exist)then
   write(6,*)trim(filesur),' does not exist'
   stop
  endif 
  ierr=nf90_open(filesur,0,ncid)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error reading sur file ',trim(filesur)
    write(6,*)trim(nf90_strerror(ierr))
  endif
  ierr=nf90_get_att(ncid,NF90_GLOBAL,'XCELL',dx)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error dx'
  endif
  ierr=nf90_get_att(ncid,NF90_GLOBAL,'YCELL',dy)
  if(ierr.ne.nf90_noerr)then
     write(6,*)'error dy'
  endif
  print *,'XCELL',dx,dy
  ierr=nf90_inq_dimid(ncid,'ROW',rowdimid)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting rowdimid sur file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_inquire_dimension(ncid,rowdimid,len=nrowin)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting nrow sur file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  if(nr.ne.nrowin)then
    write(6,*)'dimensions dont match sur file ',nr,nrowin
    stop
  endif
  ierr=nf90_inq_dimid(ncid,'COL',coldimid)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting coldimid sur file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_inquire_dimension(ncid,coldimid,len=ncolin)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting ncol sur file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  if(nc.ne.ncolin)then
    write(6,*)'dimensions dont match surfile col',nc,ncolin
    stop
  endif
! now getvariables
  ierr=nf90_inq_varid(ncid,'LAT',idvar)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting LAT id file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_get_var(ncid,idvar,lat,start,count)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting LAT surf file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_inq_varid(ncid,'LON',idvar)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting LON id file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_get_var(ncid,idvar,lon,start,count)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting LON surf file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_inq_varid(ncid,'HT',idvar)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting HT id file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_get_var(ncid,idvar,sz,start,count)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting HT surf file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_inq_varid(ncid,'MSFX2',idvar)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting MSFX2 id file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_get_var(ncid,idvar,msfx2,start,count)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting msfx2 surf file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  print *,'lat',maxval(lat),minval(lat)
  print *,'lon',maxval(lon),minval(lon)
  print *,'sz',maxval(sz),minval(sz)
  ierr=nf90_close(ncid)
  end subroutine readsurface

  subroutine readcmaqno2(file,nc,nr,nl,no2,itflag)
   use netcdf
  implicit none
  integer nc,nr,nl,ncid,idvar,rowdimid,coldimid,laydimid,ierr
  integer nrowin,ncolin,nlayin,itflag(2)
  real(kind=4) :: no2(nc,nr,nl)
  integer start(4),count(4),start3(3),count3(3),yr,jul,mon,day
  character *(*) file
  logical exist
  start=1
  start(4)=1
  count(1)=nc
  count(2)=nr
  count(3)=nl
  count(4)=1
  write(6,*)'open no2 file '
  write(6,*)trim(file)
  call flush(6)
  inquire(file=file,exist=exist)
  if(.not.exist)then
    print *,'file missing ',trim(file)
    stop
  endif
  ierr=nf90_open(file,0,ncid)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error reading no2 file ',trim(file)
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_inq_dimid(ncid,'ROW',rowdimid)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting rowdimid met file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_inquire_dimension(ncid,rowdimid,len=nrowin)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting nrow met file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  if(nr.ne.nrowin)then
    write(6,*)'dimensions dont match metfile ',nr,nrowin
    stop
  endif
  ierr=nf90_inq_dimid(ncid,'COL',coldimid)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting coldimid met file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_inquire_dimension(ncid,coldimid,len=ncolin)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting ncol met file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  if(nc.ne.ncolin)then
    write(6,*)'dimensions dont match metfile col',nc,ncolin
    stop
  endif
  ierr=nf90_inq_dimid(ncid,'LAY',laydimid)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting laydimid met file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_inquire_dimension(ncid,laydimid,len=nlayin)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting nlay met file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  if(nl.ne.nlayin)then
    write(6,*)'dimensions dont match metfile nlay',nl,nlayin
    stop
  endif
! now getvariables
  ierr=nf90_inq_varid(ncid,'NO2',idvar)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting NO2 id met file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_get_var(ncid,idvar,no2,start,count)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting NO2 met file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  write(6,*)'no2',maxval(no2),minval(no2)
  start3=1
  count3(1)=2
  count3(2)=1
  count3(3)=1
  ierr=nf90_inq_varid(ncid,'TFLAG',idvar)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting TFLAG file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  ierr=nf90_get_var(ncid,idvar,itflag,start3,count3)
  if(ierr.ne.nf90_noerr)then
    write(6,*)'error getting itflag met file '
    write(6,*)trim(nf90_strerror(ierr))
    stop
  endif
  print *,'itflagno2',itflag
  yr=itflag(1)/1000
  jul=itflag(1)-yr*1000
  
  call dayfromjul(yr,jul,mon,day)
  print *,'yr',yr,'mon',mon,'day',day
  ierr=nf90_close(ncid)

  end subroutine readcmaqno2
  subroutine dayfromjul(yr,jul,mon,day)
  integer days(12,2)
  data days/31,28,31,30,31,30,31,31,30,31,30,31, &
       31,29,31,30,31,30,31,31,30,31,30,31/ 
  integer yr,jul,mon,day,ii,im
  ileap=mod(yr,4)
  if(ileap.eq.0)then
    ii=2
  else
    ii=1
  endif
  day=jul
  do im=1,12
    if(day <= days(im,ii))then
      mon=im
      return
    endif
    day=day-days(im,ii)
  end do
  write(6,*)'error day',day,'jul',jul,'yr',yr
  stop
  end
