[38/51] [partial] incubator-sdap-nexus git commit: SDAP-1 Import all code under the SDAP SGA

[le...@apache.org]

http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/gaussInterp_f.f
----------------------------------------------------------------------
diff --git a/climatology/clim/gaussInterp_f.f b/climatology/clim/gaussInterp_f.f
new file mode 100644
index 0000000..9366870
--- /dev/null
+++ b/climatology/clim/gaussInterp_f.f
@@ -0,0 +1,219 @@
+c
+c     gaussInterp routine -- Gaussian weighted smoothing in lat, lon, and time
+c
+c Based on Ed Armstrong's routines. Designed to be called from python using f2py.
+c
+c
+c Gaussian weighting = exp( vfactor * (((x - x0)/sx)^2 + ((y - y0)/sy)^2 + ((t - t0)/st)^2 ))
+c
+c where deltas are distances in lat, lon and time and sx, sy, st are one e-folding sigmas.
+c
+c Cutoffs for neighbors allowed in the interpolation are set by distance in lat/lon (see wlat/wlon);
+c for time all epochs are included.
+
+
+      subroutine gaussInterp_f(var, mask,
+     &                         time, lat, lon,
+     &                         outlat, outlon,
+     &                         wlat, wlon,
+     &                         slat, slon, stime,
+     &                         vfactor,
+     &                         missingValue,
+     &                         verbose,
+     &                         vinterp, vweight, status,
+     &                         ntime, nlat, nlon, noutlat, noutlon)
+
+      implicit none
+      integer*4 ntime, nlat, nlon
+c                 ! prepared 3D array of variable data over time, lon, & lat
+      real*4    var(ntime, nlat, nlon)
+c                 ! variable to be interpolated, over ntime epochs
+      integer*1 mask(ntime, nlat, nlon)
+c                 ! pixel quality mask
+      integer*4 time(ntime)
+c                 ! time epochs to gaussian-interpolate over
+      real*4 lat(nlat)
+                  ! latitude corrdinate vector
+      real*4 lon(nlon)
+                  ! longitude corrdinate vector
+cf2py intent(in) var, mask, time, lat, lon        !! annotations for f2py processor
+      integer*4 noutlat, noutlon
+      real*4 outlat(noutlat), outlon(noutlon)
+c                 ! lat/lon grid to interpolate to
+cf2py intent(in) outlat, outlon
+      real*4 wlat, wlon
+c                 ! window of lat/lon neighbors to gaussian weight, expressed in delta lat (degrees)
+c                 ! if an integer, then it is the number of neighbors on EACH side
+cf2py intent(in) wlat, wlon
+      real*4 slat, slon, stime
+c                 ! sigma for gaussian downweighting with distance in lat, lon, & time
+cf2py intent(in) slat, slon, stime
+      real*4 vfactor
+c                 ! factor in front of gaussian expression
+      real*4 missingValue
+c                 ! value to mark missing values in interp result
+      integer*4 verbose
+c                 ! integer to set verbosity level
+cf2py intent(in) vfactor, missingValue
+      real*4 vinterp(noutlat, noutlon)
+c                 ! interpolated variable using gaussians, missing values not counted
+      real*4 vweight(noutlat, noutlon)
+c                 ! weight of values interpolated (can be zero), if so should become missing value
+      integer*4 status
+c                 ! negative status indicates error
+cf2py intent(out) vinterp, vweight, status
+
+      integer*4 clamp
+      integer*4 imin, imax, jmin, jmax, itmp
+      integer*4 iin, jin, kin
+      integer*4 iMidTime
+      integer*4 i, j, iwlat, iwlon
+
+      real*4 wlat2, wlon2, lat1, lon1, dlat, dlon, fac
+      real*4 varmin, varmax, val
+      real*8 midTime
+      logical*1 nLatNeighbors, nLonNeighbors
+
+      write(6, *) 'Echoing inputs ...'
+      write(6, *) 'ntime, nlat, nlon:', ntime, nlat, nlon
+      write(6, *) 'noutlat, noutlon:', noutlat, noutlon
+      write(6, *) 'wlat, wlon:', wlat, wlon
+      write(6, *) 'slat, slon, stime:', slat, slon, stime
+      write(6, *) 'vfactor:', vfactor
+      write(6, *) 'missingValue:', missingValue
+
+      status = 0
+      if (int(wlat) .eq. wlat) then
+         iwlat = int(wlat)
+         nLatNeighbors = .TRUE.
+         write(6, *) 'Using', iwlat,
+     &      'neighbors on each side in the lat direction.'
+      else
+         nLatNeighbors = .FALSE.
+      end if
+      if (int(wlon) .eq. wlon) then
+         iwlon = int(wlon)
+         nLonNeighbors = .TRUE.
+         write(6, *) 'Using', iwlon,
+     &      'neighbors on each side in the lon direction.'
+      else
+         nLonNeighbors = .FALSE.
+      end if
+
+      wlat2 = wlat / 2.
+      wlon2 = wlon / 2.
+      lat1 = lat(1)
+      lon1 = lon(1)
+      dlat = lat(2) - lat(1)
+      dlon = lon(2) - lon(1)
+      iMidTime = int(ntime/2 + 0.5)
+      midTime = time(iMidTime)
+
+      if (verbose .gt. 3) then
+          write(6, *) 'time:', time
+          write(6, *) 'lat:', lat
+          write(6, *) 'lon:', lon
+          write(6, *) 'outlat:', outlat
+          write(6, *) 'outlon:', outlon
+c          write(6, *) 'mask(iMidTime):', mask(iMidTime,:,:)
+          write(6, *) 'var(iMidTime):', var(iMidTime,:,:)
+      end if
+
+      do i = 1, noutlat
+         if (verbose .gt. 1) write(6, *) outlat(i)
+         do j = 1, noutlon
+            vinterp(i,j) = 0.0
+            vweight(i,j) = 0.0
+            if (verbose .gt. 3) then
+               write(6, *) '(i,j) = ', i, j
+               write(6, *) '(outlat,outlon) = ', outlat(i), outlon(j)
+            end if
+
+            if (nLatNeighbors) then
+               imin = clamp(i - iwlat, 1, nlat)
+               imax = clamp(i + iwlat, 1, nlat)
+            else
+               imin = clamp(int((outlat(i) - wlat2 - lat1)/dlat + 0.5),
+     &                      1, nlat)
+               imax = clamp(int((outlat(i) + wlat2 - lat1)/dlat + 0.5),
+     &                      1, nlat)
+            end if
+            if (imin .gt. imax) then
+c              ! input latitudes descending, so swap
+               itmp = imin
+               imin = imax
+               imax = itmp
+            end if
+
+            if (nLonNeighbors) then
+               jmin = clamp(j - iwlon, 1, nlon)
+               jmax = clamp(j + iwlon, 1, nlon)
+            else
+               jmin = clamp(int((outlon(j) - wlon2 - lon1)/dlon + 0.5),              
+     &                      1, nlon)     
+               jmax = clamp(int((outlon(j) + wlon2 - lon1)/dlon + 0.5),
+     &                      1, nlon)
+            end if
+
+            if (verbose .gt. 2) then
+               write(6, *) '(imin,imax,minlat,maxlat) = ',
+     &               imin, imax, lat(imin), lat(imax)
+               write(6, *) '(jmin,jmax,minlon,maxlon) = ',
+     &               jmin, jmax, lon(jmin), lon(jmax)
+            end if
+
+            if (verbose .gt. 4 .and. i .eq. noutlat/2
+     &             .and. j .eq. noutlon/2) then
+               write(6, *) 'Echoing factors ...'
+            end if
+
+            do kin = 1, ntime
+               varmin = 0.
+               varmax = 0.
+               do iin = imin, imax
+                  do jin = jmin, jmax
+                     if (mask(kin,iin,jin) .eq. 0) then
+                        fac = exp( vfactor *
+     &                         (((outlat(i) - lat(iin))/slat)**2
+     &                        + ((outlon(j) - lon(jin))/slon)**2
+     &                        + ((midTime   - time(kin))/stime)**2))
+
+                        val = var(kin,iin,jin)
+                        if (verbose .gt. 4 .and. i .eq. noutlat/2
+     &                         .and. j .eq. noutlon/2) then
+                           write(6, *) kin, iin, jin, time(kin),
+     &                           lat(iin), lon(jin), val, fac, val*fac
+                        end if
+
+                        vinterp(i,j) = vinterp(i,j) + val * fac
+                        vweight(i,j) = vweight(i,j) + fac
+                        if (verbose .gt. 3) then
+                           varmin = min(varmin, val)
+                           varmax = max(varmax, val)
+                        end if
+                     end if
+                  end do
+               end do
+               if (verbose .gt. 3) then
+                  write(6, *) '(itime, varmin, varmax) = ',
+     &                           kin, varmin, varmax
+               end if
+            end do
+            if (vweight(i,j).ne. 0.0) then
+               vinterp(i,j) = vinterp(i,j) / vweight(i,j)
+            else
+               vinterp(i,j) = missingValue
+            end if
+         end do
+      end do
+      return
+      end
+
+
+      integer*4 function clamp(i, n, m)
+      integer*4 i, n, m
+      clamp = i
+      if (i .lt. n) clamp = n
+      if (i .gt. m) clamp = m
+      end
+

http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/gaussInterp_f.mk
----------------------------------------------------------------------
diff --git a/climatology/clim/gaussInterp_f.mk b/climatology/clim/gaussInterp_f.mk
new file mode 100644
index 0000000..548ce33
--- /dev/null
+++ b/climatology/clim/gaussInterp_f.mk
@@ -0,0 +1 @@
+f2py -c -m gaussInterp_f gaussInterp_f.f

