
      PROGRAM intloc
      implicit none
      include 'netcdf.inc'



C VERSION 3 LITTLE HEADER
      CHARACTER*24 START_DATE
      CHARACTER*9 NAME

      integer yd,md,dd,hd

      CHARACTER*60 lname,nname/'intloc.nml'/,oname,fpm10,fso2
      
      
      integer nvar,numvar
      parameter (numvar=4)

      integer imx,jmx, tmx
      parameter(imx=500,jmx=1000,tmx=200)
      real lonx(jmx),latx(imx) 
      real Allvar(imx*jmx*tmx)




      integer dt,dt2
      real pm2p5hr(imx*jmx*tmx)
      real pm10hr(imx*jmx*tmx)
      real so2hr(imx*jmx*tmx)
      real no2hr(imx*jmx*tmx)


      integer times(tmx),times2(tmx),tstart
      double precision times_(tmx)
      integer intreg
      real iarr(imx,jmx)



      integer i,j,k,im,jm, time,Alltime,Alltime2

      integer maxloc
      parameter(maxloc=200000)
      real oarr(maxloc)
      integer numloc

      integer snum(maxloc), tsnum
      real slat(maxloc),slono(maxloc),slon(maxloc), tslat, tslon
      real shgt(maxloc), tshgt
      character*19 sname(maxloc), tsname


      integer maxinp
      parameter(maxinp=4)

      integer numinp(maxloc)
      integer posinpI(maxinp,maxloc)
      integer posinpJ(maxinp,maxloc)
      real valinpW(maxinp,maxloc)


      integer stat


      integer iunit/10/, ounit/11/

      real vars(maxloc,numvar)

      character*10 names(numvar)



      integer nlon,nlat


      real start_time, finish_time



      real scale, addoffset
      logical fileexist
      namelist/nlist/ lname,oname,fpm10,fso2,intreg

      print *
      print *, 'Program INTLOC_PM'
      print *
      print *,'Reads locations of the meteorological stations and'
      print *,'the prepared files with PM2.5, PM10, SO2 and NO2'
      print *,'and interpolates the data to the output locations.'
      print *, 'Produces text files'
      print*
      print *, 'Namelist file: ', nname
      print *


      call cpu_time(start_time)


C###############################################################
c Read namelist
C###############################################################



      names(1:4) =(/'PM25','PM10','SO2 ','NO2 '/)






      open(iunit,file=nname,form='formatted',status='old')
!     open(iunit,file='gfsToM3D.nml',form='formatted',status='old')


      read(iunit,nlist)
      close(iunit)



      print *, 'Stations file:  ', lname
      print *, 'Interpolation data file:    ', oname
      print *, 'pm2p5 and pm10 file:    ', fpm10
      print *, 'so2 and no2 file:    ', fso2
      print*
      if (intreg.eq.1) then
         print*,'Interpolation: inverse distance interpolation'
      else if (intreg.eq.2) then
         print*,'Interpolation: nearest point value'
      else if (intreg.eq.3) then
         print*,'Interpolation: maximal value from four nearest points'
      else
         print*,'Wrong intreg value. Intreg have to be 1,2, or 3'
         print*,'current intreg=',intreg
         stop
      endif
      print*
      print * ,'List of all variables:'
      print *,names
      print*
C###############################################################




C###############################################################
C     checking or  exist files
C###############################################################
      dt=1
      dt2=1
      INQUIRE(FILE=fpm10, EXIST=fileexist)
       if (.not. fileexist ) then
          dt=0
          print*, 'There is no file with PM10 and PM2.5!!!!!'
       endif

      INQUIRE(FILE=fso2, EXIST=fileexist)
       if (.not. fileexist ) then
          dt2=0
          print*, 'There is no file with SO2 and NO2!!!!!!'
       endif


      if (dt.eq.0 .and. dt2.eq.0) then
          stop
       endif
C###############################################################








C###############################################################
c Read output locations, heights and names
C###############################################################
      open(iunit,file=lname,form='formatted',status='old')
      numloc = 0
