      program read_nc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                               R E A D N C
c
c  this program reads netcdf dump file (-v option) into a (very large)
c  array. this program will only work on ncep time/lat/lon grid files.
c  this version reads integer*2 for slp
c
c  input to program:
c 
c  ifile               dumped format from ncdump
c  ofile               dump of data as long sequence of gravsoft grids or binary
c  n1,rmin1,rmax1      1st variable (time): no of pts, min, max
c  n2,rmin2,rmax2      2nd variable (lat):  no of pts, min, max 
c  n3,rmin3,rmax3      3rd variable (lon):  no of pts, min, max
c  imode               1: gravsoft with possible thinning,
c                      2: interpolate icesat-format data 
c  (t1,t2,ithin)       mode=1 only: minimum and maximal time wanted and thinning 
c  (file name)         mode=2 only: file name of icesat-format data
c
c  (c) rf/hsk DNSC, u of c, apr 4, 2006
c  
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc   
      implicit real*8 (a-h,o-z)
      character*72 ifile,ofile,icesatfile
      character*72 txt1,txt2
      integer*2   d(144,73,1460)
      ilondim =     144
      ilatdim =         73
      itdim =              1460
c
      offset = 119765
c
      read(*,'(a)') ifile
      read(*,'(a)') ofile
      read(*,*) nt,t1,t2
      read(*,*) nfi,rfi1,rfi2
      read(*,*) nla,rla1,rla2
      read(*,*) mode
      if (mode.eq.1) read(*,*) tmin,tmax,ithin
      if (ithin.le.1) ithin = 1
      if (mode.eq.2) read(*,'(a)') icesatfile
c
c  check times and write info
c
      if (nfi.le.1.or.nla.le.1) stop '*** nfi or nla is 1'
      dfi = (rfi2-rfi1)/(nfi-1)
      dla = (rla2-rla1)/(nla-1)
      if (nt.eq.1) stop '*** only one time epoch in file'
      dt = (t2-t1)/(nt-1)
c
      write(*,*)
      write(*,*) 'READ_NC grid sizes ',nt,nfi,nla,' total',nla*nfi*nt
      if (nt.gt.itdim) stop '*** increase nt dimension'
      if (nla.gt.ilondim) stop '*** increase nla dimension'
      if (nfi.gt.ilatdim) stop '*** increase nfi dimension'
c
      write(*,5) dt,dt/24
5     format(' data time spacing: ',f8.1,' hr, ',f8.2,' days')
      call t2ncep(1948,1,0.d0,t1948)
      call t2ncep(2003,1,0.d0,t2003)
      write(*,6) t1,t2,t1948,t2003
6     format(' ncep: start     stop         t1948        t2003'
     .,/,4f13.2)
c
      open(10,file=ifile,status='old')
      open(20,file=ofile,status='unknown')
      if (mode.eq.2) open(30,file=icesatfile,status='old')
c
c  skip netcdf header until "data" is met
c 
      read(10,'(a)') txt1
      write(*,*)
      write(*,*) txt1
10    read(10,'(a)') txt1
      if (txt1(1:5).ne.'data:') goto 10
      read(10,*)
      read(10,'(a)') txt2
      write(*,*) 'Name of data read: ',txt2(1:5),', now reading ..'
c
      read(10,*) (((d(i,j,k),i=1,nla),j=nfi,1,-1),k=1,nt) 
      write(*,*) 'first and last data: ',d(1,1,1),d(nla,1,nt)  
c      
c  output gravsoft grids
c  ---------------------
c
      if (mode.eq.1) then
        k = 0
        no = 0
        if (rla1.gt.180) rla1 = rla1-360
        if (rla2.gt.180) rla2 = rla2-360
        do 30 i1 = 1,nt
          t = t1+(i1-1)*dt
          if (t.lt.tmin.or.mod(i1,ithin).ne.0) goto 30
          if (t.gt.tmax) goto 31
c  approximate decimal years
          tyr = (t-t2003)/24/365.25 + 2003
          write(20,15) rfi1,rfi2,rla1,rla2,dfi,dla,t,tyr
          if (i1.eq.1) write(*,15) rfi1,rfi2,rla1,rla2,dfi,dla,t,tyr
15        format(/,2f10.5,2f11.5,2f9.6,f11.0,f9.2)
c
          do i2 = nfi,1,-1
            write(20,16) (d(j,i2,i1),j=1,nla)
16          format(30(/,8i9))
          enddo
          no = no+1
