      PROGRAM READICETEMP
c---------------------------------------------------------------------
c  Author:  Tim Hunter     March 2, 1994
c
c  Purpose:
c     This program reads the GLTMPICE.xxx files from Karl Schneider's
c     "Normal Water Temperature and Ice Cover of the Laurentian
c     Great Lakes: A Computer Animation, Data Base, and Analysis Tool"
c
c---------------------------------------------------------------------
c
      IMPLICIT   NONE
      INTEGER*1  DumI1, D1, M1, Buffer(26000), Dy_Mon(12)
      INTEGER*2  DumI2, RecLen, DataPoints, ImgRows, ImgCols, 
     &           CWRowSt, CWColSt, CWRowEnd, CWColEnd,
     &           Id, Flag, IErr, StMon, EndMon, BuffDat(13000), 
     &           WTemp(13000, 11), NumDays, DayCount
      INTEGER    I, J, K1, DD, MM
      LOGICAL    FindLoc, FindBth, FindTmp, FindIce
      REAL*4     DumR4, Row, Col, Lt, Ln,
     &           BMean, BStdDev, BMin, BMax, SclFact, SclSummand,
     &           TMean(12), TStdDev(12), TMin(12), TMax(12)

      CHARACTER*1  WTB(13000, 11), Ch1, Ch2, Cr
      CHARACTER*3  Bsn
      CHARACTER*5  ChMD(11)
      CHARACTER*10 Choice, MonName(12)
      CHARACTER*12 FileIn, FileOut
      CHARACTER*16 Skip16
      CHARACTER*20 Legend
      CHARACTER*30 SubT
      CHARACTER*50 Title

      EQUIVALENCE (Buffer, BuffDat)

      DATA Dy_Mon /31,28,31,30,31,30,31,31,30,31,30,31/
      DATA MonName /
     &    'January   ', 'February  ', 'March     ', 'April     ',
     &    'May       ', 'June      ', 'July      ', 'August    ',
     &    'September ', 'October   ', 'November  ', 'December  '/

      Cr = CHAR(13)

c
c Query the user to choose lake
c
      Bsn = 'xxx'
   5  Choice = '          '
      WRITE(*, 9001)
      READ(*, 9101) Choice(1:1)
      IF (Choice(1:1) .EQ. '1') Bsn = 'SUP'
      IF (Choice(1:1) .EQ. '2') Bsn = 'MIC'
      IF (Choice(1:1) .EQ. '3') Bsn = 'HUR'
      IF (Choice(1:1) .EQ. '4') Bsn = 'ERI'
      IF (Choice(1:1) .EQ. '5') Bsn = 'ONT'
      IF (Choice(1:1) .EQ. 'Q') STOP
      IF (Choice(1:1) .EQ. 'q') STOP
      IF (Bsn .EQ. 'xxx') THEN
	 WRITE(*,*) 'Invalid entry. Must be 1, 2, 3, 4, or 5.'
	 GOTO 5
      ENDIF
      WRITE(*,*)