300   if ( numloc .gt. maxloc ) then
          print *, 'Too many stations! Increase MAXLOC.'
          stop
      endif
      read(iunit,*,IOSTAT=stat) tsnum,tslat,tslon,tshgt,tsname
      if ( stat .eq. 0 ) then
          numloc = numloc + 1
          snum(numloc) = tsnum
          slat(numloc) = tslat
          slono(numloc) = tslon
          slon(numloc) = tslon
          if (slon(numloc).lt.0) slon(numloc)=slon(numloc)+360
          shgt(numloc) = tshgt
          sname(numloc) = tsname
          go to 300
      endif
      close(iunit)

      do i = 1,numloc
          call blanks_(sname(i))
      enddo

      print *, 'Number of stations:', numloc
      print *




C###############################################################
c Read dimensions
C###############################################################
      Alltime=0
      Alltime2=0
        if (dt.gt.0) then
          call  NETCDFreadDim(fpm10,'pm2p5',nlon,nlat,Alltime,4)
        endif

        if (dt2.gt.0) then
          call  NETCDFreadDim(fso2,'so2',nlon,nlat,Alltime2,5)
        endif

      print *,'Dimensions from wrfout:'
      print *,'west_east = ',nlon
      print *,'south_north = ',nlat
      print *,'Time steps for PM10 = ',Alltime
      print *,'Time steps for SO2 = ',Alltime2


      if ( nlat*nlon.gt.imx*jmx ) then
          print *, 'Big header dimensions exceed compiled-in.'
          print *, 'BIG HEADER: ',nlat,nlon
          print *, 'IMX,JMX:    ', imx,jmx
          stop
      endif






C###############################################################
c Read input coordinates
C###############################################################

      if (dt.gt.0) then 
           call  NETCDFreadRealVar(fpm10,'longitude',
     +nlon,1,1,AllVar)
           lonx(1:nlon)=AllVar(1:nlon)

           call  NETCDFreadRealVar(fpm10,'latitude',
     +     1,nlat,1,AllVar)
           latx(1:nlat)=Allvar(1:nlat)
      else if (dt2.gt.0) then
           call  NETCDFreadRealVar(fso2,'longitude',
     +      nlon,1,1,AllVar)
           lonx=AllVar(1:nlon)

           call  NETCDFreadRealVar(fso2,'latitude',
     +      1,nlat,1,AllVar)
           latx=Allvar(1:nlat)
      endif

C###############################################################










C################################################################
c For four nearest points
C################################################################
      call FillInputPositions( slat,slon,numloc,
     +                latx(1:nlat),lonx(1:nlon),nlon,nlat, maxinp,
     +                numinp,posinpI,posinpJ,
     +                valinpW)
      print *, 'Building coefficients complete.'
      print *
C################################################################





C################################################################
C     reading of time from files and comparison of time in files
C################################################################
c    from     1960-01-01 00:00:00.0 (previously from 1900-01-01 !IK2024)
      if (dt.gt.0) then
          call  NETCDFreadTimes(fpm10,Alltime,times,times_)
          tstart=times(1)
          dt=times(Alltime)-tstart+1
          write(*,*)'dt of pm10=',dt

      endif
      if (dt2.gt.0) then
          call  NETCDFreadTimes(fso2,Alltime2,times2,times_)
          tstart=times2(1)
          dt2=times2(Alltime2)-tstart+1
          write(*,*)'dt2 of s02=',dt2
      endif

      if ((times(1) .ne. times2(1)).and. (dt.gt.0).and.(dt2.gt.0)) then
        print *, 'Different start time for PM10 and for SO2!!!!!!!!'

!IK2024          call mydate('1900','01','01','00',times(1),yd,md,dd,hd)
           call mydate('1970','01','01','00',times(1),yd,md,dd,hd)
           write(start_date,'(I4,A1,I0.2,A1,I0.2,A1,I0.2,A6)')
     + yd,'-',md,'-',dd,'-',hd,':00:00'
        print *, 'for PM10:  ',start_date

!IK2024          call mydate('1900','01','01','00',times2(1),yd,md,dd,hd)
           call mydate('1970','01','01','00',times2(1),yd,md,dd,hd)
           write(start_date,'(I4,A1,I0.2,A1,I0.2,A1,I0.2,A6)')
     + yd,'-',md,'-',dd,'-',hd,':00:00'
        print *, 'for SO2:   ',start_date
      stop
      endif
