C+
C NAME:
C	sim_MOD
C PURPOSE:
C	To make a time-dependent model at even heights to use in determining projections.
C	The model gives the 3D velocities and densities in a 360 degree model from a 
C	specified Carrington rotation interval longitude location at each height at 
C	a given time in UT (iyr, doy) from the beginning time also specified in 
C	UT (iyr, doy).
C	The subroutine is inserted into the tomography mkshiftn subroutine and
C	passes those arrays required of the tomography program that are to be used
C	in the subroutine shift_MOD required to be called by the inertial MHD
C	(or other inertial) program that provides the shifts and density and velocity
C	variations required for the tomography program to work.
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	call sim_MOD(XCbe,XCtbeg,XCtend,ALng,aNday,nLng,nLat,nMap,nT,
C		      nCar,JDCar,NCoff,VV,DD,XCshift,Vratio,Dratio,RR,dRR,FALLOFF,CLRV,CLRD)
C INPUTS:
C	XCbe(2,nT)	real		Boundary values of time maps
C	XCtbeg		real		Beginning time interval
C	XCtend		real		Ending time interval
C	Alng		real		Rotations per nLng1
C	aNday		real		# of days per time interval
C	nLng		integer		# longitude points
C	nLat		integer		# latitude points
C	nMap		integer		# heights (height begins at Sun__RAu AU)
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	RR		real		Distance above sun for reference velocity map
C	dRR		real		Distance between individual maps (in AU)
C	FALLOFF		real		Power of density falloff (usually about 2.1 for the HELIOS tomography)
C	CLRV		real		Velocity map radius filter distance
C	CLRD		real		Density map radius filter distance
C OUTPUTS:
C	XCshift(nLng,nLat,nMap,nT,3) real XCshift contains the final shifts at all heights
C					in terms of fractions of a Carrington rotation
C	Vratio (nLng,nLat,nMap,nT) real	velocity ratios
C	Dratio (nLng,nLat,nMap,nT) real	density ratios
C CALLS:
C PROCEDURE:
C >	XCshift(I,J,N,K,3) is the shift needed to trace location (I,J,N,K) back 
C	to the reference surface at a given time, i.e. 
C	XC = XCfrac(I)+XCshift(I,J,N,K,3) defines the origin of point (I,J,N,K,3)
C	at the reference surface in terms of a modified Carrington variable.
C	XCshift(I,J,N,K,3) will be negative above the reference surface (N>nRR), 
C	zero at the reference surface (N=nRR) and will point to an earlier
C	or the same time. Bad grid points in VV and DD are filled with 
C	GridSphere2D
C MODIFICATION HISTORY:
C	APR-2002, B. Jackson (UCSD)
C-
	subroutine sim_MOD(XCbe,XCtbeg,XCtend,ALng,Rconst,nLng,nLat,nMap,nT,
     &			nCar,JDCar,NCoff,VV,DD,XCshift,Vratio,Dratio,RR,dRR,FALLOFF,CLRV,CLRD)

	real		VV(nLng,nLat,nT),		! Input velocity map at height RR
     &			DD(nLng,nLat,nT),		! Input density map at height RR
     &			XCshift(nLng,nLat,nMap,nT,3),	! Map of accumulated shifts at heights (from 1 RR)
     &			Vratio (nLng,nLat,nMap,nT),	! Map of accumulated velocity ratios at heights (from 1 RR)
     &			Dratio (nLng,nLat,nMap,nT),	! Map of accumulated density ratios at heights (from 1 RR)
     &			XCbe   (2,nT),			! Carrington time intervals 
     &			XLT(2),XLTM(2)			! Internal array

	parameter	(NLNGM = 73,			! Model longitude values
     &			 NLATM = 37,			! Model latitude values
     &			 NMAPM = 31,			! Model map heights
     &			 NTM   = 50)			! Model times

	parameter	(nLngt = 109,			! Tomography longitude values
     &			nLatt  = 19)			! Tomography latitude times
     
	real		V3D((NLNGM-1)*3+1,NLATM,NMAPM,NTM), ! 3D velocity from corotating tomography
     &			D3D((NLNGM-1)*3+1,NLATM,NMAPM,NTM), ! 3D density from corotating tomography
     &			R_MOD(NMAPM),			   ! Input heights for the Model
     &			DOY,DOYL,XCTXEL,FICTXL,		   ! Day of year of model time (now and last), beginning longitude (last)
     &			V3(NLNGM,NLATM,NMAPM,3),	   ! 3D inertial velocity input for shift_MOD
     &			V3L(NLNGM,NLATM,NMAPM,3),	   ! Last 3D inertial velocity input for shift_MOD
     &			D3(NLNGM,NLATM,NMAPM),		   ! 3D inertial density input for shift_MOD
     &			V3T(NLNGM,NLATM,NMAPM,3),
     &			D3T(NLNGM,NLATM,NMAPM),
     &			V3TEST(NLNGM,NLATM),
     &			D3TEST(NLNGM,NLATM),
     &			V3TMP(NLATM),
     &			D3TMP(NLATM),
     &			D3L(NLNGM,NLATM,NMAPM)		   ! Last 3D inertial density input for shift_MOD

