      SUBROUTINE XLAM(M,I,DLAT,DLONG,DATE,LQ,IFLAG,*)
C     SUBROUTINE XLAM: LAMINATION PROCEDURE FOR INVERTING X-TRACE H'(F)
C     DATA TO OBTAIN N(H) PROFILE. FOR PARABOLIC IN LOG(N) LAMINATIONS,
C     THE EQUATION USED FOR THE JTH INTERVAL IS:
C        H = H(J-1) + AA(J)*LN(N/N(J-1)) + BB(J)*LN(N/N(J-1))**2
C     WHERE AA(I),BB(I) ARE THE CORRESPONDING SLOPES. FOR LINEAR IN
C     LOG(N) LAMINATIONS,BB(I)=0. INTEGRATE BY USING GAUSSIAN QUADRATURE
C     WITH 3 COEFFICIENTS (EXCEPTION: 7 COEFFICIENTS USED WHEN Y2>0.93).
C     THE CHANGE OF VARIABLE:
C        T**2 = 1 -X/(1-YR)
C     IS USED TO KEEP THE INTEGRAND FINITE AT ALL POINTS.
C     IFLAG=1 FOR USUAL COMPUTATIONS; =2 FOR COMPUTATIONS OF FXSMAX;
C     =3 FOR COMPUTATIONS OF FXS WHICH YIELDS AA(3)=AA(2).
C     PROCESSES ONE H'(F) POINT PER CALL. RESULTS RETURNED VIA COMMON.
      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
      DIMENSION SS(10),W(10),DIFF(20)
      SAVE
      DATA SS/.77459667D0,0.D0,-.77459667D0,.94910791D0,.74153119D0,
     D         .40584515D0,0.D0,-.40584515D0,-.74153119D0,-.94910791D0/
      DATA W/.55555556D0,.88888889D0,.55555556D0,.12948497D0,.27970539D0
     D     ,.38183005D0,.41795918D0,.38183005D0,.27970539D0,.12948497D0/
      DATA I1,I2/1,2/
      IF(IFXS.EQ.1 .AND. I.EQ.4) AAI=AA(3)
      MM=M
      IQ=0
      ILIM=2
      IF (IFLAG.NE.1) ILIM=3
      FA2=1.0/(1.24D4*FA**2)
      FA1=1.0/FA
      GIRO=FH(I)
      TH(I) =TH(I-1)
      HGT=TH(I-1)
  200 HGTA = TH(I)
      YR=FH(I)*FA1
      YRM1=1.0-YR
      RYR=1.0/YRM1
      DEN(I) = (F(I)*(F(I)-FH(I)))*12400.00
      IF(DEN(I).GT.(DEN(I-1)*1.0005)) GO TO 2600
      J=I-1
      IF(IFLAG.NE.3) WRITE(6,1)I,DEN(I),J,DEN(J)
      GO TO 2500
 2600 DO 300 J=M,I
      Y1=FH(J-1)*FA1
      Y2=FH(J)*FA1
      IL=1
      IU=3
      IF (Y2.LT. 0.93) GO TO 100
      IL=4
      IU=10
  100 X1=DEN(J-1)*FA2
      X2=DEN(J)*FA2
      X1LN=DLOG(X1)
      X2LN=DLOG(X2)
C     YR=Y AT REFLECTION POINT. CHANGE FROM VARIABLE X TO T.USE B9 AND
C     C9 TO CHANGE Y WITHIN EACH LAMINATION.
      A=X1*RYR
      IF(A.LT.1.D0) GO TO 2800
      WRITE(6,3) I1,A,J,I
      GO TO 2500
 2800 T1=DSQRT(1.D0-A)
      S1 = 0.0
      S2=0.0
      B9 = DLOG(Y2/Y1)/3.0
      V=X2LN-X1LN
C     CALCULATION FOR LAST LAMINATION.
      IF (J .LT. I) GO TO 400
      T2 = 0.0
      IF (METHOD .GT. 1 .AND. I .GT. ILIM .AND. IQ .GE. 2) GO TO 500
C     LAMINATIONS LINEAR IN LOG(N).
  800 ALIM = 1.0
      RA = 1.0
      RB = 0.0
      GO TO 600
C     LAMINATIONS PARABOLIC IN LOG(N).
  500 AA(I) = AAI
      BB(I) = (TH(I)-TH(I-1)-AA(I)*V)/(V*V)
      AA(I+1)=AA(I)+(2.0)*BB(I)*V
      ALIM = AA(I+1)/AA(I)
      GO TO 700
  400 A=X2*RYR
      IF(A.LT.1.D0) GO TO 2900
      WRITE(6,3) I2,A,J,I
      GO TO 2500
C     CALCULATION FOR ALL LAMINATIONS EXCEPT THE LAST.
 2900 T2=DSQRT(1.D0-A)
      IF (METHOD .EQ. 1) GO TO 800
      ALIM = 1.0
  700 RB = BB(J)/AA(J)
      RA = 1.0+RB*V
  600 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)*YRM1
      RX=1.0/(1.0-X)
      A=X2LN-DLOG(X)
      Y = Y2/((1.0+C9*A*(1.0+RB*(2.0*V-A)))**3)
      IF (T .GT. 0.0002) GO TO 1000
