SUBROUTINE OTRACE(DLAT,DLONG,DATE,LQ)
C     SUBROUTINE OTRACE:  LAMINATION PROCEDURE FOR COMPUTING O-TRACE
C     H'(F) DATA FROM N(H) PROFILE DERIVED FROM X-TRACE.  INTEGRATE
C     USING GAUSSIAN QUADRATURE WITH 3 COEFFICIENTS (EXCEPTION:  16
C     ARE USED ON LAST TWO LAMINATIONS).  THE CHANGE OF VARIABLE
C        T**2 = 1-X
C     IS USED TO KEEP THE INTEGRAND FINITE AT ALL POINTS.
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*4 DLAT,DLONG,HGT,DATE,BN,BE,BV,BT
      COMMON /CALC/ AA(50),BB(50),F(50),HP(50),FH(50),TH(70),DEN(70)
      COMMON /MISC/ DP,FH1,ASUM,BSUM,FA,CPH,X1,X2,GMT1
      COMMON /SWJTCH/ ISAT,LOCSAT,IFXS,NCASES,ICODE(4),N,METHOD,IORD,IP
      COMMON /OCALC/ FO(50),HPO(50)
      DIMENSION SS(19),W(19)
      SAVE
      DATA SS/.77459667D0,0.D0,-.77459667D0,.98940094D0,.94457502D0,
     D        .86563120D0,.75540441D0,.61787624D0,.45801678D0,
     D        .28160355D0,.09501251D0,-.09501251D0,-.28160355D0,
     D        -.45801678D0,-.61787624D0,-.75540441D0,-.86563120D0,
     D        -.94457502D0,-.98940094D0 /
      DATA W/.55555556D0,.88888889D0,.55555556D0,.02715246D0,
     D       .06225352D0,.09515851D0,.12462897D0,.14959599D0,
     D       .16915652D0,.18260342D0,.18945061D0,.18945061D0,
     D       .18260342D0,.16915652D0,.14959599D0,.12462897D0,
     D       .09515851D0,.06225352D0,.02715246D0 /
      DTF=8.98D-3
      DO 100 I=1,N
  100 FO(I)=DTF*DSQRT(DEN(I))
      HPO(1)=0.
      FS=F(1)
C     SAVE INITIAL VALUES OF PARAMETERS.
      DLT=DLAT
      DLG=DLONG
      TH1=TH(1)
      DEN1=DEN(1)
      AA2=1./AA(2)
      FACT=DEN1*DEXP(-TH1*AA2)
      M=2
      DO 1300 I=2,N
      DP=0.
      FA=FO(I)
      IF(M.EQ.3 .OR. LOCSAT.EQ.0) GO TO 1000
      CALL FRTIME(FS,TIME)
      CALL ORBINT(1,TIME,DLAT,DLONG,HGT)
      TH(1)=HGT
      DEN(1)=FACT*DEXP(TH(1)*AA2)
      IF(TH(1).GT.TH(2)) GO TO 2800
      WRITE(6,1) TH(1),FA,DLAT,DLONG,TH(2)
      TH(1)=TH(2)
      M=3
      GO TO 1000
 2800 DO 3100 L=1,I
      CALL FIELDG(DLAT,DLONG,HGT,DATE,10,LQ,BN,BE,BV,BT)
      FH(L)=BT*2.8D-5
 3100 HGT=TH(L+1)
 1000 FA2=1.0/DEN(I)
      FA1=1.0/FA
      DO 300 J=M,I
      IL=1
      IU=3
      IF(J.LT.(I-1)) GO TO 400
      IL=4
      IU=19
  400 X1=DEN(J-1)*FA2
      X2=DEN(J)*FA2
      X1LN=DLOG(X1)
      X2LN=DLOG(X2)
      Y1=FH(J-1)*FA1
      Y2=FH(J)*FA1
C     YR=Y AT REFLECTION POINT. CHANGE FROM VARIABLE X TO T.USE B9 AND
C     C9 TO CHANGE Y WITHIN EACH LAMINATION.
      T1=DSQRT(1.D0-X1)
      T2=0.
      IF(J.NE.I) T2=DSQRT(1.D0-X2)
      S1 = 0.0
      S2=0.0
      B9 = DLOG(Y2/Y1)/3.0
      V=X2LN-X1LN
      RB=BB(J)/AA(J)
      RA = 1.0+RB*V
      C9 =(DEXP(B9)-1.0)/(V*RA)
      T1=0.5*(T1+T2)
      T2=T1-T2
C      GAUSSIAN INTEGRATION
      DO 900 K=IL,IU
      T=T2*SS(K)+T1
      X=(1.0-T**2)
      RX=1.0/(1.0-X)
      IF(CPH.GE.0.0009) GO TO 1200
C      UX IS VALUE OF INTEGRAND AT THE REFLECTION POINT
 2400 UX=T*DSQRT(RX)/X
      GO TO 1100
 1200 A=X2LN-DLOG(X)
      Y = Y2/((1.0+C9*A*(1.0+RB*(2.0*V-A)))**3)
      YL = Y * CPH
      YT2=Y**2-YL**2
      R=0.5*YT2*RX/YL
      SPSI=DSIN(DATAN(R))
      A=(1.0+SPSI)*DSQRT(1.0+R**2)
      SO=1.+YL/A
      A=X/SO
      IF (A.GE.1.0D0) GO TO 1400
      IF(A.GE.0.9999999D0) GO TO 2400
      UX=T*(1.+.5*A*(1.-SO)*(1.-(1.+X)*SPSI*RX)/SO)/(X*DSQRT(1.-A))
 1100 S1=S1+UX*W(K)
  900 S2=S2+UX*W(K)*(DLOG(X)-X1LN)
      BSUM=2.0*T2
      ASUM=BSUM*S1
      BSUM=2.0*BSUM*S2
  300 DP=DP-AA(J)*ASUM-BB(J)*BSUM
 1300 HPO(I)=DP
 1500 DLAT=DLT
      DLONG=DLG
      TH(1)=TH1
      DEN(1)=DEN1
      FH(1)=FH1
      RETURN
C     RESET IORD TO OMIT ORDINARY TRACE FROM OUTPUT.
 1400 WRITE(6,2) J,I,X,SO,A
      IORD=0
      GO TO 1500
    1 FORMAT('0SATELLITE (',F6.1,' KM) BELOW FIRST LAMINATION AT',F7.3,
     F  ' MHZ ON O-TRACE.'/' FIXED POSITION ASSUMED:  LAT =',F8.2,
     F  '  LONG =',F8.2,'  HT =',F7.1)
    2 FORMAT('0FOR OTRACE LAM',I3,'  POINT',I3,'  X =',1PD11.3,'   1-Y**
     F2/(1-X) =',D11.3/' RATIO =',D11.3,'>1. CALCULATIONS TERMINATED.')
      END