C	real		V3DL((NLNGM)*3+1,NLATM+1,NMAPM,3),			! These for NOT including 0, 360 -90 and +90
C     &			D3DL((NLNGM)*3+1,NLATM+1,NMAPM),

	real		V3DL((NLNGM-1)*3+1,NLATM,NMAPM,3),			! These for including 0, 360 -90 and +90
     &			D3DL((NLNGM-1)*3+1,NLATM,NMAPM),
     &			V3Dtmp(nLngt,nLatt,2,2,3),
     &			D3Dtmp(nLngt,nLatt,2,2)

	integer		IYR,IYRL,ICTXL			   ! Year of model time (now and last) (eg., 20XX)
C
	real		XCshiftL((NLNGM-1)*3+1,NTM,NLATM),! Latitude hifts needed for the tomography from the Model values
     &			XCshiftT((NLNGM-1)*3+1,NTM,NLATM),! Time shifts needed for the tomography from the Model values
     &			VratioM((NLNGM-1)*3+1,NTM,NLATM), ! Velocity ratios for the tomography from the Model values
     &			DratioM((NLNGM-1)*3+1,NTM,NLATM)  ! Density ratios for the tomography from the Model values

	real		XLtmpV((NLNGM-1)*3+1,NTM,NLATM), ! Scratch arrays
     & 			XLtmpD((NLNGM-1)*3+1,NTM,NLATM),
     &			XCshiftLS((NLNGM-1)*3+1,NLATM,NMAPM,NTM),
     &			XCshiftLSS((NLNGM-1)*3+1),
     &			XCTXE((NLNGM-1)*3+1,NTM),
     &			XLTXE((NLNGM-1)*3+1,NTM),
     &			V3AIXC((NLNGM-1)*3+1),
     &			D3AIXC((NLNGM-1)*3+1),
     &			XSAIXC((NLNGM-1)*3+1)

	real 		XCshift1(109,19,31,50,3),
     &			Vratio1(109,19,31,50),
     &			Dratio1(109,19,31,50),
     &			vv1(109,19,50),
     &			dd1(109,19,50)

	real*8		JDCar(nCar),JEpoch,JDXC

	integer		NDLT(2),NDLTM(2)			! Internal arrays

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

	Bad = BadR4()

	print*, 'Into sim_MOD'

		  nlnlnm = nLng*nLat*nMap*nT
C		  print *, nlnlnm, nlnlnm*3 
                  call ArrR4getminmax(nlnlnm*3,XCshift,aminm,amaxm)
C	          print *, 'aminm,amaxm XCshift', aminm,amaxm
                  call ArrR4getminmax(nlnlnm,Vratio,aminm,amaxm)
C	          print *, 'aminm,amaxm Vratio', aminm,amaxm
                  call ArrR4getminmax(nlnlnm,Dratio,aminm,amaxm)
C	          print *, 'aminm,amaxm Dratio ', aminm,amaxm

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

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

	NLNGMT3  = (NLNGM-1)*3 + 1

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

	NDLTM(1) = NLNGMT3
	NDLTM(2) = nT

	DLATM	 = 180./(NLATM - LAINC)
	DLATM2 = (1-LAINC)*DLATM/2.0
	BLAT   = DLATM2 - 90.0					! Beginning latitude value (near -90 degrees).	

	XLONM	 = 1./(NLNGM + (1 - LAINC))

	EPS = 0.00002

	Angl  = 45.			! Angle that the rotation information ends to the east of the sub-Sun point.
	Xangl = Angl/360.
	XCbeg1 = Xangl-1.0		! Carrington rotation begins about .85 to the west of the sub-Earth point

	speed = 400.			! Speed of the "average" solar wind to get to 1 AU.
	XCshif = -Sun__Spiral*(1.0-RRhght(1))/speed