C################################################################








C################################################################
C     reading of PM2.5 data and time interpolation (to 1 hour time step)
C################################################################

      if (dt.gt.0) then
          call  NETCDFreadRealVar(fpm10,'pm2p5',
     +    nlon,nlat,Alltime,AllVar)

!IK2024          call NETCDFreadAt(fpm10,'pm2p5',scale,addoffset)
          scale=1.
          addoffset=0.
          AllVar(1:nlon*nlat*Alltime)=
     +    (AllVar(1:nlon*nlat*Alltime)*scale+addoffset)*1.e6



          call timeLinearInterpolation(nlon,nlat,Alltime,
     +    Allvar(1:nlon*nlat*Alltime),times,pm2p5hr,dt)
      endif


C################################################################
C     reading of PM10 data and time interpolation (to 1 hour time step)
C################################################################

      if (dt.gt.0) then
          call  NETCDFreadRealVar(fpm10,'pm10',
     +    nlon,nlat,Alltime,AllVar)

!IK2024          call NETCDFreadAt(fpm10,'pm10',scale,addoffset)
          scale=1.
          addoffset=0.
          AllVar(1:nlon*nlat*Alltime)=
     +    (AllVar(1:nlon*nlat*Alltime)*scale+addoffset)*1.e6



          call timeLinearInterpolation(nlon,nlat,Alltime,
     +    Allvar(1:nlon*nlat*Alltime),times,pm10hr,dt)
      endif



C################################################################
C     reading of SO2 data and time interpolation (to 1 hour time step)
C################################################################

      if (dt2.gt.0) then
          call  NETCDFreadRealVar(fso2,'so2',
     +    nlon,nlat,Alltime2,AllVar)

!IK2024          call NETCDFreadAt(fso2,'so2',scale,addoffset)
          scale=1.
          addoffset=0.
          AllVar(1:nlon*nlat*Alltime2)=
     +    (AllVar(1:nlon*nlat*Alltime2)*scale+addoffset)*1.e6



          call timeLinearInterpolation(nlon,nlat,Alltime2,
     +    Allvar(1:nlon*nlat*Alltime2),times2,so2hr,dt2)
      endif



C################################################################
C     reading of NO2 data and time interpolation (to 1 hour time step)
C################################################################

      if (dt2.gt.0) then
          call  NETCDFreadRealVar(fso2,'no2',
     +    nlon,nlat,Alltime2,AllVar)

!IK2024          call NETCDFreadAt(fso2,'no2',scale,addoffset)
          scale=1.
          addoffset=0.
          AllVar(1:nlon*nlat*Alltime2)=
     +    (AllVar(1:nlon*nlat*Alltime2)*scale+addoffset)*1.e6



          call timeLinearInterpolation(nlon,nlat,Alltime2,
     +    Allvar(1:nlon*nlat*Alltime2),times2,no2hr,dt2)
      endif















C################################################################
C     write header
C################################################################
      open(ounit,file=oname,form='formatted')
C      open(ounit,file=oname,form='formatted',
C     + blocksize=1048576, buffered='yes', buffercount=1)
           
      call printhead(ounit,numvar,names)




C################################################################
C     MAIN CALCULATE
C################################################################

      do time = 1,max(dt,dt2)



C        PM2.5
          nvar=1
          name=names(nvar)
          call to2D(pm2p5hr,time,nlon,nlat,dt,iarr)
          call process(iarr,nlon,nlat,oarr,numloc,maxloc,maxinp,
     +     numinp,posinpI,posinpJ,valinpW,name,intreg)
          vars(1:numloc,nvar) = oarr(1:numloc)


C         read PM10
          nvar=2
           name=names(nvar)
          call to2D(pm10hr,time,nlon,nlat,dt,iarr)
          call process(iarr,nlon,nlat,oarr,numloc,maxloc,maxinp,
     +     numinp,posinpI,posinpJ,valinpW,name,intreg)
          vars(1:numloc,nvar) = oarr(1:numloc)