30      continue
31      write(*,*) 'number of gravsoft grids output:  ',no
      else
c
c  read in icesat-format data and output interpolated value
c  --------------------------------------------------------
c  assume all grids are equidistant
c
        open(30,file=icesatfile,status='old')
        write(*,*)
        write(*,*) 'point data read from file: ',icesatfile
        n_ice = 0
        n_skip = 0
c
40      read(30,*,end=49) iday,rlat,rlon,tsec,hg,hp
c
        if (rlon.lt.0) rlon = rlon+360
        n_ice = n_ice+1
        call t2ncep(2003,iday,tsec,tncep)   
c
c  henriettes indices are upper values
c
        ilat = (rlat-rfi1)/dfi + 2
        ilon = (rlon-rla1)/dla + 2
        itime = (tncep-t1)/dt + 2
c 
c  check for values outside interval - allow interpolation one grid 
c  element outside
c
        if (ilat.eq.1) ilat = 2
        if (ilat.eq.nfi+1) ilat = nfi
        if (ilon.eq.1) ilon = 2
        if (ilon.eq.nla+1) ilon = nla
        if (itime.eq.1) itime = 2
        if (itime.eq.nt+1) itime = nt
c
        if (ilat.lt.2.or.ilat.gt.nfi.
     .  or.ilon.lt.2.or.ilon.gt.nla.
     .  or.itime.lt.2.or.itime.gt.nt) then
          slp_ice = 9999.99
          n_skip = n_skip+1
          goto 44
        endif
c
c  slp1-slp4 are the nearest gridpoints
c  Bilinear interpolation at time t-1
c
        slp1t0=d(ilon-1,ilat-1,itime-1)
        slp2t0=d(ilon,ilat-1,itime-1)
        slp3t0=d(ilon,ilat,itime-1)
        slp4t0=d(ilon-1,ilat,itime-1)
c  Bilinear interpolation at time t
        slp1t1=d(ilon-1,ilat-1,itime)
        slp2t1=d(ilon,ilat-1,itime)
        slp3t1=d(ilon,ilat,itime)
        slp4t1=d(ilon-1,ilat,itime)
c
        t=(rlon-(rla1+(ilon-2)*dla))/dla
        u=(rlat-(rfi1+(ilat-2)*dfi))/dfi
c
        slpt0=(1-t)*(1-u)*slp1t0+t*(1-u)*slp2t0+t*u*slp3t0+
     .  (1-t)*u*slp4t0 
        slpt1=(1-t)*(1-u)*slp1t1+t*(1-u)*slp2t1+t*u*slp3t1+
     .  (1-t)*u*slp4t1 
c
c	Linear interpolation in time	
c
        r=(tncep-(t1+(itime-2)*dt))/dt
c
        slp_ice=(1-r)*slpt0+r*slpt1 + offset
c
c  output data - now make lon -180 to 180
c
44      tyr = (tncep-t2003)/24/365.25 + 2003
        if (rlon.ge.180) rlon = rlon-360
c
	pmean = 1013.3
	alfa = -13.1
	IB = alfa*((slp_ice/100) - pmean)
	hout = hp - (IB*1.E-3)
c
	write(20,45) iday,rlat,rlon,tsec,hout,slp_ice,tncep
        if (mod(n_ice,10000).eq.1.or.n_ice.lt.5)
     .  write(*, 45) iday,rlat,rlon,tsec,hout,slp_ice,tncep,tyr
45      format(i5,2f11.5,f11.2,f9.3,f9.0,f11.0,f10.3)
        goto 40
c
49      write(*,*) 'number of points: ',n_ice,', not intp.: ',n_skip
	endif
      end
c
      subroutine t2ncep(iyr,iday,tsec,tncep)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c           
c  Converts time (iyear,jd,sec) to NCEP time in hr 
c  Modified from Henriette Skourup 
c
c  rf/dnsc April 3, 2006
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
	integer iyr,iday,Y1,M,M1,DM,A,B,C,D,JD,H
      real*8 tncep,tsec
c
	Y1 = iyr-1
	M = 1
	M1 = M + 12
	DM = 1
      A = INT(Y1/100)
      B = 2 - A + INT(A/4)
	C = INT(365.25*Y1)
	D = INT(30.6001*(M1+1))
	JD = B + C + D + DM - 63 - 365 - 1
      H = JD*24
	tncep = iday*24 + tsec/3600 + H
      return
      end
      