C	print *, ' '	
C	print *, 'sim_MOD places the beginning longitude at',XCbeg1,' (west) of' 
C	print *, 'the Sun',XCshif,' rotation before the sub-Earth point at'
C	Print *, 'a lower height of',RRhght(1),' A.U.'
C	print *, ' '

	do n=1,nT						! Save all the earlier values so that they can be compared
	    do k=1,nMap
		do j=1,nLat
		    do i=1,nLng
			xcshift1(i,j,k,n,1) = xcshift(i,j,k,n,1)
			xcshift1(i,j,k,n,2) = xcshift(i,j,k,n,2)
			xcshift1(i,j,k,n,3) = xcshift(i,j,k,n,3)
			Vratio1(i,j,k,n)    = vRatio (i,j,k,n)
			Dratio1(i,j,k,n)    = dRatio (i,j,k,n)
			if (k .eq. 1) then
			    vv1(i,j,n)    = vv     (i,j,n)
			    dd1(i,j,n)    = dd     (i,j,n)

C		      if(k.eq.1.and.j.eq.10.and.N.eq.2) print *, 'sim_MOD base vv',i,j,k,N,vv1(i,j,n)

			end if
		    end do
		end do
	    end do
	end do
C
C  Model heights (begun at the same height as the tomography)
C
	do K=1,nMap			! model and tomographic heights are the same currently, but shift_MOD should not care
	  R_MOD(k) = RRhght(k)
	  nLat1 = nLat-1
	  do J=1,NLATM
	      XLA = (J - 1)*DLATM - 90.
	      XLT(2) = XLindx(XLA)
	      nLng1 = nLng-1
	      do N=1,nT			! model and tomographic times are the same currently, but shift_MOD should not care
		  XCbeg = XCbe(1,N)
		  XCend = XCbe(2,N)
		  do I=1,NLNGMT3
		      XLO = (1-I)*XLONM + XCend
		      XLT(1) = XCvar(XLO)
C	print *, 'I,J,K,NDLT(1),NDLT(2),XLT(1),XLT(2)',I,J,K,NDLT(1),NDLT(2),XLT(1),XLT(2)
		      if(k .eq. 1) then			! Make values at height 1 at the model spatial resolution
			  XLtmpV  (I,N,J) = FLINT(2,NDLT,VV(1,1,N),XLT,EPS)
			  XLtmpD  (I,N,J) = FLINT(2,NDLT,DD(1,1,N),XLT,EPS)
			  V3D(I,J,K,N)    = XLtmpV(I,N,J)
			  D3D(I,J,K,N)    = XLtmpD(I,N,J)
			  XCshiftLS(I,J,K,N) = 0.0		! Save for interpolation later
			  XCshiftL(I,N,J) = 0.0
			  XCshiftT(I,N,J) = 0.0
			  VratioM (I,N,J) = 1.0
			  DratioM (I,N,J) = 1.0

C		      if(k.eq.1.and.j.eq.19.and.N.eq.2) print *, 'sim_MOD base V3D(I,J,K,N)',i,j,k,N,V3D(I,J,K,N)

		      else				! Make values at height K at the model spatial resolution
			  XCshiftL(I,N,J) = FLINT(2,NDLT,XCshift(1,1,K,N,1),XLT,EPS)
			  XCshiftT(I,N,J) = FLINT(2,NDLT,XCshift(1,1,K,N,3),XLT,EPS)
			  VratioM (I,N,J) = FLINT(2,NDLT,Vratio(1,1,K,N),XLT,EPS)
			  DratioM (I,N,J) = FLINT(2,NDLT,Dratio(1,1,K,N),XLT,EPS)
		      end if
		  end do
	      end do
