CACK_OOPS_!@#$%^&*()*&^%$#@!_SPOO_KCAC READSRF.FOR
C+
C NAME:
C	MapReadSRFt
C PURPOSE:
C	Read file with synoptic map data into 2D array
C CALLING SEQUENCE:
C	I = MapReadSRF(iRot,nMax,nLng,nLat,ZRead,XCadj)
C INPUTS:
C	iRot	integer		Carrington rotation #
C	nLng	integer		# grid longitudes
C	nLat	integer		# grid latitudes
C	XCadj(3)real		only the input value of XCadj(1) is used
C				(see PROCEDURE)
C OUTPUTS:
C	I	integer		return status (0=Failure,1=Success)
C	ZRead(nLng,nLat)
C		real		2D synoptic map array
C	XCadj(3)real		plot margins used by MAP_PLOT
C CALLS:
C	BadR4, Say, iGetLogical, iSearch, iSetFileSpec, iGetFileSpec
C	nrInterpol, Str2Flt, ArrR4Constant, ArrR4Step, ArrR4Copy, SetRotations
C	bOpenFile, iFreeLun, MAP_CHECKDIM, GridReg2Reg, GinBar, bIsXEvent
C RESTRICTIONS:
C >	The MapRead* functions are called as external functions by MapRead
C	See MapRead for instructions on how to create additional MapRead*
C	functions for other data sets.
C PROCEDURE:
C >	The Stanford file contain mLngxM70=73x29 arrays for the radial
C	component of the magnetic field, at a 5 deg resolution, i.e.
C	covering [0,360] in longitude and [-70,70] in latitude
C >	input value XCadj(1):
C  >	XCadj(1) .ne. 0.:
C	If nLng .ne. 73, then the difference is interpreted as an extra margin
C	on both sides of the map, e.g. if nLng-mLng=16 (i.e. 80 deg in
C	longitude), then a margin of 40 degrees (extracted from preceeding and
C	following Carrington rotation) is added to either side of the map.
C	nLat corresponds to latitude range (nLat-1)/2*5*[-1,1].
C	If nLat > 29 then latitudes outside +/-70 are flagged with value BadR4().
C	If nLat < 29 only a partial maps is read.
C	This is the option used to overplot the source surface map on other
C	synoptic maps.
C  >	XCadj(1) = 0.:
C	The array ZRead(nLng,nLat) is assumed to refer to a single rotation
C	(i.e. [0,360] in longitude, [-90,+90] in latitude). The part of the
C	array above 70 and below -70 degree latitude will be flagged as
C	invalid. This option can be used to change the resolution of the input
C	maps (GridReg2Reg is used to change the resolution).
C MODIFICATION HISTORY:
C	JUL-1993, Paul Hick (UCSD)
C-
	function MapReadSRFt(iRot,iDay,nMax,nLng,nLat,ZRead,XCadj)
	parameter	(mLng   = 73,		! (73-1)*5= 360; # longitudes
     &			 mLat   = 37)		! (37-1)*5=2*90
	include 	'map_all.h'
	include		'filparts.h'
	include		'dirspec.h'
	include		'openfile.h'
C-------
C The files for the 2.5 Rsun maps are for a map with a grid of 73x29 points;
C the grid has a resolution of 5 degrees, i.e. 73 longitudes cover [0,360]
C and 29 latitudes cover [-70,+70].

	parameter	(M70	= 29,		! (29-1)*5=2*70
     &			 Xdeg	= 70.,		! Maximum lat
     &			dXdeg	= 2.*Xdeg/(M70-1))! Step size in lat
	real		X70(M70),
     &			Z70(M70),
     &			Y70(M70)	!scratch for interpolation