c
c Query the user to choose data types
c
   7  Choice = '          '
      WRITE(*, 9002)
      READ(*, 9102) Choice
      DO I=1, 9
	Ch1 = Choice(I:I)
	Ch2 = Choice(I+1:I+1)
	IF ((Ch1 .LT. '1') .OR. (Ch1 .GT. '4')) THEN
	   IF ((Ch1 .EQ. ',') .OR. (Ch1 .EQ. ' ')) THEN
	      CONTINUE
	   ELSE IF (Ch1 .EQ. 'Q' .OR. Ch1 .EQ. 'q') THEN
	      STOP
	   ELSE
	      WRITE(*, *) ' Invalid entry. Must be digits 1 to 4.'
	      GOTO 7
	   ENDIF
	ELSE
	   IF ((Ch2 .GE. '1') .AND. (Ch2 .LE. '4')) THEN
	      WRITE(*, *) ' Invalid entry. Must be digits 1 to 4.'
	      GOTO 7
	   ENDIF
	ENDIF
      ENDDO

      FindLoc = .FALSE.
      FindBth = .FALSE.
      FindTmp = .FALSE.
      FindIce = .FALSE.
      DO I = 1, 10
	Ch1 = Choice(I:I)
	IF (Ch1 .EQ. '1') FindLoc = .TRUE.
	IF (Ch1 .EQ. '2') FindBth = .TRUE.
	IF (Ch1 .EQ. '3') FindTmp = .TRUE.
	IF (Ch1 .EQ. '4') FindIce = .TRUE.
      ENDDO

      IF (FindTmp .OR. FindIce) THEN
   11    WRITE(*, 9003)     
	 READ(*, 9103, ERR=13) StMon
	 IF (StMon.LT.1 .OR. StMon.GT.12) GOTO 13
   12    WRITE(*, 9004)     
	 READ(*, 9103, ERR=14) EndMon
	 IF (EndMon.LT.1 .OR. EndMon.GT.12) GOTO 14
	 IF (EndMon.LT.StMon) GOTO 15
	 GOTO 16
   13    PRINT*, ' Invalid start month.  Must be integer 1-12.'
	 GOTO 11
   14    PRINT*, ' Invalid end month.  Must be integer 1-12.'
	 GOTO 12
   15    PRINT*, ' Invalid end month.  Must be >= start month.'
	 GOTO 12
   16    CONTINUE

	 NumDays = 0
	 DO I = StMon, EndMon
	   NumDays = NumDays + Dy_Mon(I)
	 ENDDO

      ENDIF

c
c Open the data base and output file
c
      FileIn  = 'GLTMPICE.'//Bsn
      OPEN(UNIT=11, FILE=FileIn, STATUS='OLD', RECL=256,
     &     ACCESS='DIRECT', FORM='UNFORMATTED')
      READ(11, REC=1) RecLen
      CLOSE(11)
      OPEN(UNIT=11, FILE=FileIn, STATUS='OLD', RECL=RecLen,
     &     ACCESS='DIRECT', FORM='UNFORMATTED')

      FileOut = 'DATAOUT.'//Bsn
      OPEN(UNIT=12, FILE=FileOut, STATUS='OLD', ERR=20)
      CLOSE(12, STATUS='DELETE')
  20  OPEN(UNIT=12, FILE=FileOut, STATUS='NEW')

c
c Read the header info
c
      READ(11, REC=1) RecLen, DataPoints, ImgRows,ImgCols, DumI2,
     &                DumI2, DumI2, DumI2, CWRowSt, CWColSt,
     &                CWRowEnd, CWColEnd, DumR4, DumR4, DumI2, 
     &                Title, DumI2, SubT, DumI2, Legend
      WRITE(12, 1200) SubT
c
c Read the location data
c
      IF (FindLoc) THEN
	 WRITE(*, 9049) 'Location  ', Cr
	 WRITE(*, 9049) 'Location  ', Cr
	 READ(11, REC=2) (Buffer(I), I=1, RecLen)
	 READ(11, REC=3) (Buffer(I), I=RecLen+1, DataPoints*2)

	 Id = 1      ! default value if we can't set to another lake
	 IF (Bsn .EQ. 'SUP') Id = 2
	 IF (Bsn .EQ. 'MIC') Id = 3
	 IF (Bsn .EQ. 'HUR') Id = 3
	 IF (Bsn .EQ. 'ERI') Id = 4
	 IF (Bsn .EQ. 'ONT') Id = 4
	 Flag = 2

	 OPEN(UNIT=25, FILE='SCRATCH.TMP', STATUS='OLD', ERR=30)
	 CLOSE(25, STATUS='DELETE')
  30     OPEN(UNIT=25, FILE='SCRATCH.TMP', STATUS='NEW')
	 DO I = 1, DataPoints
	   Row = INT((BuffDat(I)-1) / ImgCols)    ! rows to add to CWRowSt
	   Col = BuffDat(I) - INT((Row)*ImgCols)  ! cols to add to CWColSt
	   Row = CWRowSt + Row
	   Col = CWColSt + Col
	   CALL Mercat(Lt, Ln, Col, Row, Id, Flag, IErr)
	   WRITE(25, 1231) I, Lt, Ln
	ENDDO
	REWIND(25)
      ENDIF