C         read SO2
          nvar=3
           name=names(nvar)
          call to2D(so2hr,time,nlon,nlat,dt2,iarr)
          call process(iarr,nlon,nlat,oarr,numloc,maxloc,maxinp,
     +     numinp,posinpI,posinpJ,valinpW,name,intreg)
          vars(1:numloc,nvar) = oarr(1:numloc)


C     NO2
          nvar=4
           name=names(nvar)
          call to2D(no2hr,time,nlon,nlat,dt2,iarr)
          call process(iarr,nlon,nlat,oarr,numloc,maxloc,maxinp,
     +     numinp,posinpI,posinpJ,valinpW,name,intreg)
          vars(1:numloc,nvar) = oarr(1:numloc)







!IK2024        call mydate('1900','01','01','00',tstart+time-1,yd,md,dd,hd)
        call mydate('1970','01','01','00',tstart+time-1,yd,md,dd,hd)
        write(start_date,'(I4,A1,I0.2,A1,I0.2,A1,I0.2,A6)')
     +     yd,'-',md,'-',dd,'-',hd,':00:00'





C      Print data to loc file
C***************************************************************
          call printdata(vars,maxloc,names,
     +numvar,snum,slat,slono,shgt,sname, numloc,
     +start_date,time,ounit)

            print *, 'Done with time    ', start_date 
        enddo


C################################################################
C     final print
C################################################################
      print *
      print *, 'All done.'
      print *
      print *
      call cpu_time(finish_time)
      print *, 'Time of calculated  = ', finish_time - start_time

      END


C################################################################
C     subroutine process: from old version intloc.F
C     main subroutine for interpolate of variable
C################################################################

      subroutine process(iarr,im,jm,oarr,numloc,maxloc,maxinp,
     +  numinp,posinpI,posinpJ,valinpW,name,intreg)


      implicit none
      integer maxloc,maxinp
      integer numloc,im,jm
      real iarr(im,jm)
      real oarr(numloc)
      integer numinp(maxloc)
      integer posinpI(maxinp,maxloc)
      integer posinpJ(maxinp,maxloc)
      real valinpW(maxinp,maxloc)

      character*9 name
      integer pos, i,j,n,intreg

      real sum,suminpW

          do n = 1,numloc
              sum = 0
              suminpW = 0
             if (intreg.eq.1) then
                  do pos = 1,numinp(n)
                      i = posinpI(pos,n)
                      j = posinpJ(pos,n)

                      sum = sum + iarr(j,i) * valinpW(pos,n)
                      suminpW = suminpW + valinpW(pos,n)
                  end do
                  oarr(n) = sum / suminpW
             else if (intreg.eq.2) then
                  oarr(n)=iarr(posinpJ(1,n),posinpI(1,n))




             else if (intreg.eq.3) then
              oarr(n)=max(
     +               iarr(posinpJ(1,n),posinpI(1,n)),
     +               iarr(posinpJ(2,n),posinpI(2,n)),
     +               iarr(posinpJ(3,n),posinpI(3,n)),
     +               iarr(posinpJ(4,n),posinpI(4,n))
     +                   )
             endif

          end do

      end


C################################################################
C     subroutine to2D(pm2p5hr,time,nlon,nlat,dt,iarr)
C               get 2D array from all data
C################################################################

      subroutine to2D(var1hr,time,nlon,nlat,dt,iarr)
      integer time,nlon,nlat,dt
      real iarr(nlon,nlat)
      real var1hr(nlon,nlat,dt)
        if (time.le.dt) then
         iarr=var1hr(:,:,time)
        else
         iarr=-9999
        endif


      end




C################################################################
C     subroutine fillvar: from old version intloc.F
C################################################################

      subroutine fillvar(vars,maxloc,numvar, nvar,
     +                    oarr,numloc)
      implicit none
      integer maxloc, numvar
      real vars(maxloc,numvar)
      integer nvar
      integer numloc
      real oarr(numloc)
      integer loc
      do loc = 1,numloc
          vars(loc,nvar) = oarr(loc)
      end do
      end




C################################################################
C     subroutine printhead: from old version intloc.F
C################################################################

      subroutine printhead(ounit,numvar,names)
      implicit none
      integer ounit, i, numvar
      character*10 names(numvar)
      write(ounit,'(A20,A10,2A15,A10,A25$)')
     + ' DATETIME           ', '     SNUM ', '       SLAT    ',
     + '       SLON    ',
     + '      SHGT','     SNAME              '
      do i = 1,numvar
          write(ounit,'(A15$)') names(i)
      end do
      write(ounit,*)
      end