C      UX IS VALUE OF INTEGRAND AT THE REFLECTION POINT
 2400 UX=.707107*(2.0-YR)*RYR**2/(DSQRT(1.0+CPH**2)*DSQRT(1.0+3.0*C9*YR
     U    *ALIM*RYR))
      GO TO 1100
 1000 YL = Y * CPH
      YT2=Y**2-YL**2
      IF (CPH .LT. 0.0009) GO TO 1200
      R=0.5*YT2*RX/YL
      SPSI=DSIN(DATAN(R))
      A=(1.0+SPSI)*DSQRT(1.0+R**2)
      SX = 1.0-A*YL
      GO TO 1300
C     PROPAGATION ANGLE NEAR 90 DEGREES (TRANSVERSE PROPAGATION).
 1200 SPSI=1.0D0
      SX=1.0-RX*Y**2
 1300 A=X/SX
      IF (A.GE.1.0D0) GO TO 1400
      IF(A.GE.0.9999999D0) GO TO 2400
C     REJECTS BAD DATA CARDS
      UX=T*(1.+.5*A*(1.-SX)*(1.+(1.+X)*SPSI*RX)/SX)/(X*DSQRT(1.-A))
 1100 S1=S1+UX*W(K)
  900 S2=S2+UX*W(K)*(DLOG(X)-X1LN)
      BSUM=2.0*T2*YRM1
      ASUM=BSUM*S1
      BSUM=2.0*BSUM*S2
      IF (J .GE. I) GO TO 1500
  300 DP=DP-AA(J)*ASUM-BB(J)*BSUM
 1500 IF (IFLAG .EQ. 2) RETURN
C     ITERATE ON LAST LAMINATION USING IMPROVED TECHNIQUE.
      IF (I .GT. 2) GO TO 1600
      AA(2)=-HP(2)/ASUM
      TH(2)=TH(1)+AA(2)*V
      AAI = AA(2)
      GO TO 1800
 1600 IF (METHOD .EQ. 2 .AND. IQ .GE. 2) GO TO 1900
      AA(I)=(DP-HP(I))/ASUM
      TH(I)=TH(I-1)+AA(I)*V
      BB(I)=0.0
      HGT=TH(I)
      IF (METHOD-1) 1800,1800,2000
 1900 BB(I)=(DP-AA(I)*ASUM-HP(I))/BSUM
      TH(I)=TH(I-1)+V*(AA(I)+BB(I)*V)
 1800 HGT =TH(I)
C      CONVERGENCE TEST TO INSURE THAT SUCCESSIVE HT. ARE WITHIN 0.01 KM
      IF ((HGT-HGTA)**2 .LE. 0.0001) GO TO 2100
 2000 IQ=IQ+1
C      ITERATE ONLY FOR LAST LAMINATION (M=I),10 ITERATIONS MAX. NO.
      IF(IQ.LE.9) GO TO 3000
      WRITE(6,2) I,HGT,HGTA
      GO TO 2100
 3000 M=I
      CALL FIELDG(DLAT,DLONG,HGT,DATE,10,LQ,BN,BE,BV,BT)
      DIFF(IQ)=FH(I)-BT*2.8000E-05
      IF (IQ .GT. 1) GO TO 2200
      FH(I)=BT*2.8000E-05
      TEST=GIRO
      IF(FH(I).GT.GIRO) TEST=FH(I)
      GO TO 200
 2200 FH(I)=GIRO-DIFF(1)*(FH(I)-GIRO)/(DIFF(IQ)-DIFF(1))
      IF (FH(I) .GT. FH(I-1)) GO TO 2300
      FH(I)=FH(I-1) + 0.001
      GO TO 200
 2300 IF (FH(I) .LT. TEST) GO TO 200
      FH(I)=TEST-0.0001
      GO TO 200
 2100 AAI=AAI+2.0*BB(I)*V
      DEN(I)=(F(I)*(F(I)-FH(I)))*12400.0
      IF(DEN(I).GT.(DEN(I-1)*1.0005)) GO TO 2700
      J=I-1
      IF(IFLAG.NE.3) WRITE(6,1)I,DEN(I),J,DEN(J)
      GO TO 2500
 2700 M=MM
      RETURN
C     USE RETURN 1 IF BAD DATA FOUND.
 1400 WRITE(6,4) J,I,X,SX,A
 2500 M=MM
      IF(IFXS.EQ.0) WRITE(6,5)
      RETURN 1
    1 FORMAT(' DEN(',I2,') = ',1PD10.3,' LESS THAN DEN(',I2,') = ',
     F  1PD10.3/' DENSITIES NOT MONOTONICALLY INCREASING. PROCESSING OF
     FIONOGRAM TERMINATED.')
    2 FORMAT(' ITERATION DID NOT CONVERGE FOR I =',I3,'. LAST TWO HEIGHT
     FS:',2F8.2)
    3 FORMAT(' X',I1,'/(1-YR) =',1PD11.3,' > 1. (LAMINATION',I3,' POINT'
     F  ,I3,'). PROCESSING TERMINATED.')
    4 FORMAT(' FOR LAMINATION',I3,'  POINT',I3,'  X =',1PD11.3,'   1-Y**
     F2/(1-X) =',D11.3/' RATIO =',D11.3,'>1. CALCULATIONS TERMINATED.')
    5 FORMAT(' PROCESSING MAY BE POSSIBLE USING FXSAT SUBROUTINE (ICODE(
     F1)=2).')
      END