http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/gaussInterp_slow.py
----------------------------------------------------------------------
diff --git a/climatology/clim/gaussInterp_slow.py b/climatology/clim/gaussInterp_slow.py
new file mode 100644
index 0000000..6a7bb8f
--- /dev/null
+++ b/climatology/clim/gaussInterp_slow.py
@@ -0,0 +1,130 @@
+#
+#     gaussInterp_slow routine -- Gaussian weighted smoothing in lat, lon, and time
+#
+# Based on Ed Armstrong's routines. Pure python implementation.
+#
+#
+# Gaussian weighting = exp( vfactor * (((x - x0)/sx)^2 + ((y - y0)/sy)^2 + ((t - t0)/st)^2 ))
+#
+# where deltas are distances in lat, lon and time and sx, sy, st are one e-folding sigmas.
+#
+# Cutoffs for neighbors allowed in the interpolation are set by distance in lat/lon (see dlat/dlon);
+# for time all epochs are included.
+#
+
+import sys
+import numpy as np
+from math import exp
+from numba import jit, int32
+
+VERBOSE = 0
+
+
+def gaussInterp_slow(var,                # bundle of input arrays: masked variable, coordinates
+                varNames,                # list of names in order: primary variable, coordinates in order lat, lon, time
+                outlat, outlon,          # output lat/lon coordinate vectors
+                wlat, wlon,              # window of lat/lon neighbors to gaussian weight, expressed in delta lat (degrees)
+                slat, slon, stime,       # sigma for gaussian downweighting with distance in lat, lon (deg), & time (days)
+                vfactor=-0.6931,         # factor in front of gaussian expression
+                missingValue=-9999.,     # value to mark missing values in interp result
+                verbose=VERBOSE,         # integer to set verbosity level
+                optimization='python'):  # Mode of optimization, using 'fortran' or 'cython' or 'python'
+    '''Gaussian interpolate in lat, lon, and time to a different lat/lon grid, and over a time window to the center time.
+Bundle of arrays (var) contains a 3D masked variable and coordinate arrays for lat, lon, and time read from netdf/hdf files.
+Returns the 2D interpolated variable (masked) and a status for failures. 
+    '''
+    v = var[varNames[0]][:]
+    vmask = np.ma.getmask(v)[:]
+    vtime = var[varNames[1]][:]
+    lat = var[varNames[2]][:]
+    lon = var[varNames[3]][:]
+
+    vinterp, vweight, status = \
+         gaussInterp_(v, vmask, vtime, lat, lon,
+                      outlat, outlon, wlat, wlon, slat, slon, stime, vfactor, missingValue)
+
+    vinterp = np.ma.masked_where(vweight == 0.0, vinterp)
+    return (vinterp, vweight, status)
+
+#@jit(nopython=False)
+def gaussInterp_(var,     # variable & mask arrays with dimensions of time, lon, lat
+                 vmask,         
+                 vtime,    # coordinate vectors for inputs
+                 lat,
+                 lon,     
+                 outlat,  # coordinate vectors for grid to interpolate to
+                 outlon,
+                 wlat, wlon,           # window of lat/lon neighbors to gaussian weight, expressed in delta lat (degrees)
+                 slat, slon, stime,    # sigma for gaussian downweighting with distance in lat, lon (deg), & time (days)
+                 vfactor,              # factor in front of gaussian expression
+                 missingValue):        # value to mark missing values in interp result
+    '''Gaussian interpolate in lat, lon, and time to a different lat/lon grid, and over a time window to the center time.
+Returns the 2D interpolated variable (masked), the weight array, and a status for failures.
+    '''
+    vinterp = np.zeros( (outlat.shape[0], outlon.shape[0]), dtype=var.dtype )  # interpolated variable, missing values not counted
+    vweight = np.zeros( (outlat.shape[0], outlon.shape[0]), dtype=var.dtype )  # weight of values interpolated (can be zero)
+    status = 0     # negative status indicates error
+
+    ntime = vtime.shape[0]
+    nlat = lat.shape[0]
+    nlon = lon.shape[0]
+
+    noutlat = outlat.shape[0]
+    noutlon = outlon.shape[0]
+
+    midTime = vtime[int(ntime/2 + 0.5)]
+    wlat2 = wlat / 2.
+    wlon2 = wlon / 2.
+    lat0 = lat[0]
+    lon0 = lon[0]
+    dlat = lat[1] - lat[0]
+    dlon = lon[1] - lon[0]
+
+    for i in xrange(noutlat):
+        print >>sys.stderr, outlat[i]
+        for j in xrange(noutlon):
+           if VERBOSE: print >>sys.stderr, '\n(i,j) = %d, %d' % (i, j)
+           if VERBOSE: print >>sys.stderr, '\n(outlat,outlon) = %f, %f' % (outlat[i], outlon[j])
+
+           imin = clamp(int((outlat[i] - wlat2 - lat0)/dlat + 0.5), 0, nlat-1)
+           imax = clamp(int((outlat[i] + wlat2 - lat0)/dlat + 0.5), 0, nlat-1)
+           if imin > imax: (imin, imax) = (imax, imin)                            # input latitudes could be descending
+           if VERBOSE: print >>sys.stderr, '(imin, imax) = %d, %d' % (imin, imax)
+           if VERBOSE: print >>sys.stderr, '(minlat, maxlat) = %f, %f' % (lat[imin], lat[imax])
+           jmin = clamp(int((outlon[j] - wlon2 - lon0)/dlon + 0.5), 0, nlon-1)
+           jmax = clamp(int((outlon[j] + wlon2 - lon0)/dlon + 0.5), 0, nlon-1)
+           if VERBOSE: print >>sys.stderr, '(jmin, jmax) = %d, %d' % (jmin, jmax)
+           if VERBOSE: print >>sys.stderr, '(minlon, maxlon) = %f, %f' % (lon[jmin], lon[jmax])
+#           stencil = np.zeros( (ntime, imax-imin+1, jmax-jmin+1) )
+
+           for kin in xrange(ntime):
+               for iin in xrange(imin, imax+1):
+                   for jin in xrange(jmin, jmax+1):
+                       if not vmask[kin,iin,jin]:
+                           fac = exp( vfactor *
+                                     (((outlat[i] - lat[iin])/slat)**2
+                                    + ((outlon[j] - lon[jin])/slon)**2
+                                    + ((midTime   - vtime[kin])/stime)**2))
+#                           stencil[kin, iin-imin, jin-jmin] = fac
+                           val = var[kin, iin, jin]
+                           if VERBOSE > 1: print >>sys.stderr,  kin, iin, jin, vtime[kin], lat[iin], lon[jin], val, fac, val*fac
+
+                           vinterp[i,j] = vinterp[i,j] + val * fac
+                           vweight[i,j] = vweight[i,j] + fac
+
+
+           if vweight[i,j] != 0.0:
+               vinterp[i,j] = vinterp[i,j] / vweight[i,j]
+#               if VERBOSE > 1: print >>sys.stderr, 'stencil:\n', stencil
+#               if VERBOSE: print >>sys.stderr, 'stencil max:\n', np.max(np.max(stencil))
+#               if VERBOSE: print >>sys.stderr, 'stencil min:\n', np.min(np.min(stencil))
+           else:
+               vinterp[i,j] = missingValue
+
+    return (vinterp, vweight, status)
+
+#@jit( int32(int32,int32,int32), nopython=False)
+def clamp(i, n, m):
+    if i < n: return n
+    if i > m: return m
+    return i