C The files for the 3.25 Rsun maps are for a map with grid of 73x30 points;
C i.e. 73 longitudes covering [0,360.], and 30 latitudes equally spaced in
C sine(latitude) between -14.5/15. and +14.5/15.

	parameter	(MSIN	= 30,
     &			 Xsin	= 14.5/15.,	! Maximum sine(lat)
     &			dXsin	= 2.*Xsin/(MSIN-1))! Step size in sine(lat)
	real		Xtmp(MSIN),		! covers 30 sine(lat)
     &			Ztmp(MSIN)

	parameter	(nSrfType =3)
	integer		iSrfType	/1/,	! Default: 2.5 Rsun
     &			iDay
	save		iSrfType

	real		Radius(nSrfType)	/2.5,3.25,1.00/
	
	character	cSrf  (nSrfType)*3	/'srf','sr3','sr1'/,
     &			cDat*3,
     &			cBarStr*27	/'Radial Field Strength (muT)'/,
     &			cWildCard*8,
     &			cWildcardt*12,
     &			CmdSrf*40

	integer		nLevDef         /11/
	real		ZLevDef(nMaxLevels)/-20,-10,-5,-2,-1,0.,1,2,5,10,20,4*0./

	character	cRec*80,
     &			cEOF*7		/'_______'/,
     &			cRotWild*4,
     &			cDayWild*3

	integer		Flt2Str,		! Function
     &			Str2Str			! Function

	logical		bMargin,
     &			bIsXEvent,
     &			bOpenFile,
     &			bSinLat			/.TRUE. / !default to true 

	MapReadSRFt = 0
	Bad = BadR4()

C-----
C KLUDGE: iSrfType=2 are the 3.25 Rsun maps with MSIN=30 steps in sine latitude
C	  iSrfType=1 are the old 2.5 Rsun maps with M70=29 steps in latitude
C	  iSrfType=3 are the photospheric maps with MSIN=29 steps in latitude

c	bSinLat = iSrfType .eq. 2 	!*tjd* instead, hard wire depending on # of col.	

	cRotWild = cWildChar(:len(cRotWild))
	cDayWild = cWildChar(:len(cDayWild))
	
	cDat = cSrf(iSrfType)
c	cWildCard = cRotWild//'.'//cDat


	cWildCardt = cRotWild//'_'//cDayWild//'.'//cDat

	print *, 'wildcardt is' , cWildCardt
