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		X22 (NMAX)		! Put in specially for nrspline_old call
	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

	print *, 'In spliney 1a ', iFlag,NY,MY,NX	
	print *, 'In spliney 1b ', YNmin,YNmax,YMmin,YMmax

	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.

	print *, 'In spliney 2 - NY INDY(1) INDY(NY)', NY,INDY(1), INDY(NY)
	print *, '             - bFlag bMod         ', bFlag, bMod

	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

		if(I.eq.1000.and.KY.gt.1) then
		call arrR4getminmax(KY,ZS,amin,amax)
		print *, 'In spliney 3a KY ZS min max', KY,amin,amax
		end if

						! Calculate natural spline
	    if (KY .gt. 1) call nrSpline(KY,YS,ZS,1.0e30,1.0e30,X2)

		if(I.eq.1000.and.KY.gt.1) then
		  call arrR4getminmax(KY,YS,amin,amax)
		  call arrR4getminmax(KY,X2,xamin,xamax)
		  print *, 'In spliney 3b KY YS min max', KY,amin,amax
		  print *, '              KY X2 min max', KY,xamin,xamax
		  call nrSpline_old(KY,YS,ZS,1.0e30,1.0e30,X22)
		end if

	    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

		if(I.eq.1000.and.J.EQ.10) then
		call arrR4getminmax(MY,Ztmp,amin,amax)
		print *, 'In spliney 4 - I J MY Ztmp min max', I,J,MY,amin,amax
		print *, '             - bSpline KY         ', bSpline, KY
		end if

	    end do

		if(I.eq.1000) then
	  	call arrR4getminmax(MY,Ztmp,amin,amax)
		print *, 'In spliney 5 - MY Ztmp min max  ', MY,amin,amax
		end if

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

	end do

	  	call arrR4getminmax(-MY*NX,ZM,amin,amax)
		print *, 'In spliney 6 - -MY*NX ZM min max', MY*NX,amin,amax

	return
	end