C+
C NAME:
C	surf_MOD
C PURPOSE:
C	To make a set of boundary maps of velocity and density at the bottom of the 
C	MODEL for use by the MODEL as input for the parameters (outward velocity 
C	and density) that the tomography program can provide.  
C	The subroutine is inserted into the MODEL subroutine and is passed the MODEL
C	parameters of MODEL array and the MODEL time (IYR, DOY) and the MODEL beginning
C	longitude (ICTX, FICTX).  The program assumes that both the 0 degree and 360
C	degree ends of the MODEL array are present.  The subroutine is also passed the
C	tomography program prameters that are needed and the tomography base velocities 
C	and densities.  The subroutine outputs the results needed by the MODEL program
C	as a base velocity and density at the resolution of the MODEL program and in the
C	MODEL inertial coordinates.
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	call surf_MOD(NLNGM,NLATM,NMAPM,IYR,DOY,ICTX,FICTX,
C			XCbe,XCtbeg,XCtend,nLng,nLat,nT,nCar,JDCar,NCoff,vv,dd,
C			V3,D3)
C INPUTS:
C					  THESE ARE MODEL INPUTS
C	NLNGM			integer 	# of model longitude points (inertial frame includes 0 but not 360 deg.) 
C	NLATM			integer 	# of model latitude points (inertial frame from -90 to + 90 deg.)
C	NMAPM			integer 	# of model heights 
C	IYR			integer		Year of the current model calculation (e.g., 20XX)		(e.g., 2000)
C	DOY			real		Day of year for the current model calculation (begins at 1.0 on the first day)
C	ICTX			integer 	Beginning Carrington variable source surface at this time.	(e.g., 1965)
C	FICTX			real    	Beginning Carrington variable source surface fraction at this time. (e.g., 0.123)
C
C					  THESE ARE TOMOGRAPHY INPUTS AND THEIR VALUES
C	XCbe(2,nT)		real		Boundary values of time maps
C	XCtbeg		real		Beginning time interval
C	XCtend		real		Ending time interval
C	nLng			integer	# longitude points
C	nLat			integer	# latitude points
C	nT			integer	# time intervals
C	nCar			intger	# beginning Carrington variable Julian dates
C	JDCar			real*8	Carrington variable beginning Julian dates
C	NCoff			integer	Carrington variable offset
C	vv(nLng,nLat,nT)real		Map of velocity values at a given height RR
C	dd(nLng,nLat,nT)real		Map of density values at a given height RR
C OUTPUTS:
C					  THESE ARE OUTPUTS TO THE MODEL
C	V3(NLNGM,NLATM,NMAPM,3)	real	3D velocity on source surface (level 1) from tomography model in km/sec at ICTX + FICTX
C	D3(NLNGM,NLATM,NMAPM)   real	3D density on source surface (level 1) from tomography model in km/sec at ICTX + FICT
C FUNCTIONS/SUBROUTINES:
C PROCEDURE:
C
C MODIFICATION HISTORY:
C	May-2002, B. Jackson (UCSD)
C-
	subroutine surf_MOD(NLNGM,NLATM,NMAPM,IYR,DOY,ICTX,FICTX,
     &			XCbe,XCtbeg,XCtend,nLng,nLat,nT,nCar,JDCar,NCoff,vv,dd,
     &			V3,D3)

	real		vv(nLng,nLat,nT),		! Input velocity map at height RR
     &			dd(nLng,nLat,nT),		! Input density map at height RR
     &			XCbe      (2,nT),		! Carrington time intervals 
     &			XLT(3)				! Internal array
   
	real		V3(NLNGM,NLATM,NMAPM,3),	! 3D inertial velocity input for shift_MOD
     &			D3(NLNGM,NLATM,NMAPM)		! 3D inertial density input for shift_MOD

	integer		NDLT(3)				! Internal arrays

	real*8		JDCar(nCar)

	external	EARTH

	include		'MAPCOORDINATES.H'		! Contains executable
C			statements, so must follow last declaration statement

C	print*, 'Into surf_MOD'

	NDLT(1) = nLng
	NDLT(2) = nLat
	NDLT(3) = nT

	LOINC = 1							! MODEL longitude coordinates include 0 and 360 degrees
C	LOINC = 0

	LAINC = 1							! MODEL latitude coordinates include 0 and 180 degrees
C	LAINC = 0

	DLAT	 = 180./(NLATM - LAINC)
	DLAT2	 = (1-LAINC)*DLAT/2.0
	XLON	 = 1./(NLNGM - LOINC)
	XLONE	 = (1-LOINC)*XLON

	EPS = 0.00002

	T  = XMAP_SC_POS(EARTH,IYR,DOY,nCAR,JDCar)

	N = 1
	do iT=2,nT							! Find the nearest N
	  TM = XCtfull(iT-1) + XCbe(1,iT)
	  TP = XCtfull(iT)   + XCbe(1,iT)
	  if(T .ge. TM .and. T .lt. TP) N = iT-1
	end do

	if(T .lt. XCtfull(1) + XCbe(1,1)) then				! minimum XLT(3)
	    N = 1
	    XLT(3) = 1.
	end if
	if(T .gt. XCtfull(nT) + XCbe(1,nT)) then			! maximum XLT(3)
	    N = nT
	    XLT(3) = nT
	end if

	XLT3 = XCtvar(T - XCbe(1,N))
	if(XLT3 .ge. 1. .and. XLT3 .le. nT) XLT(3) = XLT3  

	XCTXEND  = FICTX  + ICTX  - NCoff + 1.0 + XLONE			! Ending Carrington value (beginning longitude value)
	BLAT	 = DLAT2 - 90.0						! Beginning latitude value (near -90 degrees).

	XCbeg = XCbe(1,N)						! N only used here and does not matter
	XCend = XCbe(2,N)	
	do J=1,NLATM
	    XLA = (J-1)*DLAT + BLAT
	    XLT(2) = XLindx(XLA)
	    do I=1,NLNGM
		  XLO = (1 - I)*XLON + XCTXEND
		  XLT(1) = XCvar(XLO)
		  XLT(1) = max(1.,min(1.*nLng,XLT(1)))			! limit values of longitude within the tomography range

		  V3(I,J,1,3) =  FLINT(3,NDLT,vv,XLT,EPS)		! Velocities interpolated to these long., lat and time
		  D3(I,J,1)   =  FLINT(3,NDLT,dd,XLT,EPS)		! Densities interpolated to these long., lat and time
		  if(J .eq. 19) then
C		      print *, 'surf_MOD V3 and D3',i,j,XLT(1),XLT(2),XLT(3),V3(i,j,1,3),D3(i,j,1)
		  end if
	    end do
	end do

	return
	end