c
c Read the Bathymetry data (records 4 and 5).
c Use the first column of the WTemp array in order to cut RAM requirements.
c Bathymetry is the depth (in meters here) of the water under each grid point.
c
      IF (FindBth) THEN
	WRITE(*, 9049) 'Bathymetry', Cr
	WRITE(*, 9049) 'Bathymetry', Cr
	READ(11, REC=4) DumI1, DumI1, DumI2, DumI2, DumI2, 
     &                  BMean, BStdDev, BMin, BMax, SclFact, 
     &                  SclSummand, Skip16, 
     &                  (Buffer(J), J=1, DataPoints)
	READ(11, REC=5) (Buffer(J), J=DataPoints+1, DataPoints*2)
	DO I = 1, DataPoints
	   WTemp(I,1) = BuffDat(I)
	ENDDO
      ELSE
	DO I = 1, DataPoints
	  WTemp(I,1) = -99
	ENDDO
      ENDIF

c
c Write the lat/long and bathymetry data, if necessary.  Also write the
c corresponding header info.
c
      IF (FindLoc .OR. FindBth) THEN
	 WRITE(12, 1305) RecLen, DataPoints, ImgRows, ImgCols, 
     &                   CWRowSt, CWColSt, CWRowEnd, CWColEnd
      ENDIF
      IF (FindBth) WRITE(12, 1310) BMean, BStdDev, BMin, BMax
      IF (FindLoc .AND. .NOT. FindBth) THEN
	 WRITE(12, 1201) SubT, DataPoints
	 DO I = 1, DataPoints
	   READ(25, 1231) J, Lt, Ln
	   WRITE(12, 1231) I, Lt, Ln
	 ENDDO
      ENDIF
      IF (.NOT. FindLoc .AND. FindBth) THEN
	 WRITE(12, 1202) SubT, DataPoints
	 DO I = 1, DataPoints
	   WRITE(12, 1232) I, WTemp(I,1)
	 ENDDO
      ENDIF
      IF (FindLoc .AND. FindBth) THEN
	 WRITE(12, 1203) SubT, DataPoints
	 DO I = 1, DataPoints
	   READ(25, 1231) J, Lt, Ln
	   WRITE(12, 1233) I, Lt, Ln, WTemp(I,1)
	 ENDDO
      ENDIF
      IF (FindLoc .OR. FindBth) WRITE(12, '(/)')
      IF (FindLoc) CLOSE(25, STATUS='DELETE')
      IF ((.NOT. FindTmp) .AND. (.NOT. FindIce)) PRINT*