C################################################################
C     subroutine printhead: from old version intloc.F
C                        write data to output file
C################################################################

      subroutine printdata( vars,maxloc,names,
     +numvar,snum,slat,slon,shgt,sname,numloc,start_date,time,ounit )
      implicit none
      integer maxloc,numvar,numloc,time,ounit
      real vars(maxloc,numvar)
      character*10 names(numvar)
      integer snum(numloc)
      real slat(numloc),slon(numloc)
      real shgt(numloc)
      character*24 start_date
      character*19 sname(numloc)
      character*15 str2
      character*155 str
      character*155 strall(numloc)
                  


      integer loc, var,ii

      do loc = 1,numloc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
          if ( vars(loc,1) .ne. -999999. ) then
              write(str,'(A20,I10,2G15.5,I10,A25$)')start_date(1:19),
     +          snum(loc),slat(loc),slon(loc),int(shgt(loc)),sname(loc)
              do var = 1,numvar
                  write(str2,'(G15.7$)') vars(loc,var)
                   ii=95+15*(var-1)
                   str=str(1:ii)//str2
              enddo
      write(ounit,'(A155$)') str
      write(ounit,*)
          endif
ccccccccccccccccccccccccccccccccccccccccccccccccccccccck

      enddo



      end












C################################################################
C     subroutine blanks: from old version intloc.F
C################################################################
      subroutine blanks_(name)
      implicit none
      character name(19)
      integer i
      do i = 1,19
          if ( name(i) .eq. ' ' ) name(i) = '_'
      enddo
      end

C################################################################
C     function dist2: from old version intloc.F
!C
!C exact formula for computation of distances (in km) using geographical coordinates
!C in the assumption that the earth an ideal sphere with Radius=Rearth!
C################################################################

      real function dist22(x1,y1,x2,y2)
      implicit none
      real x1,x2,y1,y2,Pi,Rearth

      Pi=3.14159265
      Rearth=6380.
      dist22=(Rearth*ACos(Cos((Pi*x1)/180.)*Cos((Pi*x2)/180.)*
     + Cos((Pi*(y1 - y2))/180.)+Sin((Pi*x1)/180.)*Sin((Pi*x2)/180.)))
      dist22=dist22*dist22

      return
      end

C################################################################
C     function dist2: from old version intloc.F
C################################################################

      real function dist2(x1,y1,x2,y2)
      implicit none
      real x1,x2,y1,y2,dy
      dy=y1-y2

      if (dy.lt.-180.) dy=dy+360.
      if (dy.gt.180.) dy=dy-360.
      dist2=(x1-x2)*(x1-x2)+dy*dy

      return
      end




C################################################################
C     subroutine FillInputPositions: from old version intloc.F
C################################################################

      subroutine FillInputPositions( slat,slon,numloc,
     +                latx,lonx,nlon,nlat, maxinp,
     +                numinp,posinpI,posinpJ,
     +                valinpW)


                                         
       implicit none
       integer numloc
       real slat(numloc),slon(numloc)
       integer nlon,nlat
       real lonx(nlon),latx(nlat)
       integer numinp(numloc),numinp1(numloc)
       integer maxinp
                                                             
       integer posinpI(maxinp,numloc)
       integer posinpJ(maxinp,numloc)
                                                                 
       real valinpW(maxinp,numloc)





      real dist2
      external dist2
                  
      integer i,j,n,k,ii,jj
      integer pos,pos1
                              
      real ilat,ilon,olat,olon
                                    
      real rbig2, rlit2, w
      real d2,r1,rtmp(4)
      integer p1

                                      