C	      call ArrR4getminmax(NLNGMT3*nT,XCshiftL(1,1,J),aminm,amaxm)
C	      print *, 'K, J, aminm,amaxm XCshiftL(1,1,J)', K, J, aminm,amaxm
C	      call ArrR4getminmax(NLNGMT3*nT,XCshiftT(1,1,J),aminm,amaxm)
C	      print *, 'K, J, aminm,amaxm XCshiftT(1,1,J)', K, J, aminm,amaxm
C	      call ArrR4getminmax(NLNGMT3*nT,VratioM(1,1,J),aminm,amaxm)
C	      print *, 'K, J, aminm,amaxm VratioM(1,1,J) ', K, J, aminm,amaxm
C	      call ArrR4getminmax(NLNGMT3*nT,DratioM(1,1,J),aminm,amaxm)
C	      print *, 'K, J, aminm,amaxm DratioM(1,1,J) ', K, J, aminm,amaxm

	      nlng1 = NLNGMT3-1
	      if(K .gt. 1) then

		  do N=1,nT
		      XCbeg = XCbe(1,N)
		      XCend = XCbe(2,N)
	              do I=1,NLNGMT3						! Determine model velocities and densities

		          dXCshiftL = XCshiftL(I,N,J) 				! Shift needed to go back to base in longitude
		          dXCshiftT = XCshiftT(I,N,J)				! Shift needed to go back to base in time

		          XLTM(1) = min(1.*NLNGMT3,XCvar(XCfull(I) + dXCshiftL)) ! Longitude on base level on scale [1,NLNGMT3]

		          XLTM(2) = max(1.,XCtvar(XCtfull(N) + dXCshiftT))	! Time on base level on scale [1,nT]

		          V3D(I,J,K,N) =  FLINT(2,NDLTM,XLtmpV(1,1,J),XLTM,EPS)*VratioM(I,N,J)
		          D3D(I,J,K,N) = (FLINT(2,NDLTM,XLtmpD(1,1,J),XLTM,EPS)*DratioM(I,N,J))*((RRhght(1)/RRhght(K))**falloff)
			  XCshiftLS(I,J,K,N) = dXCshiftL			! Save for interpolation later

		      end do
	          end do
	      end if
	  end do
	end do
C
C	do  N=1,nT
C	    nlnlnmm = NLNGMT3*NLATM*NMAPM
C	    print *, nlnlnmm
C	    call ArrR4getminmax(nlnlnmm,V3D(1,1,1,N),aminm,amaxm)
C	    print *, 'V3D N aminm,amaxm', N, aminm,amaxm
C	    call ArrR4getminmax(nlnlnmm,D3D(1,1,1,N),aminm,amaxm)
C	    print *, 'D3D N aminm,amaxm', N, aminm,amaxm
C	end do
C
C  Now into the time loop
C
	nlng1 = NLNGMT3-1
	do N=1,nT
	    XCbeg = XCbe(1,N)
	    XCend = XCbe(2,N)
	    dEarth = XCtfull(N)							! Earth longitude position difference from start
	    do k=1,NMAPM
		do j=1,NLATM					! Converts to an inertial model centered ~about Earth
		    do i=1,NLNGMT3

		      XCTXE(i,N) = XCfull(i) + dEarth + XCbeg1 + XCshif	+ 1.	! Source surface longitudes at these times
		      XLTXE(i,N) = max(1.,min(1.*NLNGMT3,XCvar(XCTXE(i,N))))	! Source surface index values at these times
		      V3AIXC(i) = FLINT(1,NLNGMT3,V3D(1,j,k,N),XLTXE(i,N),EPS)	! Velocities interpolated to these long. & times
		      D3AIXC(i) = FLINT(1,NLNGMT3,D3D(1,j,k,N),XLTXE(i,N),EPS)	! Densities interpolated to these long. & times
		      XSAIXC(i) = FLINT(1,NLNGMT3,XCshiftLS(1,j,k,N),XLTXE(i,N),EPS) ! Shifts interpolated to these long. & times

