use netcdf
implicit none
integer ncidres,nciddiff,ierr,idno2res,idno2diff,dim3d(4),ncidnew,idno2,idno2outdiff
integer nc,nr,nlay,latdimid,londimid,levdimid,timeid,timedimid,i,j,k
integer xlatid,xlonid,dim2d(2),lenpath,sdate,stime
character *256 cmaqres,cmaqdiff,sigf00,gsidate*10,cmaqdiffdir,gsidir
logical exist
real*4, allocatable,dimension(:,:,:,:) :: no2new,no2old,no2diff,no2newest,no2diff2
real*4, allocatable,dimension(:,:,:) :: alat,alon
real*4 maxabs,no2max
cmaqres=' '
call getenv('CMAQRESTART',cmaqres)
write(6,*)'cmaqres',trim(cmaqres)
gsidate=' '
call getenv('GSIDATE',gsidate)
write(6,*)'gsidate',trim(gsidate)
call flush(6)
if(gsidate.eq.' ')stop
if(cmaqres.eq.' ')stop
cmaqdiffdir=' '
call getenv('CMAQGSIDIFFDIR',cmaqdiffdir)
write(6,*)'cmaqdiffdir',trim(cmaqdiffdir)
if(cmaqdiffdir.eq.' ')stop
inquire(file=cmaqres,exist=exist)
if(.not.exist)then
   write(6,*)'restart file missing'
   stop
endif
lenpath=len_trim(cmaqdiffdir)
if (cmaqdiffdir(lenpath:lenpath).eq.'/')then
  cmaqdiff=trim(cmaqdiffdir)//'incno2.'//trim(gsidate)
else
  cmaqdiff=trim(cmaqdiffdir)//'/incno2.'//trim(gsidate)