c
c Find the corresponding output points
c
      do n = 1, numloc



       i=CEILING((slat(n) - latx(1))/(latx(2)-latx(1)))
       j=CEILING((slon(n) - lonx(1))/(lonx(2)-lonx(1)))



       posinpI(1,n) = i
       posinpJ(1,n) = j

       posinpI(2,n) = i+1
       posinpJ(2,n) = j

       posinpI(3,n) = i
       posinpJ(3,n) = j+1


       posinpI(4,n) = i+1
       posinpJ(4,n) = j+1


       rtmp(1)=dist2(slat(n),slon(n),latx(i),lonx(j))
       rtmp(2)=dist2(slat(n),slon(n),latx(i+1),lonx(j))
       rtmp(3)=dist2(slat(n),slon(n),latx(i),lonx(j+1))
       rtmp(4)=dist2(slat(n),slon(n),latx(i+1),lonx(j+1))



      if (j.eq.nlon) then
       posinpJ(3,n) = 1
       posinpJ(4,n) = 1
      end if

       numinp(n) = 4

!C    Sorting of interpolation points according distance
       do j=1,4
       r1=1000
       do i=j,4
       if ((rtmp(i).le.r1)) then
       r1=rtmp(i)
       rtmp(i)=rtmp(j)
       rtmp(j)=r1
      
       p1= posinpI(i,n)
       posinpI(i,n)= posinpI(j,n)
       posinpI(j,n)=p1
      
       p1= posinpJ(i,n)
       posinpJ(i,n)= posinpJ(j,n)
       posinpJ(j,n)=p1

       endif
       enddo




          enddo










c       if (mod(n,20).eq.0) write(*,'(A)',advance='no') '*'



      enddo
c      stop
        


c
c Build coefficients
c
      do n = 1, numloc
        olat = slat(n)
        olon = slon(n)
          do pos = 1, numinp(n)
            i = posinpI(pos,n)
            j = posinpJ(pos,n)
            ilat = lonx(j)
            ilon = latx(i)
            rlit2 = dist2(ilat,ilon,olat,olon)+1.e-6
            w = 1 / rlit2
            valinpW(pos,n) = w
           end do
        end do
      end



C################################################################
C     subroutine NETCDFreadDim: find dimension array for variable
C                               from name NetCDF file and name var
C################################################################
!IK2024      subroutine NETCDFreadDim(namef,namevar,nlat,nlon,Alltime)
      subroutine NETCDFreadDim(namef,namevar,nlat,nlon,Alltime,ind_time)
      implicit none
      include 'netcdf.inc'
      character*(*) namef,namevar
      character*(20) dname
      integer retval,ncid,varid
      integer nlat,nlon,Alltime,did(5)
      integer ind_time !IK2024  time index is different for single level
      ! and multi-level fields
            
      if(ind_time>5)then
              write(*,*)'ind_time=',ind_time,'is too large'
              pause
      endif

      did=0
                                                                   
                                                        
C     Open the file..
      retval = nf_open(namef, NF_NOWRITE, ncid)
      if (retval .ne. nf_noerr) then
      print *,namef
      call handle_err(retval)
      endif
                                                          
C
C     Get the varid of the data variable, based on its name.
      retval = nf_inq_varid(ncid, namevar, varid)
      if (retval .ne. nf_noerr) call handle_err(retval)
            
C     Get the dimension id of the data variable.
      retval = nf_inq_vardimid(ncid,  varid,did)
      if (retval .ne. nf_noerr) call handle_err(retval)
                        
C     Get the dimension wrfout
      retval = nf_inq_dimlen(ncid, did(1), nlat)
      if (retval .ne. nf_noerr) call handle_err(retval)
      retval = nf_inq_dimlen(ncid, did(2), nlon)
      if (retval .ne. nf_noerr) call handle_err(retval)
      !IK2024 retval = nf_inq_dimlen(ncid, did(3), Alltime)
      retval = nf_inq_dimlen(ncid, did(ind_time), Alltime)
      if (retval .ne. nf_noerr) call handle_err(retval)
                                                



                              
C     Close the file, freeing all resources.
      retval = nf_close(ncid)
      if (retval .ne. nf_noerr) call handle_err(retval)
            
            
      end
                  




C################################################################
C     subroutine NETCDFreadAt: 
C                              
C################################################################
      subroutine NETCDFreadAt(namef,namevar,scale,addoffset)
      implicit none
      include 'netcdf.inc'
      character*(*) namef
      character*(*) namevar
      real scale,addoffset
      integer ncid,varid,retval,did