C		      if(k.eq.15.and.j.eq.19.and.N.eq.2) print *, 'sim_MOD base V3AIXC(i), XSAIXC(i)',i,j,k,N,V3AIXC(i),XSAIXC(i)
C		      if(k.eq.15.and.j.eq.10.and.N.eq.2) print *, 'sim_MOD',i,j,k,N,dEarth,XCTXE(i,N),XLTXE(i,N),V3AIXC(i),D3AIXC(i),XSAIXC(i)
C		      if(k.eq.15.and.j.eq.10.and.N.eq.32) print *,'sim_MOD',i,j,k,N,dEarth,XCTXE(i,N),XLTXE(i,N),V3AIXC(i),D3AIXC(i),XSAIXC(i)
		    end do
		    do i=1,NLNGM
		      V3(i,j,k,1) = 0.
		      V3(i,j,k,2) = 0.
		      V3(i,j,k,3) = Bad
		      D3(i,j,k)   = Bad		      
		    end do
		    do i=1,NLNGMT3-1
		      AIXCM = XCvar(XSAIXC(i)   + XCTXE(i,N))		! At higher heights shift right to get back to inertial
		      AIXCP = XCvar(XSAIXC(i+1) + XCTXE(i+1,N))		! P & M Indicies at these heights, longitudes and times
		      V3M   = V3AIXC(i)
		      V3P   = V3AIXC(i+1)
		      D3M   = D3AIXC(i)
		      D3P   = D3AIXC(i+1)
		      do II=1,NLNGMT3					! Step through each index at the location of the above
			AINDEX = II				
			if (AINDEX .ge. AIXCM .and. AINDEX .lt. AIXCP) then
			  if(i .ge. NLNGM .and. i .le. NLNGMT3-NLNGM+1) then	! the index wanted is in the interval center 
			    V3(i-NLNGM+1,j,k,3) = V3M + (AINDEX - AIXCM)*(V3P - V3M)/(AIXCP - AIXCM)
			    D3(i-NLNGM+1,j,k)   = D3M + (AINDEX - AIXCM)*(D3P - D3M)/(AIXCP - AIXCM)
C		      if(k.eq.15.and.j.eq.19.and.N.eq.2) print *, 'sim_MOD base V3(II-NLNGM+1,j,k,3)',i,II-NLNGM+1,j,k,N,V3(II-NLNGM+1,j,k,3)
C		      if(k.eq.15.and.j.eq.10.and.N.eq.2) print *,'sim_MOD',i,II,j,k,N,V3(II-NLNGM+1,j,k,3),D3(II-NLNGM+1,j,k),V3P,V3M,AIXCP,AIXCM
C		      if(k.eq.15.and.j.eq.10.and.N.eq.32) print *,'sim_MOD',i,II,j,k,N,V3(II-NLNGM+1,j,k,3),D3(II-NLNGM+1,j,k),V3P,V3M,AIXCP,AIXCM
			  end if
			end if
		      end do
		    end do
		end do
                call ArrR4getminmax(NLNGM*NLATM,V3(1,1,k,3),aminV,amaxV)
C	        print *, 'sim_MOD V3(1,1,k,3), N, k, aminV,amaxV ',N, k,aminV,amaxV
	        if(amaxV .eq. Bad) print *, 'Bad sim_MOD V3(1,1,k) max,   N, k, aminD,amaxD ',N, k,aminV,amaxV
		if(amaxV .ne. Bad) then
		    call GridSphere2D(1.,NLNGM,NLATM,1,V3(1,1,k,3),90./(nLat-1),4,0.0,90.0)	! Fill any holes
		end if
                call ArrR4getminmax(NLNGM*NLATM,D3(1,1,k),aminD,amaxD)
C	        print *, 'sim_MOD D3(1,1,k),   N, k, aminD,amaxD ',N, k,aminD,amaxD
	        if(amaxD .eq. Bad) print *, 'Bad sim_MOD D3(1,1,k) max,   N, k, aminD,amaxD ',N, k,aminD,amaxD
		if(amaxD .ne. Bad) then
		    call GridSphere2D(1.,NLNGM,NLATM,1,D3(1,1,k),90./(nLat-1),4,0.0,90.0)	! Fill any holes
		end if

                call ArrR4getminmax(NLNGM*NLATM,V3(1,1,k,3),aminm,amaxm)
	        if(aminm .eq. Bad) then
		    print *, 'sim_MOD - Still Bad values in k =',k,' V3 3'
	        end if
                call ArrR4getminmax(NLNGM*NLATM,D3(1,1,k),aminm,amaxm)
	        if(aminm .eq. Bad) then
		    print *, 'sim_MOD - Still Bad values in k =',k,' D3'
	        end if

	    end do
