C+
C NAME:
C	SplineY
C PURPOSE:
C	Interpolation of 2D arrays in second (vertical) dimension.
C CATEGORY:
C	Interpolation using natural spline
C CALLING SEQUENCE:
	subroutine SplineY(iFlag,NY,YN,MY,YM,NX,ZN,ZM)
C INPUTS:
C	iFlag		integer		1st bit set: flag bad point in output ZN if
C					    nearest grid point in ZM is bad.
C
C					2nd bit set: if the range of YM extends
C					    outside the range of YN then y-coordinates outside
C					    the range of YN are mapped back inside
C					    by a mod ( max(YN)-min(YN)) operation. 
C
C	NY		integer		second (vertical) dimension of input ZN
C	YN(NY)		real		vertical coordinates for input ZN
C	MY		integer		second (vertical) dimension of output ZM
C	YM(MY)		real		vertical coordinates for output ZM
C	ZN(NX,NY)	real		fnc-values
C OUTPUTS:
C	ZM(NX,MY)	real		interpolated fnc-values
C CALLS:
C	ArrR4GetMinMax, IndexR4, BadR4, nrSplInt, nrSpline, Say
C SEE ALSO:
C	SplineX, SplineGridX, SplineGridY
C RESTRICTIONS:
C	The in- and output vertical coordinate ranges must be non-zero
C PROCEDURE:
C >	1st bit of iFlag is set: for each grid point of the output grid the
C	nearest grid point of the output grid is determined. If this nearest
C	point was flagged (contained value BadR4()) then the value of the
C	output grid is also flagged.
C	!! `Nearest' means the nearest grid point in the
C	   vertical direction only (along a single column);
C >	For each column in ZN a spline is calculated, using all valid fnc-values
C	in that column. YN provides the ordinates for the spline
C >	If the column dimension is not changed (NY=MY), SplineY and has the effect of
C	filling in the invalid function values (if 1st bit of iFlag is set)
C	In this case the same array can be used for both ZN and ZM.
C MODIFICATION HISTORY:
C	JUL-1993, Paul Hick (UCSD/CASS)
C	APR-1994, Paul Hick (UCSD/CASS)
C	    added argument bFlag to control flagging of function values in the output array
C	MAY-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C	    Generalized version of old SplineGridY. Instead of working with evenly spaced
C	    ordinate values, now arbitrary ordinate values for in- and output arrays
C	    are allowed.
C	    Changed logical argument bFlag to integer argument iFlag. The original bFlag
C	    functionality is now controlled by the first bit in iFlag. The second bit
C	    is now used to contol the mod( max(YN)-min(YN) ) behavior.
C-
	    integer	iFlag
	    integer	NY
	    real	YN(NY)
	    integer	MY
	    real	YM(MY)
	    integer	NX
	    real	ZN(NX,NY)		! Input  data array
	    real	ZM(NX,MY)		! Output data array

	character	cSay*7			/'SplineY'/

	logical		bFlag
	logical		bMod
	logical		bSpline

	parameter	(NMAX = 1000)

	real		YS  (NMAX)		! Scratch arrays for splines
	real		ZS  (NMAX)
	real		X2  (NMAX)
	real		Ztmp(NMAX)
	integer		INDY(NMAX)


	if (max(NY,MY) .gt. NMAX) call Say(cSay,'E','Not enough','scratch space. Increase parameter NMAX')


	call ArrR4GetMinMax(NY,YN,YNmin,YNmax)
	call ArrR4GetMinMax(MY,YM,YMmin,YMmax)

	if (YNmin .eq. YNmax .or. YMmin .eq. YMmax)
     &		call Say(cSay,'E','zero range','in- and output range must be non-zero')

	call IndexR4(1,NY,1,NY,YN,INDY)		! Sort into monotonic increasing order

	Bad = BadR4()

	bFlag = iand(iFlag,1) .ne. 0		! Test 1st bit
	bMod  = iand(iFlag,2) .ne. 0		! Test 2nd bit

	YNrange = YNmax-YNmin

	bSpline = .TRUE.

	do I=1,NX				! Spline each column separately
	    KY = 0

	    do J=1,NY				! Pick up all valid grid points
		JS = INDY(J)			! Process in sorted order
		if (ZN(I,JS) .ne. Bad) then
		    KY = KY+1
		    YS(KY) = YN(JS)		! YS monotonic increasing
		    ZS(KY) = ZN(I,JS)
		end if
	    end do
						! Calculate natural spline
	    if (KY .gt. 1) call nrSpline(KY,YS,ZS,1.0e30,1.0e30,X2)

	    do J=1,MY				! Calculate new grid fnc-values

		YY = YM(J)
		if (bMod) YY = YNmin+mod( mod(YY-YNmin,YNrange)+YNrange, YNrange )

		if (bFlag) then
		    jNear = 1			! Nearest input grid point
		    dNear = abs(YY-YN(1))	! Minimum distance

		    do N=2,NY
			dY = abs(YY-YN(N))
			if (dY .lt. dNear) then
			    dNear = dY
			    jNear = N
			end if
		    end do

		    bSpline = ZN(I,jNear) .ne. Bad ! Nearest grid point flagged?
		end if

		if (KY .gt. 1 .and. bSpline) then
		    call nrSplInt(KY,YS,ZS,X2, YY,Ztmp(J))
		else
		    Ztmp(J) = Bad
		end if

	    end do

	    do J=1,MY				! Copy result into ZM array
		ZM(I,J) = Ztmp(J)
	    end do

	end do

	return
	end