C     Open the file..
      retval = nf_open(namef , NF_NOWRITE, ncid)
      if (retval .ne. nf_noerr) then 
      print *,namef
      call handle_err(retval)
      endif

C     Get the varid of the data variable, based on its name.
      retval = nf_inq_varid(ncid, namevar, varid)
      if (retval .ne. nf_noerr) call handle_err(retval)

            
C     Get the attribute value.
      retval = nf_get_att_real(ncid,varid,'scale_factor',scale)
      if (retval .ne. nf_noerr) call handle_err(retval)
      retval = nf_get_att_real(ncid,varid,'add_offset',addoffset)
      if (retval .ne. nf_noerr) call handle_err(retval)
      

      
C     Close the file, freeing all resources.
      retval = nf_close(ncid)
      if (retval .ne. nf_noerr) call handle_err(retval)
      end




C################################################################
C     subroutine NETCDFreadRealVar: get 2D array {arr(x,y,t)} from
C                                   name NetCDF file, name var and
C                                   position in list (t)
C################################################################
      subroutine NETCDFreadRealVar(namef,namevar,
     +iLen,jLen,tLen,arrray)
      implicit none
      include 'netcdf.inc'
      integer iLen,jLen,tLen,retval,ncid,varid,ii,jj,imax,jmax
      character*(*) namef
      character*(*) namevar
      real arrray(ilen,jLen,tLen)




C     Open the file. NF_NOWRITE tells netCDF we want read-only access to
C     the file.
      retval = nf_open(namef, NF_NOWRITE, ncid)
      if (retval .ne. nf_noerr) call handle_err(retval)
      

C     Get the varid of the data variable, based on its name.
      retval = nf_inq_varid(ncid, namevar, varid)





C     Read the data.
      retval = nf_get_var_real(ncid, varid, arrray)
      if (retval .ne. nf_noerr) call handle_err(retval)      
             

                  


C     Close the file, freeing all resources.
210      retval = nf_close(ncid)
      if (retval .ne. nf_noerr) call handle_err(retval)



      end





C################################################################
C     subroutine timeLinearInterpolation(nlon,nlat,Alltime,Allvar)
C            interpolating any data to time_step = 1 hour
C################################################################
      subroutine timeLinearInterpolation(nlon,nlat,
     +ntime,Allvar,times,var1hr,dt)
      integer  nlon,nlat,ntime,dt
      integer times(ntime)
      real Allvar(nlon,nlat,ntime)
      real var1hr(nlon,nlat,dt)

      real y2(nlon,nlat),y1(nlon,nlat)
      integer t2,t1

      integer t,k,step



      do step=1,ntime-1
        t1=times(step)
        t2=times(step+1)
        y1=Allvar(:,:,step)
        y2=Allvar(:,:,step+1)
        do t=0, t2-t1
            k=1+t1-times(1)+t
            var1hr(:,:,k)=y1+(y2-y1)*t/(t2-t1)
        enddo
      enddo

      end








      subroutine NETCDFreadTimes(namef,tLen,times,times_)
      implicit none
      include 'netcdf.inc'
      integer tLen,tid,retval,ncid,varid
      character*(*) namef
      integer times(tLen),i
      double precision times_(tLen)
      
      times_=0

C     Open the file. NF_NOWRITE tells netCDF we want read-only access to
C     the file.
      retval = nf_open(namef, NF_NOWRITE, ncid)
      if (retval .ne. nf_noerr) call handle_err(retval)
      

C     Get the varid of the data variable, based on its name.
!IK2024      retval = nf_inq_varid(ncid, 'time', varid)
       retval = nf_inq_varid(ncid, 'valid_time', varid) !IK2024


C     Read the data.
!IK2024      retval = nf_get_var_int(ncid, varid, times)
      retval = nf_get_var_double(ncid, varid, times_)
      do i=1,tLen
        times(i)=IDINT(times_(i)/3600.)
        write(*,*)'times(',i,')=',times(i)
      enddo

      if (retval .ne. nf_noerr) call handle_err(retval)
             

                  


C     Close the file, freeing all resources.
      retval = nf_close(ncid)
      if (retval .ne. nf_noerr) call handle_err(retval)



      end