C
C  Determine model times to use in the subroutine shift_MOD
C
	    if(LOINC .eq. 0 .or. LAINC .eq. 0) then	! Change the model according to how the real model is to be tested
	      ISTART = 1 + (LOINC-1)
	      NLATEND = NLATM + (LAINC-1)
	      do i=ISTART,NLNGM
		do j=1,NLATM
		    V3TMP(j) = V3(i,j,1,3)
		    D3TMP(j) = D3(i,j,1)
		end do
		do j=1,NLATEND 
		    XLA = (J-1)*DLATM + BLAT
		    XLTT = XLindx(XLA)
		    V3(i,j,1,3) =  FLINT(1,NLATM,V3TMP,XLTT,EPS)	! Velocities interpolated to these long., lat and time
		    D3(i,j,1)   =  FLINT(1,NLATM,D3TMP,XLTT,EPS)	! Densities interpolated to these long., lat and time	
		end do
	      end do
	    end if

	    if(amaxV .ne. Bad .and. amaxD .ne. Bad) then
	      if(N .eq. 1) then
		IYRL   = 0
		DOYL   = 0.
		XCTXEL = 0.
		FICTXL = 0.
		ICTXL  = 0
	      end if
	      XCVarT = dEarth + XCbeg					! Time of measurement in terms of a Carrington variable
	      I = XCVarT
	      JDXC = JDCar(I) + (JDCar(I+1) - JDCar(I))*(XCVarT - I)
	      call Julian(1,IYR,DOY,JDXC,JEpoch)			! Time of measurement in year and day

	      ICTX  = XCTXE(NLNGMT3,N) + 1.0				! Carrington variable beginning at surface (a 360 deg CR)
	      FICTX = XCTXE(NLNGMT3,N) + 1.0 - ICTX			! Fraction of Carrington variable beginning at surface
	      ICTX = ICTX + NCoff					! Total Carrington variable

C	      print *, ' Before shift_MOD XCVarT, IYR, DOY, IYRL, DOYL, XCTXE, ICTX, FICTX, NCoff', XCVarT, IYR, DOY, IYRL, DOYL, XCTXE(NLNGMT3,N),ICTX, FICTX, NCoff 

	      NLNGMM = NLNGM - (1 - LOINC)
	      NLATEND = NLATM + (LAINC-1)

	      call shift_MOD(NLNGMM,NLATEND,NMAPM,NTM,V3,D3,V3L,D3L,R_MOD,IYR,DOY,IYRL,DOYL,ICTX,FICTX,ICTXL,FICTXL,
     &			     V3D,D3D,V3DL,D3DL,V3DTMP,D3DTMP,
     &			     nCar,JDCar,NCoff,xcbe,xctbeg,xctend,aLng,nLng,nLat,nMap,nT,rr,drr,Rconst,
     &			     cLrV,cLrD,vv,dd,xcshift,vRatio,dRatio)

	      call surf_MOD(NLNGM,NLATM,NMAPM,IYR,DOY,ICTX,FICTX,
     &			xcbe,xctbeg,xctend,nLng,nLat,nT,nCar,JDCar,NCoff,vv,dd,
     &			V3T,D3T)
	    end if


C The below are all tests to see that the outputs of shift_MOD are OK.
C
C  This first write checks the output of surf_MOD to determine if surf_MOD is providing the surface boundary answers to
C  the MODEL program.  The check is done by comparing the MODEL surface inputs (V3 and D3) provided by sim_MOD, with the 
C  output values of surf_MOD (V3T and D3T), which are inertial-frame (not important) MODEL-resolution inputs that are 
C  derived from the tomography program surface.  The output compared is only a timy portion of the sim_SURF output over 
C  longitude at a MODEL mid-latitude.  Every model longitude at that latitude is printed out at specific MODEL times (n).  
C  Tomographic resolution often gives the surf-MOD output a granular and sometimes non-varying look.  Portions of longitude 
C  that are not deconvolved by the tomography program remain the same in the output from surf_MOD.
C
		do i=1,NLNGM
		    do j=1,NLATM
			if(n.eq.2.and.j.eq.19) then
			    print *, 'sim-surf_MOD V3, D3',i,j,n,V3(i,j,1,3),V3T(i,j,1,3),D3(i,j,1),D3T(i,j,1)
			end if
			if(n.eq.20.and.j.eq.19) then
			    print *, 'sim-surf_MOD V3, D3',i,j,n,V3(i,j,1,3),V3T(i,j,1,3),D3(i,j,1),D3T(i,j,1)
			end if
		    end do
		end do
	end do