c	print *, 'cFile is', cFile
					! Check existence $SRF logical
	if (iGetLogical(cEnvi//cDat,cFile) .ne. 0) then
	    iDir = iSearch(11,cFile,cFile)! Ensures that cFile is a fully
	    iDir = iSetFileSpec(cFile)	! .. qualified directory name
	    iDir = iGetFileSpec(0,FIL__DIRECTORY,cFile)
	else
	    call SetRotations(1,cWildCard,I,iDir)
	    iDir = iGetLogical(cEnvi//cDat,cFile)
	end if				! cFile(:iDir) is the data directory
					! Set up x-coordinate arrays for
	if (bSinLat) then		! .. the spline converting from sine(lat) to lat
	    call ArrR4Step(-Xdeg,dXdeg,M70 ,X70 )
	    call ArrR4Step(-Xsin,dXsin,MSIN,Xtmp)
	    do I=1,MSIN			! Convert to degrees
		Xtmp(I) = asind(Xtmp(I))
	    end do
	end if

	bMargin = XCadj(1) .ne. 0

	if (bMargin) then	! Overplot mode: check for extra longitudinal
	    iLng = nLng		! .. and latitudinal margins
	    iLat = nLat
	else			! Read 1 rotation into mLngxmLat array
	    iLng = mLng		! .. (top and bottom 4 rows=20 deg are flagged)
	    iLat = mLat
	end if

	LngOff = (iLng-mLng)/2	! Longitudinal margin (=0 for bMargin=.FALSE.)
	LatOff = (iLat-M70 )/2	!  Latitudinal margin (=4 for bMargin=.FALSE.)

	iRotOff = 0
	if (LngOff .gt. 0) iRotOff = 1		! Non-zero longitude margins:
						! .. read neighbouring rotations

	call ArrR4Constant(LatOff*iLng,Bad,ZFile)
	call ArrR4Constant(LatOff*iLng,Bad,ZFile( (iLat-LatOff)*iLng+1 ) )

	do Joff=-iRotOff,iRotOff		! If iRotOff=0 (e.g. bMargin=.FALSE.)
	    cFile(iDir+1:) = cWildCardt		! .. the loop is executed only once
	    I = index(cFile,cRotWild)
	    write (cFile(I+5:I+7),'(I3.3)') iDay
	    write (cFile(I:I+3),'(I4.4)') iRot+Joff
	    iRecl = 0

	    if (.not. bOpenFile(OPN__TRYINPUT+OPN__ONEPASS+OPN__READONLY+OPN__TEXT,iU,cFile,iRecl)) then
		if (Joff .eq. 0) return		! Main rotation missing

	    else
		inquire (iU,name=cFile)

		call Say('ReadSRF','I','Reading ',cFile)
		read (iU,'(A)',iostat=iEOF) cRec
		do while (iEOF .eq. 0 .and. cRec(:2) .ne. 'CT')	! Read past header records
		    read (iU,'(A)',iostat=iEOF) cRec
		end do

		if (iEOF .ne. 0 .or. cRec(:7) .eq. cEOF) call Say('ReadSRF','E',cFile,'is an empty file')

		read (cRec(3: 6),'(I4)') J	! Carrington rotation number
		read (cRec(8:10),'(I3)') I	! Longitude of first entry

		if (J .ne. iRot+Joff .and. .not. (J .eq. iRot+Joff+1 .and. I .eq. 360))
     &			call Say('ReadSRF','E','BAD FILE','rotation # in file does not match file name')
		if (I .ne. 360) call Say('ReadSRF','W','INCOMPLETE FILE','rotation does not start at 360 deg')

C-------	For the first Carrington rotation read records until start 
C  longitude is located
C  Successive longitude bins are labeled e.g. 1670:360, 1670:355, 1670:350, 
C  etc. The last two are 1670:005 and 1671:360. The last longitude bin is
C  translated to 1670:000

		do while (iEOF .eq. 0 .and. cRec(:7) .ne. cEOF)

		    read (cRec(3: 6),'(I4)') J	! Rotation number
		        
		    read (cRec(8:10),'(I3)') I	! Longitude

		    if (J .eq. iRot+Joff+1 .and. I .eq. 360) then I = 0  ! Last bin
				
		    L1 = 1+LngOff+I/5-Joff*(mLng-1)
	
		    if (L1 .ge. 1 .and. L1 .le. iLng) then
  
			N = 8
			call Str2Flt(cRec(11:),N,Ztmp)
			if (N .ne. 6) then
			 print *, 'ReadSRF: ouch-1'
			  print *, 'problem, N =',N
			  print *, 'while L1 =', L1, ' I =',I
			 end if
			
			I = N
			do J=1,3
			    read (iU,'(A)') cRec
			    N = 8
			    call Str2Flt(cRec,N,Ztmp(I+1))
			    if ((J .eq. 3) .and. (N .eq. 7)) then
			    	 bSinLat = .FALSE.
			    	 print *, 'in mapread file, bsinlat set to false' 
			    end if
			    I = I+N
		
			end do

C I is total # of numbers read for each longitude (29 or 30)

			do J=1,I/2      ! Reverse order
			    ZZtmp = Ztmp(I+1-J)
			    Ztmp(I+1-J) = Ztmp(J)
			    Ztmp(J) = ZZtmp
			end do

c			I = M70
c			if (bSinLat) I = MSIN

c			do J=1,N	! Move to top of Ztmp and reverse order
c			    Ztmp(I+1-J) = Ztmp(J)
c			end do		! Fill in rest of latitudes in Ztmp
c			read (iU,*) (Ztmp(J),J=I- 6,I-13,-1)
c			read (iU,*) (Ztmp(J),J=I-14,I-21,-1)
c			read (iU,*) (Ztmp(J),J=I-22,   1,-1)	! 7 or 8 numbers

			if (bSinLat) then	! Convert to M70 equal steps in latitude	
			
			    call nrInterpol(MSIN,Xtmp,Ztmp,M70,X70,Z70,Y70, 4)
			else			! Already in M70 equal steps in latitude

			    call ArrR4Copy(M70,Ztmp,Z70)
			end if
				
			do J=1,M70
			    L2 = J+LatOff
			    if (1 .le. L2 .and. L2 .le. iLat) ZFile((L2-1)*iLng+L1) = Z70(J)
			end do
		    else
			read (iU,'(A)') cRec
			read (iU,'(A)') cRec
			read (iU,'(A)') cRec

		    end if

		    read (iU,'(A)',iostat=iEOF) cRec

		end do
		
		
		
		if (Joff .le. 0 .and. L1 .ne. 1+LngOff-Joff*(mLng-1)) then	

     			call Say('ReadSRF','W','INCOMPLETE FILE','rotation does not end at 0 deg')

		end if
		iU = iFreeLun(iU)
	    end if
	
	end do

c	print *, 'Zfile(1), Zfile(2)',ZFile(1), ZFile(2)
c	print *, 'Zfile(M70*73+3),Zfile(M70*73+4),Zfile(M70*73+5)',ZFile(M70*73+3),ZFile(M70*73+4),ZFile(M70*73+5)

	if (bMargin) then
	    call ArrR4Copy(iLng*iLat,ZFile,ZRead)
	else
	    if (MapCheckDim(nMax,mLng,mLat,nLng,nLat,.TRUE.) .eq. 0) return

C-------	Average in case the spatial resolution has to be changed
C If nLng=mLng and nLat=mLat and DIST (1st argument)=1., then GridReg2Reg
C has the effect of copying ZFile into ZRead

	    call GridReg2Reg(-1.,mLng,mLat,ZFile,nLng,nLat,ZRead)
	end if

c	print *, 'Zfile(1), Zfile(mLng*mLat), Zfile(mLng*mLat+1)',Zfile(1), Zfile(mLng*mLat), Zfile(mLng*mLat+1)

	XCadj(1) = LngOff*5./360.
	XCadj(2) = -XCadj(1)
	XCadj(3) =  XCadj(1)

	MapReadSRFt = 1

c	print *, 'Zread in MapReadFILE'
c	print *, 'mLng,mLat,nLng 19 =,nLat10 =',mLng,mLat,nLng,nLat
c	print *, 'Zread(1), Zread(230), Zread(29*70)',Zread(1), Zread(230), Zread(29*70)
	
	return


c ;TJD commented out all below
c	entry iGetMapState_SRF(nLevOut,ZLevOut,RadiusOut,cArrOut,cBarOut,iExtraOut,cExtraOut)
c
c	iGetMapState_SRF = 0

c	nLevOut = nLevDef
c	call ArrR4Copy(nMaxLevels,ZLevDef,ZLevOut)

c	RadiusOut = Radius(iSrfType)

c	cRotWild = cWildChar(:len(cRotWild))
c	cArrOut = cSrf(iSrfType)//cRotWild//'.'//cRotWild(:2)
c	cBarOut = cBarStr
c
c	iExtraOut = 0
c	cExtraOut(iExtraOut) = ' '
c
c	return
c
c	entry iGetMapWildCard_SRF(cWildOut)
c
c	iGetMapWildCard_SRF = 0
c
c	cRotWild = cWildChar(:len(cRotWild))
c	cWildCard = cRotWild//'.'//cSrf(iSrfType)
c	cWildOut  = cWildCard
c
c	return
c
c	entry iCMD_SRF(ID)
c
c	iCMD_SRF = 0
c
c	iSrfOld = iSrfType
c
c	I = 0
c	do J=1,nSrfType
c	    if (J .ne. iSrfType) then
c		I = I+Int2Str(J		 ,CmdSrf(I+1:))
c		I = I+Str2Str('-'	 ,CmdSrf(I+1:))
c		I = I+Flt2Str(Radius(J),2,CmdSrf(I+1:))
c		I = I+Str2Str(' Rsun,'	 ,CmdSrf(I+1:))
c	    end if
c	end do
c	CmdSrf(I:) = ' '		! Erase trailing comma
c
c	I = -4
c	call GinBar(ID,I,30.,0.,40.,CmdSrf)
c	if (bIsXEvent(I)) then
c	    continue
c	else if (1 .le. I .and. I .le. nSrfType) then
c	    iCMD_SRF = 3
c	    iSrfType = I
c	end if
c
c	return

	end