endif
gsidir=' '
call getenv('GSIDIR',gsidir)
write(6,*)'gsidir',trim(gsidir)
if(gsidir.eq.' ' )stop
! update cmaqrestart file with no2 and make no2 diff netcdf file
write(6,*)'open cmaq sigf00 for update'
write(6,*)trim(GSIDIR)//trim(gsidate)//'/sigf00'
call flush(6)
ierr=nf90_open(trim(GSIDIR)//trim(gsidate)//'/sigf00',0,ncidnew)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'error oopen sigf00'
  write(6,*)trim(nf90_strerror(ierr))
  stop
endif
ierr=nf90_inq_varid(ncidnew,'NO2',idno2)
write(6,*)'open ncidres ',trim(cmaqres)
call flush(6)

ierr=nf90_open(cmaqres,NF90_WRITE,ncidres)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'error open restart file'
  write(6,*)trim(nf90_strerror(ierr))
endif
ierr=nf90_inq_varid(ncidres,'NO2',idno2res)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'error inq varid ncidres',ncidres,' NO2'
  write(6,*)trim(nf90_strerror(ierr))
  stop
endif
ierr=nf90_inq_dimid(ncidnew,'south-north',latdimid)
if(ierr.ne.NF90_NOERR)then
   write(6,*)'errror latdimid'
   call flush(6)
  write(6,*)trim(nf90_strerror(ierr))
  stop
endif
ierr=nf90_inq_dimid(ncidnew,'west-east',londimid)
ierr=nf90_inq_dimid(ncidnew,'lev',levdimid)
ierr=nf90_inquire_dimension(ncidnew,latdimid,len=nr)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'error latdim'
  write(6,*)trim(nf90_strerror(ierr))
endif
ierr=nf90_inquire_dimension(ncidnew,londimid,len=nc)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'error londim'
  write(6,*)trim(nf90_strerror(ierr))
endif
ierr=nf90_inquire_dimension(ncidnew,levdimid,len=nlay)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'error levdim'
  write(6,*)trim(nf90_strerror(ierr))
endif
write(6,*)'nc',nc,nr,nlay
call flush(6)
allocate (no2new(nc,nr,nlay,1),no2old(nc,nr,nlay,1),no2diff(nc,nr,nlay,1))
allocate (no2newest(nc,nr,nlay,1),no2diff2(nc,nr,nlay,1))
allocate (alon(nc,nr,1),alat(nc,nr,1))
ierr=nf90_inq_varid(ncidnew,'NO2',idno2)
ierr=nf90_get_var(ncidnew,idno2,no2new)
write(6,*)'find no2new ',maxval(no2new),minval(no2new)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'cant get no2new'
  call flush(6)
  stop
endif
ierr=nf90_inq_varid(ncidnew,'XLAT',xlatid)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'cant get xlatid'
  write(6,*)trim(nf90_strerror(ierr))
  stop
endif
ierr=nf90_inq_varid(ncidnew,'XLONG',xlonid)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'cant get xlonid'
  write(6,*)trim(nf90_strerror(ierr))
  stop
endif
  
ierr=nf90_get_var(ncidnew,xlatid,alat)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'cant get alat'
  write(6,*)trim(nf90_strerror(ierr))
  stop
endif
ierr=nf90_get_var(ncidnew,xlonid,alon)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'cant get alon'
  write(6,*)trim(nf90_strerror(ierr))
  stop
endif
ierr=nf90_get_var(ncidres,idno2res,no2old)
no2max=maxval(no2old)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'cant get no2old'
  write(6,*)trim(nf90_strerror(ierr))
  stop
endif
ierr=nf90_get_att(ncidres,NF90_GLOBAL,'SDATE',sdate)
ierr=nf90_get_att(ncidres,NF90_GLOBAL,'STIME',stime)
if(no2max>4.)then
! ppbV input convert to ppmV
  no2old=no2old*.001
  write(6,*)'update ppbV input'
else if(no2max<1.e-6)then
! ppv input convert to ppmV
  no2old=no2old*1.e6
  write(6,*)'update ppV input'
endif
no2diff=no2new-no2old
write(6,*)'no2diff',maxval(no2diff),minval(no2diff),' no2new',maxval(no2new)
if (maxval(abs(no2diff)) > 0.0 ) then
write(6,*)'cmaqdiff',trim(cmaqdiff)
ierr=nf90_create(cmaqdiff,NF90_CLOBBER,nciddiff)
ierr=nf90_def_dim(nciddiff,'south-north',nr,latdimid)
ierr=nf90_def_dim(nciddiff,'west-east',nc,londimid)
ierr=nf90_def_dim(nciddiff,'lev',nlay,levdimid)
ierr=nf90_def_dim(nciddiff,'time',1,timedimid)
dim3d(1)=londimid
dim3d(2)=latdimid
dim3d(3)=levdimid
dim3d(4)=timedimid
ierr=nf90_def_var(nciddiff,'INCNO2',NF90_FLOAT,dim3d,idno2diff)
ierr=nf90_put_att(nciddiff,idno2diff,'coordinates','XLONG XLAT')
ierr=nf90_put_att(nciddiff,idno2diff,'units','ppmV')
ierr=nf90_def_var(nciddiff,'NO2',NF90_FLOAT,dim3d,idno2outdiff)
ierr=nf90_put_att(nciddiff,idno2res,'coordinates','XLONG XLAT')
ierr=nf90_put_att(nciddiff,idno2res,'units','ppmV')
dim2d(1)=londimid
dim2d(2)=latdimid
ierr=nf90_def_var(nciddiff,'XLAT',NF90_FLOAT,dim2d,xlatid)
ierr=nf90_def_var(nciddiff,'XLONG',NF90_FLOAT,dim2d,xlonid)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'cant def xlon'
  write(6,*)'cant get alon'
endif
ierr=nf90_put_att(nciddiff,xlatid,'coordinates','XLONG XLAT')
ierr=nf90_put_att(nciddiff,xlatid,'units','degree_north')
ierr=nf90_put_att(nciddiff,xlonid,'coordinates','XLONG XLAT')
ierr=nf90_put_att(nciddiff,xlonid,'units','degree_east')
ierr=nf90_put_att(nciddiff,NF90_GLOBAL,'SDATE',sdate)
ierr=nf90_put_att(nciddiff,NF90_GLOBAL,'STIME',stime)
ierr=nf90_enddef(nciddiff)
ierr=nf90_put_var(nciddiff,xlatid,alat)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'cant put alat '
  write(6,*)trim(nf90_strerror(ierr))
  stop
!else
!  write(6,*)'alat',maxval(alat),minval(alat)
!  call flush(6)
endif
ierr=nf90_put_var(nciddiff,xlonid,alon)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'cant put alon '
  write(6,*)trim(nf90_strerror(ierr))
  stop
!else
!  write(6,*)'alat',maxval(alat),minval(alat)
!  call flush(6)
endif
ierr=nf90_put_var(nciddiff,idno2diff,no2diff)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'cant put no2diff '
  write(6,*)trim(nf90_strerror(ierr))
  stop
else
  write(6,*)'no2diff',maxval(no2diff),minval(no2diff)
  call flush(6)
endif
ierr=nf90_put_var(nciddiff,idno2outdiff,no2new)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'cant put no2new '
  write(6,*)trim(nf90_strerror(ierr))
  stop
else
  write(6,*)'no2new',maxval(no2new),minval(no2new)
  call flush(6)
endif
ierr=nf90_put_var(ncidres,idno2res,no2new)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'cant put no2new '
  write(6,*)trim(nf90_strerror(ierr))
  stop
endif
ierr=nf90_close(ncidres)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'cant close cmaqres'
  write(6,*)trim(nf90_strerror(ierr))
  stop
endif
ierr=nf90_close(nciddiff)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'cant close cmaqdiff'
  write(6,*)trim(nf90_strerror(ierr))
  stop
endif
! check that did change
ierr=nf90_open(cmaqres,0,ncidres)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'error open restart file'
  write(6,*)trim(nf90_strerror(ierr))
  stop
endif
ierr=nf90_inq_varid(ncidres,'NO2',idno2res)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'error get id restart file'
  write(6,*)trim(nf90_strerror(ierr))
  stop
endif
ierr=nf90_get_var(ncidres,idno2res,no2newest)
if(ierr.ne.NF90_NOERR)then
  write(6,*)'cant get idno2 newest'
  write(6,*)trim(nf90_strerror(ierr))
  stop
endif
no2diff2=no2newest-no2new
maxabs=maxval(abs(no2diff2))
if (maxabs>1.e-20) then
  write(6,*)'error reading no2newest',maxabs
  do k=1,nlay
  do j=1,nr
  do i=1,nc
     if(no2diff2(i,j,k,1)>1.e-10)then
       write(6,*)i,j,k,no2newest(i,j,k,1),no2new(i,j,k,1)
     endif
  end do
  end do
  end do
endif
ierr=nf90_close(ncidres)
endif
stop
end