c
c Read/write the watertemp data
c
      IF (FindTmp) THEN
	 WRITE(12, 1205) SubT, DataPoints
	 WRITE(12, 1205) SubT, DataPoints
	 Id = 0    ! day index (will go from 1 to 11 only)
	 DD = 1    ! first day
	 MM = 1    ! first month
	 DayCount = 0
	 DO 100, I = 1, 365
	   IF (MM .GT. EndMon) GOTO 100
	   IF (MM .LT. StMon) THEN
	      DD = DD + 1
	      IF (DD .GT. Dy_Mon(MM)) THEN
		 DD = 1
		 MM = MM + 1
	      ENDIF
	      GOTO 100
	   ENDIF
	   IF (MM .EQ. EndMon) THEN
	      J = DD + Id
	      IF (J .GT. Dy_Mon(EndMon)) GOTO 100
	   ENDIF
	   Id = Id + 1
	   DayCount = DayCount + 1
	   WRITE(*, 9051) 'Temperature', MonName(MM), 
     &                    DayCount*100./NumDays, Cr
	   IF (Id .LE. 11) THEN
	      READ(11, REC=I+5)  D1, M1, DumI2, DumI2, DumI2, 
     &            TMean(Id), TStdDev(Id), TMin(Id), TMax(Id), 
     &            SclFact, SclSummand, Skip16,
     &            (WTB(J,Id), J=1, DataPoints)
	   ELSE
	      DO Id = 1, 11
		DO J = 1, DataPoints
		  K1 = ICHAR(WTB(J,Id))
		  IF (K1 .GT. 10) THEN
		    WTemp(J,Id) = NINT( (K1-SclSummand)/SclFact*100. )
		  ELSE
		    WTemp(J,Id) = -9999
		  ENDIF
		ENDDO
	      ENDDO
	      CALL GetDayMon(DD, MM, 11, ChMD)
	      WRITE(12, 1215) (ChMD(J), J=1,11)
	      WRITE(12, 1245) 'Mean ', (TMean(J), J=1,11)
	      WRITE(12, 1245) 'StdDv', (TStdDev(J), J=1,11)
	      WRITE(12, 1245) 'Min  ', (TMin(J), J=1,11)
	      WRITE(12, 1245) 'Max  ', (TMax(J), J=1,11)
	      DO J = 1, DataPoints
		WRITE(12, 1235) J, (WTemp(J,Id)/100., Id=1,11)
	      ENDDO
	      WRITE(12, *)
	      Id = 1
	      READ(11, REC=I+5)  D1, M1, DumI2, DumI2, DumI2, 
     &            TMean(Id), TStdDev(Id), TMin(Id), TMax(Id), 
     &            SclFact, SclSummand, Skip16,
     &            (WTB(J,Id), J=1, DataPoints)
	   ENDIF
 100     CONTINUE
	 DO I = 1, Id
	    DO J = 1, DataPoints
	      K1 = ICHAR(WTB(J,I))
	      IF (K1 .GT. 10) THEN
		WTemp(J,I) = NINT( (K1-SclSummand)/SclFact*100. )
	      ELSE
		WTemp(J,I) = -9999
	      ENDIF
	    ENDDO
	 ENDDO
	 CALL GetDayMon(DD, MM, 11, ChMD)
	 WRITE(12, 1215) (ChMD(J), J=1,Id)
	 WRITE(12, 1245) 'Mean ', (TMean(J), J=1,Id)
	 WRITE(12, 1245) 'StdDv', (TStdDev(J), J=1,Id)
	 WRITE(12, 1245) 'Min  ', (TMin(J), J=1,Id)
	 WRITE(12, 1245) 'Max  ', (TMax(J), J=1,Id)
	 DO J = 1, DataPoints
	   WRITE(12, 1235) J, (WTemp(J,I)/100., I=1,Id)
	 ENDDO
	 WRITE(12, '(/)')
	 IF (.NOT. FindIce) PRINT*
      ENDIF

c
c Read/write the icecover data (reuse the WTEMP array)
c
      IF (FindIce) THEN
	 WRITE(12, 1206) SubT, DataPoints
	 WRITE(12, 1206) SubT, DataPoints
	 Id = 0    ! day index (will go from 1 to 11 only)
	 DD = 1    ! first day
	 MM = 1    ! first month
	 DayCount = 0
	 DO 200, I = 1, 365
	   IF (MM .GT. EndMon) GOTO 200
	   IF (MM .LT. StMon) THEN
	      DD = DD + 1
	      IF (DD .GT. Dy_Mon(MM)) THEN
		 DD = 1
		 MM = MM + 1
	      ENDIF
	      GOTO 200
	   ENDIF
	   IF (MM .EQ. EndMon) THEN
	      J = DD + Id
	      IF (J .GT. Dy_Mon(EndMon)) GOTO 200
	   ENDIF
	   Id = Id + 1
	   DayCount = DayCount + 1
	   WRITE(*, 9051) 'Ice Cover', MonName(MM), 
     &                    DayCount*100./NumDays, Cr
	   IF (Id .LE. 11) THEN
	      READ(11, REC=I+5)  D1, M1, DumI2, DumI2, DumI2, 
     &            DumR4, DumR4, DumR4, DumR4, SclFact, SclSummand,
     &            Skip16, (WTB(J,Id), J=1, DataPoints)
	   ELSE
	      DO Id = 1, 11
		DO J = 1, DataPoints
		  K1 = ICHAR(WTB(J,Id))
		  IF (K1 .LE. 10 .AND. K1 .GE. 1) THEN
		    WTemp(J,Id) = (0-(K1-11)) * 10
		  ELSE
		    WTemp(J,Id) = 0
		  ENDIF
		ENDDO
	      ENDDO
	      CALL GetDayMon(DD, MM, 11, ChMD)
	      WRITE(12, 1215) (ChMD(J), J=1,11)
	      DO J = 1, DataPoints
		WRITE(12, 1236) J, (WTemp(J,Id), Id=1,11)
	      ENDDO
	      WRITE(12, *)
	      Id = 1
	      READ(11, REC=I+5)  D1, M1, DumI2, DumI2, DumI2, 
     &            DumR4, DumR4, DumR4, DumR4, SclFact, SclSummand,
     &            Skip16, (WTB(J,Id), J=1, DataPoints)
	   ENDIF
 200     CONTINUE
	 DO I = 1, Id
	    DO J = 1, DataPoints
	      K1 = ICHAR(WTB(J,I))
	      IF (K1 .LE. 10 .AND. K1 .GE. 1) THEN
		WTemp(J,I) = (0-(K1-11)) * 10
	      ELSE
		WTemp(J,I) = 0
	      ENDIF
	    ENDDO
	 ENDDO
	 CALL GetDayMon(DD, MM, 11, ChMD)
	 WRITE(12, 1215) (ChMD(J), J=1,Id)
	 DO J = 1, DataPoints
	   WRITE(12, 1236) J, (WTemp(J,I), I=1,Id)
	 ENDDO
	 PRINT* 
      ENDIF

