Document title: Fortran subroutine OAREAD to be used in program OAWATSREAD
Project: DE
NDADS Datatype: WATS_2S_ASCII
EID: SOFTWARE
Super-EID: SOFTWARE
There may be other documents also identified by this super-EID.
NDADS filename: wats_oaread_de.for
TRF entry: b47609.for in NSSDC's controlled digital document library, Dec. 1998.
Document text follows:
----------------------
C---------------------------------------------------------
C
C To read the 4 minute data in SYS$OA:OABASE.DAT
C
C OAREAD Simulation for single record return
C A. E. Hedin 11/20/84
C End of day correction 3/8/85
C Modified for 4 minute data set 7/5/85 (S.Cohn)
C
C Input parameters:
C
C IDATE - Date of data (YYDDD).
C ITIME - Time of data (Millisec).
C MAP - Array containing number of required parameters
C plus list of parameters (as on Sigma 9)
C 13 Orbit number
C 32 Altitude
C 33 Latitude
C 34 Longitude
C 37 Local Solar Time
C 38 Local Magnetic Time
C 39 L Shell
C 40 Invariant Latitude
C 76 Solar Zenith Angle
C
C Output parameters:
C
C IERR - Error return from file open
C ierr = -20 means gap in data prevented interpolation (may not occur)
C ierr = -21 means read error from OABASE.DAT (probably
C trying to access non-existant record)
C DATA - Array of data according to MAP
C
C I/O Unit 79 for input.
C
C--------------------------------------------------------------------
SUBROUTINE OAREAD(IDATE,ITIME,IERR,MAP,DATA)
PARAMETER LATITUDE = 2 ! D(2,*) ARE THE LATITUDES.
INTEGER*2 IDATA2(8)
INTEGER MAP(1), MO(76), ITIM(4), JPTS(4), ZPTS(4), IPTS(4)
INTEGER ITPOLE, J, K, L, NALT, NSZA, NIL
REAL DATA(1), D(9,4), ADD(9), DELTA, DT, DIFFL, DATAIL
LOGICAL FIRST, TIMLOG
DATA MO/12*0,9,0,17*0,1,2,3,0,0,5,8,7,6,35*0,4/
DATA IDATEL/0/, ITIM1,ITIM2/-100000,-100000/
DATA ITFDIF/60000/, ERROUT/-999./, ERRIN/9000./, IDITO/-99999/
DATA ADD/999.,999.,90.,999.,6.,999.,999.,6.,999./
DATA FIRST/.TRUE./, IPTS,JPTS,ZPTS/12*0/, DATAIL/-999./
IERR=0
JUMPKEY=0
C Open file for new date
IF(FIRST) THEN
FIRST=.FALSE.
OPEN(UNIT=79,NAME='OABASE.DAT',TYPE='OLD',
& ACCESS='KEYED',KEY=(1:4:INTEGER),READONLY,
& RECL=6,FORM='UNFORMATTED',ORGANIZATION='INDEXED',
& SHARED)
ENDIF
C Get FOUR records around requested time
ITIMEB = ITIME
ID = (IDATE-80000)
IT = (INT(ITIMEB/240000) - 1)*240 !the -1 for new quadratic fit
IF((IT*1000+240000) .EQ. ITIME) THEN
TIMLOG=.TRUE.
ELSE
TIMLOG=.FALSE.
ENDIF
IGAPSIZE=120000
IBIGDIFFO=240001
IBIGDIFFE=240001
IF(IT .LT. 0) THEN !Previous day
IT = 86400 + IT
ITIMEB = ITIMEB + 86400000
IDM = IFIX(ID*.001)*1000
IF(ID .EQ. IDM+1) THEN
ID = IDM + 365 - 1000
ELSE
ID = ID - 1
ENDIF
ENDIF
IDIT=ID*100000+IT
IF(IDIT .EQ. IDITO) GOTO 133 !Same times
IDITO = IDIT
C Input records contain: 1)altitude, 2) latitude, 3) longitude,
C 4) Solar Zenith Angle, 5) Local solar time, 6) Invariant lat,
C 7) L shell, 8) Local magnetic time, 9) Orbit number.
NALT=0 !counter for number of OK alt points
NSZA=0 !counter for number of OK SZA points
NIL=0 !counter for number of OK Inv_Lat points
NEAREST=0 !index of nearest point for ORBIT choice
READ(79,ERR=89,KEYEQ=IDIT) IDIT1, (IDATA2(K),K=1,8), IORB
ITIME1=IDIT1/100000
ITIM(1)=(IDIT1-ITIME1*100000) * 1000
IF(ABS(ITIMEB-ITIM(1)) .LE. igapsize) NEAREST=1
DO L = 1, 8
D(L,1) = FLOAT(IDATA2(L))*.1
ENDDO
D(9,1) = FLOAT(IORB)*.1
IF(D(1,1) .LT. ERRIN) THEN !alt
NALT = NALT + 1
JPTS(NALT)=1
ENDIF
IF(D(4,1) .LT. ERRIN) THEN
NSZA = NSZA + 1
ZPTS(NSZA)=1
ENDIF
IF(D(6,1) .LT. ERRIN) THEN !IL
C Adjustment for bad data due to pole crossing
IF((D(6,1) .EQ. -999.) .AND.
& (ABS(D(LATITUDE,1)) .GT. 45)) THEN
D(6,1)=87.5
ENDIF
IF(D(6,1) .NE. -999.) THEN !not bad data
NIL = NIL + 1
IPTS(NIL)=1
ENDIF
ENDIF
c Now sequential read
READ(79) IDIT2, (IDATA2(K), K = 1, 8), IORB
ITIME2=IDIT2/100000
ITIM(2)=(IDIT2-ITIME2*100000) * 1000
IF(ITIME1 .EQ. ITIME2) THEN
CONTINUE
ELSEIF(ITIME1 .EQ. ITIME2 - 1) THEN
ITIM(2) = ITIM(2) + 86400000
JUMPKEY=1
ELSEIF(((ITIME1+636) .EQ. ITIME2) .AND.
& mod(ITIME2,1000) .EQ. 1) THEN !Yr Jump
ITIM(2) = ITIM(2) + 86400000
JUMPKEY=1
ELSE
IOS=-20
GOTO 90
ENDIF
IF(ABS(ITIMEB-ITIM(2)) .LE. IGAPSIZE) NEAREST=2
IF((ITIM(2) - ITIM(1)) .GT. IBIGDIFFE) THEN
IOS=-20
GOTO 90
ENDIF
DO L = 1, 8
D(L,2) = FLOAT(IDATA2(L))*.1
ENDDO
D(9,2) = FLOAT(IORB)*.1
IF(D(1,2) .LT. ERRIN) THEN
NALT = NALT + 1
JPTS(NALT)=2
ENDIF
IF(D(4,2) .LT. ERRIN) THEN
NSZA = NSZA + 1
ZPTS(NSZA)=2
ENDIF
IF(D(6,2) .LT. ERRIN) THEN !IL
C Adjustment for bad data due to pole crossing
IF((D(6,2) .EQ. -999.) .AND.
& (ABS(D(LATITUDE,2)) .GT. 45)) THEN
D(6,2)=87.5
ENDIF
IF(D(6,2) .NE. -999.) THEN !not bad data
NIL = NIL + 1
IPTS(NIL)=2
ENDIF
ENDIF
C Third record
READ(79) IDIT3, (IDATA2(K), K = 1, 8), IORB
ITIME3=IDIT3/100000
ITIM(3)=(IDIT3-ITIME3*100000) * 1000
IF(ITIME2 .EQ. ITIME3) THEN
ITIM(3) = ITIM(3) + 86400000 * JUMPKEY
ELSEIF(ITIME2 .EQ. ITIME3 - 1) THEN
ITIM(3) = ITIM(3) + 86400000
JUMPKEY = 1
ELSEIF(((ITIME2+636) .EQ. ITIME3) .AND.
& MOD(ITIME3,1000) .EQ. 1) THEN !Yr Jump
ITIM(3) = ITIM(3) + 86400000
JUMPKEY = 1
ELSE
IOS=-20
GOTO 90
ENDIF
IF(ABS(ITIMEB-ITIM(3)) .LE. IGAPSIZE) NEAREST=3
IF((ITIM(3) - ITIM(2)) .GT. IBIGDIFFO) THEN
IOS=-20
GOTO 90
ENDIF
DO L = 1, 8
D(L,3) = FLOAT(IDATA2(L))*.1
ENDDO
D(9,3) = FLOAT(IORB)*.1
IF(D(1,3) .LT. ERRIN) THEN
NALT = NALT + 1
JPTS(NALT)=3
ENDIF
IF(D(4,3) .LT. ERRIN) THEN
NSZA = NSZA + 1
ZPTS(NSZA)=3
ENDIF
IF(D(6,3) .LT. ERRIN) THEN !IL
C Adjustment for bad data due to pole crossing
IF((D(6,3) .EQ. -999.) .AND.
& (ABS(D(LATITUDE,3)) .GT. 45)) THEN
D(6,3)=87.5
ENDIF
IF(D(6,3) .NE. -999.) THEN !not bad data
NIL = NIL + 1
IPTS(NIL)=3
ENDIF
ENDIF
C Fourth record
READ(79) IDIT4, (IDATA2(K), K = 1, 8), IORB
ITIME4=IDIT4/100000
ITIM(4)=(IDIT4-ITIME4*100000) * 1000
IF(ITIME3 .EQ. ITIME4) THEN
ITIM(4) = ITIM(4) + 86400000 * JUMPKEY
ELSEIF(ITIME3 .EQ. ITIME4 - 1) THEN
ITIM(4) = ITIM(4) + 86400000
ELSEIF(((ITIME3+636) .EQ. ITIME4) .AND.
& MOD(ITIME4,1000) .EQ. 1) THEN !Yr Jump
ITIM(4) = ITIM(4) + 86400000
ELSE
IOS=-20
GOTO 90
ENDIF
IF(ABS(ITIMEB-ITIM(4)) .LE. IGAPSIZE) NEAREST=4
IF((ITIM(4) - ITIM(3)) .GT. IBIGDIFFE) THEN
IOS=-20
GOTO 90
ENDIF
DO L = 1, 8
D(L,4) = FLOAT(IDATA2(L))*.1
ENDDO
D(9,4) = FLOAT(IORB)*.1
IF(D(1,4) .LT. ERRIN) THEN
NALT = NALT + 1
JPTS(NALT)=4
ENDIF
IF(D(4,4) .LT. ERRIN) THEN
NSZA = NSZA + 1
ZPTS(NSZA)=4
ENDIF
IF(D(6,4) .LT. ERRIN) THEN !IL
C Adjustment for bad data due to pole crossing
IF((D(6,4) .EQ. -999.) .AND.
& (ABS(D(LATITUDE,4)) .GT. 45)) THEN
D(6,4)=87.5
ENDIF
IF(D(6,4) .NE. -999.) THEN !not bad data
NIL = NIL + 1
IPTS(NIL)=4
ENDIF
ENDIF
IF(NEAREST .EQ. 0) THEN
IOS=-30
GOTO 90
ENDIF
133 DT=FLOAT(ITIM(3)-ITIM(2)) !the size of the middle time interval
C Find time of pole crossing
ITPOLE=-1
D5=ABS(D(3,3)-D(3,2)) !absolute difference in longitude
IF((D5 .GT. ADD(3) .AND. D5 .LT. 3.*ADD(3)) .OR. !>90 long. diff.
& ABS(D(2,2)) .EQ. 90. .OR. ABS(D(2,3)) .EQ. 90.) THEN !at the pole
IF(D(2,2) .LT. ERRIN .AND. D(2,3) .LT. ERRIN) THEN
XPOLE=SIGN(90.,D(2,2))
ITPOLE=DT*(XPOLE-D(2,2))/(2.*XPOLE-D(2,3)-D(2,2))+ITIM(2)
ENDIF
ENDIF
C Interpolate data
DELTA=FLOAT(ITIMEB-ITIM(2))/DT !fraction of time interval to ITIMEB
DO 20 J=1,MAP(1)
K=MO(MAP(J+1))
IF(TIMLOG .OR. K .EQ. 9) THEN !exact time or orbit number
IF(K .EQ. 4) THEN !exact time SZA, convert radians to degrees
DATA(J)=D(K,NEAREST)
ELSE
DATA(J)=D(K,NEAREST) !otherwise we interpolate below
ENDIF
ELSEIF(K .EQ. 6) THEN !find IL
IF(D(K,2) .LT. ERRIN .AND. D(K,3) .LT. ERRIN .AND.
& D(K,2) .NE. -999 .AND. D(K,3) .NE. -999) THEN !General interp
IF(NIL .GT. 2) THEN !Lagrange's quadratic interpolation
DATAIL=F_LAGRANGE(ITIM,D,6,ITIMEB,NIL,IPTS)
ELSE !linear interpolation
C Adjust for crossing a pole and any subsequent wraparound
IF(ITPOLE .GT. -1) THEN !pole crossing
DATAIL=D(K,2)+DELTA*(180.-D(K,2)-D(K,3))
IF(DATAIL .GT. 90) DATAIL=180.-DATAIL !wraparound
ELSE
DATAIL=D(K,2)+DELTA*(XU-D(K,2))
ENDIF
ENDIF
ELSE !Fill data input
DATAIL=AMIN1(D(K,2),D(K,3)) !will probably equal -999 here
ENDIF
DATA(J)=DATAIL
ELSEIF(K .EQ. 7) THEN !L Shell
IF(DATAIL .NE. -999) THEN !requires that Inlat be earlier in MAP
DATA(J)=(1/COS(DATAIL*0.01745329))**2
ELSE !a bad data point in the data base or Inlat not done 1st
DATA(J)=DATAIL
ENDIF
ELSEIF(K .GT. 0 .AND.
& D(K,2) .LT. ERRIN .AND. D(K,3) .LT. ERRIN) THEN !no cRaZy data
XU=D(K,3) !use this for the end point for linear interpolation
C here do altitude fit (K=1)
IF(K .EQ. 1) THEN !altitude interpolation
IF(NALT .GT. 2) THEN !Lagrange's quadratic interp
DATA(J)=F_LAGRANGE(ITIM,D,K,ITIMEB,NALT,JPTS)
ELSEIF(NALT .EQ. 2) THEN !linear interpolation
DATA(J)=DELTA*(XU-D(K,2))+D(K,2)
ELSE !error default value
DATA(J)=-999
ENDIF
ELSEIF(K .EQ. 2) THEN !linear interpolation of LATITUDE
IF(ITPOLE .GT. -1) THEN !polar Xing
DATA(J)=D(K,2)+DELTA*(2.*XPOLE-XU-D(K,2))
IF(ABS(DATA(J)) .GT. 90.) DATA(J)=2.*XPOLE-DATA(J)
ELSE !not a polar Xing
DATA(J)=D(K,2)+DELTA*(XU-D(K,2))
ENDIF
ELSEIF(K .EQ. 3 .AND. ITPOLE .EQ. -1) THEN !LONG., no pole Xing
IF(XU*D(K,2) .LT. 0.) THEN !crosses zero degrees long.?
DIFFL=ABS(XU-D(K,2))
IF(DIFFL .GT. 180.) THEN !probably Xes +/-180 long. instead
IF(D(K,2) .LT. 0.) THEN
DATA(J)=D(K,2)-DELTA*(360-DIFFL) !use adjusted DIFFL
IF(DATA(J) .LT. -180.) DATA(J)=DATA(J)+360.
ELSE !D(K,2)>0.
DATA(J)=D(K,2)+DELTA*(360-DIFFL) !use adjusted DIFFL
IF(DATA(J) .GT. 180.) DATA(J)=DATA(J)-360.
ENDIF
ELSE !zero degree crossing (Xing)
DATA(J)=D(K,2)+DELTA*(XU-D(K,2))
ENDIF
ELSE !treat as normal
DATA(J)=D(K,2)+DELTA*(XU-D(K,2))
ENDIF
ELSEIF(K .EQ. 4) THEN !SZA
IF(NSZA .GT. 2) THEN !quadratic interpolation
DATA(J)=F_LAGRANGE(ITIM,D,K,ITIMEB,NSZA,ZPTS)
ELSEIF(NSZA .EQ. 2) THEN !linear interpolation
DATA(J)=D(K,2)+DELTA*(D(K,3)-D(K,2))
ELSE !error default value
DATA(J)=-999
ENDIF
C Here get LST when not crossing a pole or LMT
ELSEIF((K .EQ. 5 .AND. ITPOLE .EQ. -1) .OR. K .EQ. 8) THEN
DIFFL=ABS(XU-D(K,2))
IF(DIFFL .GT. 12.) THEN !probably crosses 24/0 hours
IF(D(K,2) .LT. 12.) THEN
DATA(J)=D(K,2)-DELTA*(24.-DIFFL) !use adjusted DIFFL
IF(DATA(J) .LT. 0.) DATA(J)=DATA(J)+24.
ELSE !D(K,2)>12.
DATA(J)=D(K,2)+DELTA*(24.-DIFFL) !use adjusted DIFFL
IF(DATA(J) .GT. 24.) DATA(J)=DATA(J)-24.
ENDIF
ELSE !doesn't cross LS or LM midnight
DATA(J)=D(K,2)+DELTA*(XU-D(K,2))
ENDIF
C At a polar crossing event, we presume that LONG and LST don't
C change slowly enough to interpolate and will have the same value as the
C recorded data (D(K,i)) before or after the pole crossing, as determined
C by comparing the time we're interpolating (ITIMEB) vs. the time of
C crossing (ITPOLE)
ELSEIF((K .EQ. 3 .OR. K .EQ. 5) .AND.
& ITPOLE .GT. -1) THEN
IF(ITIMEB .GE. ITPOLE) THEN !if near but past pole
DATA(J)=D(K,3)
ELSE !not past pole yet
DATA(J)=D(K,2)
ENDIF
ELSE
DATA(J)=ERROUT
ENDIF !of "IF(K .EQ. 1..."
ELSE
DATA(J)=ERROUT
ENDIF !of "IF(TIMLOG..."
20 ENDDO
RETURN
C Error exit
89 IOS=-21
90 IERR=IOS
DO J=1,MAP(1)
DATA(J)=ERROUT
ENDDO
RETURN
END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Function to interpolate ALTITUDE from 3 or 4 points.
C For use with OAREAD4
FUNCTION F_LAGRANGE(ITIM,D,K,ITIMEB,NUM,JPTS)
DIMENSION ITIM(4),D(9,4),JPTS(4)
REAL*8 Li,A
F_LAGRANGE = 0.
T0 = FLOAT(ITIMEB)
DO I=1,NUM
A = D(K,JPTS(I))
Ti = FLOAT(ITIM(JPTS(I)))
Li = 1.
DO J=1,NUM
T = FLOAT(ITIM(JPTS(J)))
IF( I.ne.J ) Li = Li * (T0-T)/(Ti-T)
ENDDO
F_LAGRANGE = F_LAGRANGE + A*Li
ENDDO
RETURN
END