C################################################################
C     subroutine handle_err: for identefication error in libriary
C                            "netcdf.inc"
C################################################################

      subroutine handle_err(errcode)
      implicit none
      include 'netcdf.inc'
      integer errcode
                            
      print *, 'Error: ', nf_strerror(errcode)
      stop
      end







C#############################################################

      subroutine bred(a,im,jm)
      integer im,jm
      real a(im,jm)
      print *, a
      end




! subroutine mydate takes input arguments: cry, crm, crd, crh, chh - year,
! month, day, hour, offset hour (which is to be added to date) and returns.
! resulting ry, rm, rd, rh
! Author: A. Kuschan, adapted by I. Kovalets
      subroutine mydate(cry,crm,crd,crh,hh,ry,rm,rd,rh)
      implicit none
!C      external iargc
!C      integer iargc
            
      character*4 cry
      character*2 crm
      character*2 crd
      character*2 crh
                                    

                                          
      integer ry,rm,rd,rh,hh
                                                
      integer mo(12)/31,28,31,30,31,30,31,31,30,31,30,31/
                                                      
      logical leap/.false./
                                                            
      integer i,shift
                                                                  


      read(cry,*) ry
      read(crm,*) rm
      read(crd,*) rd
      read(crh,*) rh

                              

      if ( ry.lt.1900 .or. ry.gt.2200 ) then
      write(0,*) 'Cannot deal with year:', ry
      call exit
      endif
                                
      if ( rm.lt.1 .or. rm.gt.12 ) then
      write(0,*) 'Cannot deal with month:', rm
      call exit
      endif

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!c
!c find if this is a leap year
!c
      leap=.false. !IK2024
      if ( mod(ry,100) .ne. 0 ) then
      if ( mod(ry,4).eq.0 ) leap = .true.
      else
      if ( mod(ry,400).eq.0 ) leap = .true.
      endif
      if ( leap ) then

      mo(2) = 29
      else
      mo(2) = 28
      endif
                                                                            
!c
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      if ( rd.lt.1 .or. rd.gt.mo(rm) ) then
      write(0,*) 'Cannot deal with day ', rd, ' of month ', rm
!c         print *, mo(rm)
          call exit
      endif

      if ( rh.lt.0 .or. rh.gt.24 ) then
      write(0,*) 'Cannot deal with hour ', rh
      call exit
      endif
                                                                
!c
!c cycle through hours
!c
      if ( hh .ge. 0 )  then
        do i = 1,hh
                
         rh = rh + 1
        if ( rh.eq.24 ) then
          rh = 0
                                                              
          rd = rd + 1
          if ( rd.gt.mo(rm) ) then
          rd = 1

          rm = rm + 1
          if ( rm.gt.12 ) then
          rm = 1

        ry = ry + 1

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!c
!c find if this is a leap year
!c
      leap = .false.
      if ( mod(ry,100) .ne. 0 ) then
           if ( mod(ry,4).eq.0 ) leap = .true.
          else
            if ( mod(ry,400).eq.0 ) leap = .true.
          endif
        if ( leap ) then
         mo(2) = 29
        else
         mo(2) = 28
        endif
!c
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

       endif
        endif
          endif
!c             print *, ry,rm,rd,rh
        enddo
         else

          do i = -1,hh,-1
      
        rh = rh - 1
        if ( rh.eq.-1 ) then
        rh = 23
                                                        
        rd = rd - 1
        if ( rd.eq.0 ) then

!c
!c find previous date
!c
      if ( rm .eq. 1 ) then
        rm = 12
                                                
      ry = ry - 1

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!c
!c find if this is a leap year
!c
       if ( mod(ry,100) .ne. 0 ) then
        if ( mod(ry,4).eq.0 ) leap = .true.
          else
        if ( mod(ry,400).eq.0 ) leap = .true.
          endif
        if ( leap ) then
          mo(2) = 29
        else
         mo(2) = 28
        endif

!c
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

          else
        rm = rm - 1
          endif
                                                                      
        rd = mo(rm)
          endif
        endif
!c             print *, ry, rm, rd, rh
        enddo
        endif

!C      print '(I4,3I3.2)', ry, rm, rd, rh

      end
                                                                                                                                                  