c
c Close the files
c
      CLOSE(11)
      CLOSE(12)
      PRINT*, ' Finished!                                        '

c
c  Eliminate spurious compiler warnings
c
      Skip16 = Skip16
      DumI1 = DumI1
      DumI2 = DumI2
      DumR4 = DumR4
      D1 = D1
      M1 = M1

c
c  FORMATs
c
 1200 FORMAT('Data for ', A30)
 1201 FORMAT('The following location data is for ', A30, /,
     &       'There are ', I5, ' gridpoints.', /,
     &       'GPNum   Lat   Lon')
 1202 FORMAT('The following bathymetry data is for ', A30, /,
     &       'There are ', I5, ' gridpoints.', /,
     &       'GPNum  Bth(meters)')
 1203 FORMAT('The following location and bathymetry data is for ',
     &       A30, /, 'There are ', I5, ' gridpoints.', /,
     &       'GPNum   Lat   Lon  Bth(meters)')
 1205 FORMAT('The following water temp data (deg C) is for ', A30, /,
     &       'There are ', I5, ' gridpoints.',
     &       '   -99 indicates ice covered.')
 1206 FORMAT('The following ice cover data (%) is for ', A30, /,
     &       'There are ', I5, ' gridpoints.')

 1215 FORMAT('GPNum', 1X, 11(1X, A5))

 1231 FORMAT(I5, 2F6.2)
 1232 FORMAT(I5, I5)
 1233 FORMAT(I5, 2F6.2, I5)
 1235 FORMAT(I5, 1X, 11F6.2)
 1236 FORMAT(I5, 1X, 11I6)
 1245 FORMAT(A5, 1X, 11F6.2)

 1305 FORMAT(' Record Length=', I6, 3X, 'DataPoints   =', I6, /,
     &       ' Image Rows   =', I6, 3X, 'Image Columns=', I6, /,
     &       ' CoastWatch Start Row  =', I4, 3X, 
     &        'CoastWatch Start Column=', I4, /,
     &       ' CoastWatch End Row    =', I4, 3X, 
     &        'CoastWatch End Column=', I4)
 1310 FORMAT(' Mean Bathymetry (meters) = ', F7.2, 3X,
     &       ' Standard Dev. of Bathymetry = ', F7.2, /
     &       ' Minimum Bathymetry (m)   = ', F7.2, 3X,
     &       ' Maximum Bathymetry (m)   = ', F7.2)

 9001 FORMAT(' Which lake do you wish to process?', /,
     &       ' (1)Superior (2)Michigan (3)Huron',
     &       ' (4)Erie (5)Ontario  (Q)Quit', /,
     &       ' Enter the number [1-5,Q]: ', $)
 9002 FORMAT(' Which data do you wish to extract from the data base?',/
     &       ' (1)Locations (2)Bathymetry (3)Watertemp',
     &       ' (4)Ice Cover', /,
     &       ' Enter number(s) separated by commas ',
     &       '(or Q to quit): ', $)
 9003 FORMAT(' Please enter the starting month for temp/ice data: ', $)
 9004 FORMAT(' Please enter the ending month for temp/ice data: ', $)
 9049 FORMAT(' Now extracting ', A10, ' Data        ', A1, $)
 9050 FORMAT(' Now extracting ', A10, ' Data        ', I5, '% ', A1, $)
 9051 FORMAT(' Now extracting ', A12, ' Data for ', A10,I5, '% ',A1, $)

 9101 FORMAT(A1)
 9102 FORMAT(A10)
 9103 FORMAT(I4)

      END