http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/interp.f
----------------------------------------------------------------------
diff --git a/climatology/clim/interp.f b/climatology/clim/interp.f
new file mode 100644
index 0000000..8977798
--- /dev/null
+++ b/climatology/clim/interp.f
@@ -0,0 +1,302 @@
+#include "interp.h"
+c $Revision: 2.12 $
+c subprogram interp.f
+c reads hdf file, extract data, call bisum(), calculates
+c and writes out interpolated maps
+
+
+      subroutine interp( temp, sst_inter, sst_seq, sl, 
+     &		  x, y, z, f,
+     &            xx, yy, zz, vsum, ixmin, ixmax,
+     &		  iymin, iymax, izmin, izmax, ndatamax,
+     &		  x_length, y_length, imax, jmax, kmax )
+
+      integer x_length,y_length
+      integer dummy
+
+c-- dimension the arrays
+      real*4 temp(x_length,y_length)
+      real*4 sl(imax,jmax,kmax)
+      real*4 sst_inter(kmax), sst_seq(imax)
+      real*4 x(ndatamax),y(ndatamax),z(ndatamax),f(ndatamax)
+      real*4 xx(imax),yy(jmax),zz(kmax)
+      real*4 vsum(2,imax,jmax,kmax)
+      real*4 ULlon, ULlat, LRlon, LRlat
+      integer*4 ixmin(ndatamax),ixmax(ndatamax)
+      integer*4 iymin(ndatamax),iymax(ndatamax)
+      integer*4 izmin(ndatamax),izmax(ndatamax)
+
+c-- hdf data array must be hardwired, no dync alloc :(
+      byte in_data( XSIZE, YSIZE )
+
+      character*80 plabel
+      character*50 infile
+      character*30 outfile
+      character*150 datafile
+      character*30 basename
+      character*35 tmpfile
+      character*10 mapformat
+      character*15 outformat
+      character*3 cday
+      character*4 cyr
+      character*200 zcatcmd
+
+      common /parms/ imin,jmin,kmin,
+     &               xwin2,ywin2,zwin2,
+     &               xh,yh,zh,
+     &               dx,dy,dz,
+     &               xlo,xhi,ylo,yhi,zlo,
+     &               dxd,dyd
+
+      common /fileio/ infile,outfile,istart,
+     &	              iend,mapformat,outformat	
+
+      common /hdfindex/ icolLeft,icolRight,irowUpper,irowLower
+
+      common /basefile/ basename
+
+c-- open output file
+      open(unit=20,status='new',form=outformat,file=outfile,
+     &     err=1000,iostat=ierr )
+      plabel='sst interpolation'
+
+c-- initialize vsum to '0's
+      do k1=1,2
+        do k2=1,kmax
+          do k3=1,jmax
+            do k4=1,imax
+               vsum(k1,k4,k3,k2)=0.
+            enddo
+          enddo
+        enddo
+      enddo
+
+c-- initialize arrays of long, lat and time intervals
+      do ii=imin,imax
+          xx(ii)=xlo+float(ii-1)*dx  ! Long to interpolate to
+      enddo
+      do jj=jmin,jmax
+          yy(jj)=ylo+float(jj-1)*dy  ! Lat to interpolate to
+      enddo
+      do kk=kmin,kmax
+          zz(kk)=zlo+float(kk-1)*dz  ! Time to interpolate to
+      enddo
+
+      do n=1,ndatamax
+          x(n)=0.0
+          y(n)=0.0
+          z(n)=0.0
+          f(n)=0.0
+	  ixmin(n) = 0
+	  ixmax(n) = 0
+	  iymin(n) = 0
+	  iymax(n) = 0
+	  izmin(n) = 0
+	  izmax(n) = 0
+      enddo
+
+c-- slope (cal) and offset to convert DN to SST
+      cal= 0.15
+      offset = 3.0
+
+c-- Open input file list
+      open(unit=15,status='old',file=infile,access='direct',
+     &  recl=RECSIZE,form='formatted',err=1500,iostat=ierr)
+
+c--  read a infile one record at a time and process . . .
+      do k=istart,iend
+         ipts=0
+         read( 15,'(a)',rec=k ) datafile
+c	 print *, datafile
+
+        icompress = 0 
+c--  hack to get trailing 'Z' (if compressed)
+c	do islash = 150, 1, -1
+c	    if ( datafile(islash:islash) .eq. 'Z' ) then
+c		icompress = 1 
+c	        goto 101	
+c	    endif
+c	enddo
+c101     continue
+
+c-- 12/27/99: C call for basename implemented
+c-- variable basename returned via common statement
+        call getbase( datafile ) 
+	print *,'\n file: ',  basename
+
+c--  cyr and cday parsed from basename 
+c--  (e.g., parsed from '1988332h54dd-gdm.hdf')
+         cday=basename( 5:7 )
+         cyr=basename( 1:4 )
+
+c-- if unix compressed, zcat to a temporary filename
+        if( basename(22:22) .eq. 'Z' ) then
+	    icompress = 1
+	    tmpfile = '/tmp/' // basename(1:20)
+	    zcatcmd = 'zcat ' // datafile( 1:len(datafile)-1 ) 
+     &                        // ' > ' // tmpfile
+	    call system( zcatcmd )
+	    datafile = tmpfile
+	endif
+
+c--  convert iday character to integer
+         read(cday,'(i3)') iday
+         read(cyr,'(i4)') iyr
+c         write(6,*)'cday = ',cday
+c         write(6,*)'cyr = ',cyr
+          
+c***> HDF call: d8gimg() 
+         retn=d8gimg(datafile,in_data,x_length,
+     &               y_length,dummy)
+
+c--  read hdf DN values and convert to SST
+         do  i=icolLeft, icolRight
+            do j=irowUpper, irowLower
+                xlat=-89.75+dyd*float(j-1)
+
+c-- center output in Pacific Ocean 
+c                if( i .lt. x_length ) ix = i-(x_length/2)
+c                if( i .le. (x_length/2) ) ix =  i+(x_length/2)
+
+c-- center output in Atlantic (default)
+                ix = i
+
+c--  convert signed byte to float
+c                if ( in_data(i,j) .lt. 0 .and. abs(xlat) .lt. 70. ) then
+                if ( in_data(i,j).lt. 0 ) then
+                    temp(ix,j)=float( in_data(i,j) ) + 256
+                else
+                    temp(ix,j)=float( in_data(i,j) )
+                endif
+
+
+
+C MULTIPLY THE PATHFINDER DIGITAL NUMBER BY THE CALIBRATION NUMBER (0.15)
+C AND ADD THE OFFSET (-3.0) TO GET DEGREES CELSIUS
+
+                 if ( temp(ix,j).gt.0 ) then
+                     ipts=ipts+1
+                     f(ipts)=( cal*temp(ix,j) ) - offset
+                     x(ipts) = ( dxd*float(ix-1) ) - 180. + dxd/2 
+                     y(ipts) = ( 90. - (dyd*float(j-1)) ) - dyd/2
+                     z(ipts)=float(iday)+float(iyr-1985)*365.
+                     if(iyr.gt.1988) z(ipts)=z(ipts)+1
+                     if(iyr.gt.1992) z(ipts)=z(ipts)+1
+                     if(iyr.gt.1996) z(ipts)=z(ipts)+1
+                     if(iyr.gt.2000) z(ipts)=z(ipts)+1
+                     if(iyr.gt.2004) z(ipts)=z(ipts)+1
+                 endif
+            enddo
+         enddo
+
+         nfnew=ipts
+         print *, ' no of pts in file ',' = ', nfnew
+
+c-- calculate interpolation weights and vsum 
+c-- arrays passed directly... a common statement 
+c-- does not seem to work.
+        call  binsum( nfnew, x, y, z, f, 
+     &                xx, yy, zz, vsum,
+     &                ixmin, ixmax, iymin,
+     &                iymax, izmin, izmax, 
+     &                imax, jmax, kmax, ndatamax )
+
+        if( icompress .eq. 1 ) call system( 'rm -f ' // tmpfile )
+c-- ..... read next hdf file from infile
+      enddo
+
+c-- all input files processed; calculate interpolated SST   
+      do kk=1,kmax
+        do jj=1,jmax
+          do ii=1,imax
+            sl(ii,jj,kk)=0.
+            if (vsum(2,ii,jj,kk).gt.0) then
+                sl(ii,jj,kk)=vsum(1,ii,jj,kk)/vsum(2,ii,jj,kk)
+            endif
+          enddo
+        enddo
+      enddo
+
+c-- write output as map "interleaved" or map "sequential"
+c-- "interleaved" is the original implementation
+c-- both formats preceded by header 
+
+c-- "interleaved" refers to each row in the output
+c-- file representing a unique lon, lat position with columns 
+c-- of associated sst values (maps) 
+c-- i.e. row one: lon1 lat1 sst1 sst2 sst3....lastsst
+c--      row two: lon2 lat1 sst1 sst2 sst3....lastsst
+
+c-- "sequential" refers to each map written as rows and
+c-- columns of lat, lon with the array element representing the
+c-- sst at that geo-position for that map.
+c-- each map is written sequentially in the file 
+
+c-- geo-positions of UL and LR corners
+          ULlon = -(180 - ( ((icolLeft-1) * dxd) + dxd/2 ))  
+	  ULlat = (90 - ( ((irowUpper-1) * dyd) + dyd/2 ))
+          LRlon = -(180 - ( (icolRight * dxd) - dxd/2 ))  
+	  LRlat = (90 - ( (irowLower * dyd) - dyd/2 ))
+
+c-- version number, "f" refers to fortran version
+      version = 'f2.9' 
+
+c-- write the 3 record header
+      if( outformat .eq. 'formatted' ) then
+          write(20,'(a)'), plabel
+          write(20,*) imax,jmax,kmax
+          write(20,*) dx,dy,dz 
+	  write(20,*) ULlon,ULlat,LRlon,LRlat
+      elseif( outformat .eq. 'unformatted' ) then
+          write(20) imax,jmax,kmax
+          write(20) dx,dy,dz
+	  write(20) ULlon,ULlat,LRlon,LRlat
+      endif
+
+      if( mapformat .eq. 'interleave' ) then
+        print *, '\n map output is interleave'
+        do jj=1,jmax
+          do ii=1,imax
+            do kk=1,kmax
+                sst_inter(kk)=sl(ii,jj,kk)
+            enddo
+            if( outformat .eq. 'formatted' ) then
+		write(20,*) ii,jj,sst_inter
+            else 
+	        write(20) ( sst_inter(i), i=1,kmax ) 
+            endif
+          enddo
+        enddo
+
+      else if( mapformat .eq. 'sequential' ) then
+        print *, '\n map output is sequential'
+        do kk=1,kmax
+          do jj=1,jmax
+            do ii=1,imax
+               sst_seq(ii)=sl(ii,jj,kk)
+             enddo
+             if( outformat .eq. 'formatted' ) then
+		 write(20,*) jj,kk,sst_seq
+             else
+		 write(20) ( sst_seq(i), i=1,imax )
+             endif
+          enddo
+        enddo
+      endif
+
+      print *,  '\ndone. . . '
+      close( 15,status='keep',err=2000,iostat=ierr )
+      close( 20,status='keep',err=2500,iostat=ierr )
+      stop
+
+ 1000 print *, 'Error opening output: ', outfile, 'error num: ', ierr
+      goto 102
+ 1500 print *, 'Error opening input: ', infile, 'error num: ', ierr
+      goto 102
+ 2000 print *, 'Error closing input: ', infile, 'error num: ', ierr
+      goto 102
+ 2500 print *, 'Error closing output: ', outfile, 'error num: ', ierr
+      goto 102
+ 102  continue
+
+      end

http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/jobClimatology2.py
----------------------------------------------------------------------
diff --git a/climatology/clim/jobClimatology2.py b/climatology/clim/jobClimatology2.py
new file mode 100644
index 0000000..4b3933e
--- /dev/null
+++ b/climatology/clim/jobClimatology2.py
@@ -0,0 +1,22 @@
+#
+# jobClimatology2.py -- mrjob class for computing climatologies by gaussian interpolation
+#
+
+from mrjob.job import MRJob
+from climatology2 import climByAveragingPeriods
+
+
+class MRClimatologyByGaussInterp(MRJob):
+
+    def mapper(self, _, line):
+        yield "chars", len(line)
+        yield "words", len(line.split())
+        yield "lines", 1
+
+    def reducer(self, key, values):
+        yield key, sum(values)
+
+
+if __name__ == '__main__':
+    MRWordFrequencyCount.run()
+

http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/jobTest.py
----------------------------------------------------------------------
diff --git a/climatology/clim/jobTest.py b/climatology/clim/jobTest.py
new file mode 100644
index 0000000..5b1b0f7
--- /dev/null
+++ b/climatology/clim/jobTest.py
@@ -0,0 +1,18 @@
+
+from mrjob.job import MRJob
+
+
+class MRWordFrequencyCount(MRJob):
+
+    def mapper(self, _, line):
+        yield "chars", len(line)
+        yield "words", len(line.split())
+        yield "lines", 1
+
+    def reducer(self, key, values):
+        yield key, sum(values)
+
+
+if __name__ == '__main__':
+    MRWordFrequencyCount.run()
+

http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/C/README
----------------------------------------------------------------------
diff --git a/climatology/clim/orig/C/README b/climatology/clim/orig/C/README
new file mode 100644
index 0000000..7e50ce0
--- /dev/null
+++ b/climatology/clim/orig/C/README
@@ -0,0 +1,6 @@
+2 Mar 00
+
+`test1-100day' was produced by gaussinterp (C version) using /sst/vol7/PATHFINDER/avhrrlists/test.list
+and /sst/vol7/PATHFINDER/parmfiles/interp.test.parms. First 100 files in test.list were processed.
+
+

http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/C/binsum.c
----------------------------------------------------------------------
diff --git a/climatology/clim/orig/C/binsum.c b/climatology/clim/orig/C/binsum.c
new file mode 100644
index 0000000..724245c
--- /dev/null
+++ b/climatology/clim/orig/C/binsum.c
@@ -0,0 +1,125 @@
+/* $Revision: 1.6 $ */
+/*  calculate interpolation weights and vsum */
+
+#include <math.h>
+#include "interp.h"
+
+#define min(x,y)  ( (x) < (y) ? (x) : (y) )
+#define max(x,y)  ( (x) > (y) ? (x) : (y) )
+#define XFAC -0.6931
+
+ binsum( int nf ) {
+
+  register int i, ii, jj, kk, n;
+  int lon_180_index;  
+  float fac;
+  float long_diff_deg, long_diff_rad;
+  float west_long_endpoint, east_long_endpoint;
+
+  /* index in array xx (longitude interp  bins) for long. 180 
+     in extended portion of array beyond imax */
+  if( all_zones )
+      lon_180_index = parms.imax + long_extend / 2;
+  
+  /* determine spatial/temporal interpolation bins for every point */
+  for( n=0; n < nf; n++ ) {
+
+      /* longitude bins ranges */
+      /* Check for the discontinuities along long. 180. If crossing long. 180
+         set ixmin and ixmax to reflect indicies in "extended" portion of 
+	 longitude bins in array xx */
+      
+      west_long_endpoint = x[n] - parms.xwin2;
+      east_long_endpoint = x[n] + parms.xwin2;  
+
+      ixmin[n] = (int) floorf( (x[n]-parms.xwin2-parms.xlo)/parms.dx );
+      ixmax[n] = (int) floorf( (x[n]+parms.xwin2-parms.xlo)/parms.dx ); 
+
+      /* crossing west of 180 in search and global interp */
+      if( west_long_endpoint < -180. && east_long_endpoint > -180. &&
+          all_zones ) { 
+         
+          ixmin[n] = lon_180_index - abs( ixmin[n] );
+          ixmax[n] = lon_180_index + ixmax[n];
+/* printf( " ixmin, ixmax %d %d in extended array\n", ixmin[n], ixmax[n] ); */
+      }
+
+      /* crossing east of 180 in search and global interp */
+      else if( west_long_endpoint < 180. && east_long_endpoint > 180. &&
+	       all_zones ) {
+        
+         ixmin[n] = lon_180_index - ( parms.imax - ixmin[n] );
+         ixmax[n] = lon_180_index + abs( parms.imax - ixmax[n] );
+/* printf( "ixmin, ixmax %d %d in extended array\n", ixmin[n], ixmax[n] ); */
+      }
+      
+      /* no crossing of 180 in search */
+      else {         
+	 ixmin[n] = max( (int) floorf( (x[n]-parms.xwin2-parms.xlo)/parms.dx ), 0 );
+         ixmax[n] = min( (int) floorf( (x[n]+parms.xwin2-parms.xlo)/parms.dx ), parms.imax - 1 );
+     }
+
+      /* lat bin ranges */
+      iymin[n] = max( (int) floorf( (y[n]-parms.ywin2-parms.ylo)/parms.dy ), 0 );
+      iymax[n] = min( (int) floorf( (y[n]+parms.ywin2-parms.ylo)/parms.dy ), parms.jmax - 1 ); 
+   
+      /* temporal bin ranges */
+      izmin[n] = max( (int) floorf( (z[n]-parms.zwin2-parms.zlo)/parms.dz ), 0 );
+      izmax[n] = min( (int) floorf( (z[n]+parms.zwin2-parms.zlo)/parms.dz ), parms.kmax - 1 ); 
+
+  }
+
+ 
+  /* weight the SST and sum */
+  for( n=0; n < nf; n++ ) {
+    for( kk=izmin[n]; kk <= izmax[n]; kk++ ) {
+      for( jj=iymin[n]; jj <= iymax[n]; jj++ ) {
+        for( i=ixmin[n]; i <= ixmax[n]; i++ ) {
+             
+  	/* Derive the guassian weights. Because
+	   longitudes can cross 180, e.g., the ranges
+	   can be 175 to -175, use the acos(cos( long_diff ))) */
+            long_diff_deg = x[n] - xx[i];
+            long_diff_rad = acos( cos(deg2rad(long_diff_deg)) );
+
+	    fac = exp( XFAC * ( pow( (rad2deg(long_diff_rad)
+                                      /parms.xh),2 ) 
+	                    +   pow( ((y[n]-yy[jj])/parms.yh),2 )
+    	                    +   pow( ((z[n]-zz[kk])/parms.zh),2 ) ));
+
+            /* if crossing long. 180 must drop weighted SST 
+               in correct longitude bin between 0 and imax */        
+
+	    /* range crossing east of 180 */  
+            if( i >= lon_180_index && all_zones ) {
+                    ii =  i - lon_180_index; 
+             
+/* printf(" orig bin is %d, correct bin is %d\n", i, ii ); */
+            }
+            /* range crossing west of 180 */
+            else if( (i < lon_180_index && i >= parms.imax)  && all_zones ) {
+                    ii = i - long_extend / 2;
+
+/* printf(" orig bin is %d, correct bin is %d\n", i, ii ); */
+            }
+            /* no crossing */
+            else {
+	        ii = i;
+            }
+                 
+            vsum[0][ii][jj][kk] += f[n]*fac;
+            vsum[1][ii][jj][kk] += fac;
+        }
+      }
+    }
+  }
+ }
+
+
+
+
+
+
+
+
+

http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/C/clouderosion.c
----------------------------------------------------------------------
diff --git a/climatology/clim/orig/C/clouderosion.c b/climatology/clim/orig/C/clouderosion.c
new file mode 100644
index 0000000..6dec692
--- /dev/null
+++ b/climatology/clim/orig/C/clouderosion.c
@@ -0,0 +1,33 @@
+/* $Revision:$ */
+/* Apply the Casey cloud erosion to the satellite data.
+   This algorithm flags as cloudy all data with one 
+   nearest neighbor of a Pathfinder flagged cloud pixel. */
+
+#include "interp.h"
+#include "hdf.h"
+
+clouderosion( float in_data, int i, int j ) {
+    int row, col;
+
+    /* for each image pixel loop over
+       nearest neighbors; if any neightbor is
+       a cloud, set temp[j][i] to cloud value (0)
+    */
+    for( col = i-1; col <= i+1; col++ ) {
+        for( row = j-1; row <= j+1; row++ ) {
+	    printf( "col is %d, row is %d\n", col, row );
+	    /* make sure we are searching within image array bounds */
+	    if( col >= 0 && col < parms.x_length &&  
+		row >= 0 && row < parms.y_length ) {
+	    printf( "valid col is %d, valid row is %d\n\n", col, row );
+		if( in_data[col][row] == 0 ) {
+		    temp[j][i] = 0.;
+		    printf( "in_data[%d][%d] is %c\n", col, row, in_data[col][row] );
+		    printf( " setting temp[%d][%d] to zero\n\n\n", j,i );
+		    break;
+                }
+            }
+         }
+    }
+
+}

http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/C/gaussinterp.readme
----------------------------------------------------------------------
diff --git a/climatology/clim/orig/C/gaussinterp.readme b/climatology/clim/orig/C/gaussinterp.readme
new file mode 100644
index 0000000..4a0e017
--- /dev/null
+++ b/climatology/clim/orig/C/gaussinterp.readme
@@ -0,0 +1,159 @@
+05 Oct 00
+ed armstrong	NASA/JPL/Caltech
+
+`gaussinterp' is a program designed to spatially and temporally
+interpolate SST pathfinder data using gaussian methodology. This
+description is for the C version. To run the program from the command
+line:
+
+% gaussinterp <infile> <beg_rec> <parmfile> <outfile> 
+
+where
+	infile = file containing list of files to process 
+	(ascii; fixed length records; temporally sorted from 
+	old to recent)
+	beg_rec = first record (file) in infile processed
+	parmfile = input parameter file
+	outfile = interpolated output file (ascii or binary)
+
+
+`infile' must be ascii with a fixed record length of 150 bytes. Use
+`pad.pl' to convert variable or other fixed record length files to
+150.  `infile' records should be in monotonically ascending in time
+(year and julian day). beg_rec is the record number of the temporal
+start of the interpolation. The day for the file corresponding to this
+record number is also the first map center date. The next map center is
+'first_map_center' + time_step (parameter dz; see parameter file
+description) and so forth.  `outfile' can be either an ascii or binary
+file, sequential or interleaved format; see below for more
+information.  `parmfile' is an ascii file read by the program that
+contains many interpolation and output parameters.
+
+
+** Parameter file description **
+
+~~~~ example of parameter file: ~~~~~
+200000
+720 360
+10
+3.0 3.0 20.0
+1.0 1.0 10.0
+1.0 1.0  8.0
+-180. 180. 
+-90. 90.
+sequential
+unformatted
+1
+1
+6
+~~~~ end example ~~~~
+
+record 1: ndatamax = maximum points in each file
+       2: xlength, y_length = HDF data array: cols and rows
+       3: kmax = number of "maps" or time "steps"  
+       4: xwin2, ywin2, zwin2 = longitude, latitude "search" 
+	  range (decimal degrees), time search range (days)
+       5: xh, yh, zh = e folding scaling for long, lat (decimal
+		       degrees) and time (days) 
+       6: dx, dy, dz = longitude, latitude interpolation 
+	  	       step (decimal degrees), and 
+		       time step (days)
+       7: xlo, xhi = minimum, maximum longitude (-180 to 180 decimal degrees)
+       8: ylo, yhi = minimum, maximum latitude (-90 to 90 decimal degrees)
+       9: 'interleave' or 'sequential' map format 
+      10: 'formatted' (ascii) or 'unformatted' (binary) output file
+      11: cloud erosion flag: 1=cloud erosion; 0=no cloud erosion
+      12: all pixel data flag: 1="all pixel data"; 0="best pixel" data
+      13: quality flag value: quality flag value (0-7) to use for "all pixel" data; 
+			    ignored if "best pixel" data used
+
+
+Some notes on the parameter file: 
+
+	* ndatamax should be x_length*y_length but in reality is much
+	less since much of the data are missing due to clouds.
+
+	* the user can subset into the image array by choosing proper
+	xlo, xhi, ylo, yhi. xlo must always be less than xhi; the same
+	for ylo and yhi.
+
+	* interleave format refers to each row in the output file
+	representing a unique lon, lat position with columns of
+	associated sst values (maps) 
+	i.e. row one: lon1 lat1 sst1 sst2 sst3....lastsst 
+	     row two: lon2 lat1 sst1 sst2 sst3....lastsst
+
+	* sequential format refers to each map written as rows and
+	columns of lat, lon with the array element representing the sst
+	at that geo-position for that map. Each map is written
+	sequentially in the output file.
+
+	* generally unformatted (binary) output should be chosen as
+	this results in a smaller output file size. Using binary format
+	each pixel is written as a 1 byte char (0-255).
+
+	* cloud erosion refers to the Casey cloud erosion algorithm
+	used in the NSIPP (Casey) Pathfinder climatology. 
+
+
+** Header description **
+
+For each output format a header is written:
+     	ASCII format:
+	record 1: some text
+	record 2: imax, jmax, kmax
+		imax = width of output in pixels
+		jmax = height of output in pixels 
+		kmax = number of "maps" or images
+        record 3: xwin2, ywin2, zwin2
+		spatial and temporal search "radii" 
+		in decimal degrees and days
+        record 4: xh, yh, zh
+		e folding scales in decimal degrees
+	record 5: dx,dy,dz
+		interpolation intervals
+		in decimal degrees
+	record 6: ULx, ULy, LRx, LRy
+		upper left and lower right longitude
+		and latitudes in decimal degrees
+        record 7: date.startdate, date.enddate
+		first and last map center dates in
+		YYYYDDD (e.g., 1995010)
+
+        BINARY (C) [word is 4 bytes, types are indicated]: 
+		word 1 [char]: magic number
+		words 2-4 [int]: imax, jmax, kmax
+		words 5-7 [float]: xwin2, ywin2, zwin2
+		words 8-10 [float]: xh, yh, zh
+		words 11-13 [float]: dx,dy,dz
+		words 14-17 [float]: ULx, ULy, LRx, LRy
+		word 18 [int]: first map center date
+		word 19 [int]: last map center date
+		[total header size is 76 bytes]
+
+** Additional information **
+
+- The parameter file can be edited to change/update all entries but
+care must be taken to avoid typos. The program is compiled for a fixed
+record length for the input datafile. These can be modified in
+`interp.h' and the program recompiled.
+
+- In the case of global interpolation where xlo=-180 and xhi=180, the
+program will correctly interpolate across longitude 180 East (all_zones
+is true). That is, there is no "cutoff" in the interpolation if the
+search window for a particular satellite pixel crosses this longitude.
+In all other cases this is not so.  For example, if xlo=-180 and
+xhi=-150, interpolations in the vicinity of longitude -180 are "cutoff"
+at this boundary, i.e., points west of -180 do not contribute to the
+interpolation.  Since this is a regional subset, the interpolations are
+also "cutoff" in the vicinity of longitude -150 (data east of -150 do
+not contribute). By extension this is also true for a subset in
+latitude.
+
+- The first line in the usage dump (type gaussinterp with no arguments)
+will specify if the executable is from Fortran or C code.
+
+- Sample parameter files can be found in
+/sst/vol7/PATHFINDER/interp/parmfiles .  Datafile lists can be found in
+/sst/vol7/PATHFINDER/interp/avhrrlists .
+

http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/C/gaussinterp_C_code.tar
----------------------------------------------------------------------
diff --git a/climatology/clim/orig/C/gaussinterp_C_code.tar b/climatology/clim/orig/C/gaussinterp_C_code.tar
new file mode 100644
index 0000000..5db4b0e
Binary files /dev/null and b/climatology/clim/orig/C/gaussinterp_C_code.tar differ

http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/C/interp.c
----------------------------------------------------------------------
diff --git a/climatology/clim/orig/C/interp.c b/climatology/clim/orig/C/interp.c
new file mode 100644
index 0000000..c057875
--- /dev/null
+++ b/climatology/clim/orig/C/interp.c
@@ -0,0 +1,448 @@
+/* $Revision: 1.20 $ */
+
+/* SYNOPSIS
+   called by setupinterp(). This subroutine reads each
+   HDF file, passes data to binsum() for interpolation
+   calculation and writes output image(s).
+*/
+
+#include <stdio.h>
+#include <string.h>
+#include <libgen.h>
+#include <hdf.h>
+#include "interp.h"
+
+/* C binary file id = 3890345459 (0xE7E1F5F3) */
+#define MAGIC 3890345459L
+
+ interp() {
+  uint8 in_data[parms.y_length][parms.x_length];   
+  uint8 quality_data[parms.y_length][parms.x_length];   
+  uint8 aerosol_corr[parms.y_length][parms.x_length];   
+
+  intn retn, haspal;
+  int32 width, height;
+
+  float ULlon, ULlat, LRlon, LRlat;
+  float cal, offset, xlat;
+  int k1,k2,k3,k4,ii,jj,kk,n,i,j,k;
+  int ipts, iyr, iday;
+  int jj_reverse;
+  int row, col;
+  unsigned long int magic_num; 
+
+  FILE *fid_infile, *fid_outfile, *fid_datafile;
+  FILE *fid_aerosol_list, *fid_aerosol_file;
+
+  char plabel[40];
+  char datafile[RECSIZE];
+  char aerosol_file[200];
+  char *basefile;
+  char cday[3];
+  char cyr[4];
+  char tmpstr[22]; 
+  char aerosol_tmpstr[30];
+  char tmpfile[35];
+  char aerosol_tmpfile[45];
+  char zcatcmd[200];
+  char rmcmd[45];
+  char compress[1];
+  char version[4];
+  char sstDN;
+  /* short int sstDN; */
+
+  /* open output file */
+  fid_outfile = fopen( fileio.outfile, "w" );
+
+
+  /* Set arrays of interpolation bins.
+     Long and lat bin centers set to 1/2 
+     interpolation step */
+
+  /*  Long to interpolate to:
+      1st for loop sets interpolation bins
+      between the ranges of xlo to xhi.
+      If all_zones is true, 2nd and 3rd for 
+      loops set bins for those "extended" interpolations 
+      that must cross longitude 180. */
+  for( ii=parms.imin; ii < parms.imax; ii++ )
+      xx[ii] = parms.xlo + (float)(ii)*parms.dx + parms.dx/2;
+
+  if( all_zones ) {
+      /* west of long. 180 */
+      for( ii=0; ii < long_extend/2; ii++ )
+          xx[ii+parms.imax] = 180 - (long_extend/2 * parms.dx) + (float)(ii)*parms.dx 
+                              + parms.dx/2;
+
+      /* east of long. 180; sets center point */
+      for( ii=0; ii <= long_extend/2; ii++ )
+         xx[ii+parms.imax+long_extend/2] = -180 + (float)(ii)*parms.dx 
+                                           + parms.dx/2;
+  }
+
+  /* for( ii = 0; ii <=  parms.imax + long_extend; ii++) 
+    printf( " ii: %d xx[ii] = %f \n", ii, xx[ii] ); */
+ 
+  
+  /* Lat to interpolate to */
+  for( jj=parms.jmin; jj < parms.jmax; jj++ )
+      yy[jj] = parms.ylo + (float)(jj)*parms.dy + parms.dy/2;
+
+  /* Time to interpolate to */ 
+  for( kk=parms.kmin; kk < parms.kmax; kk++ )
+      zz[kk] = parms.zlo + (float)(kk)*parms.dz;	
+
+  /* slope (cal) and offset to convert DN to SST */
+  cal = 0.15;
+  offset = 3.0;
+
+  /* open input data file list */
+  if( (fid_infile = fopen( fileio.infile, "r" )) == NULL ){
+      printf( "\nCan not open %s\n", fileio.infile );
+      exit(1);
+  }
+
+  /* position to proper record in file */
+  fseek( fid_infile, ( (fileio.istart-1)*RECSIZE ), 0 );
+
+
+
+  /* open input aerosol file list if necessary */
+  if(aerosol) {
+      if( (fid_aerosol_list = fopen( fileio.aerosol_list, "r" )) == NULL ){
+          printf( "\nCan not open %s\n", fileio.aerosol_list );
+          exit(1);
+      }
+      /* position to proper record in file */
+      fseek( fid_aerosol_list, ( (fileio.istart-1)*RECSIZE ), 0 ); 
+
+  }
+
+
+  /* read infile one record at a time and process . . . */
+  for( k=fileio.istart; k <= fileio.iend; k++ ){
+
+      ipts = 0; 
+      strcpy( compress, " " );
+      strcpy( tmpstr, "                      " );
+
+      fgets( datafile, sizeof(datafile)+1, fid_infile ); 
+
+      /* strip off newline if necessary */
+      if( datafile[ strlen(datafile) ] == '\n' ){ 
+          datafile[ strlen(datafile) ] = '\0';
+      }
+      /* printf("datafile is %s \n", datafile ); */
+
+      /* check to make sure datafile actually exists 
+         if not ..... quit program  */
+      fid_datafile = fopen( datafile, "r" ); 
+      if( fid_datafile == NULL ) {
+	  printf( "could not open:\n" );
+	  printf( "%s \n", datafile );
+	  printf( ". . . quitting program \n" );
+	  exit(2);
+      }
+      fclose( fid_datafile ); 
+
+      basefile = basename( &datafile[0] );
+      printf( "\nfile: %s\n", basefile );
+
+ /* e.g, parse year and jul day from 1985004h54dd-gdm.hdf */
+      strcpy( tmpstr, basefile ); 
+      strncpy( cyr, &tmpstr[0], 4 ); 
+      iyr = atoi( cyr );
+      strncpy( cday, &tmpstr[4], 3 );
+      iday = atoi( cday );
+   /* printf(" iday is %d\n", iday ); 
+      printf("iyr is %d\n", iyr ); 
+   */
+
+ /* if unix compressed, zcat to a temporary filename */
+     strncpy( compress, &tmpstr[21], 1 );
+     if( strcmp( compress, "Z" ) == 0 ) {
+         strcpy( tmpstr, basefile );
+         strcpy( tmpfile, "/tmp/" ); 
+	 strncat( tmpfile, &tmpstr[0], 20 ); 
+
+	 strcpy( zcatcmd, "zcat ");
+	 strcat( zcatcmd, datafile );
+	 strcat( zcatcmd, " > " );
+	 strcat( zcatcmd, tmpfile );
+
+         (void)  system( zcatcmd );
+	 strcpy( datafile,  tmpfile );
+     }
+
+     /* do the same for the aerosol file if necessary */
+     if( aerosol ) {
+      strcpy( compress, " " );
+      strcpy( aerosol_tmpstr, "                             " );
+
+      fgets( aerosol_file, sizeof(aerosol_file)+1, fid_aerosol_list ); 
+      /* printf("aerosol_file is %s \n", aerosol_file ); */
+
+      /* strip off newline if necessary  */
+      if( aerosol_file[ strlen(aerosol_file) ] == '\n' ){ 
+          aerosol_file[ strlen(aerosol_file) ] = '\0';
+      }
+
+
+      basefile = basename( &aerosol_file[0] );
+      strcpy( aerosol_tmpstr, basefile ); 
+      printf( "aerosol file: %s\n", basefile ); 
+
+
+     /* if unix compressed, zcat to a temporary filename  */
+     strncpy( compress, &aerosol_tmpstr[27], 1 );
+     if( strcmp( compress, "Z" ) == 0 ) {
+         strcpy( aerosol_tmpstr, basefile );
+         strcpy( aerosol_tmpfile, "/tmp/" ); 
+	 strncat( aerosol_tmpfile, &aerosol_tmpstr[0], 26 ); 
+
+	 strcpy( zcatcmd, "zcat ");
+	 strcat( zcatcmd, aerosol_file );
+	 strcat( zcatcmd, " > " );
+	 strcat( zcatcmd, aerosol_tmpfile );
+
+         (void)  system( zcatcmd );
+	 strcpy( aerosol_file,  aerosol_tmpfile );
+     }
+     }
+
+         /* read the aerosol file into an array */
+	 if( aerosol ) { 
+	     fid_aerosol_file = fopen( aerosol_file, "r" );
+             fread( aerosol_corr, sizeof(uint8), parms.y_length * parms.x_length, fid_aerosol_file );
+	     fclose( fid_aerosol_file );
+	 }
+
+         /* read HDF SST data dimensions */
+	 retn = DFR8getdims( datafile, &width, &height, &haspal );
+	 /* printf( "HDF width and height: %d %d \n", width, height ); */
+
+	 if( width != parms.x_length || height != parms.y_length ) {
+	     printf( "\nbounds mismatch after reading HDF file dimensions !!\n" );
+	     printf( "width: %d, parms.x_length: %d \n", width, parms.x_length ); 
+	     printf( "height: %d, parms.y_length: %d \n", height, parms.y_length ); 
+	     printf( "skipping %s . . . \n\n", datafile );
+	     continue;
+         }
+
+ /* ------------------------------------------------- */
+ /* Read the SST image of the HDF file.
+    HDF call: DFR8getimage()  */
+ /* ------------------------------------------------- */
+
+	 retn=DFR8getimage(datafile, (VOIDP)in_data, parms.x_length, 
+			   parms.y_length, NULL);
+
+ /* ------------------------------------------------- */
+ /* If "all pixel" data choosen read the quality  
+    flag image */
+ /* ------------------------------------------------- */
+	
+	if( sst_type )  
+	   retn=DFR8getimage(datafile, (VOIDP)quality_data, 
+			   parms.x_length, parms.y_length, NULL);
+
+
+
+ /* ------------------------------------------------- */
+ /* Read hdf file DN values and convert to SST.
+    If "all pixel" data read only those pixels 
+    with a quality flag greater than or equal to 
+    quality_flag */
+ /* ------------------------------------------------- */
+
+        for( i=hdf.icolLeft; i <= hdf.icolRight; i++ ) {
+	  for( j=hdf.irowUpper; j <= hdf.irowLower; j++ ) {
+
+              /* only conider aerosol correction greater than 5 cnt */
+	      if( aerosol && aerosol_corr[j][i] > 5  && in_data[j][i] > 0 )
+		  in_data[j][i] = in_data[j][i] + aerosol_corr[j][i];
+
+	      /* perform Casey cloud erosion on non-zero SST data */
+	      /* if nearest neighbor is a cloud, set pixel to cloud (0) */
+	      if( in_data[j][i] > 0 && cloud_erosion ) {
+                  for( col = i-1; col <= i+1; col++ ) {
+                     for( row = j-1; row <= j+1; row++ ) { 
+	             /* make sure we are searching within 
+		     image array bounds */
+	              if( col >= 0 && col < parms.x_length &&  
+		          row >= 0 && row < parms.y_length && 
+		          in_data[row][col] == 0 ) {
+		            in_data[j][i] = 0;
+			    goto end_erosion;
+                      }
+                    }
+                  }
+	      }
+	      end_erosion:
+
+
+              /* if all_pixel (sst_type = 1) filter for quality level, 
+		 or if best_sst accept data values greater than 0  */
+		  
+		  /*printf( "quality flag is %d \n", quality_flag);
+		  printf (" sst_type is %d \n", sst_type);
+		  */
+
+	      if( (sst_type && quality_data[j][i] >= quality_flag) || 
+		  (!sst_type && in_data[j][i] > 0) ) {
+		  /* printf (" quality_data = %d  \n", quality_data[j][i]); */
+
+                  f[ipts] = cal * (float)in_data[j][i]  - offset;
+                  x[ipts] = -180. + parms.dxd * (float)(i) + parms.dxd/2;
+                  y[ipts] = 90. - parms.dyd * (float)(j) - parms.dyd/2;
+                  z[ipts] = (float)( iday+( (iyr-1985)*365 ) );
+
+                  if(iyr > 1988) z[ipts]=z[ipts]+1;
+                  if(iyr > 1992) z[ipts]=z[ipts]+1;
+                  if(iyr > 1996) z[ipts]=z[ipts]+1;
+                  if(iyr > 2000) z[ipts]=z[ipts]+1;
+                  if(iyr > 2004) z[ipts]=z[ipts]+1;
+                  if(iyr > 2008) z[ipts]=z[ipts]+1;
+		  ipts++;
+	      }
+	  }
+        } 
+
+ /* calculate interpolation weights and vsum  */
+
+	  printf(" num of pts in file: %d \n", ipts );
+	  binsum( ipts );
+
+        if( strcmp( compress, "Z" ) == 0 ) {
+	    strcpy( rmcmd, "rm -f " );
+	    strcat( rmcmd, tmpfile );
+	    system( rmcmd ); 
+	    if( aerosol ) {
+	        strcpy( rmcmd, "rm -f " );
+	        strcat( rmcmd, aerosol_tmpfile );
+	        system( rmcmd ); 
+            }
+        }
+ /* ..... read next hdf file from infile */
+
+  }
+  fclose( fid_infile );
+  if( aerosol )
+      fclose( fid_aerosol_list );
+
+ /* all input files processed; calculate interpolated SST maps */
+      for( kk=0; kk < parms.kmax; kk++ ) {
+        for( jj=0; jj < parms.jmax; jj++ ) {
+          for( ii=0; ii < parms.imax; ii++ ) {
+            if ( vsum[1][ii][jj][kk] > 0 ) {  
+
+		/* reverse jj indicies to flip N to S */
+		jj_reverse = parms.jmax - jj - 1;
+                sl[ii][jj_reverse][kk] = 
+		vsum[0][ii][jj][kk] / vsum[1][ii][jj][kk];
+            }
+          }
+        }
+      }
+
+ /* write output as map "interleaved" or map "sequential".
+    both formats preceded by header.
+    see gaussinterp.readme for more info.
+ */
+
+ /* geo-positions of UL and LR corners */
+      ULlon = parms.xlo + parms.dx/2;
+      ULlat = parms.yhi - parms.dy/2;
+      LRlon = parms.xhi - parms.dx/2;
+      LRlat = parms.ylo + parms.dy/2;
+
+ /* version number, "c" refers to C version */
+      strcpy( version, "c1.3" );
+      strcpy( plabel, "sst interpolation");
+
+ /*  write the header */
+      if( strcmp( fileio.outformat, "formatted" ) == 0 ) { 
+	  fprintf( fid_outfile, "%s -- %s\n ", version, plabel );
+	  fprintf( fid_outfile, "%d %d %d\n", 
+		   parms.imax,parms.jmax,parms.kmax );
+	  fprintf( fid_outfile, "%f %f %f\n", 
+		   parms.xwin2,parms.ywin2,parms.zwin2 );
+	  fprintf( fid_outfile, "%f %f %f\n", 
+		   parms.xh,parms.yh,parms.zh );
+	  fprintf( fid_outfile, "%f %f %f\n", 
+		   parms.dx,parms.dy,parms.dz );
+	  fprintf( fid_outfile, "%f %f %f %f\n", 
+	  	   ULlon,ULlat,LRlon,LRlat );
+          fprintf( fid_outfile, "%ld %ld\n",
+                   date.startdate, date.enddate );
+      }
+      else if( strcmp( fileio.outformat,  "unformatted" ) == 0 ) { 
+	  magic_num = MAGIC;
+	  fwrite( &magic_num, sizeof(magic_num),1,fid_outfile ); 
+	  fwrite( &parms.imax, sizeof(parms.imax),1,fid_outfile );
+	  fwrite( &parms.jmax, sizeof(parms.jmax),1,fid_outfile );
+	  fwrite( &parms.kmax, sizeof(parms.kmax),1,fid_outfile );
+	  fwrite( &parms.xwin2, sizeof(parms.xwin2),1,fid_outfile );
+	  fwrite( &parms.ywin2, sizeof(parms.ywin2),1,fid_outfile );
+	  fwrite( &parms.zwin2, sizeof(parms.zwin2),1,fid_outfile );
+	  fwrite( &parms.xh, sizeof(parms.xh),1,fid_outfile );
+	  fwrite( &parms.yh, sizeof(parms.yh),1,fid_outfile );
+	  fwrite( &parms.zh, sizeof(parms.zh),1,fid_outfile );
+	  fwrite( &parms.dx, sizeof(parms.dx),1,fid_outfile );
+	  fwrite( &parms.dy, sizeof(parms.dy),1,fid_outfile );
+	  fwrite( &parms.dz, sizeof(parms.dz),1,fid_outfile );
+	  fwrite( &ULlon, sizeof(ULlon),1,fid_outfile );
+	  fwrite( &ULlat, sizeof(ULlat),1,fid_outfile );
+	  fwrite( &LRlon, sizeof(LRlon),1,fid_outfile );
+	  fwrite( &LRlat, sizeof(LRlat),1,fid_outfile );
+          fwrite( &date.startdate, sizeof(date.startdate),1,fid_outfile);
+	  fwrite( &date.enddate, sizeof(date.enddate),1,fid_outfile);
+	  /* binary header is 76 bytes */
+      }
+
+ /*  write the maps */
+      if( strcmp( fileio.mapformat, "interleave" ) == 0 ) {
+        printf( "\n map output is interleave\n" );
+        for( jj=0; jj < parms.jmax; jj++ ) {
+          for( ii=0; ii < parms.imax; ii++ ) {
+            for( kk=0; kk < parms.kmax; kk++ ) {
+                sst_inter[kk] = sl[ii][jj][kk]; 
+
+		/* scale sst back to 0-255 DNs */
+		if( sst_inter[kk] != 0 ) 
+		    sstDN = (char) nint( (sst_inter[kk] + offset) / cal );
+                else
+		    sstDN = 0;	/* O's in array sst_inter are missing data */ 
+
+	        if( strcmp( fileio.outformat, "formatted" ) == 0 ) 
+		    fprintf( fid_outfile, "%d ", sstDN ); 
+                else 
+		    fwrite( &sstDN, sizeof(char),1,fid_outfile ); 
+
+             }
+          }
+        }
+      } else if( strcmp( fileio.mapformat, "sequential" ) == 0 ) {
+        printf( "\n map output is sequential\n" );
+        for( kk=0; kk < parms.kmax; kk++ ) {
+          for( jj=0; jj < parms.jmax; jj++ ) {
+            for( ii=0; ii < parms.imax; ii++ ) {
+               sst_seq[ii] = sl[ii][jj][kk]; 
+
+		/* scale sst back to 0-255 DNs */
+	       if( sst_seq[ii] != 0 ) 
+	           sstDN = (char) nint( (sst_seq[ii] + offset) / cal );
+               else
+		   sstDN = 0;	/* O's in array sst_seq are missing data */ 
+
+               if( strcmp( fileio.outformat, "formatted" ) == 0 ) 
+	           fprintf( fid_outfile, "%d ", sstDN ); 
+               else 
+	     	   fwrite( &sstDN, sizeof(char),1,fid_outfile );
+
+	    }
+          }
+        }
+      }
+  fclose( fid_outfile );
+}

http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/C/makefile
----------------------------------------------------------------------
diff --git a/climatology/clim/orig/C/makefile b/climatology/clim/orig/C/makefile
new file mode 100644
index 0000000..ab004b3
--- /dev/null
+++ b/climatology/clim/orig/C/makefile
@@ -0,0 +1,33 @@
+# makefile for guassinterp (C version). 
+
+# C compiler 
+CC  =  cc
+#HDF = /usr/local/hdf
+HDF = /usr/local/hdf4.1r3-n32
+
+CFLAGS =  -n32 -I$(HDF)/include  -O -w
+#LIBS    = -L$(HDF)/lib -lmfhdf -ldf -ljpeg -lz -lgen -lm 
+LIBS    = -L$(HDF)/lib  -ldf -ljpeg -lz -lgen -lm 
+BIN = /sst/vol7/PATHFINDER/bin/sgi
+PGM = gaussinterp
+
+SRC = setupinterp.c interp.c binsum.c 
+OBJ = setupinterp.o interp.o binsum.o 
+
+$(PGM):	$(OBJ)
+	$(CC) $(CFLAGS) -o $(PGM) $(OBJ) $(LIBS) 
+
+setupinterp.o: 
+	$(CC) $(CFLAGS) $(DEFINE) -c setupinterp.c
+interp.o:
+	$(CC) $(CFLAGS) $(DEFINE) -c interp.c
+binsum.o:
+	$(CC) $(CFLAGS) $(DEFINE) -c binsum.c
+
+
+install: $(PGM)
+	 cp  $(PGM) $(BIN)
+
+clean:
+	 rm -f $(OBJ) $(PGM)
+

http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/C/setupinterp.c
----------------------------------------------------------------------
diff --git a/climatology/clim/orig/C/setupinterp.c b/climatology/clim/orig/C/setupinterp.c
new file mode 100644
index 0000000..9ca67b3
--- /dev/null
+++ b/climatology/clim/orig/C/setupinterp.c
@@ -0,0 +1,431 @@
+/* NAME gaussinterp (C version)
+   16 Jan 00    earmstrong      NASA/JPL
+   $Revision: 1.16 $
+
+  DESCRIPTION
+  Gaussian interpolation program modified to map SST (HDF) data.
+  This program is a C port of the fortran version of gaussinterp.
+
+  SYNOPSIS
+  setupinterp() allocates space for arrays, calls interp()
+  which does the interpolation
+
+  USAGE
+  % gaussinterp <infile> <beg_rec> <parmfile> <outfile> */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <libgen.h>
+#include <math.h>
+#include "interp.h"
+/* see interp.h for global structs and arrays */
+      
+main( int argc, char *argv[] ) {
+
+      int iday, iyr, calloc_bytes, rec_step;
+      long max_infile_size;
+
+      int value;
+      FILE *fid_infile, *fid_parm;
+
+      char datafile[RECSIZE];
+      char *basefile;
+      char *parmfile;
+      char *startrec, *endrec;
+      char cday[3], cyr[4], cdate[7];
+      char *cstartdate;
+      char *cenddate;
+      char tmpline[80];
+
+/* check command line arguments */
+      if( argc != 5 && argc != 6 ) {
+	  usage();
+	  exit(1);
+      }
+      fileio.infile = argv[1];
+      startrec = argv[2];
+      parmfile = argv[3];
+      fileio.outfile = argv[4];
+      if( argc == 6 )
+          fileio.aerosol_list = argv[5];  
+
+      /* open parameter file and read line by line */
+      if( (fid_parm = fopen( parmfile, "r" )) == NULL ){
+	  printf( "\nCan not open %s\n", parmfile ); 
+	  exit(1);
+      }
+
+      fgets( tmpline, sizeof(tmpline), fid_parm );
+      sscanf( tmpline, "%d", &parms.ndatamax ); 
+      fgets( tmpline, sizeof(tmpline), fid_parm );
+      sscanf( tmpline, "%d %d", &parms.x_length, &parms.y_length ); 
+      fgets( tmpline, sizeof(tmpline), fid_parm );
+      sscanf( tmpline, "%d", &parms.kmax ); 
+      fgets( tmpline, sizeof(tmpline), fid_parm );
+      sscanf( tmpline, "%f %f %f", &parms.xwin2, &parms.ywin2, &parms.zwin2 ); 
+      fgets( tmpline, sizeof(tmpline), fid_parm );
+      sscanf( tmpline, "%f %f %f", &parms.xh, &parms.yh, &parms.zh ); 
+      fgets( tmpline, sizeof(tmpline), fid_parm );
+      sscanf( tmpline, "%f %f %f", &parms.dx, &parms.dy, &parms.dz ); 
+      fgets( tmpline, sizeof(tmpline), fid_parm );
+      sscanf( tmpline, "%f %f", &parms.xlo, &parms.xhi ); 
+      fgets( tmpline, sizeof(tmpline), fid_parm );
+      sscanf( tmpline, "%f %f", &parms.ylo, &parms.yhi ); 
+      fgets( fileio.mapformat, sizeof(fileio.mapformat), fid_parm );
+      fgets( fileio.outformat, sizeof(fileio.outformat), fid_parm );
+      fgets( tmpline, sizeof(tmpline), fid_parm );
+      sscanf( tmpline, "%d", &cloud_erosion );
+      fgets( tmpline, sizeof(tmpline), fid_parm );
+      sscanf( tmpline, "%d", &sst_type );
+      fgets( tmpline, sizeof(tmpline), fid_parm );
+      sscanf( tmpline, "%d", &quality_flag );
+      fgets( tmpline, sizeof(tmpline), fid_parm );
+      sscanf( tmpline, "%d", &aerosol );
+
+      fileio.mapformat[ strlen(fileio.mapformat)-1 ] = '\0';
+      fileio.outformat[ strlen(fileio.outformat)-1 ] = '\0';
+
+      fclose( fid_parm );
+      
+      if( strcmp(fileio.mapformat, "interleave") != 0  &&   
+          strcmp(fileio.mapformat, "sequential") != 0 ) {
+	   printf( "check parameter file for interleave or sequential maps !\n" );
+	   exit(1);
+      } 
+      if( strcmp(fileio.outformat, "formatted") != 0  &&   
+          strcmp(fileio.outformat, "unformatted") != 0 ) {
+	   printf( "check parameter file for formatted or unformatted output !\n" );
+	   exit(1);
+      }
+
+    /* Calculate number of interpolation "steps" based on
+       long and lat ranges divided by interpolation interval.
+       Origin is at lon=-180, lat=90 */
+      parms.imax = floorf((parms.xhi + 180) / parms.dx) - floorf((parms.xlo + 180) / parms.dx );
+      parms.jmax = floorf((parms.yhi - 90) / parms.dy) - floorf((parms.ylo - 90) / parms.dy );
+      printf( "imax, jmax, kmax: %d %d %d \n", 
+		parms.imax, parms.jmax, parms.kmax );
+
+      parms.imin = 0;
+      parms.jmin = 0;
+      parms.kmin = 0;
+
+    /* calculate degrees per grid point */
+      parms.dxd =  360. / (float)parms.x_length;
+      parms.dyd =  180. / (float)parms.y_length;
+
+    /* find columns and rows for subsetting into hdf image array 
+       remember @ hdfarray(0,0) -->  lon=-180, lat=90 */
+      hdf.icolLeft  = floorf( (parms.xlo + 180.) / parms.dxd );
+      hdf.icolRight = floorf( ((parms.xhi + 180.) / parms.dxd) - 1. );
+      hdf.irowUpper = floorf( fabsf(parms.yhi - 90.) / parms.dyd );
+      hdf.irowLower = floorf( (fabsf(parms.ylo - 90.) / parms.dyd) - 1. );
+
+      /* printf( "array corners: %d %d %d %d \n", hdf.icolLeft, hdf.irowUpper, 
+		  hdf.icolRight ,hdf.irowLower); */
+
+   /* set global interpolation flag if necessary */
+      if( parms.xlo == -180. && parms.xhi == 180. ) {
+          all_zones = 1;  /* true */
+          printf( "all_zones is true: interpolating across -180 West\n" );
+      }
+
+      if( cloud_erosion )
+	  printf( "performing cloud erosion filtering on SST data \n" );
+      if( sst_type )
+	  printf( "working with 'all pixel' data: quality flag is %d \n", quality_flag );
+      if( aerosol )
+	  printf( "performing aerosol correction \n" );
+
+   /* allocate space for arrays  */
+      calloc_bytes = 0;
+      if( (temp = fmal2d(  parms.x_length, parms.y_length )) == NULL ) 
+	  printf( "not enough memory for array temp\n" ), exit(-1);
+      else calloc_bytes = parms.x_length * parms.y_length * sizeof(float);  
+
+      if( (in_data = cmal2d( parms.x_length, parms.y_length )) == NULL ) 
+	  printf( "not enough memory for array in_data\n" ), exit(-1);
+      else calloc_bytes += parms.x_length * parms.y_length * sizeof(char);  
+
+      if( (quality_data = cmal2d( parms.x_length, parms.y_length )) == NULL ) 
+	  printf( "not enough memory for array quality_data\n" ), exit(-1);
+      else calloc_bytes += parms.x_length * parms.y_length * sizeof(char);  
+
+       if( (sl = fmal3d( parms.imax, parms.jmax, parms.kmax )) == NULL )
+	  printf( "not enough memory for array sl\n" ), exit(-1);
+      else calloc_bytes += parms.imax * parms.jmax * parms.kmax * sizeof(float);  
+
+      if( (vsum = fmal4d( 2, parms.imax, parms.jmax, parms.kmax )) == NULL )
+	  printf( "not enough memory for array vsum\n" ), exit(-1);
+      else calloc_bytes += parms.imax * parms.jmax * parms.kmax * 2 * sizeof(float);  
+    
+      if( (sst_seq = (float *)calloc( parms.imax, sizeof(float) )) == NULL ) 
+	  printf( "not enough memory for array sst_seq\n" ), exit(-1);
+      else calloc_bytes += parms.imax * sizeof(float);  
+
+      if( (sst_inter = (float *)calloc( parms.kmax, sizeof(float) )) == NULL )
+	  printf( "not enough memory for array sst_inter\n" ), exit(-1);
+      else calloc_bytes += parms.kmax * sizeof(float);  
+
+      if( (x = (float *)calloc( parms.ndatamax, sizeof(float) )) == NULL ) 
+	  printf( "not enough memory for array x\n" ), exit(-1);
+      else calloc_bytes += parms.ndatamax * sizeof(float);  
+
+      if( (y = (float *)calloc( parms.ndatamax, sizeof(float) )) == NULL )
+	  printf( "not enough memory for array y\n" ), exit(-1);
+      else calloc_bytes += parms.ndatamax * sizeof(float);  
+
+      if( (z = (float *)calloc( parms.ndatamax, sizeof(float) )) == NULL )
+	  printf( "not enough memory for array z\n" ), exit(-1);
+      else calloc_bytes += parms.ndatamax * sizeof(float);  
+
+      if( (f = (float *)calloc( parms.ndatamax, sizeof(float) )) == NULL )
+	  printf( "not enough memory for array f\n" ), exit(-1);
+      else calloc_bytes += parms.ndatamax * sizeof(float);  
+
+      /* Extend array xx by long_extend to handle crossing
+	 longitude 180 during global interpolation.
+         Extend by 4*search radius (plus 1 in xx dimension
+         for -180 West crossing point */
+      if( all_zones ) {
+          long_extend = nint( 4 * parms.xwin2 / parms.dx );
+          /* make sure long_extend is even for equal number of
+	     bins on either side of -180 W  */
+          if( long_extend % 2 != 0 )
+              long_extend++;
+          if( (xx = (float *)calloc( parms.imax + long_extend + 1, sizeof(float) )) == NULL )
+	      printf( "not enough memory for array xx\n" ), exit(-1);
+          else calloc_bytes += (parms.imax + long_extend + 1) * sizeof(float);
+      } else {
+             if( (xx = (float *)calloc( parms.imax, sizeof(float) )) == NULL )
+	         printf( "not enough memory for array xx\n" ), exit(-1);
+             else calloc_bytes += parms.imax * sizeof(float);
+      }  
+ 
+      if( (yy = (float *)calloc( parms.jmax, sizeof(float) )) == NULL)
+	  printf( "not enough memory for array yy\n" ), exit(-1);
+      else calloc_bytes += parms.jmax * sizeof(float);  
+
+      if( (zz = (float *)calloc( parms.kmax, sizeof(float) )) == NULL )
+	  printf( "not enough memory for array zz\n" ), exit(-1);
+      else calloc_bytes += parms.kmax * sizeof(float);  
+
+      if( (ixmin = (int *)calloc( parms.ndatamax, sizeof(int) )) == NULL ) 
+	  printf( "not enough memory for array ixmin\n" ), exit(-1);
+      else calloc_bytes += parms.ndatamax * sizeof(int);  
+
+      if( (ixmax = (int *)calloc( parms.ndatamax, sizeof(int) )) == NULL ) 
+	  printf( "not enough memory for array ixmax\n" ), exit(-1);
+      else calloc_bytes += parms.ndatamax * sizeof(int);  
+
+      if( (iymin = (int *)calloc( parms.ndatamax, sizeof(int) )) == NULL ) 
+	  printf( "not enough memory for array iymin\n" ), exit(-1);
+      else calloc_bytes += parms.ndatamax * sizeof(int);  
+
+      if( (iymax = (int *)calloc( parms.ndatamax, sizeof(int) )) == NULL ) 
+	  printf( "not enough memory for array iymax\n" ), exit(-1);
+      else calloc_bytes += parms.ndatamax * sizeof(int);  
+
+      if( (izmin = (int *)calloc( parms.ndatamax, sizeof(int) )) == NULL ) 
+	  printf( "not enough memory for array izmin\n" ), exit(-1);
+      else calloc_bytes += parms.ndatamax * sizeof(int);  
+
+      if( (izmax = (int *)calloc( parms.ndatamax, sizeof(int) )) == NULL ) 
+	  printf( "not enough memory for array izmax\n" ), exit(-1);
+      else calloc_bytes += parms.ndatamax * sizeof(int);  
+
+
+      printf( "done calloc arrays\n" );  
+      printf( "total bytes used: %d (%5.1f Mb)\n\n", calloc_bytes, calloc_bytes/1e6 ); 
+
+
+      /* Determine zlo which corresponds to record number in infile for
+	 first map center and number of days since 1 Jan 1985 */
+      fileio.istart = atoi( startrec );
+ 
+      /* Open input file list and read first record */
+      if( (fid_infile = fopen( fileio.infile, "r" )) == NULL ){
+	  printf( "\nCan not open %s\n", fileio.infile ); 
+	  exit(1);
+      }
+
+      /* position to proper record in file */
+      fseek( fid_infile, (fileio.istart-1)*RECSIZE, SEEK_SET );
+      fgets( datafile, sizeof(datafile), fid_infile );
+      fclose( fid_infile );
+      basefile = basename( &datafile[0] );
+
+      strcpy( cday, "   " );  /* initialize strings */
+      strcpy( cyr, "    " ); 
+      /* e.g, parse year and jul day from 1985004h54dd-gdm.hdf */
+      strncpy( cyr, &basefile[0], 4 );
+      iyr = atoi( cyr );
+      strncpy( cday, &basefile[4], 3 );
+      iday = atoi( cday );
+
+      /* set parms.zlo: days since 1 Jan 1985 */
+      parms.zlo = (float)iday+(float)(iyr-1985)*365.;
+      if(iyr > 1988) parms.zlo = parms.zlo + 1.;
+      if(iyr > 1992) parms.zlo = parms.zlo + 1.;
+      if(iyr > 1996) parms.zlo = parms.zlo + 1.;
+      if(iyr > 2000) parms.zlo = parms.zlo + 1.;
+      if(iyr > 2004) parms.zlo = parms.zlo + 1.;
+      if(iyr > 2008) parms.zlo = parms.zlo + 1.;
+      printf( " zlo is  %f\n", parms.zlo );
+      printf( " first map centered on day: %d in year: %d\n", iday, iyr );
+
+      /* Set the first center dates */
+      cstartdate = strcat( cdate, cyr );
+      cstartdate = strcat( cdate, cday );
+      date.startdate = atoi( cstartdate );
+      printf( " start date is %d \n", date.startdate );
+
+      /* reopen file to determine last map center date */
+      fopen( fileio.infile, "r" );
+      fileio.iend = fileio.istart + (parms.kmax-1) * parms.dz;
+      fseek(fid_infile, (fileio.iend-1)*RECSIZE, SEEK_SET );
+      fgets( datafile, sizeof(datafile), fid_infile );
+      fclose( fid_infile );
+      basefile = basename( &datafile[0] );
+  
+      /* parse year and jul day  */
+      strcpy( cdate, "" );
+      strncpy( cyr, &basefile[0], 4 );
+      strncpy( cday, &basefile[4], 3 );
+      cenddate = strcat( cdate, cyr );
+      cenddate = strcat( cdate, cday );
+      date.enddate = atoi( cenddate );
+      printf( " end date is %d\n", date.enddate );
+
+      /* Now determine temporal interpolation (record) ranges */
+      /* set istart and iend so the first and last maps will span a complete
+	 temporal interpolation range if possible. e.g., if first map is centered
+	 on 1 Jan 1986 (this is date of first record to read from the 
+	 input datafile [fileio.istart]), with a search radius of 5 days (zwin2),
+         set istart back five days to 27 Dec 1985 and iend to 6 Jan 1986 */
+
+      /* set initial step back and forward size to temporal search radius */ 
+      rec_step = (int) parms.zwin2;  
+            
+      /* reopen input datafile */
+      if( (fid_infile = fopen( fileio.infile, "r" )) == NULL ){
+	  printf( "\nCan not open %s\n", fileio.infile ); 
+	  exit(1);
+      }
+
+      /* determine max bytes in infile */  
+      fseek( fid_infile, 0L, SEEK_END );
+      max_infile_size = ftell( fid_infile );
+
+      /* step back into infile; use return of fseek() to determine if 
+         file pointer is set before beginning of file  */
+      while( fseek(fid_infile, (fileio.istart-1-rec_step)*RECSIZE, SEEK_SET) != 0 ) 
+          rec_step--;
+      fileio.istart -= rec_step;
+
+      /* step forward into infile */
+      rec_step = parms.zwin2; 
+      fseek( fid_infile, (fileio.iend-1+rec_step)*RECSIZE, SEEK_SET );
+
+      /* if the file pointer is set beyond the EOF: step back */
+      while( ftell(fid_infile) >= max_infile_size ) {
+        rec_step--;  
+	fseek( fid_infile, (fileio.iend-1+rec_step)*RECSIZE, SEEK_SET );
+      }
+      fileio.iend += rec_step;
+ 
+      printf( " records: istart %d and iend %d \n\n", fileio.istart, fileio.iend );
+      fclose( fid_infile );           
+    
+      /* call interp() */
+      interp();
+      
+      free(temp), free(sl), free(vsum), free(sst_seq), free(sst_inter);
+      free(x), free(y), free(z), free(f), free(xx), free(yy), free(zz);
+      free(ixmin), free(ixmax), free(iymin), free(iymax);
+      free(izmin), free(izmax), free(in_data), free(quality_data);
+ }
+
+/* functions to allocate space for float arrays */
+/* 2d dimension */
+float **fmal2d( int ix, int iy ){
+
+ float **ival;
+ int i;
+
+ ival = (float **)calloc( ix, sizeof(float *) );
+ for( i=0; i < ix; i++ )
+     ival[i] = (float *)calloc( iy, sizeof(float) );
+
+ return(ival);
+}
+
+ /* 3d dimension */
+ float ***fmal3d( int ix, int iy, int iz ){
+
+ float ***ival;
+ int i, j;
+
+ ival = (float ***)calloc( ix, sizeof(float **) );
+ for( i=0; i < ix; i++ )
+     ival[i] = (float **)calloc( iy, sizeof(float *) );
+
+ for( i=0; i < ix; i++ )
+     for( j=0; j < iy; j++ )
+         ival[i][j] = (float *)calloc( iz, sizeof(float) );
+ 
+ return(ival);
+ }
+
+/* 4d dimension */
+float ****fmal4d( int ix, int iy, int iz, int ia ){
+     
+    float ****ival;
+    int i, j, k;
+
+    ival = (float ****)calloc( ix, sizeof(float ***) );
+    for( i=0; i < ix; i++ )
+	ival[i] = (float ***)calloc( iy, sizeof(float **) );
+    for( i=0; i < ix; i++ )
+	for( j=0; j < iy; j++)
+	    ival[i][j] = (float **)calloc( iz, sizeof(float *) );   
+    for( i=0; i < ix; i++ )
+	for( j=0; j < iy; j++)
+	    for( k=0; k < iz; k++)
+	        ival[i][j][k] = (float *)calloc( ia, sizeof(float) );   
+
+    return(ival);
+}
+
+/* function to allocate space for char array */
+/* 2d dimension */
+char **cmal2d( int ix, int iy ) {
+
+ char **ival;
+ int i;
+
+ ival = (char **)calloc( ix, sizeof(char *) );
+ for( i=0; i < ix; i++ )
+     ival[i] = (char *)calloc( iy, sizeof(char) );
+
+ return(ival);
+}
+
+
+
+usage() {
+	  printf( "	~~ gaussian interpolation of pathfinder SST data (C version) ~~\n\n " );
+	  printf(  "USAGE:  gaussinterp <infile> <beg_rec> <parmfile> <outfile>\n" );
+	  printf(  " -- infile: file containing list to process\n" );
+	  printf(  " -- beg_rec: first record (file) in infile processed\n" );
+	  printf(  " -- parmfile: input parameter file (e.g., interp.parms)\n" );
+	  printf(  " -- outfile: interpolated output file (ascii or binary)\n" );
+	  printf(  "\n e.g.,  gaussinterp list.dat 100 interp.parms output.dat\n" );
+	  printf(  "[process records of list.dat starting from record 100; interpolated \
+result is output.dat]\n" );
+	  printf(  "note: see interp.parms for example format of parameter file\n" );
+	  printf(  "      see gaussinterp.readme for more general information\n" );
+	  return(0);
+}

http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/Fortran/armstrong_interp_code.tar
----------------------------------------------------------------------
diff --git a/climatology/clim/orig/Fortran/armstrong_interp_code.tar b/climatology/clim/orig/Fortran/armstrong_interp_code.tar
new file mode 100755
index 0000000..4e125b3
Binary files /dev/null and b/climatology/clim/orig/Fortran/armstrong_interp_code.tar differ

http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/Fortran/binsum.f
----------------------------------------------------------------------
diff --git a/climatology/clim/orig/Fortran/binsum.f b/climatology/clim/orig/Fortran/binsum.f
new file mode 100644
index 0000000..96b65af
--- /dev/null
+++ b/climatology/clim/orig/Fortran/binsum.f
@@ -0,0 +1,64 @@
+c $Revision: 2.4 $
+c     subprogram binsum
+
+c-- calculate interpolation weights and vsum
+
+      subroutine binsum( nf, x, y, z, f, 
+     &			 xx, yy, zz, vsum,
+     &			 ixmin, ixmax, iymin,
+     &			 iymax, izmin, izmax, 
+     &			 imax, jmax, kmax, ndatamax )
+
+      real*4 x(ndatamax),y(ndatamax),z(ndatamax),f(ndatamax)
+      real*4 xx(imax),yy(jmax),zz(kmax)
+      real*4 vsum(2,imax,jmax,kmax)
+      integer*4 ixmin(ndatamax),ixmax(ndatamax)
+      integer*4 iymin(ndatamax),iymax(ndatamax)
+      integer*4 izmin(ndatamax),izmax(ndatamax)
+
+      common /parms/ imin,jmin,kmin,
+     &               xwin2,ywin2,zwin2,
+     &               xh,yh,zh,
+     &               dx,dy,dz,
+     &               xlo,xhi,ylo,yhi,zlo,
+     &               dxd,dyd
+
+      data xfac/-0.6931/
+
+      do n=1,nf
+        ixmin(n)=max(nint((x(n)-xwin2-xlo+0.5*dx)/dx),1)
+        ixmax(n)=min(nint((x(n)+xwin2-xlo+0.5*dx)/dx),imax)
+        iymin(n)=max(nint((y(n)-ywin2-ylo+0.5*dy)/dy),1)
+        iymax(n)=min(nint((y(n)+ywin2-ylo+0.5*dy)/dy),jmax)
+        izmin(n)=max(nint((z(n)-zwin2-zlo+0.5*dz)/dz),1)
+        izmax(n)=min(nint((z(n)+zwin2-zlo+0.5*dz)/dz),kmax)
+c        print *, x(n),y(n),z(n), f(n)
+c	print *,' ixmin, ixmax', ixmin(n), ixmax(n)
+c	print *,' iymin, iymax', iymin(n), iymax(n)
+c	print *,' izmin, izmax', izmin(n), izmax(n)
+      enddo
+
+
+
+      do n=1,nf 
+        do kk=izmin(n),izmax(n)
+          do jj=iymin(n),iymax(n)
+            do ii=ixmin(n),ixmax(n)
+             
+c- - this is the algorithm coded for weights
+
+                fac=exp( xfac*(((x(n)-xx(ii))/xh)**2
+     &                       + ((y(n)-yy(jj))/yh)**2
+     &			     + ((z(n)-zz(kk))/zh)**2) )
+
+c            print *, 'x, xx,  y, yy,  z, zz, fac f',
+c     &        x(n), xx(ii),  y(n), yy(jj), z(n), zz(kk), fac, f(n)
+
+                vsum(1,ii,jj,kk)=vsum(1,ii,jj,kk)+f(n)*fac
+                vsum(2,ii,jj,kk)=vsum(2,ii,jj,kk)+fac
+            enddo
+          enddo
+        enddo
+      enddo
+      return
+      end

http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/Fortran/interp.f
----------------------------------------------------------------------
diff --git a/climatology/clim/orig/Fortran/interp.f b/climatology/clim/orig/Fortran/interp.f
new file mode 100644
index 0000000..8977798
--- /dev/null
+++ b/climatology/clim/orig/Fortran/interp.f
@@ -0,0 +1,302 @@
+#include "interp.h"
+c $Revision: 2.12 $
+c subprogram interp.f
+c reads hdf file, extract data, call bisum(), calculates
+c and writes out interpolated maps
+
+
+      subroutine interp( temp, sst_inter, sst_seq, sl, 
+     &		  x, y, z, f,
+     &            xx, yy, zz, vsum, ixmin, ixmax,
+     &		  iymin, iymax, izmin, izmax, ndatamax,
+     &		  x_length, y_length, imax, jmax, kmax )
+
+      integer x_length,y_length
+      integer dummy
+
+c-- dimension the arrays
+      real*4 temp(x_length,y_length)
+      real*4 sl(imax,jmax,kmax)
+      real*4 sst_inter(kmax), sst_seq(imax)
+      real*4 x(ndatamax),y(ndatamax),z(ndatamax),f(ndatamax)
+      real*4 xx(imax),yy(jmax),zz(kmax)
+      real*4 vsum(2,imax,jmax,kmax)
+      real*4 ULlon, ULlat, LRlon, LRlat
+      integer*4 ixmin(ndatamax),ixmax(ndatamax)
+      integer*4 iymin(ndatamax),iymax(ndatamax)
+      integer*4 izmin(ndatamax),izmax(ndatamax)
+
+c-- hdf data array must be hardwired, no dync alloc :(
+      byte in_data( XSIZE, YSIZE )
+
+      character*80 plabel
+      character*50 infile
+      character*30 outfile
+      character*150 datafile
+      character*30 basename
+      character*35 tmpfile
+      character*10 mapformat
+      character*15 outformat
+      character*3 cday
+      character*4 cyr
+      character*200 zcatcmd
+
+      common /parms/ imin,jmin,kmin,
+     &               xwin2,ywin2,zwin2,
+     &               xh,yh,zh,
+     &               dx,dy,dz,
+     &               xlo,xhi,ylo,yhi,zlo,
+     &               dxd,dyd
+
+      common /fileio/ infile,outfile,istart,
+     &	              iend,mapformat,outformat	
+
+      common /hdfindex/ icolLeft,icolRight,irowUpper,irowLower
+
+      common /basefile/ basename
+
+c-- open output file
+      open(unit=20,status='new',form=outformat,file=outfile,
+     &     err=1000,iostat=ierr )
+      plabel='sst interpolation'
+
+c-- initialize vsum to '0's
+      do k1=1,2
+        do k2=1,kmax
+          do k3=1,jmax
+            do k4=1,imax
+               vsum(k1,k4,k3,k2)=0.
+            enddo
+          enddo
+        enddo
+      enddo
+
+c-- initialize arrays of long, lat and time intervals
+      do ii=imin,imax
+          xx(ii)=xlo+float(ii-1)*dx  ! Long to interpolate to
+      enddo
+      do jj=jmin,jmax
+          yy(jj)=ylo+float(jj-1)*dy  ! Lat to interpolate to
+      enddo
+      do kk=kmin,kmax
+          zz(kk)=zlo+float(kk-1)*dz  ! Time to interpolate to
+      enddo
+
+      do n=1,ndatamax
+          x(n)=0.0
+          y(n)=0.0
+          z(n)=0.0
+          f(n)=0.0
+	  ixmin(n) = 0
+	  ixmax(n) = 0
+	  iymin(n) = 0
+	  iymax(n) = 0
+	  izmin(n) = 0
+	  izmax(n) = 0
+      enddo
+
+c-- slope (cal) and offset to convert DN to SST
+      cal= 0.15
+      offset = 3.0
+
+c-- Open input file list
+      open(unit=15,status='old',file=infile,access='direct',
+     &  recl=RECSIZE,form='formatted',err=1500,iostat=ierr)
+
+c--  read a infile one record at a time and process . . .
+      do k=istart,iend
+         ipts=0
+         read( 15,'(a)',rec=k ) datafile
+c	 print *, datafile
+
+        icompress = 0 
+c--  hack to get trailing 'Z' (if compressed)
+c	do islash = 150, 1, -1
+c	    if ( datafile(islash:islash) .eq. 'Z' ) then
+c		icompress = 1 
+c	        goto 101	
+c	    endif
+c	enddo
+c101     continue
+
+c-- 12/27/99: C call for basename implemented
+c-- variable basename returned via common statement
+        call getbase( datafile ) 
+	print *,'\n file: ',  basename
+
+c--  cyr and cday parsed from basename 
+c--  (e.g., parsed from '1988332h54dd-gdm.hdf')
+         cday=basename( 5:7 )
+         cyr=basename( 1:4 )
+
+c-- if unix compressed, zcat to a temporary filename
+        if( basename(22:22) .eq. 'Z' ) then
+	    icompress = 1
+	    tmpfile = '/tmp/' // basename(1:20)
+	    zcatcmd = 'zcat ' // datafile( 1:len(datafile)-1 ) 
+     &                        // ' > ' // tmpfile
+	    call system( zcatcmd )
+	    datafile = tmpfile
+	endif
+
+c--  convert iday character to integer
+         read(cday,'(i3)') iday
+         read(cyr,'(i4)') iyr
+c         write(6,*)'cday = ',cday
+c         write(6,*)'cyr = ',cyr
+          
+c***> HDF call: d8gimg() 
+         retn=d8gimg(datafile,in_data,x_length,
+     &               y_length,dummy)
+
+c--  read hdf DN values and convert to SST
+         do  i=icolLeft, icolRight
+            do j=irowUpper, irowLower
+                xlat=-89.75+dyd*float(j-1)
+
+c-- center output in Pacific Ocean 
+c                if( i .lt. x_length ) ix = i-(x_length/2)
+c                if( i .le. (x_length/2) ) ix =  i+(x_length/2)
+
+c-- center output in Atlantic (default)
+                ix = i
+
+c--  convert signed byte to float
+c                if ( in_data(i,j) .lt. 0 .and. abs(xlat) .lt. 70. ) then
+                if ( in_data(i,j).lt. 0 ) then
+                    temp(ix,j)=float( in_data(i,j) ) + 256
+                else
+                    temp(ix,j)=float( in_data(i,j) )
+                endif
+
+
+
+C MULTIPLY THE PATHFINDER DIGITAL NUMBER BY THE CALIBRATION NUMBER (0.15)
+C AND ADD THE OFFSET (-3.0) TO GET DEGREES CELSIUS
+
+                 if ( temp(ix,j).gt.0 ) then
+                     ipts=ipts+1
+                     f(ipts)=( cal*temp(ix,j) ) - offset
+                     x(ipts) = ( dxd*float(ix-1) ) - 180. + dxd/2 
+                     y(ipts) = ( 90. - (dyd*float(j-1)) ) - dyd/2
+                     z(ipts)=float(iday)+float(iyr-1985)*365.
+                     if(iyr.gt.1988) z(ipts)=z(ipts)+1
+                     if(iyr.gt.1992) z(ipts)=z(ipts)+1
+                     if(iyr.gt.1996) z(ipts)=z(ipts)+1
+                     if(iyr.gt.2000) z(ipts)=z(ipts)+1
+                     if(iyr.gt.2004) z(ipts)=z(ipts)+1
+                 endif
+            enddo
+         enddo
+
+         nfnew=ipts
+         print *, ' no of pts in file ',' = ', nfnew
+
+c-- calculate interpolation weights and vsum 
+c-- arrays passed directly... a common statement 
+c-- does not seem to work.
+        call  binsum( nfnew, x, y, z, f, 
+     &                xx, yy, zz, vsum,
+     &                ixmin, ixmax, iymin,
+     &                iymax, izmin, izmax, 
+     &                imax, jmax, kmax, ndatamax )
+
+        if( icompress .eq. 1 ) call system( 'rm -f ' // tmpfile )
+c-- ..... read next hdf file from infile
+      enddo
+
+c-- all input files processed; calculate interpolated SST   
+      do kk=1,kmax
+        do jj=1,jmax
+          do ii=1,imax
+            sl(ii,jj,kk)=0.
+            if (vsum(2,ii,jj,kk).gt.0) then
+                sl(ii,jj,kk)=vsum(1,ii,jj,kk)/vsum(2,ii,jj,kk)
+            endif
+          enddo
+        enddo
+      enddo
+
+c-- write output as map "interleaved" or map "sequential"
+c-- "interleaved" is the original implementation
+c-- both formats preceded by header 
+
+c-- "interleaved" refers to each row in the output
+c-- file representing a unique lon, lat position with columns 
+c-- of associated sst values (maps) 
+c-- i.e. row one: lon1 lat1 sst1 sst2 sst3....lastsst
+c--      row two: lon2 lat1 sst1 sst2 sst3....lastsst
+
+c-- "sequential" refers to each map written as rows and
+c-- columns of lat, lon with the array element representing the
+c-- sst at that geo-position for that map.
+c-- each map is written sequentially in the file 
+
+c-- geo-positions of UL and LR corners
+          ULlon = -(180 - ( ((icolLeft-1) * dxd) + dxd/2 ))  
+	  ULlat = (90 - ( ((irowUpper-1) * dyd) + dyd/2 ))
+          LRlon = -(180 - ( (icolRight * dxd) - dxd/2 ))  
+	  LRlat = (90 - ( (irowLower * dyd) - dyd/2 ))
+
+c-- version number, "f" refers to fortran version
+      version = 'f2.9' 
+
+c-- write the 3 record header
+      if( outformat .eq. 'formatted' ) then
+          write(20,'(a)'), plabel
+          write(20,*) imax,jmax,kmax
+          write(20,*) dx,dy,dz 
+	  write(20,*) ULlon,ULlat,LRlon,LRlat
+      elseif( outformat .eq. 'unformatted' ) then
+          write(20) imax,jmax,kmax
+          write(20) dx,dy,dz
+	  write(20) ULlon,ULlat,LRlon,LRlat
+      endif
+
+      if( mapformat .eq. 'interleave' ) then
+        print *, '\n map output is interleave'
+        do jj=1,jmax
+          do ii=1,imax
+            do kk=1,kmax
+                sst_inter(kk)=sl(ii,jj,kk)
+            enddo
+            if( outformat .eq. 'formatted' ) then
+		write(20,*) ii,jj,sst_inter
+            else 
+	        write(20) ( sst_inter(i), i=1,kmax ) 
+            endif
+          enddo
+        enddo
+
+      else if( mapformat .eq. 'sequential' ) then
+        print *, '\n map output is sequential'
+        do kk=1,kmax
+          do jj=1,jmax
+            do ii=1,imax
+               sst_seq(ii)=sl(ii,jj,kk)
+             enddo
+             if( outformat .eq. 'formatted' ) then
+		 write(20,*) jj,kk,sst_seq
+             else
+		 write(20) ( sst_seq(i), i=1,imax )
+             endif
+          enddo
+        enddo
+      endif
+
+      print *,  '\ndone. . . '
+      close( 15,status='keep',err=2000,iostat=ierr )
+      close( 20,status='keep',err=2500,iostat=ierr )
+      stop
+
+ 1000 print *, 'Error opening output: ', outfile, 'error num: ', ierr
+      goto 102
+ 1500 print *, 'Error opening input: ', infile, 'error num: ', ierr
+      goto 102
+ 2000 print *, 'Error closing input: ', infile, 'error num: ', ierr
+      goto 102
+ 2500 print *, 'Error closing output: ', outfile, 'error num: ', ierr
+      goto 102
+ 102  continue
+
+      end

http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/Fortran/makefile
----------------------------------------------------------------------
diff --git a/climatology/clim/orig/Fortran/makefile b/climatology/clim/orig/Fortran/makefile
new file mode 100644
index 0000000..a624096
--- /dev/null
+++ b/climatology/clim/orig/Fortran/makefile
@@ -0,0 +1,46 @@
+# makefile for guassinterp. 
+# choose f77 or f90 compiler 
+# and appropriate flags 
+
+# f77 compiler 
+FC  = f77
+FFLAGS = -g    
+#FFLAGS = -O -n32   
+DEFINE =  -DLANG_F77
+HDF = /usr/local/hdf
+#HDF = /usr/local/hdf4.1r3-n32
+
+# f90 compiler
+#FC = f90
+#FFLAGS = -n32 -bytereclen -cpp -extend_source
+#DEFINE =  -DLANG_F90
+
+PGM = gaussinterp
+LIBS    = -L$(HDF)/lib -ldf -ljpeg -lz -lgen
+INCLUDE = -Wf,-I/usr/local/include/hdf
+BIN = /sst/vol7/PATHFINDER/bin/sgi
+
+SRC = setupinterp.f interp.f binsum.f getbase.c passbase.f
+OBJ = setupinterp.o interp.o binsum.o getbase.o passbase.o
+
+$(PGM):	$(OBJ)
+	$(FC) $(FFLAGS) -o $(PGM) $(OBJ) $(LIBS) 
+
+setupinterp.o: 
+	$(FC) $(FFLAGS) $(DEFINE) -c setupinterp.f
+interp.o:
+#	$(FC) $(DEFINE) -c interp.f
+	$(FC) $(FFLAGS) $(DEFINE) -c interp.f
+binsum.o:
+	$(FC) $(FFLAGS) $(DEFINE) -c binsum.f
+passbase.o:
+	$(FC) $(FFLAGS) $(DEFINE) -c passbase.f
+getbase.o:
+	$(CC) $(FFLAGS) -c getbase.c
+
+install: $(PGM)
+	 cp  $(PGM) $(BIN)
+
+clean:
+	 rm -f $(OBJ) $(PGM)
+

http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/Fortran/passbase.f
----------------------------------------------------------------------
diff --git a/climatology/clim/orig/Fortran/passbase.f b/climatology/clim/orig/Fortran/passbase.f
new file mode 100644
index 0000000..79e2d8e
--- /dev/null
+++ b/climatology/clim/orig/Fortran/passbase.f
@@ -0,0 +1,9 @@
+      subroutine passbase( filename )
+
+      character*25 filename
+      character*25 basename
+      common /basefile/ basename
+
+      basename = filename
+      end
+