C  This write checks the sim_MOD input V3 and D3 at the model base to determine if sim_MOD is providing the  
C  surface boundary values that are commensurate with tomography program values at the base of the MODEL to shift_MOD.
C  The check is done by comparing the MODEL surface inputs (V3 and D3) provided by sim_MOD, with model-resolution  
C  values interpolated from the tomography program (V3TEST and D3TEST) which are located at the MODEL inertial-frame 
C  inputs.  The output compared is only a timy portion of the sim_MOD input over longitude at a MODEL mid-latitude.  
C  Model longitudes at that latitude are printed out at very selected model times (N).  Tomographic resolution may 
C  give the surf-MOD output a granular and sometimes non-varying look.  If the tomographic output can not resolve the
C  determine the input to the MODEL, then it will give a set of same values at one end of the MODEL provided to 
C  shift_MOD.    

	print *, ' '

	nLng1 = nLng-1
	nLat1 = nLat-1
	do N=1,nT
	    XCbeg = XCbe(1,N)
	    XCend = XCbe(2,N)
	    dEarth = XCtfull(N)	
	    do J=1,NLATM
	      XLA = (J - 1)*DLATM + BLAT
	      XLT(2) = XLindx(XLA)
	      do I=1,NLNGM
		XLO = dEarth + XCbeg1 + XCshif + 1. + XCbeg
		XLT(1) = XCvar(XLO + (1. - I)*XLONM)
		V3TEST(I,J) =	FLINT(2,NDLT,VV(1,1,1),XLT,EPS)		! Velocities interpolated to these long. & lats.
		D3TEST(I,J) =	FLINT(2,NDLT,DD(1,1,1),XLT,EPS)		! Densities interpolated to these long. & lats.
	        if(n.eq.2.and.j.eq.19) then
		  print *, 'sim_MOD', i,j,n,XLT(1),v3(i,j,1,3),v3test(i,j),d3(i,j,1),d3test(i,j)
		end if
	        if(n.eq.20.and.j.eq.19) then
		  print *, 'sim_MOD', i,j,n,XLT(1),v3(i,j,1,3),v3test(i,j),d3(i,j,1),d3test(i,j)
		end if
	      end do
	    end do

C  This write checks the tomographic program values of xcshift and vratio and dratio (sim_MOD 1) with the shift_MOD 
C  outputs (sim_MOD  ) of xcshift and vratio and dratio.  The output compared is only a timy portion of the shift_MOD 
C  output over longitude at a tomography mid-latitude (10) and tomography height (15 - about 1.5 AU).  Tomography 
C  longitudes at that latitude are printed out at a few selected model times (N).

	    print *, ' '

	    do k=1,nMap
		do j=1,nLat
		    do i=1,nLng
			if(n.eq.2.and.k.eq.15.and.j.eq.10) then
C			    print *, 'sim_MOD 1-blank',i,j,k,n,vv1(i,j,n),vv (i,j,n),dd1(i,j,n),dd (i,j,n)
			end if
			if(n.eq.2.and.k.eq.15.and.j.eq.10) then
			    print *, 'sim_MOD 1',i,j,k,n,xcshift1(i,j,k,n,1),xcshift1(i,j,k,n,2),xcshift1(i,j,k,n,3),Vratio1(i,j,k,n),Dratio1(i,j,k,n)
			    print *, 'sim_MOD  ',i,j,k,n,xcshift (i,j,k,n,1),xcshift (i,j,k,n,2),xcshift (i,j,k,n,3),Vratio (i,j,k,n),Dratio (i,j,k,n)
			end if
			if(n.eq.20.and.k.eq.15.and.j.eq.10) then
			    print *, 'sim_MOD 1',i,j,k,n,xcshift1(i,j,k,n,1),xcshift1(i,j,k,n,2),xcshift1(i,j,k,n,3),Vratio1(i,j,k,n),Dratio1(i,j,k,n)
			    print *, 'sim_MOD  ',i,j,k,n,xcshift (i,j,k,n,1),xcshift (i,j,k,n,2),xcshift (i,j,k,n,3),Vratio (i,j,k,n),Dratio (i,j,k,n)
			end if
		    end do
		end do
	    end do
	end do

	return
	end