C-------------------------------------------------------------------------
      SUBROUTINE GetDayMon(DD, MM, Cnt, ChMD)
      IMPLICIT  NONE
      INTEGER   DD, MM, Cnt, I, D, M
      INTEGER   Dy_Mon(12) /31,28,31,30,31,30,31,31,30,31,30,31/

      CHARACTER*5  ChMD(11)

c
c Fill the ChMD array
c
      DO I = 1, Cnt
	ChMD(I) = '  /  '
	D = DD + I - 1
	M = MM
	IF (D .GT. Dy_Mon(M)) THEN
	   D = D - Dy_Mon(M)
	   M = M + 1
	ENDIF
	IF (M .EQ. 13) M = 1
	WRITE(ChMD(I)(1:2), '(I2)') M
	WRITE(ChMD(I)(4:5), '(I2)') D
	IF (D .LT. 10) ChMD(I)(4:4) = '0'
      ENDDO

c
c Update the DD and MM variables to the day after the last day above
c
      D = D+1
      IF (D .GT. Dy_Mon(M)) THEN
	 D = D - Dy_Mon(M)
	 M = M + 1
      ENDIF
      IF (M .EQ. 13) M = 1
      DD = D
      MM = M

      END

C-------------------------------------------------------------------------
c
	subroutine mercat (plat,plon,px,py,id,flag,ierr)
c
c calculates geogr. coordinates from row and col and vice versa
c
c input parameters
c latitude, longitude if flag = 1
c row,col if flag = 2
c id = 1= 5 Lakes; 2= Superior, 3= Michigan-Huron 4= Erie-Ontario
c flag= Input =1 longitude and latitude to x,y
c              2 x,y to flogitude and latitude
c ierr = 1 wrong id
c x,y coordinates start with 1,1
c
c subroutine based on C Routine for SSTMAP
c translated to FORTRAN by K. Schneider; Date: 07/20/92
c
	implicit integer*2 (a-z)
c
	real*4 xlat,xlon,resl(4),x,y, plat,plon,px,py
	integer*2 xoff(4),yoff(4)
	
c
	real*8 radpdg,rpd2,pi4,earthr,offset
	real*8 scale,argt,ycor,xcor,row,col
c
c offset values and resolution from header IMGMAP
c
	data xoff /2854,5711,5440,5199/
	data yoff /5512,11149,11386,11386/
	data resl /3.6,1.8,1.8,1.8/
c
	IF (Flag .EQ. 1) THEN
	  xlat=plat
	  xlon=plon
	ELSE
	  x=px
	  y=py
	ENDIF

	radpdg=0.01745329252
	rpd2=0.0087266462
	pi4=0.785398163
	earthr=6371.2
	offset=4.14159203
c
	ierr=0
	if (id.lt.1.or.id.gt.4) then
	  ierr=1
	  return
	endif
c
	scale=earthr/resl(id)
c
	if (flag.eq.1) then
	  xcor=xlon*radpdg
	  if (xlat.gt.85.0) xlat=85.0
	  if (xlat.lt.-85.0) xlat=-85.0
	  argt=pi4+abs(xlat)*rpd2
	  ycor=log(tan(argt))
	  if (xlat.gt.0.0) ycor=-ycor
	  row=ycor+offset
	  col=xcor
	  y=row*scale
	  x=col*scale
	  x=xoff(id)-x
	  y=y-yoff(id)
	else
	  hemi=-1
	  x=xoff(id)-x
	  y=yoff(id)+y
	  col=x
	  row=y
	  col=col/scale
	  row=row/scale
	  xcor=col
	  ycor=row-offset
	  if (ycor.lt.0.0) then
	    hemi=1
	    ycor=-ycor
	  endif
	  argt=atan(exp(ycor))
	  xlat=(argt-pi4)/rpd2
	  if (hemi.eq.-1) xlat=-xlat
	  xlon=xcor/radpdg
	endif

	IF (Flag .EQ. 1) THEN
	  px=x
	  py=y
	ELSE
	  plat=xlat
	  plon=xlon
	ENDIF

	return
	end

