      PROGRAM THUNDER
      COMMON RE,RO,ZM,SIGM,PI,RAD,SRPI,DPHI,JMX,DTHETA,IMX,ABC1,SR2PI,
     1ABC2,THET(73),THETA(73),PHII(37),PHI(37),DTHET,DPH,OROG(37,73),
     2FS(37,73),Z1,SIGSL,SIGO,ALP,SIG,DRL,SIG1,AC1(37,37),BC1(37,37),
     3AS1(37,37),BS1(37,37),AC2(37,37),BC2(37,37),AS2(37,37),BS2(37,37),
     4GL(4000),POT(37,73),
     5GCNK(37,37),GSNK(37,37),TES(26011),
     6  GOO,F(37,73),
     7  LMX,ZSIG(37,73),
     8  UBC(37,73),USNK(37,37),UCNK(37,37),UOO,
     9  AC3(37,37),AS3(37,37),BC3(37,37),BS3(37,37),SBQ,DPHINF,
     a  mtns(72,36), xmtns(73,37), cmtns(73,37) 
      COMMON/PHIFT/ PHINF
      COMMON/POSNEG/CURPOS,CURNEG
      COMMON / ANGLE / ANGX,ANGY
      COMMON /CURT/ TCUR,SUMCUP,SUMCDN
      COMMON/XYXP/ITILT
      COMMON/TRANS/ DSIGZ
      COMMON/CLDD/CLOUDL(37,73),CLOUDH(37,73)
      COMMON/CLDDD/DRR(37,73),CBB,ICLDS,ACLR,CSH
      DIMENSION CP(37),P(37),W(97)
      COMMON/GGGG/FSLS(37,73),SIGSLO
      DIMENSION UBCC(37,73)
      DIMENSION YPLT(73),ZPLT(73)
      COMMON/GRDDCC/GRDCUR
      DIMENSION UBCD(37,73)
      DIMENSION POTO(37,73),POTP(37,73),POTT(37,73)
      DIMENSION CCKI(37,73),SSKI(37,73),GZIM(37),GZIP(37),AQ1(37),AQ2(37
     .),SPOT(37,5),S(5,37),POTU(37,73)
      DIMENSION GZII(37), GZIMM(37), GZIPP(37)
      DIMENSION AQ3(37), AQ4(37)
      DIMENSION XPLT(73), WRKK(5512)
	dimension scld(37,73)

c size of main array = 37 X 73
	data ArrSz /2701/	
	data IXSZ /73/
	data IYSZ /37/
      	character*50 title, xlabel, ylabel
      character*24 histvol
      common/input/ histvol,mtime(3),dsigo,iubcd,ilwrg,idelt

      ALPNK(N,K)=SQRT(FLOAT((N-K+1)*(N+K)))
      INDEX(N,K,J,JM)=((N-1)*N/2+K-1)*JM+J
C;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
C Variable definitions (defined from getinp.f):
C_____________________________________________________________________________
C     Z1 IS THE HEIGHT OF THE SIGMA 1 SURFACE
C     SIGSL IS THE SEA LEVEL CONDUCTIVITY
C     SIGO IS THE SEA LEVEL CONDUCTIVITY THAT DEFINES THE UPPER
C            EXPONENTIAL CONDUCTIVITY IN REGIONS 1 AND 2
C     ZM IS THE HEIGHT OF THE MAGNETOSPHERE IN METERS
C     ALP DEFINES THE EXPONENTIAL FALL OFF OF CONDUCTIVITY IN
C            REGIONS 1 AND 2
C     DSIGO GIVES THE LATITUDINAL VARIATION OF CONDUCTIVITY ALONG
C            THE SIGMA1 SURFACE
C     IUBCD IS THE PARAMETER THAT SPECIFIES WHETHER THE IONOSPHERIC
C            POTENTIAL IS COUPLED INTO MODEL CALCULATIONS
C            IUBCD = 0 FOR NO IONOSPHERIC POTENTIAL
C            IUBCD = 1 INCLUDES COUPLING WITH IONOSPHERIC POTENTIAL
C     ITILT IS THE PARAMETER THAT GOVERNS GEOMAGNETIC AND GEOGRAPHIC
C            COORDINATES
C            ITILT = 0 GEOGRAPHIC POLE IS COINCIDENT WITH
C                      GEOMAGNETIC POLE
C            ITILT = 1 GEOGRAPHIC AND GEOMAGNETIC POLES ARE
C                      SEPARATED
C     IDELT IS THE PARAMETER CONTROLLING THE PLOTS BELOW THE
C            SIGMA1 SURFACE
C            IDELT = 0 PLOT FULL SOLUTION WITH NEGATIVE REGIONS
C            IDELT = 1 PLOT ONLY POSITIVE REGIONS GIVING SURFACE DETAILS
C     ICLDS IS THE PARAMETER CONTROLLING THE PRESENCE OF CLOUDS
C            ICLDS = 0 NO CLOUDS IN THE CALCULATION
C            ICLDS = 1 CLOUDS INCLUDED GCM OUTPUT
C     ILWRG IS THE PARAMETER CONTROLLING THE SOLUTION TECHNIQUE IN
C            REGION 0
C            ILWRG = 0 ANALYTIC SOLUTION USING EXPONENTIAL FORM
C            ILWRG = 1 TWO POINT BOUNDARY VALUE NUMERICAL SOLUTION
C
      call getinp
c
c     Z1=8.E+3
      Z11=Z1
c     SIGSL=1.14E-14
c     SIGO=1.14E-13
c     ZM=105.E+3
c     ALP=1./12.E+3
c     DSIGO=0.4
      SIG1=SIGO*EXP(2.*ALP*Z1)
      ZHTO=12.E+3
c     IUBCD=1
      RE=6371.E+3
c     ITILT=1
      RO=RE
      PI=3.1415926
      CB=2.E+3
      CBB=CB
      PPI=PI
      ACLR=1.
      CSH=2.5E+3
      ANGX=315.
      ANGY=45.
      RAD=57.295
      SRPI=SQRT(PI)
      DPHI=5.
      JMX=37
      JMX1=JMX-1
      DTHETA=5.
c     IDELT=1
      IMX=73
      ISM=2
      ISMC=37
      ISM=3
      ISMC=37*38/2
      ABC1=1./SRPI
      SR2PI=SQRT(2.*PI)
      ABC2=1./SR2PI
      DTHET=DTHETA/RAD
      DPH=DPHI/RAD
      DO 1 I=1,IMX
      THET(I)=0.+(I-1)*DTHETA
      THETA(I)=THET(I)/RAD
  1   CONTINUE
      DO 2 J=1,JMX
      PHII(J)=0.+(J-1)*DPHI
      PHI(J)=PHII(J)/RAD
  2   CONTINUE

c Initialize access to NCAR Graphics package.
	call opngks

c Select two-point BC numerical solution for low level potential field.
c     ILWRG= 0

c No clouds right now.
c     ICLDS= 0

      IF(ICLDS.GT.0) then
	open (unit=7, file= "cloud.data")
	READ (7,1000) CLOUDL
        READ (7,1000) CLOUDH
 1000   FORMAT(16F5.2)
	close (unit= 7)
        CALL EZCNTR(CLOUDL,37,73)
        CALL EZSRFC(CLOUDL,37,73,ANGX,ANGY,WRKK)
        CALL EZCNTR(CLOUDH,37,73)
        CALL EZSRFC(CLOUDH,37,73,ANGX,ANGY,WRKK)
      endif

      DO 20 N=1,JMX
        DO 20 K=1,N
          NN=N-1
          KK=K-1
c	  Ix= jmx/(jmx - n)
          CALL LFK(NN,KK,CP,W)
          CALL LFP(NN,KK,JMX,CP,P)
          ABC=ABC1
          IF(KK.EQ.0) ABC=ABC2
          DO 20 J=1,JMX
            NKJ=INDEX(N,K,J,JMX)
            TES(NKJ)=P(J)*ABC
  20  CONTINUE

      CALL STORMS

      DO 816 K=1,JMX
      DO 816 I=1,IMX
        CCKI(K,I)=COS((K-1)*THETA(I))
        SSKI(K,I)=SIN((K-1)*THETA(I))
  816 CONTINUE

      DO 817 N=1,JMX
        NN=N-1
        GZI=SQRT(1.+NN*(NN+1)/(ALP*ALP*RO*RO) )
        GZII(N)=GZI
        GZIM(N)=-0.5*(GZI+1.)
        GZIP(N)=0.5*(GZI-1.)
        GZIMM(N)=-0.5*(GZI+3.)
        GZIPP(N)=0.5*(GZI-3.)
  817 CONTINUE

c*** POT is a temporary array used for storage of orographic data.
	open(unit= 3, file= "orog.data")
      	READ (3,2001, END= 700) ((POT(J,I),I=1,73),J=1,37)
 2001 	FORMAT(18F4.0)
  700   continue
	close(unit= 3)

      DO 35 I=1,IMX
      DO 35 J=1,JMX
        JJ=JMX-J+1
        FS(J,I)=POT(JJ,I)
  35  CONTINUE


      DO 350 I=1,IMX
      DO 350 J=1,JMX
        IF(ITILT.EQ.1) GO TO 4334
        OROG(J,I)=FS(J,I)
        GO TO 4335
 4334   CONTINUE
c;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
C  Geomagnetic --> geographics coordinate xform.
c;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
        PL=PHII(J)
        TL=THET(I)
        CALL GMTOG(1,PL,TL,1,PLL,TLL)
        OROG(J,I)=OROGRY(PLL,TLL)
 4335   CONTINUE
  350 CONTINUE

      SUD=RE*RE*DTHETA*DPHI/(RAD*RAD)
      SUMA=0.
      SUMB=0.
      SUMC=0.
      SUMD=0.
      DO 2002 J=1,JMX
      DO 2002 I=1,IMX
C     DSISL=0.15
C     FSLS(J,I)=1.-DSISL/2.*(1.+COS(2.*(PHII(J)-30.)/60.*PI/2.))
C     IF(PHII(J).LT.30..OR.PHII(J).GT.150.) FSLS(J,I)=1.-DSISL
      FS(J,I)=1.+DSIGO/2.*(1.+COS(2.*(PHII(J)-30.)/60.*PI/2.))
      DDSIGO=0.
      IF(PHII(J).LT.30..OR.PHII(J).GT.150.) FS(J,I)=1.+DSIGO+DDSIGO
C     SIGSL=SIGSLO*FSLS(J,I)
      CONS1=ALOG(SIG1/SIGSL)
      CONS2=ALOG(SIG1/(SIGO*FS(J,I)))
      BET=CONS2/(CONS1*2.*ALP)
      Z1S=CONS2/(2.*ALP)
      SIGG=SIGSL*EXP(CONS1*OROG(J,I)/Z1S)
      DR=BET*(1./SIGG-1./SIG1)
      IF(ICLDS.EQ.0) GO TO 714
      CL=CLOUDH(J,I)+CLOUDL(J,I)
      CLO=CL
      CLO=AMIN1(1.0,CLO)
      CLO=AMAX1(CLO,0.)
      ORO=OROG(J,I)
      DRR(J,I)=DRLL(Z1S,ORO,SIG1,SIGSL,CLO,CB,PPI)
      DR=DRR(J,I)
  714 CONTINUE
      F(J,I)=1./(1.+2.*ALP*SIG1*DR)
      DRTOP=(EXP(-2.*ALP*Z1S)-EXP(-2.*ALP*ZM))/(SIGO*FS(J,I)*2.*ALP)
      Q=(90.-PHII(J))/RAD
      SUMA=SUMA+(1./(DR+DRTOP))*SUD*COS(Q)
      SUMC=SUMC+SUD*COS(Q)
      SUMD=SUMD+(1./(DR+DRTOP))
 2002 CONTINUE

c;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
c  TIE-GCM POTENTIAL IS READ HERE INTO UBC ARRAY
c   "fipd_R()" accesses the TIGCM potential data at runtime.
c;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
      call fipd_R(histvol,mtime,UBC)

      SUMB=4.*PI*RE*RE
      SUMA=SUMB/SUMA
      SBQ=SUMA/SUMB
      SBQQ=1./(SUMD*SUMB)

      WRITE(6,3114) SBQ,SUMA,SUMB,SUMC,SUMD,SBQQ
 3114 FORMAT(2X,19HGLOBAL RESISTANCE= ,E15.4,2X,5E15.4)
	print *,"plotting orography"
       	title= "Orography (M)"
	xlabel= "GM LONGITUDE"
	ylabel= "GM LATITUDE"
        call mbcntr (OROG, title, xlabel, ylabel, IXSZ, IYSZ)
	call scaleit (0.02, IXSZ, IYSZ, orog, scld, Zmin, Zmax)
	call mbsrfc(scld, title, xlabel, ylabel, Zmin, Zmax, IXSZ, IYSZ)

      title= "F Parameter"
      call mbcntr (F, title, xlabel, ylabel, IXSZ, IYSZ)

	title= "TIE-GCM Ionospheric Potential (kV)"
	xlabel= "Geographic Longitude"
	ylabel= "Geographic Latitude"

	factor= 1.0e-3
       	call scaleit(factor, IXSZ, IYSZ, UBC, SCLD, Zmin, Zmax)
	call mbcntr (SCLD, title, xlabel, ylabel, IXSZ, IYSZ)
	
c;;;;;;;;;;;;;;;;;;; TEMPORARY termination point :::::::::::::::::::::::
c	call clsgks
c	stop 0

	  title= "DRR parameter"
	  xlabel= "GM LONGITUDE"
	  ylabel= "GM LATITUDE"
	  call mbcntr(DRR, title, xlabel, ylabel, IXSZ, IYSZ)
C         call scaleit (1.e-3, IXSZ, IYSZ, DRR, IXSZ, IYSZ, scld, Zmin, Zmax)
C         call mbsrfc (scld, title, xlabel, ylabel, Zmin, Zmax, IXSZ, IYSZ)
C
C   COMPUTE COEFFICIENTS
C
      CALL SPHCOF(JMX,F,GCNK,GSNK)
      CALL SPHCOF(JMX,UBC,UCNK,USNK)

      DO 402 N=1,JMX
      DO 402 K=1,N
        ABC=ABC1
        IF(K.EQ.1) ABC=ABC2
        GCNK(K,N)=GCNK(K,N)/ABC
        GSNK(K,N)=GSNK(K,N)/ABC
        UCNK(K,N)=UCNK(K,N)/ABC
        USNK(K,N)=USNK(K,N)/ABC
        IF(K.EQ.1.OR.K.EQ.JMX) GCNK(K,N)=0.5*GCNK(K,N)
        IF(K.EQ.1.OR.K.EQ.JMX) UCNK(K,N)=UCNK(K,N)*0.5
  402 CONTINUE

      GOO=GCNK(1,1)*ABC2/SQRT(2.)
      UOO=UCNK(1,1)*ABC2/SQRT(2.)

      DO 2118 J=1,JMX
      DO 2118 I=1,IMX
        F(J,I)= 0.0
        UBC(J,I)= 0.0

       DO 2118 N=1,JMX
        NN=N-1
       DO 2118 K=1,N
        KK=K-1
        NKJ=INDEX(N,K,J,JMX)
        TS=TES(NKJ)
        F(J,I)= F(J,I) + GCNK(K,N)*TS*CCKI(K,I) + 
     &       GSNK(K,N)*TS*SSKI(K,I)
         UBC(J,I)= UBC(J,I) + UCNK(K,N)*TS*CCKI(K,I) + 
     &       USNK(K,N)*TS*SSKI(K,I)
 2118 CONTINUE

      title= "Reconstructed F Parameter"
      call mbcntr(F, title, xlabel, ylabel, IXSZ, IYSZ)

	title= "Reconstructed Ionospheric Potential  (kV)"
	factor= 1.e-3
      	call scaleit(factor, IXSZ, IYSZ, UBC, SCLD, Zmin, Zmax)
      	call mbcntr (SCLD, title, xlabel, ylabel, IXSZ, IYSZ)

      IF(IUBCD.EQ.1) GO TO 3388
      DO 784 J=1,JMX
      DO 784 I=1,IMX
      PF=PHI(J)
      TF=THETA(I)
      UBC(J,I)=0.
  784 CONTINUE
      DO 7844 N=1,JMX
      DO 7844 K=1,JMX
      UCNK(K,N)=0.
      USNK(K,N)=0.
 7844 CONTINUE
 3388 CONTINUE
      PX=90.
      PY=180.
      CALL COND(ZM,PX,PY)
      SIGM=SIG
      CALL COND(ZHTO,PX,PY)
      SIGMS=SIG
      DO 40 N=1,JMX
      DO 40 K=1,N
      CALL GNK(N,K)
  40  CONTINUE
      ZHT=Z1
      CALL COND(ZHT,PX,PY)
      SIGMA=SIG
      DO 8188 N=1,JMX
      AQ1(N)=SIGMA**GZIM(N)
      AQ2(N)=SIGMA**GZIP(N)
      AQ3(N)=((-1.-GZII(N))/2.)*SIGMA**GZIMM(N)
      AQ4(N)=((-1.+GZII(N))/2.)*SIGMA**GZIPP(N)
 8188 CONTINUE
      DO 3131 I=1,IMX
      DO 3131 J=1,JMX
      POTU(J,I)=0.
 3131 CONTINUE
      DO 3222 J=1,JMX
      DO 3111 I=1,IMX
      DO 3333 N=1,JMX
      SPOT(N,1)=0.
      DO 3444 K=1,N
      AC=AC3(N,K)
      BC=BC3(N,K)
      AS=AS3(N,K)
      BS=BS3(N,K)
      NKJ=INDEX(N,K,J,JMX)
      TS=TES(NKJ)
      HW=(AC*AQ1(N) + BC*AQ2(N))*TS*CCKI(K,I) + 
     &   (AS*AQ1(N)+BS*AQ2(N))*TS*SSKI(K,I)
      POTU(J,I)=POTU(J,I)+HW
 3444 CONTINUE
 3333 CONTINUE
 3111 CONTINUE
 3222 CONTINUE
	title= "POTENTIAL Z1 SURFACE, 8KM"
	xlabel= "GM LONGITUDE"
	ylabel= "GM LATITUDE"
      	call mbcntr (POTU, title, xlabel, ylabel, IXSZ, IYSZ)

      DO 1111 J=1,JMX
       DO 1111 I=1,IMX
        ZSIG(J,I)= ALOG(SIGMA/(SIGO*FS(J,I)))/(2.*ALP)
        POT(J,I)= POTU(J,I)
 1111 CONTINUE

      SUM=0.

      DO 7 J=1,JMX
       PL= PHII(J)
       WRITE(6,892) PHII(J)
       DO 77 I=1,IMX
        OR= OROG(J,I)
        TL= THET(I)
        CALL COND(OR,PL,TL)
        Z1= ZSIG(J,I)
C       SIGSL= SIGSLO*FSLS(J,I)
        TOBET= 2.*ALP*(ALOG(SIG1/SIGSL))/(ALOG(SIG1/(SIGO*FS(J,I))))
        POTU(J,I)= POTU(J,I)*TOBET*EXP(-TOBET*OROG(J,I))/(EXP(-TOBET*
     &    OROG(J,I))-EXP(-TOBET*Z1))
         Q= (90.0-PHII(J))/RAD
         UBC(J,I)= SIG*POTU(J,I)
         SUM= SUM + UBC(J,I)*SUD*COS(Q)
  77   CONTINUE
       WRITE(6,100) (POTU(J,I),I=1,IMX)
       WRITE(6,100) (UBC(J,I),I=1,IMX)
   7  CONTINUE

	  title= "ELECTRIC FIELD ON Z1 SURFACE @ 8KM"
	  xlabel= " GM Longitude"
	  ylabel= " GM Latitude"
	  call mbcntr (POTU, title, xlabel, ylabel, IXSZ, IYSZ)

	  title= "CURRENT ON Z1 SURFACE @ 8KM"
	  xlabel= " GM LONGITUDE"
	  ylabel= " GM LATITUDE"
	  call mbcntr (UBC, title, xlabel, ylabel, IXSZ, IYSZ)


      DPHINF=-SUM*SBQ
      PHINFSV=PHINF
      PHINF=0.
      if (ILWRG.GT.0) then
	CALL LOWLAY
        DPHINF=-GRDCUR*SBQ
      endif
      PHINF=PHINFSV
      DPHINFS=DPHINF
      WRITE(6,8) DPHINF,SUM,SBQ
   8  FORMAT(2X,*DPHINF=*,E12.4,2X,*TCURENT=*,E12.4,2X,*T RESISTANCE=*,
     &   E12.4)
      COO=1./SQRT(4.*PI)
      BC3(1,1)=BC3(1,1)+DPHINF/COO
      Z1=Z11

C Start the ZHT parameter loop (# 30)
      DO 30 IZ=1,5
      ZHT= (25.0 + (IZ-1)*25.0)*1.e+3
      IF(IZ.EQ.4) ZHT=ZM
      IF(IZ.EQ.5) ZHT=Z1
      CALL COND(ZHT,PX,PY)
      SIGMA=SIG

      DO 818 N=1,JMX
       AQ1(N)=SIGMA**GZIM(N)
       AQ2(N)=SIGMA**GZIP(N)
       AQ3(N)=((-1.-GZII(N))/2.)*SIGMA**GZIMM(N)
       AQ4(N)=((-1.+GZII(N))/2.)*SIGMA**GZIPP(N)
  818 CONTINUE

      DO 31 I=1,IMX
      DO 31 J=1,JMX
      POT(J,I)=0.
      F(J,I)=0.
      POTT(J,I)=0.
      POTP(J,I)=0.
      POTU(J,I)=0.
  31  CONTINUE

      DO 322 J=1,JMX
      DO 311 I=1,IMX
      DO 33 N=1,JMX
      NN=N-1
      SPOT(N,1)=0.
      SPOT(N,2)=0.
      SPOT(N,3)=0.
      SPOT(N,4)=0.

      DO 34 K=1,N
      KK=K-1
      IF(SIGMA.GT.SIGMS) GO TO 444
      AC=AC1(N,K)
      BC=BC1(N,K)
      AS=AS1(N,K)
      BS=BS1(N,K)
      GO TO 445
  444 CONTINUE
      AC=AC2(N,K)
      BC=BC2(N,K)
      AS=AS2(N,K)
      BS=BS2(N,K)
  445 CONTINUE
      NKJ=INDEX(N,K,J,JMX)
      TS=TES(NKJ)
      HH=(AC*AQ1(N)+BC*AQ2(N))*TS*CCKI(K,I)+(AS*AQ1(N)+BS*AQ2(N))*TS*
     1  SSKI(K,I)
      HI=(AC*AQ3(N)+BC*AQ4(N))*TS*CCKI(K,I)+(AS*AQ3(N)+BS*AQ4(N))*TS*
     1  SSKI(K,I)
      HP=(AC3(N,K)*AQ1(N)+BC3(N,K)*AQ2(N))*TS*CCKI(K,I)+(AS3(N,K)*AQ1(N)
     1   +BS3(N,K)*AQ2(N))*TS*SSKI(K,I)
      POTU(J,I)=POTU(J,I)+HP

      IF(N.EQ.1) GO TO 58
      IF(K.EQ.1) GO TO 59
      IF(K.EQ.N) GO TO 60
      GO TO 61
  58  CONTINUE
      DTS=0.
      GO TO 62
  59  CONTINUE
      NKJ=INDEX(N,2,J,JMX)
      TS1=TES(NKJ)
      DTS=SQRT(FLOAT(NN*(NN+1)))*TS1
      GO TO 62
  60  CONTINUE
      NKJ=INDEX(N,NN,J,JMX)
      TS2=TES(NKJ)
      T2=ALPNK(NN,KK)
      DTS=0.5*T2*TS2
      GO TO 62
  61  CONTINUE
      KP1=K+1
      T3=ALPNK(NN,KK)
      T4=ALPNK(NN,K)
      NKJ=INDEX(N,KK,J,JMX)
      TS3=TES(NKJ)
      NKJ=INDEX(N,KP1,J,JMX)
      TS4=TES(NKJ)
      DTS=0.5*T3*TS3-0.5*T4*TS4
  62  CONTINUE
      HK=(AC*AQ1(N)+BC*AQ2(N))*DTS*CCKI(K,I)+(AS*AQ1(N)+BS*AQ2(N))*DTS
     .   *SSKI(K,I)
      HK=HK/RE
      HJ=(AC*AQ1(N)+BC*AQ2(N))*TS*(-KK*SSKI(K,I))+(AS*AQ1(N)+BS*AQ2(N))
     .   *TS*(KK*CCKI(K,I))
      IF(J.EQ.1.OR.J.EQ.JMX) GO TO 447
      HJ=HJ/(RE*COS(PI/2.-PHI(J)))
      GO TO 448
  447 HJ=0.
  448 CONTINUE
      HI=HI*DSIGZ
      SPOT(N,1)=SPOT(N,1)+HH
      SPOT(N,2)=SPOT(N,2)+HI
      SPOT(N,3)=SPOT(N,3)+HJ
      SPOT(N,4)=SPOT(N,4)+HK
  34  CONTINUE
  33  CONTINUE
      DO 8181 IFZ=1,4
C   SUCCESSIVE CESARO SUMMATIONS START HERE
C   FOR ISM=2 THIS CODE SHOULD CORRESPOND TO ONE ITERATION OF CESARO
C   SUMMATION
C   STORE PARTIAL SUMS OF THE ORIGINAL SERIES
      S(1,1)=SPOT(1,IFZ)
      DO 822 K=2,JMX
      S(1,K)=S(1,K-1)+SPOT(K,IFZ)
  822 CONTINUE
C   S(1,N) WOULD BE THE VALUES WITHOUT CESARO SUMMATION
C   OUTER CESARO ITERATION LOOP
      DO 819 IS=2,ISM
      S(IS,1)=S(IS-1,1)
C   COMPUTE PARTIAL SUMS OF PREVIOS SEQUENCE BUT DO NOT DIVIDE YET
      DO 820 K=2,JMX
      S(IS,K)=S(IS,K-1)+S(IS-1,K)
  820 CONTINUE
  819 CONTINUE
      IF(IFZ.EQ.1) POT(J,I)=S(ISM,JMX)/ISMC
      IF(IFZ.EQ.2) F(J,I)=S(ISM,JMX)/ISMC
      IF(IFZ.EQ.3) POTT(J,I)=-S(ISM,JMX)*1.E+3/ISMC
      IF(IFZ.EQ.4) POTP(J,I)=-S(ISM,JMX)*1.E+3/ISMC
 8181 CONTINUE
  311 CONTINUE
  322 CONTINUE
      IF(IZ.EQ.5) PHINF=0.
      IF(IZ.EQ.5) DPHINF=0.
      DO 411 J=1,JMX
      DO 411 I=1,IMX
      POT(J,I)=POT(J,I)+POTU(J,I)
      POT(J,I)=POT(J,I)-PHINF-DPHINF
  411 CONTINUE
      IF(IZ.EQ.2) GO TO 1212
      IF(IZ.EQ.3) GO TO 1212
      DO 43 J=1,JMX
      WRITE(6,101) PHII(J),ZHT,SIGMA,PHINF
  101 FORMAT(1X,4E20.4)
  100 FORMAT(1X,12E11.3)
      WRITE(6,100) (POT(J,I),I=1,IMX)
  43  CONTINUE
 1212 CONTINUE

	write (title, '("POTENTIAL ON ZHT SURFACE @",f5.1," KM")') 
     &   ZHT*1.e-3 
	xlabel= " GM LONGITUDE"
	ylabel= " GM LATITUDE"
        call mbcntr (POT, title, xlabel, ylabel, IXSZ, IYSZ)

	call scaleit (1.0, IXSZ, IYSZ, POT, scld, Zmin, Zmax)
  	call mbsrfc(scld, title, xlabel, ylabel, Zmin, Zmax, IXSZ, IYSZ)

	write (title, '("ELECTRIC FIELD ON ZHT SURFACE @",f5.1," KM")')
     &   ZHT*1.e-3
	xlabel= " GM LONGITUDE"
	ylabel= " GM LATITUDE"
        call mbcntr (F, title, xlabel, ylabel, IXSZ, IYSZ)

	write (title, '("CURRENT ON ZHT SURFACE @",f5.1," KM")') ZHT*1.e-3
	xlabel= " GM LONGITUDE"
	ylabel= " GM LATITUDE"
        call mbcntr (POTT, title, xlabel, ylabel, IXSZ, IYSZ)

       write (title, 
     &  '("HORIZONTAL ELECTRIC FIELD ZHT SURFACE @",f5.1," KM")')
     &  ZHT*1.e-3
	xlabel= " GM LONGITUDE"
	ylabel= " GM LATITUDE"
        call mbcntr (POTP, title, xlabel, ylabel, IXSZ, IYSZ)
	call mbvec (POTP, POTT, title, xlabel, ylabel)

      DO 111 J=1,JMX
       DO 111 I=1,IMX
        ZSIG(J,I)= ALOG(SIGMA/(SIGO*FS(J,I)))/(2.*ALP)
  111 CONTINUE

      write (title, '("HEIGHT OF CONDUCTIVITY SURFACE @",f5.1," KM")')
     &  ZHT*1.e-3
	xlabel= " GM LONGITUDE"
	ylabel= " GM LATITUDE"
        call mbcntr (ZSIG, title, xlabel, ylabel, IXSZ, IYSZ)

      IF(IZ.GT.1) GO TO 30
      DO 489 J=1,JMX
        DO 490 I=1,IMX
          XPLT(I)= POT(J,I)
          YPLT(I)= F(J,I)
  490   CONTINUE
C     CALL EZXY(THET,XPLT,IMX, "POTEN")
  489 CONTINUE
  30  CONTINUE

      WRITE(6,788)
  788 FORMAT(12X, "PHINF",10X, "TCURRENT", 7X, "CURRENT UP",5X,
     1 "CURRENT DOWN")
      WRITE(6,777) PHINFSV,TCUR,SUMCUP,SUMCDN,DPHINFS
  777 FORMAT(5X,5E15.5)
C     PHINF=0.
      IF(ILWRG.GT.0) GO TO 4762
      DO 888 IZO=1,4
      ZHT=4.E+3-(IZO-1)*1.E+3
      DO 889 J=1,JMX
      DO 889 I=1,IMX
      Z1=ZSIG(J,I)
C     SIGSL=SIGSLO*FSLS(J,I)
      TO BET=2.*ALP*(ALOG(SIG1/SIGSL))/(ALOG(SIG1/(SIGO*FS(J,I))))
      POTO(J,I)=(POT(J,I)+PHINF)*(EXP(-TOBET*OROG(J,I))-EXP(-TOBET*ZHT))
     .  /(EXP(-TOBET*OROG(J,I))-EXP(-TOBET*Z1))
      WWW=ZHT-OROG(J,I)
      IF(WWW.LE.0.) POTO(J,I)=0.
      IF(IDELT.EQ.1.AND.POTO(J,I).LT.0.) POTO(J,I)=0.
  889 CONTINUE

       write (title, '("POTENTIAL ON ZHT SURFACE @",f5.1," KM")')
     &  ZHT*1.e-3
	xlabel= " GM LONGITUDE"
	ylabel= " GM LATITUDE"
        call mbcntr (POTO, title, xlabel, ylabel, IXSZ, IYSZ)
	call scaleit (1.0, IXSZ, IYSZ, POTO, scld, Zmin, Zmax)
  	call mbsrfc (scld, title, xlabel, ylabel, Zmin, Zmax, IXSZ, IYSZ)
  888 CONTINUE

      DO 8890 J=1,JMX
       DO 890 I=1,IMX
        Z1= ZSIG(J,I)
C       SIGSL= SIGSLO*FSLS(J,I)
        TOBET= 2.0*ALP*(ALOG(SIG1/SIGSL))/(ALOG(SIG1/(SIGO*FS(J,I))))
        POTO(J,I)= (POT(J,I)+PHINF)*TOBET*EXP(-TOBET*OROG(J,I))/
     &   (EXP(-TOBET*OROG(J,I))-EXP(-TOBET*Z1))
        XPLT(I)= POTO(J,I)
        IF (IDELT.EQ.1.AND.POTO(J,I).LT.0.) POTO(J,I)= 0.0
  890  CONTINUE
       IF (IZO.GT.1) GO TO 8890
C      CALL EZXY(THET,XPLT,IMX,5HPOTEN)
 8890 CONTINUE

	title= "GROUND ELECTRIC FIELD (V/M)"
	xlabel= " GM LONGITUDE"
	ylabel= " GM LATITUDE"
        call mbcntr (POTO, title, xlabel, ylabel, IXSZ, IYSZ)
        call scaleit (1.0, IXSZ, IYSZ, POTO, scld, Zmin, Zmax)
	call mbsrfc (scld, title, xlabel, ylabel, Zmin, Zmax, IXSZ, IYSZ)


      DO 891 J=1,JMX
       WRITE(6,892) PHII(J)
  892  FORMAT(5X,E20.5)
       WRITE(6,100) (POTO(J,I),I=1,IMX)
  891 CONTINUE

      SUME=0.
      SUMH=0.
      SUMI=0.
      SUMG=0.
      DO 8911 J=1,JMX
      PL=PHII(J)
      DO 8911 I=1,IMX
      OR=OROG(J,I)
      TL=THET(I)
      CALL COND(OR,PL,TL)
      POT(J,I)=POTO(J,I)*SIG
      Q=(90.-PHII(J))/RAD
      SUMF=POT(J,I)*SUD*COS(Q)
      SUMJ=SUD*COS(Q)
      IF(SUMF.LT.0.) GO TO 4616
      SUME=SUME+SUMF
      SUMH=SUMH+SUMJ
      GO TO 4617
 4616 SUMG=SUMG+SUMF
      SUMI=SUMI+SUMJ
 4617 CONTINUE
 8911 CONTINUE
      WRITE(6,3183) SUME,SUMG,SUMH,SUMI
 3183 FORMAT(2X,*TOTAL CURRENT =*,4E20.5)

	title= "GROUND CURRENT"
	xlabel= " GM Longitude"
	ylabel= " GM Latitude"
        call mbcntr (POT, title, xlabel, ylabel, IXSZ, IYSZ)
c*     CALL EZCNTR(POT,JMX,IMX)
	call scaleit (1.e12, IXSZ, IYSZ, POT, scld, Zmin, Zmax)
	call mbsrfc (scld, title, xlabel, ylabel, Zmin, Zmax, IXSZ, IYSZ)
c*     CALL EZSRFC(POT,37,73,ANGX,ANGY,WRKK)

 4762 CONTINUE
      IF(ILWRG.GT.0) CALL LOWLAY

c Close gks and the "workstation" device.
	call clsgks

      STOP 0
      END


      SUBROUTINE COND (ZHE,PHISE,THETSE)
      COMMON RE,RO,ZM,SIGM,PI,RAD,SRPI,DPHI,JMX,DTHETA,IMX,ABC1,SR2PI,
     1ABC2,THET(73),THETA(73),PHII(37),PHI(37),DTHET,DPH,OROG(37,73),
     2FS(37,73),Z1,SIGSL,SIGO,ALP,SIG,DRL,SIG1,AC1(37,37),BC1(37,37),
     3AS1(37,37),BS1(37,37),AC2(37,37),BC2(37,37),AS2(37,37),BS2(37,37),
     4GL(4000),POT(37,73),
     5GCNK(37,37),GSNK(37,37),TES(26011)
     6  ,GOO,F(37,73)
     7  ,LMX,ZSIG(37,73)
     8   ,UBC(37,73),USNK(37,37),UCNK(37,37),UOO
     9   ,AC3(37,37),AS3(37,37),BC3(37,37),BS3(37,37),SBQ,DPHINF
      COMMON/TRANS/ DSIGZ
      COMMON/GGGG/FSLS(37,73),SIGSLO
      COMMON/TEST/IL,JL,FSM,ORO
      COMMON/CLDD/CLOUDL(37,73),CLOUDH(37,73)
      COMMON/CLDDD/DRR(37,73),CBB,ICLDS,ACLR,CSH
C     SET PARAMETERS INTERNALLY
      ZH=ZHE
      PHIS=PHISE
      THETS=THETSE
      JL=PHIS/DPHI+1
      IL=THETS/DTHETA+1
      JL=MAX0(JL,1)
      IL=MAX0(IL,1)
      JL1=JL+1
      IL1=IL+1
      IF(JL1.GT.JMX) JL1=1
      IF(IL1.GT.IMX) IL1=1
      AHW=THETS-THET(IL)
C     AVOID ROUNDOFF POSSIBILITY IN NEXT TEST
      IF(AHW.LT.0.) THETS=360.+THETS
      FS1=FS(JL,IL)+(FS(JL,IL1)-FS(JL,IL))/DTHET*(THETS/RAD-THETA(IL))
      FS2=FS(JL1,IL)+(FS(JL1,IL1)-FS(JL1,IL))/DTHET*(THETS/RAD-THETA(IL)
     .)
      FSM=FS1+(FS2-FS1)/DPH*(PHIS/RAD-PHI(JL))
      SIG=SIGO*FSM*EXP(2.*ALP*ZH)
      DSIGZ=2.*ALP*SIG
C     FSL1=FSLS(JL,IL)+(FSLS(JL,IL1)-FSLS(JL,IL))/DTHET*(THETS/RAD-THETA
C    1(IL))
C     FSL2=FSLS(JL1,IL)+(FSLS(JL1,IL1)-FSLS(JL1,IL))/DTHET*(THETS/RAD-
C    1THETA(IL))
C     FSLM=FSL1+(FSL2-FSL1)/DPH*(PHIS/RAD-PHI(JL))
C     SIGSL=SIGSLO*FSLM
      IF(ICLDS.EQ.0) GO TO 888
      CLH1=CLOUDH(JL,IL)+(CLOUDH(JL,IL1)-CLOUDH(JL,IL))/DTHET*(THETS/RAD
     .   -THETA(IL))
      CLL1=CLOUDL(JL,IL)+(CLOUDL(JL,IL1)-CLOUDL(JL,IL))/DTHET*(THETS/RAD
     .   -THETA(IL))
      CLH2=CLOUDH(JL1,IL)+(CLOUDH(JL1,IL1)-CLOUDH(JL1,IL))/DTHET*(THETS
     .   /RAD-THETA(IL))
      CLL2=CLOUDL(JL1,IL)+(CLOUDL(JL1,IL1)-CLOUDL(JL1,IL))/DTHET*(THETS
     .   /RAD-THETA(IL))
      CLH=CLH1+(CLH2-CLH1)/DPH*(PHIS/RAD-PHI(JL))
      CLL=CLL1+(CLL2-CLL1)/DPH*(PHIS/RAD-PHI(JL))
      CL=CLH+CLL
      CLO=CL
C     CLO=CL-0.2
C     IF(PHIS.LT.30..OR.PHIS.GT.150.) CLO=1.5*CL-0.2
      CLO=AMAX1(CLO,0.)
C     CLO=AMIN1(0.95,CLO)
      CLO=AMIN1(1.0,CLO)
  888 CONTINUE
      IF(ICLDS.EQ.0) CLO=0.
      ORO1=OROG(JL,IL)+(OROG(JL,IL1)-OROG(JL,IL))/DTHET*(THETS/RAD-THETA
     .(IL))
      ORO2=OROG(JL1,IL)+(OROG(JL1,IL1)-OROG(JL1,IL))/DTHET*(THETS/RAD-
     .THETA(IL))
      ORO=ORO1+(ORO2-ORO1)/DPH*(PHIS/RAD-PHI(JL))
      C1=ALOG(SIG1/SIGSL)
      C2=ALOG(SIG1/(SIGO*FSM))
      BET1=C2/(C1*2.*ALP)
      Z1S=C2/(2.*ALP)
      SI1=SIGSL*EXP(C1*ORO/Z1S)
      DRL=BET1*(1./SI1-1./SIG1)
      IF(ICLDS.EQ.0.AND.SIG.GE.SIG1) RETURN
      IF(ICLDS.EQ.0) GO TO 889
      DRL1=DRR(JL,IL)+(DRR(JL,IL1)-DRR(JL,IL))/DTHET*(THETS/RAD-THETA(IL
     .))
      DRL2=DRR(JL1,IL)+(DRR(JL1,IL1)-DRR(JL1,IL))/DTHET*(THETS/RAD-
     .THETA(IL))
      DRL=DRL1+(DRL2-DRL1)/DPH*(PHIS/RAD-PHI(JL))
  889 CONTINUE
      IF(SIG.GE.SIG1) RETURN
C     FC=CLO*ACLR*SIN(PI*(ZH-CBB)/CSH)
      FC=0.9*CLO
      IF(ZH.LE.CBB) FC=0.
      IF(ZH.GE.CBB+CSH) FC=0.
      CONS=ALOG(SIG1/SIGSL)
      Z1S=ALOG(SIG1/(SIGO*FSM))/(2.*ALP)
      SIG=SIGSL*EXP(CONS*ZH/Z1S)
      SIGA=SIG
      SIG=(1.-FC)*SIG
C     DFC=ACLR*CLO*COS(PI*(ZH-CBB)/CSH)*PI/CSH
      DFC=0.
      IF(ZH.LE.CBB) DFC=0.
      IF(ZH.GE.CBB+CSH) DFC=0.
      DSIGZ=CONS/Z1S*SIGA-DFC*SIGA
      RETURN
      END
      FUNCTION DRLL(Z1S,ORO,SIG1,SIGSL,CLO,CB,PPI)
      COMMON/CLDDD/DRR(37,73),CBB,ICLDS,ACLR,CSH
      NPTS=50
      DINT=(Z1S-ORO)/NPTS
      PI=PPI
C     FC=CLO*ACLR*SIN(PI*(ORO-CBB)/CSH)
      FC=0.9*CLO
      IF(ORO.LE.CBB) FC=0.
      IF(ORO.GE.CBB+CSH) FC=0.
      CONS=ALOG(SIG1/SIGSL)
      SG1=SIGSL*EXP(CONS*ORO/Z1S)*(1.-FC)
      SUM=0.
      ZHT=ORO
      DO 1 I=2,NPTS
      ZHT=ZHT+DINT
C     FC=CLO*ACLR*SIN(PI*(ZHT-CBB)/CSH)
      FC=0.9*CLO
      IF(ZHT.LE.CBB) FC=0.
      IF(ZHT.GE.CBB+CSH) FC=0.
      SG2=SIGSL*EXP(CONS*ZHT/Z1S)*(1.-FC)
      X=1./SG1+1./SG2
      SUM=SUM+X*DINT*0.5
      SG1=SG2
    1 CONTINUE
      DRLL=SUM
      RETURN
      END
      SUBROUTINE GNK(N,K)
      COMMON RE,RO,ZM,SIGM,PI,RAD,SRPI,DPHI,JMX,DTHETA,IMX,ABC1,SR2PI,
     1ABC2,THET(73),THETA(73),PHII(37),PHI(37),DTHET,DPH,OROG(37,73),
     2FS(37,73),Z1,SIGSL,SIGO,ALP,SIG,DRL,SIG1,AC1(37,37),BC1(37,37),
     3AS1(37,37),BS1(37,37),AC2(37,37),BC2(37,37),AS2(37,37),BS2(37,37),
     4GL(4000),POT(37,73),
     5GCNK(37,37),GSNK(37,37),TES(26011)
     6  ,GOO,F(37,73)
     7  ,LMX,ZSIG(37,73)
     8   ,UBC(37,73),USNK(37,37),UCNK(37,37),UOO
     9   ,AC3(37,37),AS3(37,37),BC3(37,37),BS3(37,37),SBQ,DPHINF
      COMMON/STRMS/CURL(4000),ZL(4000),PHIL(4000),THETL(4000)
      COMMON/PHIFT/ PHINF
      COMMON /CURT/ TCUR,SUMCUP,SUMCDN
      DIMENSION CP(37),W(97)
      IF(N.GT.1) GO TO 501
      AS1(N,K)=0.
      BS1(N,K)=0.
      SUM1=0.
      SUM2=0.
      SUM3=0.
      COO=1./SQRT(4.*PI)
      DO 301 L=1,LMX
      ZLL=ZL(L)
      PHILL=PHIL(L)
      THETLL=THETL(L)
      CALL COND(ZLL,PHILL,THETLL)
      QCOO=CURL(L)*COO/(SIG*RO*RO)
      SUM1=SUM1+QCOO
      GL(L)=1./(1.+2.*ALP*SIG1*DRL)
      SUM2=SUM2+QCOO*(SIG/SIG1-GL(L))
      SUM3=SUM3+QCOO*SIG
  301 CONTINUE
      AC1(N,K)=-SUM3/(2.*ALP)
      PHINF=SUM2/(2.*ALP*GOO)
      BC1(N,K)=SUM1/(2.*ALP)+PHINF/COO
      BS2(N,K)=0.
      AC2(N,K)=0.
      AS2(N,K)=0.
      BC2(N,K)=PHINF/COO
      AC3(N,K)=0.
      AS3(N,K)=0.
      BS3(N,K)=0.
      BC3(N,K)=UCNK(K,N)
      GO TO 500
  501 CONTINUE
      NN=N-1
      KK=K-1
      GZI=SQRT(1.+NN*(NN+1)/(ALP*ALP*RO*RO))
      GZIP=(1.+GZI)/2.
      GZIM=(1.-GZI)/2.
      SUM1=0.
      SUM2=0.
      SUM3=0.
      SUM4=0.
      SUM5=0.
      SUM6=0.
      SUMCUP=0.
      SUMCDN=0.
      DO 302 L=1,LMX
      ZLL=ZL(L)
      PHILL=PHIL(L)
      THETLL=THETL(L)
      CALL COND(ZLL,PHILL,THETLL)
      PHILR=PHILL/RAD
      CALL LFK(NN,KK,CP,W)
      CALL LFPT(NN,KK,PHILR,CP,PP)
      ABC=ABC1
      IF(KK.EQ.0) ABC=ABC2
      PP=PP*ABC
      QCNK=CURL(L)*PP*COS(KK*THETLL/RAD)/(SIG*RO*RO)
      QSNK=CURL(L)*PP*SIN(KK*THETLL/RAD)/(SIG*RO*RO)
      GW1L=SIG**GZIP/GZI-GL(L)*SIG1**GZIP
      SUM1=SUM1+QCNK*GW1L
      SUM2=SUM2+QSNK*GW1L
      CW1=SIG**GZIP
      CW2=GL(L)*SIG1**GZIP
      SUM3=SUM3+QCNK*(SIG**GZIM)/GZI
      SUM4=SUM4+QCNK*CW2
      SUM5=SUM5+QSNK*CW2
      SUM6=SUM6+QSNK*(SIG**GZIM)/GZI
      SUMCUP=SUMCUP+CURL(L)*(1.-SIG1*GL(L)/SIG)
      SUMCDN=SUMCDN+CURL(L)*SIG1*GL(L)/SIG
  302 CONTINUE
      AC2(N,K)=-PHINF*SIG1**GZIP*GCNK(K,N)+SUM1/(2.*ALP)
      AS2(N,K)=-PHINF*SIG1**GZIP*GSNK(K,N)+SUM2/(2.*ALP)
      NTEST=NN-KK
      NST=NTEST-2*(NTEST/2)
      IF(NST.EQ.0) GO TO 12
      GO TO 13
  12  CONTINUE
      ETAP=4.*ALP*ALP*RO*RO/(NN*(NN+1))/(SIGM**GZI)
      BC2(N,K)=ETAP*AC2(N,K)
      BS2(N,K)=ETAP*AS2(N,K)
      GO TO 14
  13  CONTINUE
      ETAM=-(1./(SIGM**GZI))
      BC2(N,K)=ETAM*AC2(N,K)
      BS2(N,K)=ETAM*AS2(N,K)
  14  CONTINUE
      AC1(N,K)=-PHINF*SIG1**GZIP*GCNK(K,N)-SUM4/(2.*ALP)
      AS1(N,K)=-PHINF*SIG1**GZIP*GSNK(K,N)-SUM5/(2.*ALP)
      BC1(N,K)=BC2(N,K)+SUM3/(2.*ALP)
      BS1(N,K)=BS2(N,K)+SUM6/(2.*ALP)
      DZIP=NN*(NN+1)/(4.*ALP*ALP*RO*RO)
      DZIPP=1.+2.*DZIP
      AC3(N,K)=-UCNK(K,N)*SIG1**DZIPP*GCNK(K,N)*SIGM**GZIM
      AS3(N,K)=-USNK(K,N)*SIG1**DZIPP*GSNK(K,N)*SIGM**GZIM
      BC3(N,K)=UCNK(K,N)*SIGM**GZIM
      BS3(N,K)=USNK(K,N)*SIGM**GZIM
  500 CONTINUE
      RETURN
      END
      SUBROUTINE STORMS
      COMMON RE,RO,ZM,SIGM,PI,RAD,SRPI,DPHI,JMX,DTHETA,IMX,ABC1,SR2PI,
     1ABC2,THET(73),THETA(73),PHII(37),PHI(37),DTHET,DPH,OROG(37,73),
     2FS(37,73),Z1,SIGSL,SIGO,ALP,SIG,DRL,SIG1,AC1(37,37),BC1(37,37),
     3AS1(37,37),BS1(37,37),AC2(37,37),BC2(37,37),AS2(37,37),BS2(37,37),
     4GL(4000),POT(37,73),
     5GCNK(37,37),GSNK(37,37),TES(26011)
     6  ,GOO,F(37,73)
     7  ,LMX,ZSIG(37,73)
     8   ,UBC(37,73),USNK(37,37),UCNK(37,37),UOO
     9   ,AC3(37,37),AS3(37,37),BC3(37,37),BS3(37,37),SBQ,DPHINF
      COMMON/STRMS/CURL(4000),ZL(4000),PHIL(4000),THETL(4000)
      COMMON/TEST/IL,JL,FSM,ORO
      COMMON /CURT/ TCUR,SUMCUP,SUMCDN
      COMMON/XYXP/ITILT
      LMX=2000
      LMXL=LMX/2
      ALAT1=90.
      DLAT1=10.
      ALONG1=20.
      DLONG1=20.
      ALAT2=80.
      DLAT2=20.
      ALONG2=110.
      DLONG2=30.
      ALAT3=65.
      DLAT3=15.
      DLAT3=10.
      ALONG3=270.
      DLONG3=30.
      DLONG3=10.
      ALAT4=120.
      DLAT4=10.
      ALONG4=320.
      DLONG4=20.
      ALAT5=45.
      DLAT5=5.
      ALONG5=7.5
      DLONG5=7.5
      IF(ITILT.EQ.0) GO TO 49
      A=ALAT1
      B=ALONG1
      CALL GTOGM(1,A,B,1,C,D)
      ALAT1=C
      ALONG1=D
      A=ALAT2
      B=ALONG2
      CALL GTOGM(1,A,B,1,C,D)
      ALAT2=C
      ALONG2=D
      A=ALAT3
      B=ALONG3
      CALL GTOGM(1,A,B,1,C,D)
      ALAT3=C
      ALONG3=D
      A=ALAT4
      B=ALONG4
      CALL GTOGM(1,A,B,1,C,D)
      ALAT4=C
      ALONG4=D
      A=ALAT5
      B=ALONG5
      CALL GTOGM(1,A,B,1,C,D)
      ALAT5=C
      ALONG5=D
  49  CONTINUE
      CUR1=4.06
      CURO=-4.06
      DZ1=2.E+3
      ZSH1=10.E+3
      DHTT=3.E+3
      ZSH2=ZSH1-DHTT
      ATN=1.0
      DO 1 L=1,300
      XP=RANF()
      YP=RANF()
      PHIL(L)=-DLAT1+2.*DLAT1*XP+ALAT1
      THETL(L)=-DLONG1+2.*DLONG1*YP+ALONG1
      IF(THETL(L).LT.0.) THETL(L)=360.+THETL(L)
      ZP=RANF()
      ZL(L)=-DZ1+2.*DZ1*ZP+ZSH1
      CURL(L)=CUR1
  1   CONTINUE
      DO 2 L=301,350
      XP=RANF()
      YP=RANF()
      PHIL(L)=-DLAT2+2.*DLAT2*XP+ALAT2
      THETL(L)=-DLONG2+2.*DLONG2*YP+ALONG2
      IF(THETL(L).LT.0.) THETL(L)=360.+THETL(L)
      ZP=RANF()
      ZL(L)=-DZ1+2.*DZ1*ZP+ZSH1
      CURL(L)=CUR1
  2   CONTINUE
      DO 3 L=351,850
      CURL(L)=CUR1
      XP=RANF()
      YP=RANF()
      ZP=RANF()
      ZL(L)=-DZ1+2.*DZ1*ZP+ZSH1
      THETL(L)=-DLONG3+2.*DLONG3*YP+ALONG3
      IF(THETL(L).LT.0.) THETL(L)=360.+THETL(L)
      PHIL(L)=-DLAT3+2.*DLAT3*XP+ALAT3
  3   CONTINUE
      DO 4 L=851,990
      CURL(L)=CUR1
      XP=RANF()
      YP=RANF()
      ZP=RANF()
      ZL(L)=-DZ1+2.*DZ1*ZP+ZSH1
      PHIL(L)=-DLAT4+2.*DLAT4*XP+ALAT4
      THETL(L)=-DLONG4+2.*DLONG4*YP+ALONG4
      IF(THETL(L).LT.0.) THETL(L)=360.+THETL(L)
   4  CONTINUE
      DO 5 L=991,1000
      CURL(L)=CUR1
      XP=RANF()
      YP=RANF()
      ZP=RANF()
      ZL(L)=-DZ1+2.*DZ1*ZP+ZSH1
      PHIL(L)=-DLAT5+2.*DLAT5*XP+ALAT5
      THETL(L)=-DLONG5+2.*DLONG5*YP+ALONG5
      IF(THETL(L).LT.0.) THETL(L)=360.+THETL(L)
   5  CONTINUE
      DO 11 L=1001,1300
      LL=L-LMXL
      PHIL(L)=PHIL(LL)
      THETL(L)=THETL(LL)
      CURL(L)=CURO*ATN
      ZL(L)=ZL(LL)-DHTT
  11  CONTINUE
      DO 22 L=1301,1350
      LL=L-LMXL
      PHIL(L)=PHIL(LL)
      THETL(L)=THETL(LL)
      CURL(L)=CURO*ATN
      ZL(L)=ZL(LL)-DHTT
  22  CONTINUE
      DO 33 L=1351,1850
      LL=L-LMXL
      PHIL(L)=PHIL(LL)
      THETL(L)=THETL(LL)
      CURL(L)=CURO*ATN
      ZL(L)=ZL(LL)-DHTT
  33  CONTINUE
      DO 44 L=1851,1990
      LL=L-LMXL
      PHIL(L)=PHIL(LL)
      THETL(L)=THETL(LL)
      CURL(L)=CURO*ATN
      ZL(L)=ZL(LL)-DHTT
  44  CONTINUE
      DO 55 L=1991,2000
      LL=L-LMXL
      PHIL(L)=PHIL(LL)
      THETL(L)=THETL(LL)
      ZL(L)=ZL(LL)-DHTT
      CURL(L)=CURO*ATN
  55  CONTINUE
      TCUR=0.
      DO 1000 L=1,LMX
      TCUR=TCUR+CURL(L)
 1000 CONTINUE
      RETURN
      END

      SUBROUTINE EFLD(CT,CL,SL,POT,UE,UM)
C CT=COS(THETA), CL=COS(PHI), AND SL=SIN(PHI),  WHERE THETA AND PHI ARE DEFINED
C IN RICHMOND (JGR 81, P. 1447, 1976).  POT IS THE ELECTROSTATIC POTENTIAL IN
C VOLTS, UE IS U(SUBSCRIPT PHI), AND UM IS U(SUBSCRIPT PERPENDICULAR) AS DEFINED
C BY RICHMOND.  THIS SUBROUTINE RUNS FASTER IF CALLS ARE ARRANGED SUCH THAT CT
C IS USUALLY THE SAME IN SUCCESSIVE CALLS.
      COMMON/COEF/Q(7),RO(6),R(6,6),A(27),SQ2,CTO
      DIMENSION P(7),FC(6),FS(6),FCPOC(6),FSPOC(6)
      DATA IFRST/1/
      IF (IFRST.EQ.1) CALL CNST
      IFRST = 0
      IF (CT.EQ.CTO) GO TO 90
      CTO = CT
      X = 1. - CT*CT
      IF (X.LT.0.) X = 0.
      ST = SQRT(X)
      DO 30 M=1,6
      FC(M) = 0.
      FS(M) = 0.
      FCPOC(M) = 0.
   30 FSPOC(M) = 0.
      FO = 0.
      FOPOC = 0.
      P(1) = CT/RO(1)
      P(2) = (CT*P(1) - RO(1))/RO(2)
      DO 50 N=3,6
   50 P(N) = (CT*P(N-1) - RO(N-1)*P(N-2))/RO(N)
      I = 0
      DO 54 N=2,6,2
      I = I + 1
      FO = FO + A(I)*P(N)
      IF (ABS(CT).LT.1.E-5) GO TO 53
      IF (    ST .LT.1.E-5) GO TO 52
      PPOC = (FLOAT(N)*P(N) - FLOAT(2*N+1)*RO(N)*P(N-1)/CT)/ST
      GO TO 54
   52 PPOC = 0.
      GO TO 54
   53 PPOC = FLOAT(N*(N+1))*P(N)
   54 FOPOC = FOPOC + A(I)*PPOC
C FROM HERE TO STATEMENT 80, P IS P(N,M)/SIN(THETA)
      P(1) = Q(1)
      X = 0.
      DO 80 M=1,6
      MP = M + 1
      IF (MP.EQ.7) GO TO 61
      DO 60 N=MP,6
      IF (N.GT.MP) X = P(N-2)
   60 P(N) = (CT*P(N-1) - R(N-1,M)*X)/R(N,M)
   61 DO 64 N=M,6,2
      I = I + 2
      FC(M) = FC(M) + A(I-1)*P(N)
      FS(M) = FS(M) + A(I)*P(N)
      IF (ABS(CT).LT.1.E-5) GO TO 63
      PPOC = FLOAT(N)*P(N)
      IF(N.GT.M) PPOC = PPOC - FLOAT(2*N+1)*R(N,M)*P(N-1)/CT
      GO TO 635
   63 PPOC = FLOAT(N*(N+1) - M*M)*P(N)
  635 FCPOC(M) = FCPOC(M) + A(I-1)*PPOC
   64 FSPOC(M) = FSPOC(M) + A(I)*PPOC
   80 P(MP) = Q(MP)*P(M)*ST
   90 CM = SQ2*CL
      SM = SQ2*SL
      POT = FO
      UE = -FOPOC
      UM = 0.
      DO 100 M=1,6
      POT = POT + ST*(FC(M)*CM + FS(M)*SM)
      UE = UE - FCPOC(M)*CM - FSPOC(M) *SM
      UM = UM + FLOAT(M)*(FC(M)*SM - FS(M)*CM)
      X = CM*CL - SM*SL
      SM = CM*SL + SM*CL
  100 CM = X
      UE = UE/346.892
      UM = UM/(346.892*SQRT(.25 + .75*CT*CT))
      RETURN
      END
      SUBROUTINE CNST
      REAL MS,NS
      COMMON/COEF/Q(7),RO(6),R(6,6),A(27),SQ2,CTO
      SQ2 = SQRT(2.E0)
      CTO = 2.
      Q(7) = 0.
      DO 5 M=1,6
      Q(M) = SQRT(1. + 1.E0/FLOAT(2*M))
      MS = M*M
      RO(M) = SQRT(MS/(4.E0*MS - 1.))
      DO 5 N=M,6
      NS = N*N
    5 R(N,M) = SQRT((NS-MS)/(4.E0*NS - 1.))
      READ 7, (A(I),I=1,27)
    7 FORMAT (8F10.0)
      RETURN
      END

      FUNCTION POTOP(COLAT,XLONG)
C
C     THIS FUNCTION RETURNS AN ANALYTIC UPPER BOUNDARY CONDITION GIVEN
C     GEOMAGNETIC CO-LATTITUDE, LONGITUDE IN RADIANS; COLAT,XLONG, 
C     RESPECTIVELY.  IT IS ASSUMED THAT COLAT LIES BETWEEN 0 AND 180 DEGREES 
C     (i.e. 0 to PIE RADIANS) AND THAT XLONG LIES BETWEEN 0 AND 360 DEGREES. 
C
      DATA KSW / 0 /
      IF (KSW.NE.0) GO TO 1
      KSW=1
C     SET INTERNAL CONSTANTS INDEPENDENT OF XLAT,XLONG ONCE ONLY
      PII=4.0*ATAN(1.0)
      TWOPII=PII+PII
      T0=15.*PII/180.
      P0=PII/9.
      P1=3.*PII/4.
      ST0=SIN(T0)
      SP0=SIN(P0)
      SP1 = SIN(P1)
    1 CONTINUE
C     SET COORDINATES INTERNALLY
      T=COLAT
      P=XLONG
C     SET EXPONENT IN ANALUTIC FORMULA
      IF(T.GT.T0) KEXP=-3
      IF(T.LE.T0 .OR. T.GE.PII-T0) KEXP=1
C     SET MULTIPLICATIVE CONSTANT
      IF(P.GE.0..AND.P.LE.PII) PSI0=-0.4E+5
      IF(P.GE.PII.AND.P.LE.PII+PII) PSI0=0.3E+5
      ST = SIN(T)
      SP=SIN(P)
      STK = PSI0*(ST/ST0)**KEXP
      IF( P.GE.0.0 .AND. P.LE.P0) POTOP   =STK*SP/SP0
      IF(P.GE.TWOPII-P0 .AND. P.LE.TWOPII)POTOP   =STK*SP/SIN(TWOPII-P0)
      IF((P.GT.P0.AND.P.LE.P1).OR.(P.GT.TWOPII-P1.AND.P.LE.TWOPII-P0))
     1POTOP = STK
      IF(P.GT.P1.AND.P.LE.PII) POTOP   =STK*SP/SP1
      IF(P.LE.TWOPII-P1.AND.P.GT.PII) POTOP   =STK*SP/SIN(TWOPII-P1)
      RETURN
      END

      FUNCTION OROGRY(P,T)
      COMMON RE,RO,ZM,SIGM,PI,RAD,SRPI,DPHI,JMX,DTHETA,IMX,ABC1,SR2PI,
     1ABC2,THET(73),THETA(73),PHII(37),PHI(37),DTHET,DPH,OROG(37,73),
     2FS(37,73),Z1,SIGSL,SIGO,ALP,SIG,DRL,SIG1,AC1(37,37),BC1(37,37),
     3AS1(37,37),BS1(37,37),AC2(37,37),BC2(37,37),AS2(37,37),BS2(37,37),
     4GL(4000),POT(37,73),
     5GCNK(37,37),GSNK(37,37),TES(26011)
     6  ,GOO,F(37,73)
     7  ,LMX,ZSIG(37,73)
     8   ,UBC(37,73),USNK(37,37),UCNK(37,37),UOO
     9   ,AC3(37,37),AS3(37,37),BC3(37,37),BS3(37,37),SBQ,DPHINF
      PP=P
      TT=T
      JL=PP/DPHI+1
      IL=TT/DTHETA+1
      JL1=JL+1
      IL1=IL+1
      IF(JL1.GT.JMX) JL1=1
      IF(IL1.GT.IMX) IL1=1
      AHW=TT-THET(IL)
      IF(AHW.LT.0.) TT=360.+TT
      FS1=FS(JL,IL)+(FS(JL,IL1)-FS(JL,IL))/DTHET*(TT/RAD-THETA(IL))
      FS2=FS(JL1,IL)+(FS(JL1,IL1)-FS(JL1,IL))/DTHET*(TT/RAD-THETA(IL))
      FSM=FS1+(FS2-FS1)/DPH*(PP/RAD-PHI(JL))
      OROGRY=FSM
      RETURN
      END
      SUBROUTINE GMTOG(I1,DTHETA,DPHI,I2,GTHETA,GPHI)
C     THIS SUBROUTINE TRANSFORMS GEOMAGNETIC TO GEOGRAPHIC COORDINATES.
C     THE CORRESPONDING THETAS AND PHIS ARE THE COLATITUDE AND LONGITUDE EAST
C     RESPECTIVELY. THE INPUT IS IN DEGREES OR RADIANS FOR I1 OF 1 OR 2 RESPECT
C     IVELY. THE CORRESPONDING CHOICE FOR I2 GOVERNS THE OUTPUT.
      CHI=11.65036*1.74532925E-2
      PD=290.3202 *1.74532925E-2
      SINCHI=SIN(CHI)
      COSCHI=COS(CHI)
      A=DTHETA
      B=DPHI
      GOTO(1,2),I1
    1 A=A*1.74532925E-2
      B=B*1.74532925E-2
  2   SINTH=SIN(A)
      COSTH=COS(A)
      SINPHI=SIN(B)
      COSPHI=COS(B)
      GTHETA=ACOS(-SINTH*COSPHI*SINCHI+COSTH*COSCHI)
      Y=SINTH*SINPHI
      X=SINTH*COSPHI*COSCHI+COSTH*SINCHI
      GPHI=ATAN2(Y,X)+6.28318531+PD
      GPHI=AMOD(GPHI,6.28318531)
      GOTO(3,4),I2
    3 GTHETA=GTHETA/1.74532925E-2
      GPHI=GPHI/1.74532925E-2
    4 CONTINUE
      RETURN
      END
      SUBROUTINE GTOGM(I1,DTHETA,DPHI,I2,GTHETA,GPHI)
      DIMENSION ACC(2)
      PHI=90.-DTHETA
      ALAM=DPHI
      RAD=1./1.74532925E-2
      PHIP=11.65036/RAD
      ALAMP=290.3202/RAD
      I=1
      DLAM=5.
      AS=SIN(PHIP)
      AT=COS(PHIP)
      PHIR=PHI/RAD
      A5=ALAM-DLAM
      IF(A5.LT.0.) A5=360.+A5
      A6=SIN(PHIR)*AT+COS(PHIR)*AS*COS(A5/RAD-ALAMP)
      A7=SQRT(1.-A6*A6)
      IF(A7.EQ.0.) A7=1.E-6
      A8=COS(PHIR)*SIN(A5/RAD-ALAMP)/A7
      ACC(1)=A8
      ALAMR=ALAM/RAD
      AB=SIN(PHIR)*AT+COS(PHIR)*AS*COS(ALAMR-ALAMP)
      IF(AB.GT.1.) AB=1.
      IF(AB.LT.-1.) AB=-1.
      GPHI=ASIN(AB)
      DQ=SQRT(1.-AB*AB)
      IF(DQ.EQ.0.) DQ=1.E-6
      AC=COS(PHIR)*SIN(ALAMR-ALAMP)/DQ
      IF(AC.GT.1.) AC=1.
      IF(AC.LT.-1.) AC=-1.
      ACC(I+1)=AC
      DERV=ACC(I+1)-ACC(I)
      GLAM=ASIN(AC)
       GPHIR=GPHI
      GPHI=GPHI*RAD
      GLAM=GLAM*RAD
      GLAN=GLAM
      PHIA=PHI
      IF(PHIA.LT.-78.3) GO TO 88
      IF(PHIA.GT.78.3.AND.AC.LE.0..AND.DERV.GT.0.) GLAM=180.-GLAM
      IF(PHIA.GT.78.3.AND.AC.GE.0..AND.DERV.GE.0.) GLAM=180.-GLAM
      IF(AC.LE.0..AND.DERV.LE.0.) GLAM=180.-GLAM
      IF(AC.GE.0..AND.DERV.LE.0.) GLAM=180.-GLAM
      IF(PHI.GT.78.3) GO TO 77
      IF(AC.GE.0..AND.DERV.GE.0.) GLAM=GLAM
      IF(AC.LE.0..AND.DERV.GT.0.) GLAM=360.+GLAM
      GO TO 77
  88  CONTINUE
      GLAM=GLAM
  77  CONTINUE
      IF(AC.LT.-.95) GO TO 3
      IF(AC.GT..95) GO TO 3
      GO TO 4
  3   B5=ALAM+DLAM
      B6=SIN(PHIR)*AT+COS(PHIR)*AS*COS(B5/RAD-ALAMP)
      B7=SQRT(1.-B6*B6)
      IF(B7.EQ.0.) B7=1.E-6
      B8=COS(PHIR)*SIN(B5/RAD-ALAMP)/B7
      ACP1=B8
      IF(ACP1.GT.0..AND.ACP1.GE.ACC(I)) GLAM=GLAN
      IF(ACP1.GT.0..AND.ACP1.LT.ACC(I)) GLAM=180.-GLAN
      IF(ACP1.LT.0..AND.ACP1.GE.ACC(I)) GLAM=360.+GLAN
      IF(GLAM.LT.0.) GLAM=360.+GLAM
      IF(ACP1.LT.0..AND.ACP1.LT.ACC(I)) GLAM=180.-GLAN
  4   CONTINUE
      IF(GLAM.LT.1.E-6) GLAM=0.
      IF(GLAM.LT.0.) GLAM=360.+GLAM
      GPHI=90.-GPHI
      IF(GPHI.LT.1.E-3) GPHI=0.
      GTHETA=GPHI
      GPHI=GLAM
      IF(GTHETA.EQ.0..OR.GTHETA.EQ.180.) GPHI=0.
      RETURN
      END
      SUBROUTINE WGHT(N,M,L,CP,WT)
C     DIMENSION CP(1),WT(1),CW(LAT)
      DIMENSION CP(1),WT(1),CW(37)
      DO 23 I=1,L
      CW(I)=0.
23    CONTINUE
      NMOD=MOD(N,2)
      MMOD=MOD(M,2)
      IF(NMOD)1,1,2
1     IF(MMOD)3,3,4
C
C     N EVEN, M EVEN
C
3     KDO=N/2+1
      CP(1)=.5*CP(1)
      DO 5 IP1=1,L,2
      I=IP1-1
      SUM=0.
      DO 6 KP1=1,KDO
      K=KP1-1
      T1=FLOAT(1-(K+K+I)**2)
      T2=FLOAT(1-(K+K-I)**2)
      IF(T1)7,6,7
7     IF(T2)8,6,8
8     SUM=SUM+CP(KP1)*(T1+T2)/(T1*T2)
6     CONTINUE
      CW(IP1)=SUM
5     CONTINUE
      CALL COST(L,CW,WT)
      WT(1)=.5*WT(1)
      WT(L)=.5*WT(L)
      RETURN
C
C     N EVEN, M ODD
C
4     KDO=N/2
      DO 9 IP1=3,L,2
      I=IP1-1
      SUM=0.
      DO 10 K=1,KDO
      T1=FLOAT(1-(K+K+I)**2)
      T2=FLOAT(1-(K+K-I)**2)
      IF(T1)11,10,11
11    IF(T2)12,10,12
12    SUM=SUM+CP(K)*(T1-T2)/(T1*T2)
10    CONTINUE
      CW(IP1)=SUM
9     CONTINUE
      CALL SINT(L,CW,WT)
      RETURN
2     IF(MMOD)13,13,14
C
C     N ODD, M EVEN
C
13    KDO=(N+1)/2
      DO 15 IP1=2,L,2
      I=IP1-1
      SUM=0.
      DO 16 K=1,KDO
      T1=FLOAT(1-(K+K-1+I)**2)
      T2=FLOAT(1-(K+K-1-I)**2)
      IF(T1)17,16,17
17    IF(T2)18,16,18
18    SUM=SUM+CP(K)*(T1+T2)/(T1*T2)
16    CONTINUE
      CW(IP1)=SUM
15    CONTINUE
      CALL COST(L,CW,WT)
      WT(1)=.5*WT(1)
      WT(L)=.5*WT(L)
      RETURN
C
C     N ODD, M ODD
C
14    KDO=(N+1)/2
      DO 19 IP1=2,L,2
      I=IP1-1
      SUM=0.
      DO 20 K=1,KDO
      T1=FLOAT(1-(K+K-1+I)**2)
      T2=FLOAT(1-(K+K-1-I)**2)
      IF(T1)21,20,21
21    IF(T2)22,20,22
22    SUM=SUM+CP(K)*(T1-T2)/(T1*T2)
20    CONTINUE
      CW(IP1)=SUM
19    CONTINUE
      CALL SINT(L,CW,WT)
      RETURN
      END
      SUBROUTINE SPHCOF(L,G,A,B)
C     COMPUTE THE SPHERICAL HARMONIC COEFFICIENTS
C     DIMENSION G(LAT,LONG),WT(LAT),A(LAT,LAT),B(LAT,LAT),CP(LAT),
C     W(2.5*LAT+4.),WB(LONG)
      DIMENSION G(37,73),WT(37),A(37,37),B(37,37),CP(37),WB(73),W(97)
      NTOT=L+L-2
      FN=FLOAT(NTOT)
      LM1=L-1
      DO 60 I=2,LM1
      DO 61 J=1,NTOT
      W(J)=G(I,J)
61    CONTINUE
      CALL PF(NTOT,W,WB)
      DO 62 J=1,NTOT
      G(I,J)=WB(J)
62    CONTINUE
60    CONTINUE
      G(1,NTOT)=2.*G(1,NTOT)
      G(L,NTOT)=2.*G(L,NTOT)
      NTM=NTOT-1
      DO 63 J=1,NTM
      G(1,J)=0.
      G(L,J)=0.
63    CONTINUE
      DO 100 NP1=1,L
      N=NP1-1
      DO 100 MP1=1,NP1
      M=MP1-1
      MA=MP1
      MB=MP1+L-1
      CALL LFK(N,M,CP,W)
      CALL WGHT(N,M,L,CP,WT)
      IF(M)50,50,51
50    SUMA=0.
      DO 52 I=1,L
      SUMA=SUMA+WT(I)*G(I,NTOT)
52    CONTINUE
      A(MP1,NP1)=SUMA
      GO TO 100
  51  IF(MP1-L) 522,53,53
53    SUMA=0.
      DO 54 I=1,L
      SUMA=SUMA+WT(I)*G(I,NTOT-1)
54    CONTINUE
      A(MP1,NP1)=SUMA
      GO TO 100
  522 SUMA=0.
      SUMB=0.
      DO 102 I=1,L
      SUMA=SUMA+WT(I)*G(I,2*MP1-3)
      SUMB=SUMB+WT(I)*G(I,2*MP1-2)
102   CONTINUE
      A(MP1,NP1)=SUMA
      B(MP1,NP1)=SUMB
100   CONTINUE
C
C     FOR M=0 SET B TO ZERO
C
      DO 101 NP1=1,L
      B(1,NP1)=0.
101   CONTINUE
C
C     FOR N=L-1 EVEN AND M ODD  SET A AND B EQUAL TO ZERO
C
      DO 104 MP1=2,L,2
      A(MP1,L)=0.
      B(MP1,L)=0.
104   CONTINUE
      B(L,L)=0.
      RETURN
      END
      SUBROUTINE INVSPH(L,G,A,B)
C
C     COMPUTE THE FUNCTION FROM ITS SPECTRAL COEFFICIENTS
C
C     DIMENSION G(LAT,LONG),A(LAT,LAT),B(LAT,LAT),CP(LAT),P(LAT),
C      W(2.5*LAT+4.),WB(LONG)
      DIMENSION G(37,73),A(37,37),B(37,37),CP(37),P(37),WB(73),W(97)
      B(L,L)=0.
      DO 70 N=1,L
      B(1,N)=0.
70    CONTINUE
      DO 80 M=2,L,2
      A(M,L)=0.
      B(M,L)=0.
80    CONTINUE
      NTOT=L+L-2
      DO 81 J=1,NTOT
      DO 81 I=1,L
      G(I,J)=0.
  81  CONTINUE
      DO 100 NP1=1,L
      N=NP1-1
      DO 100 I=1,L
      CALL LFMA(N,L,I,P,W)
      IF(NP1.EQ.1) GO TO 89
      DO 90 MP1=2,NP1
      G(I,2*MP1-3)=G(I,2*MP1-3)+A(MP1,NP1)*P(MP1)
      G(I,2*MP1-2)=G(I,2*MP1-2)+B(MP1,NP1)*P(MP1)
90    CONTINUE
89    G(I,NTOT)=G(I,NTOT)+A(1,NP1)*P(1)
100   CONTINUE
      DO 50 I=1,L
      DO 51 J=1,NTOT
      W(J)=G(I,J)
51    CONTINUE
      CALL PB(NTOT,W,WB)
      DO 52 J=1,NTOT
      G(I,J)=WB(J)
52    CONTINUE
50    CONTINUE
      RETURN
      END
      SUBROUTINE PF(N,Y,YB)
C     N MUST BE DIVISIBLE BY 2
      DIMENSION Y(1),YB(1)
      PI=3.14159265358979
      TSN=2./N
      NM2=N-2
      DT=PI/N
      SUM1=0.
      SUM2=0.
      DO 1 I=2,N,2
      SUM1=SUM1+Y(I)
      SUM2=SUM2+Y(I-1)
1     CONTINUE
      YB(N-1)=TSN*(SUM1-SUM2)
      YB(N)=TSN*(SUM1+SUM2)
      DO 4 NU=2,NM2,2
      SUM1=0.
      SUM2=0.
      DO 3 I=1,N
      ARG=NU*I*DT
      SUM1=SUM1+Y(I)*COS(ARG)
      SUM2=SUM2+Y(I)*SIN(ARG)
3     CONTINUE
      YB(NU-1)=TSN*SUM1
      YB(NU)=TSN*SUM2
4     CONTINUE
      RETURN
      END
      SUBROUTINE PB(N,XB,X)
C     N MUST BE DIVISIBLE BY 2
      DIMENSION XB(1),X(1)
      PI=3.14159265358979
      DT=PI/N
      NM2=N-2
      SGN=1.
      DO 1 I=1,N
      SGN=-SGN
      SUM=.5*(XB(N)+SGN*XB(N-1))
      DO 2 NU=2,NM2,2
      ARG=NU*I*DT
      SUM=SUM+XB(NU-1)*COS(ARG)+XB(NU)*SIN(ARG)
2     CONTINUE
      X(I)=SUM
1     CONTINUE
      RETURN
      END
      SUBROUTINE SINT(NN,Y,YB)
      DIMENSION Y(1),YB(1)
      N=NN-2
      TSN=2./(N+1)
      CALL SINB(N,Y(2),YB(2))
      DO 1 NU=1,N
      YB(NU+1)=TSN*YB(NU+1)
1     CONTINUE
      YB(1)=0.
      YB(NN)=0.
      RETURN
      END
      SUBROUTINE SINB(N,XB,X)
      DIMENSION XB(1),X(1)
      PI=3.14159265358979
      DT=PI/(N+1)
      DO 1 I=1,N
      SUM=0.
      DO 2 NU=1,N
      ARG=I*NU*DT
      SUM=SUM+SIN(ARG)*XB(NU)
2     CONTINUE
      X(I)=SUM
1     CONTINUE
      RETURN
      END
      SUBROUTINE COST(N,Y,YB)
      DIMENSION Y(1),YB(1)
      TSN=2./(N-1)
      CALL COSB(N,Y,YB)
      DO 1 I=1,N
      YB(I)=TSN*YB(I)
1     CONTINUE
      RETURN
      END
      SUBROUTINE COSB(N,XB,X)
      DIMENSION XB(1),X(1)
      PI=3.14159265358979
      NM2=N-2
      DT=PI/(N-1)
      SGN=-1.
      DO 1 I=1,N
      SGN=-SGN
      SUM=.5*(XB(1)+SGN*XB(N))
      DO 2 NU=1,NM2
      ARG=(I-1)*NU*DT
      SUM=SUM+XB(NU+1)*COS(ARG)
2     CONTINUE
      X(I)=SUM
1     CONTINUE
      RETURN
      END
      FUNCTION FAC(L)
      FAC=1.
      LH=L
1     IF(LH)3,3,2
2     FAC=LH*FAC
      LH=LH-1
      GO TO 1
3     RETURN
      END
      FUNCTION  GAM(L)
      GAM=1.
      FL=L
1     IF(FL) 3,3,2
2     GAM=FL*GAM
      FL=FL-2
      GO TO 1
3     RETURN
      END
      SUBROUTINE LFK (N,M,CP,W)                                                 
C                                                                               
C     FILE ALFPAC DOCUMENTATION ,DEC 1976 , VERSION 0                           
C     PAUL N SWARZTRAUBER , NATIONAL CENTER FOR ATMOSPHERIC RESEARCH            
C     BOULDER COLORADO 80303 , WHICH IS SPONSORED BY THE NATIONAL               
C     SCIENCE FOUNDATION                                                        
C                                                                               
C     THIS FILE CONTAINS THE DOCUMENTATION FOR FILE ALFPAC                      
C     WHICH CONTAINS SEVERAL PROGRAMS FOR COMPUTING THE NORMALIZED              
C     ASSOCIATED LEGENDRE FUNCTIONS PBAR(N,M,THETA) WHERE THETA                 
C     IS COLATITUDE.  THEY ARE NORMALIZED SO THAT THE INTEGRAL OF               
C     PBAR(N,M,THETA)**2*SIN(THETA) ON THE INTERVAL THETA=0 TO THETA=PI         
C     IS EQUAL TO 1.  THE RELATION BETWEEN PBAR(N,M,THETA) AND THE              
C     ASSOCIATED LEGENDRE FUNCTION P(N,M,THETA) IS GIVEN BY -                   
C                                                                               
C     PBAR(N,M,THETA)=SQRT((N+.5)*FACTORIAL(N-M)/FACTORIAL(N+M))                
C                                                                               
C                             *P(N,M,THETA)                                     
C                                                                               
C     FOR MORE DETAIL SEE BELOUSOV - TABLES OF NORMALIZED ASSOCIATED            
C     LEGENDRE POLYNOMIALS, MACMILLAN,1962                                      
C                                                                               
C     A BRIEF DISCRIPTION OF EACH SUBROUTINE FOLLOWS WITH THE DETAILS           
C     OF THE ARGUMENTS GIVEN FURTHER BELOW                                      
C                                                                               
C     SUBROUTINE LFK                                                            
C                                                                               
C     SUBROUTINE LFK COMPUTES THE FOURIER COEFFICIENTS IN THE                   
C     TRIGNOMETRIC SERIES REPRESENTATION OF PBAR(N,M,THETA)                     
C                                                                               
C     SUBROUTINE LFP                                                            
C                                                                               
C     SUBROUTINE LFP USES THE COEFFICIENTS COMPUTED BY SUBROUTINE               
C     LFK AND TABULATES PBAR(N,M,THETA) AT THE COLATITUDES                      
C     THETA(I)=(I-1)*PI/(L-1) I=1,...L.                                         
C                                                                               
C     SUBROUTINE LFPT                                                           
C                                                                               
C     SUBROUTINE LFPT USES THE COEFFICIENTS COMPUTED BY SUBROUTINE              
C     LFK TO COMPUTE PBAR(N,M,THETA) AT A FIXED COLATITUDE THETA.               
C                                                                               
C     SUBROUTINE LFMA                                                           
C                                                                               
C     GIVEN N,I AND L THEN SUBROUTINE LFM TABULATES PBAR(N,M,THETA)             
C     FOR M=0,...,N AT COLATITUDE THETA=(I-1)*PI/(L-1)                          
C                                                                               
C                                                                               
C     SUBROUTINE LFMB                                                           
C                                                                               
C     LFMB IS THE SAME AS LFMA EXCEPT THAT IT WILL BE MORE EFFICIENT            
C     IN CERTIAN APPLICATIONS. SEE THE DISCRIPTIONS OF LFMA AND LFMB            
C     GIVEN BELOW                                                               
C                                                                               
C     SUBROUTINE LFNA                                                           
C                                                                               
C     GIVEN M, I AND L THEN SUBROUTINE LFNA TABULATES PBAR(N,M,THETA)           
C     FOR N=M,...,L-1 AT THETA=(I-1)*PI/(L-1)                                   
C                                                                               
C     SUBROUTINE LFNB                                                           
C                                                                               
C     LFNB IS THE SAME AS LFNA EXCEPT THAT IT WILL BE MORE EFFICIENT            
C     IN CERTIAN APPLICATIONS. SEE THE DISCRIPTIONS OF LFNA AND LFNB            
C     GIVEN BELOW                                                               
C                                                                               
C     SUBROUTINE BELSOV                                                         
C                                                                               
C     GIVEN N,L AND I, SUBROUTINE BELSOV USES THE METHOD OF BELOUSOV            
C     TO TABULATE PBAR(N,M,THETA) FOR M=0,...,N AT THETA=(I-1)*PI/(L-1).        
C                                                                               
C                                                                               
C     THE ACCURACY OF THE SUBROUTINES WAS CHECKED USING DOUBLE PRECISION        
C     VERSIONS OF EACH PROGRAM WHICH WITH THE EXCEPTION OF BELSOV ARE           
C     ALSO ON THIS FILE. THESE PROGRAMS CAN BE ACCESSED BY PRECEEDING           
C     THE NAME WITH D.  I.E. DLFK CALLS THE DOUBLE VERSION OF LFK.              
C                                                                               
C                 TABLE OF MAXIMUM ABSOLUTE ERRORS                              
C                 ON THE CONTROL DATA 7600 FOR                                  
C                 N=0,...,L  M=0,...,N AND I=1,...,L                            
C                                                                               
C                 THE 7600 HAS A MACHINE EPSILON OF                             
C                 2.**(-48)=3.6E-15                                             
C                                                                               
C                                                                               
C              LFK      LFMA     LFMB     LFNA     LFNB    BELSOV               
C                                                                               
C                                                                               
C     L= 51  7.5E-12  7.5E-12  2.8E-12  3.6E-11  2.2E-11  9.6E-13               
C                                                                               
C     L=101  8.0E-11  2.2E-11  7.8E-12  2.9E-10  1.2E-10  3.4E-12               
C                                                                               
C     L=201  6.3E-10  4.8E-11  1.6E-11  3.3E-09  1.2E-09  7.8E-12               
C                                                                               
C                                                                               
C                 TIMES IN MILLISECONDS FOR N=0,...,L                           
C                 M=0,...,N AND I=1,...,L ON THE                                
C                 CONTROL DATA 7600.                                            
C                                                                               
C                                                                               
C              LFK      LFMA     LFMB     LFNA     LFNB    BELSOV               
C                                                                               
C     L= 51    264       722     1279      870     1246     1752                
C                                                                               
C     L=101   1787      4789     8118     5761     8837    12391                
C                                                                               
C     L=201  12773     33965    56123    40993    65499    90998                
C                                                                               
C                                                                               
C     SUBROUTINE LFP TAKES APPROXIMATELY L TIMES SUBROUTINE LFK.  WHEN          
C     THE FAST FOURIER TRANSFORM IS INCORPORATED INTO LFP IT WILL TAKE          
C     ABOUT LOG2(L) TIMES SUBROUTINE LFK.                                       
C                                                                               
C                                                                               
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *         
C                                                                               
C     SUBROUTINE LFK (N,M,CP,W)                                                 
C                                                                               
C     SUBROUTINE LFK COMPUTES THE FOURIER COEFFICIENTS IN THE                   
C     TRIGNOMETRIC SERIES REPRESENTATION OF PBAR(N,M,THETA).  LFK IS            
C     INTENDED TO BE FOLLOWED BY A CALL OF SUBROUTINE LFP WHICH                 
C     TABULATES THE FUNCTION. THE FORM OF THE TRIGNOMETRIC SERIES               
C     IS GIVEN BELOW IN THE DOCUMENTATION FOR SUBROUTINE LFP                    
C                                                                               
C     IN ADDITION TO COMPUTING THE ASSOCIATED LEGENDRE FUNCTIONS                
C     THESE COEFFICIENTS CAN BE USED TO OBTAIN ACCURATE VALUES                  
C     FOR EXPRESSIONS INVOLVING DERIVATIVES OR INTEGRALS OF THE                 
C     FUNCTIONS                                                                 
C                                                                               
C     INPUT PARAMETERS                                                          
C                                                                               
C     N      THE DEGREE (SUBSCRIPT) OF PBAR(N,M,THETA)                          
C                                                                               
C     M      THE ORDER (SUPERSCRIPT) OF PBAR(N,M,THETA)                         
C                                                                               
C     NOTE-  LFK CAN BE CALLED IN ANY ORDER OF N AND M, HOWEVER                 
C            IT IS SOMEWHAT MORE EFFICIENT IF M CHANGES MORE                    
C            RAPIDLY THAN N                                                     
C                                                                               
C     W      A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST                    
C            INT(2.5*N)+4. IF M=0 OR M=1 THIS ARRAY IS NOT USED AND             
C            W MAY BE A DUMMY VARIABLE.                                         
C                                                                               
C     OUTPUT PARAMETERS                                                         
C                                                                               
C     CP     AN ARRAY OF LENGTH AT LEAST INT(N/2)+1 LOCATIONS WHICH             
C            CONTAINS THE COEFFICIENTS                                          
C                                                                               
C     W      CONTAINS RESULTS WHICH MUST NOT BE DESTROYED IF LFK                
C            WILL BE CALLED AGAIN WITH N UNCHANGED. THIS DOES NOT APPLY         
C            IF M=0 OR M=1 IN WHICH CASE W MAY BE A DUMMY VARIABLE.             
C                                                                               
      DIMENSION       CP(1)      ,W(1)                                          
C                                                                               
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR         
C                                                                               
      CALL Q8QST4 (4HXLIB,6HALFPAC,3HLFK,10HVERSION  1)                         
      IDW = N/2+1                                                               
      IW1 = N+2                                                                 
      IW2 = IW1+IDW                                                             
      IW3 = IW2+IDW                                                             
      CALL LFK1 (N,M,CP,W,W(IW1),W(IW2),W(IW3))                                 
      RETURN                                                                    
      END                                                                       
      SUBROUTINE LFP (N,M,L,CP,P)                                               
C                                                                               
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *         
C                                                                               
C     SUBROUTINE LFP (N,M,L,CP,P)                                               
C                                                                               
C     SUBROUTINE LFP USES THE COEFFICIENTS CP GENERATED BY SUBROUTINE           
C     LFK AND TABULATES PBAR(N,M,THETA) AT THETA=(I-1)*PI/(L-1) FOR             
C     I=1,...,L.                                                                
C                                                                               
C     INPUT PARAMETERS                                                          
C                                                                               
C     N      THE DEGREE (SUBSCRIPT) OF PBAR(N,M,THETA)                          
C                                                                               
C     M      THE ORDER (SUPERSCRIPT) OF PBAR(N,M,THETA)                         
C                                                                               
C     CP     AN ARRAY OF LENGTH AT LEAST INT(N/2)+1 LOCATIONS WHICH             
C            CONTAINS THE COEFFICIENTS COMPUTED BY SUBROUTINE LFK.              
C                                                                               
C     L      THE NUMBER OF COLATITUDES THETA(I)=(I-1)*PI/(L-1)                  
C            FOR I=1,...,L.                                                     
C                                                                               
C     OUTPUT PARAMETER                                                          
C                                                                               
C     P      AN ARRAY OF LENGTH L WHICH CONTAINS A TABULATION OF                
C            PBAR(N,M,THETA).   THE TABULATION IS COMPUTED FROM ONE             
C            OF THE FOLLOWING SERIES DEPENDING ON THE PARITY OF N AND M         
C                                                                               
C        FOR N EVEN,M EVEN AND FOR I=1,...L                                     
C                                                                               
C            P(I)=.5*CP(1) PLUS THE SUM FROM K=1 TO K=N/2 OF                    
C                                                                               
C                 CP(K)*COS(2*K*THETA(I))                                       
C                                                                               
C        FOR N EVEN,M ODD AND FOR I=1,...L                                      
C                                                                               
C            P(I)=THE SUM FROM K=1 TO K=N/2 OF                                  
C                                                                               
C                 CP(K)*SIN(2*K*THETA(I))                                       
C                                                                               
C        FOR N ODD,M EVEN AND FOR I=1,...L                                      
C                                                                               
C            P(I)=THE SUM FROM K=1 TO K=(N+1)/2 OF                              
C                                                                               
C                 CP(K)*COS((2*K-1)*THETA(I))                                   
C                                                                               
C        FOR N ODD,M ODD AND FOR I=1,...L                                       
C                                                                               
C            P(I)=THE SUM FROM K=1 TO K=(N+1)/2 OF                              
C                                                                               
C                 CP(K)*SIN((2*K-1)*THETA(I))                                   
C                                                                               
C                                                                               
      DIMENSION       CP(1)      ,P(1)                                          
C                                                                               
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR         
C                                                                               
      CALL Q8QST4 (4HXLIB,6HALFPAC,3HLFP,10HVERSION  1)                         
      IF (N)  10, 10, 40                                                        
   10 IF (M)  20, 20, 40                                                        
   20 DO  30 I=1,L                                                              
         P(I) = .707106781186548                                                
   30 CONTINUE                                                                  
      GO TO 240                                                                 
   40 LS2 = (L+1)/2                                                             
      LM1 = L-1                                                                 
      NP1 = N+1                                                                 
      PI = 3.14159265358979                                                     
      DT = PI/FLOAT(LM1)                                                        
      NMOD = MOD(N,2)                                                           
      MMOD = MOD(M,2)                                                           
      IF (NMOD)  50, 50,120                                                     
   50 IF (MMOD)  60, 60, 90                                                     
   60 KDO = N/2+1                                                               
      DO  80 I=1,LS2                                                            
         TH = FLOAT(I-1)*DT                                                     
         CDT = COS(TH+TH)                                                       
         SDT = SIN(TH+TH)                                                       
         CT = 1.                                                                
         ST = 0.                                                                
         SUM = .5*CP(1)                                                         
         DO  70 KP1=2,KDO                                                       
C                                                                               
C     K=KP1-1                                                                   
C                                                                               
            CTH = CDT*CT-SDT*ST                                                 
            ST = SDT*CT+CDT*ST                                                  
            CT = CTH                                                            
            SUM = SUM+CP(KP1)*CT                                                
C                                                                               
C     SUM=SUM+CP(KP1)*COS((K+K)*TH)                                             
C                                                                               
   70    CONTINUE                                                               
         P(I) = SUM                                                             
   80 CONTINUE                                                                  
      GO TO 190                                                                 
   90 KDO = N/2                                                                 
      DO 110 I=1,LS2                                                            
         TH = FLOAT(I-1)*DT                                                     
         CDT = COS(TH+TH)                                                       
         SDT = SIN(TH+TH)                                                       
         CT = 1.                                                                
         ST = 0.                                                                
         SUM = 0.                                                               
         DO 100 K=1,KDO                                                         
            CTH = CDT*CT-SDT*ST                                                 
            ST = SDT*CT+CDT*ST                                                  
            CT = CTH                                                            
            SUM = SUM+CP(K)*ST                                                  
C                                                                               
C     SUM=SUM+CP(K)*SIN((K+K)*TH)                                               
C                                                                               
  100    CONTINUE                                                               
         P(I) = SUM                                                             
  110 CONTINUE                                                                  
      GO TO 190                                                                 
  120 KDO = (N+1)/2                                                             
      IF (MMOD) 130,130,160                                                     
  130 DO 150 I=1,LS2                                                            
         TH = FLOAT(I-1)*DT                                                     
         CDT = COS(TH+TH)                                                       
         SDT = SIN(TH+TH)                                                       
         CT = COS(TH)                                                           
         ST = -SIN(TH)                                                          
         SUM = 0.                                                               
         DO 140 K=1,KDO                                                         
            CTH = CDT*CT-SDT*ST                                                 
            ST = SDT*CT+CDT*ST                                                  
            CT = CTH                                                            
            SUM = SUM+CP(K)*CT                                                  
C                                                                               
C     SUM=SUM+CP(K)*COS((K+K-1)*TH)                                             
C                                                                               
  140    CONTINUE                                                               
         P(I) = SUM                                                             
  150 CONTINUE                                                                  
      GO TO 190                                                                 
  160 DO 180 I=1,LS2                                                            
         TH = FLOAT(I-1)*DT                                                     
         CDT = COS(TH+TH)                                                       
         SDT = SIN(TH+TH)                                                       
         CT = COS(TH)                                                           
         ST = -SIN(TH)                                                          
         SUM = 0.                                                               
         DO 170 K=1,KDO                                                         
            CTH = CDT*CT-SDT*ST                                                 
            ST = SDT*CT+CDT*ST                                                  
            CT = CTH                                                            
            SUM = SUM+CP(K)*ST                                                  
C                                                                               
C     SUM=SUM+CP(K)*SIN((K+K-1)*TH)                                             
C                                                                               
  170    CONTINUE                                                               
         P(I) = SUM                                                             
  180 CONTINUE                                                                  
  190 IF (MOD(N-M,2)) 220,200,220                                               
  200 DO 210 I=1,LS2                                                            
         LMI = L-I                                                              
         P(LMI+1) = P(I)                                                        
  210 CONTINUE                                                                  
      GO TO 240                                                                 
  220 DO 230 I=1,LS2                                                            
         LMI = L-I                                                              
         P(LMI+1) = -P(I)                                                       
  230 CONTINUE                                                                  
  240 RETURN                                                                    
      END                                                                       
      SUBROUTINE LFPT (N,M,TH,CP,P)                                             
C                                                                               
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *         
C                                                                               
C     SUBROUTINE LFPT (N,M,THETA,CP,P)                                          
C                                                                               
C     SUBROUTINE LFPT USES THE COEFFICIENTS CP COMPUTED BY SUBROUTINE           
C     LFK TO COMPUTE PBAR(N,M,THETA) AT A SINGLE COLATITUDE THETA.              
C     IF A TABULATION ON A COLATITUDE GRID IS DESIRED THEN SUBROUTINE           
C     LFP SHOULD BE USED                                                        
C                                                                               
C     INPUT PARAMETERS                                                          
C                                                                               
C     N      THE DEGREE (SUBSCRIPT) OF PBAR(N,M,THETA)                          
C                                                                               
C     M      THE ORDER (SUPERSCRIPT) OF PBAR(N,M,THETA)                         
C                                                                               
C     CP     AN ARRAY OF LENGTH AT LEAST INT(N/2)+1 LOCATIONS WHICH             
C            CONTAINS THE COEFFICIENTS COMPUTED BY SUBROUTINE LFK.              
C                                                                               
C     TH     THE COLATITUDE IN RADIANS                                          
C                                                                               
C     OUTPUT PARAMETER                                                          
C                                                                               
C     P      PBAR(N,M,THETA) COMPUTED FROM THE APPROPRIATE SERIES               
C            GIVEN IN THE DOCUMENTATION FOR SUBROUTINE LFP                      
C                                                                               
      DIMENSION       CP(1)                                                     
C                                                                               
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR         
C                                                                               
      CALL Q8QST4 (4HXLIB,6HALFPAC,4HLFPT,10HVERSION  1)                        
      IF (N)  10, 10, 30                                                        
   10 IF (M)  20, 20, 30                                                        
   20 P = SQRT(.5)                                                              
      GO TO 140                                                                 
   30 NP1 = N+1                                                                 
      NMOD = MOD(N,2)                                                           
      MMOD = MOD(M,2)                                                           
      IF (NMOD)  40, 40, 90                                                     
   40 IF (MMOD)  50, 50, 70                                                     
   50 KDO = N/2+1                                                               
      CDT = COS(TH+TH)                                                          
      SDT = SIN(TH+TH)                                                          
      CT = 1.                                                                   
      ST = 0.                                                                   
      SUM = .5*CP(1)                                                            
      DO  60 KP1=2,KDO                                                          
         CTH = CDT*CT-SDT*ST                                                    
         ST = SDT*CT+CDT*ST                                                     
         CT = CTH                                                               
         SUM = SUM+CP(KP1)*CT                                                   
   60 CONTINUE                                                                  
      P = SUM                                                                   
      GO TO 140                                                                 
   70 KDO = N/2                                                                 
      CDT = COS(TH+TH)                                                          
      SDT = SIN(TH+TH)                                                          
      CT = 1.                                                                   
      ST = 0.                                                                   
      SUM = 0.                                                                  
      DO  80 K=1,KDO                                                            
         CTH = CDT*CT-SDT*ST                                                    
         ST = SDT*CT+CDT*ST                                                     
         CT = CTH                                                               
         SUM = SUM+CP(K)*ST                                                     
   80 CONTINUE                                                                  
      P = SUM                                                                   
      GO TO 140                                                                 
   90 KDO = (N+1)/2                                                             
      IF (MMOD) 100,100,120                                                     
  100 CDT = COS(TH+TH)                                                          
      SDT = SIN(TH+TH)                                                          
      CT = COS(TH)                                                              
      ST = -SIN(TH)                                                             
      SUM = 0.                                                                  
      DO 110 K=1,KDO                                                            
         CTH = CDT*CT-SDT*ST                                                    
         ST = SDT*CT+CDT*ST                                                     
         CT = CTH                                                               
         SUM = SUM+CP(K)*CT                                                     
  110 CONTINUE                                                                  
      P = SUM                                                                   
      GO TO 140                                                                 
  120 CDT = COS(TH+TH)                                                          
      SDT = SIN(TH+TH)                                                          
      CT = COS(TH)                                                              
      ST = -SIN(TH)                                                             
      SUM = 0.                                                                  
      DO 130 K=1,KDO                                                            
         CTH = CDT*CT-SDT*ST                                                    
         ST = SDT*CT+CDT*ST                                                     
         CT = CTH                                                               
         SUM = SUM+CP(K)*ST                                                     
  130 CONTINUE                                                                  
      P = SUM                                                                   
  140 RETURN                                                                    
      END                                                                       
      SUBROUTINE LFMA (N,L,I,P,W)                                               
C                                                                               
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *         
C                                                                               
C     SUBROUTINE LFMA (N,L,I,P,W)                                               
C                                                                               
C     GIVEN N,I AND L THEN SUBROUTINE LFMA TABULATES PBAR(N,M,THETA)            
C     FOR M=0,...,N AT THETA=(I-1)*PI/(L-1)                                     
C                                                                               
C     NOTE-  THE REASON FOR SPECIFYING L AND I RATHER THAN JUST THETA           
C            IS TO ENABLE LFMA TO MAKE CERTAIN INITIAL CALCULATIONS             
C            WHICH DO NOT HAVE TO BE REPEATED WHEN I CHANGES.                   
C                                                                               
C                                                                               
C     INPUT PARAMETERS                                                          
C                                                                               
C     N      THE DEGREE (SUBSCRIPT) OF PBAR(N,M,THETA)                          
C                                                                               
C     L      THE NUMBER OF COLATITUDES INCLUDING THE POLES                      
C                                                                               
C     I      THE INDEX OF THE COLATITUDE THETA=(I-1)*PI/(L-1)                   
C                                                                               
C                                                                               
C     NOTE-  LFMA CAN BE CALLED IN ANY ORDER OF N,L AND I HOWEVER IT            
C            IS MORE EFFICIENT IF I CHANGES MORE RAPIDLY THAN EITHER N          
C            OR L.                                                              
C                                                                               
C     W      A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 4*L                
C                                                                               
C     OUTPUT PARAMETERS                                                         
C                                                                               
C     W      CONTAINS RESULTS WHICH MUST NOT BE DESTROYED IF LFMA               
C            WILL BE CALLED AGAIN WITH N AND L UNCHANGED                        
C                                                                               
C     P      AN ARRAY OF LENGTH AT LEAST N+1 WHICH CONTAINS                     
C            PBAR(N,M,THETA) STORED AS P(M+1)=PBAR(N,M,THETA) M=0,...,N.        
C                                                                               
      DIMENSION       P(1)       ,W(1)                                          
C                                                                               
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR         
C                                                                               
      CALL Q8QST4 (4HXLIB,6HALFPAC,4HLFMA,10HVERSION  1)                        
      IW1 = L+1                                                                 
      IW2 = IW1+L                                                               
      IW3 = IW2+L                                                               
      CALL LFMA1 (N,L,I,P,W,W(IW1),W(IW2),W(IW3))                               
      RETURN                                                                    
      END                                                                       
      SUBROUTINE LFMB (N,L,I,P,W)                                               
C                                                                               
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *         
C                                                                               
C     SUBROUTINE LFMB (N,L,I,P,W)                                               
C                                                                               
C     GIVEN N,I AND L THEN SUBROUTINE LFMA TABULATES PBAR(N,M,THETA)            
C     FOR M=0,...,N AT THETA=(I-1)*PI/(L-1)                                     
C                                                                               
C     NOTE-  THE REASON FOR SPECIFYING L AND I RATHER THAN JUST THETA           
C            IS TO ENABLE LFMB TO MAKE CERTAIN INITIAL CALCULATIONS             
C            WHICH DO NOT HAVE TO BE REPEATED WHEN I CHANGES.                   
C                                                                               
C                                                                               
C     INPUT PARAMETERS                                                          
C                                                                               
C     N      THE DEGREE (SUBSCRIPT) OF PBAR(N,M,THETA)                          
C                                                                               
C     NOTE-  LFMB CAN BE CALLED IN ANY ORDER OF N,L AND I HOWEVER IT            
C            IS MORE EFFICIENT IF N CHANGES MORE RAPIDLY THAN EITHER I          
C            OR L. THIS NOTE AND THE DISCRIPTION OF THE OUTPUT PARAMETER        
C            W ARE THE ONLY DIFFERENCES BETWEEN SUBROUTINES LFMA                
C            AND LFMB.                                                          
C                                                                               
C     L      THE NUMBER OF COLATITUDES INCLUDING THE POLES                      
C                                                                               
C     I      THE INDEX OF THE COLATITUDE THETA=(I-1)*PI/(L-1)                   
C                                                                               
C                                                                               
C     W      A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 4*L                
C                                                                               
C     OUTPUT PARAMETERS                                                         
C                                                                               
C     W      CONTAINS RESULTS WHICH MUST NOT BE DESTROYED IF LFMB               
C            WILL BE CALLED AGAIN WITH I AND L UNCHANGED                        
C                                                                               
C     P      AN ARRAY OF LENGTH AT LEAST N+1 WHICH CONTAINS                     
C            PBAR(N,M,THETA) STORED AS P(M+1)=PBAR(N,M,THETA) M=0,...,N.        
C                                                                               
      DIMENSION       P(1)       ,W(1)                                          
C                                                                               
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR         
C                                                                               
      CALL Q8QST4 (4HXLIB,6HALFPAC,4HLFMB,10HVERSION  1)                        
      IW1 = L+1                                                                 
      IW2 = IW1+L                                                               
      IW3 = IW2+L                                                               
      CALL LFMB1 (L,N,I,P,W,W,W(IW1),W(IW2),W(IW2),W(IW3))                      
      RETURN                                                                    
      END                                                                       
      SUBROUTINE LFNA (M,L,I,P,W)                                               
C                                                                               
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *         
C                                                                               
C     SUBROUTINE LFNA (M,L,I,P,W)                                               
C                                                                               
C     GIVEN M, I AND L THEN SUBROUTINE LFNA TABULATES PBAR(N,M,THETA)           
C     FOR N=M,...,L-1 AT THETA=(I-1)*PI/(L-1)                                   
C                                                                               
C                                                                               
C     INPUT PARAMETERS                                                          
C                                                                               
C     M      THE ORDER (SUPERSCRIPT) OF PBAR(N,M,THETA)                         
C                                                                               
C     L      THE MAXIMUM VALUE OF N (PLUS ONE).  WHEN APPROXIMATING             
C            FUNCTIONS ON THE SURFACE OF THE SPHERE L IS EQUAL TO               
C            THE NUMBER OF COLATITUDES INCLUDING THE POLES.                     
C                                                                               
C     I      THE INDEX OF THE COLATITUDE THETA=(I-1)*PI/(L-1)                   
C                                                                               
C     NOTE-  LFNA CAN BE CALLED IN ANY ORDER OF M,L OR I HOWEVER IT             
C            IS MORE EFFICIENT IF I CHANGES MORE RAPIDLY THAN EITHER            
C            M OR L.                                                            
C                                                                               
C     W      A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 4*L                
C                                                                               
C     OUTPUT PARAMETERS                                                         
C                                                                               
C     W      CONTAINS RESULTS WHICH MUST NOT BE DESTROYED IF LFNA               
C            WILL BE CALLED AGAIN WITH L AND M UNCHANGED.                       
C                                                                               
C     P      AN ARRAY OF LENGTH AT LEAST L WHICH CONTAINS                       
C            PBAR(N,M,THETA) STORED AS P(N+1)=PBAR(N,M,THETA) FOR               
C            N=M,...,L-1.                                                       
C                                                                               
C                                                                               
      DIMENSION       P(1)       ,W(1)                                          
C                                                                               
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR         
C                                                                               
      CALL Q8QST4 (4HXLIB,6HALFPAC,4HLFNA,10HVERSION  1)                        
      IW1 = L+1                                                                 
      IW2 = IW1+L                                                               
      IW3 = IW2+L                                                               
      CALL LFNA1 (M,L,I,P,W,W(IW1),W(IW2),W(IW3))                               
      RETURN                                                                    
      END                                                                       
      SUBROUTINE LFNB (M,L,I,P,W)                                               
C                                                                               
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *         
C                                                                               
C     SUBROUTINE LFNB (M,L,I,P,W)                                               
C                                                                               
C     GIVEN M, I AND L THEN SUBROUTINE LFNB TABULATES PBAR(N,M,THETA)           
C     FOR N=M,...,L-1 AT THETA=(I-1)*PI/(L-1)                                   
C                                                                               
C                                                                               
C     INPUT PARAMETERS                                                          
C                                                                               
C     M      THE ORDER (SUPERSCRIPT) OF PBAR(N,M,THETA)                         
C                                                                               
C     NOTE-  LFNB CAN BE CALLED IN ANY ORDER OF M,L OR I HOWEVER IT             
C            IS MORE EFFICIENT IF M CHANGES MORE RAPIDLY THAN EITHER            
C            I OR L. THIS NOTE AND THE DISCRIPTION OF THE OUTPUT                
C            PARAMETER W ARE THE ONLY DIFFERENCES BETWEEN LFNA AND LFNB         
C                                                                               
C     L      THE MAXIMUM VALUE OF N (PLUS ONE).  WHEN APPROXIMATING             
C            FUNCTIONS ON THE SURFACE OF THE SPHERE L IS EQUAL TO               
C            THE NUMBER OF COLATITUDES INCLUDING THE POLES.                     
C                                                                               
C     I      THE INDEX OF THE COLATITUDE THETA=(I-1)*PI/(L-1)                   
C                                                                               
C     W      A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 4*L                
C                                                                               
C     OUTPUT PARAMETERS                                                         
C                                                                               
C     W      CONTAINS RESULTS WHICH MUST NOT BE DESTROYED IF LFNB               
C            WILL BE CALLED AGAIN WITH L AND I UNCHANGED.                       
C                                                                               
C     P      AN ARRAY OF LENGTH AT LEAST L WHICH CONTAINS                       
C            PBAR(N,M,THETA) STORED AS P(N+1)=PBAR(N,M,THETA) FOR               
C            N=M,...,L-1.                                                       
C                                                                               
      DIMENSION       P(1)       ,W(1)                                          
C                                                                               
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR         
C                                                                               
      CALL Q8QST4 (4HXLIB,6HALFPAC,4HLFNB,10HVERSION  1)                        
      IW1 = L+1                                                                 
      IW2 = IW1+L                                                               
      IW3 = IW2+L                                                               
      CALL LFNB1 (M,L,I,P,W,W(IW1),W(IW2),W(IW3))                               
      RETURN                                                                    
      END                                                                       
      SUBROUTINE BELSOV (L,N,I,P,W)                                             
C                                                                               
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *         
C                                                                               
C     SUBROUTINE BELSOV(L,N,I,P,W)                                              
C                                                                               
C     GIVEN N,L AND I, SUBROUTINE BELSOV USES THE METHOD OF BELOUSOV            
C     DESCRIBED IN (TABLES OF NORMALIZED ASSOCIATED LEGENDRE FUNCTIONS)         
C     TO TABULATE PBAR(N,M,THETA) FOR M=0,...,N AND THETA=(I-1)*PI/(L-1)        
C     IT MUST BE CALLED STARTING WITH N=0 AND PROCEEDING THROUGH N=L-1          
C     WITH N INCREMENTED BY ONE BETWEEN CALLS.                                  
C                                                                               
C                                                                               
C     INPUT PARAMETERS                                                          
C                                                                               
C     L      THE MAXIMUM VALUE OF N (PLUS ONE). WHEN APPROXIMATING              
C            FUNCTIONS ON THE SURFACE OF THE SPHERE L IS EQUAL TO               
C            THE NUMBER OF COLATITUDES INCLUDING THE POLES.                     
C                                                                               
C     N      THE DEGREE (SUBSCRIPT) OF THE LEGENDRE FUNCTION. BELSOV            
C            MUST BE CALLED STARTING WITH N=0 AND PROCEEDING THROUGH            
C            N=L-1 WITH N INCREMENTED BY ONE BETWEEN CALLS.                     
C                                                                               
C     I      THE INDEX OF THE COLATITUDE THETA=(I-1)*PI/(L-1)                   
C                                                                               
C     W      A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 4*L                
C                                                                               
C     OUTPUT PARAMETERS                                                         
C                                                                               
C     P      AFTER EACH CALL OF BELSOV, P(M+1) CONTAINS PBAR(N,M,THETA)         
C            FOR M=0,...,N. THUS P SHOULD BE DIMENSIONED AT LEAST L.            
C                                                                               
C     W      CONTAINS RESULTS WHICH MUST NOT BE DESTROYED BETWEEN               
C            CALLS WHERE ONLY N CHANGES                                         
C                                                                               
      DIMENSION       P(1)       ,W(1)                                          
C                                                                               
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR         
C                                                                               
      CALL Q8QST4 (4HXLIB,6HALFPAC,6HBELSOV,10HVERSION  1)                      
      IW1 = L+1                                                                 
      IW2 = IW1+L                                                               
      IW3 = IW2+L                                                               
      CALL BEL1 (L,N,I,P,W,W,W(IW1),W(IW2),W(IW2),W(IW3))                       
      RETURN                                                                    
      END                                                                       
      SUBROUTINE LFK1 (N,M,CP,D,A,B,C)                                          
      DIMENSION       CP(1)      ,D(1)       ,A(1)       ,B(1)       ,          
     1                C(1)                                                      
      DATA NLAST /-1/                                                           
      IF (M-1)  10, 10,130                                                      
C                                                                               
C     COEFFICIENTS FOR M=0                                                      
C                                                                               
   10 NMOD = MOD(N,2)                                                           
      FN = FLOAT(N)                                                             
      IF (NMOD)  20, 20, 30                                                     
   20 K = -1                                                                    
      CP(1) = 2.*SQRT(FN+.5)*QP(N-1)**2                                         
      GO TO  40                                                                 
   30 K = 0                                                                     
      CP(1) = 2.*SQRT(FN+.5)*(FN+1.)/FN*QP(N)**2                                
   40 KS = 1                                                                    
      FK = K                                                                    
      FNPK = FN+FK                                                              
      FNMK = FN-FK                                                              
   50 K = K+2                                                                   
      IF (N-K)  70, 70, 60                                                      
   60 KS = KS+1                                                                 
      FNPK = FNPK+2.                                                            
      FNMK = FNMK-2.                                                            
      CP(KS) = FNPK*(FNMK+1.)/FNMK*CP(KS-1)/(FNPK+1.)                           
      GO TO  50                                                                 
   70 IF (M) 340,340, 80                                                        
C                                                                               
C     COEFFICIENTS FOR M=1                                                      
C                                                                               
   80 ENM = 1./SQRT(FN*(FN+1.))                                                 
      IF (NMOD)  90, 90,110                                                     
   90 KSM = KS-1                                                                
      TK = 0.                                                                   
      DO 100 K=1,KSM                                                            
         TK = TK+2.                                                             
         CP(K) = ENM*TK*CP(K+1)                                                 
  100 CONTINUE                                                                  
      GO TO 340                                                                 
  110 TK = -1.                                                                  
      DO 120 K=1,KS                                                             
         TK = TK+2.                                                             
         CP(K) = ENM*TK*CP(K)                                                   
  120 CONTINUE                                                                  
      GO TO 340                                                                 
  130 NP1 = N+1                                                                 
      MSQ = M*M                                                                 
      TMSQ = MSQ+MSQ                                                            
      MMOD = MOD(M,2)                                                           
      IF (N-NLAST) 140,250,140                                                  
  140 NLAST = N                                                                 
      NNP1 = N*(N+1)                                                            
      NMOD = MOD(N,2)                                                           
      FN = FLOAT(N)                                                             
      FNNP1 = FN*(FN+1.)                                                        
      FNNP12 = FNNP1-2.                                                         
      FNNP16 = FNNP1-6.                                                         
      HFNM2 = FNNP12/2.                                                         
      TN = FLOAT(N+N)                                                           
      IF (NMOD) 200,150,200                                                     
C                                                                               
C     STARTING COEFFICIENTS FOR N EVEN                                          
C                                                                               
  150 D(1) = 2.*SQRT(FN+.5)*QP(N-1)**2                                          
      IP = 0                                                                    
      MH = 0                                                                    
      FM = 0.                                                                   
      FNMM = FN-FM                                                              
      FNPM = FN+FM                                                              
  160 MH = MH+1                                                                 
      IP = 1-IP                                                                 
      IF (N-MH) 250,170,170                                                     
  170 FM = FM+1.                                                                
      FNMM = FNMM-1.                                                            
      FNPM = FNPM+1.                                                            
      ALPHA = SQRT(FNPM*(FNMM+1.))                                              
      IF (IP) 180,180,190                                                       
  180 D(MH+1) = HFNM2/(FM-1.)*D(MH)/ALPHA                                       
      GO TO 160                                                                 
  190 D(MH+1) = ALPHA/HFNM2*FM*D(MH)                                            
      GO TO 160                                                                 
C                                                                               
C     STARTING COEFFICIENTS FOR N ODD                                           
C                                                                               
  200 D(1) = 2.*SQRT(FN+.5)*(FN+1.)/FN*QP(N)**2                                 
      IP = 0                                                                    
      MH = 0                                                                    
      FM = 0.                                                                   
      FNMM = FN-FM                                                              
      FNPM = FN+FM                                                              
  210 MH = MH+1                                                                 
      IP = 1-IP                                                                 
      IF (N-MH) 250,220,220                                                     
  220 FM = FM+1.                                                                
      FNMM = FNMM-1.                                                            
      FNPM = FNPM+1.                                                            
      ALPHA = SQRT(FNPM*(FNMM+1.))                                              
      IF (IP) 230,230,240                                                       
  230 D(MH+1) = ALPHA/(FM-1.)*D(MH)                                             
      GO TO 210                                                                 
  240 D(MH+1) = FM/ALPHA*D(MH)                                                  
      GO TO 210                                                                 
  250 IF (NMOD) 260,260,300                                                     
C                                                                               
C     N IS EVEN                                                                 
C                                                                               
  260 NH = N/2                                                                  
      NHDO = NH+1                                                               
      TI = -2.                                                                  
      DO 270 IP1=1,NHDO                                                         
         TI = TI+2.                                                             
         A(IP1) = (TI-1.)*(TI-2.)-FNNP1                                         
         B(IP1) = 2.*(FNNP1-TI*TI-TMSQ)                                         
         C(IP1) = (TI+1.)*(TI+2.)-FNNP1                                         
  270 CONTINUE                                                                  
      IF (MMOD) 280,280,290                                                     
C                                                                               
C     N EVEN, M EVEN                                                            
C                                                                               
  280 B(1) = FNNP1-TMSQ                                                         
      C(1) = 2.-FNNP1                                                           
      CP(1) = D(M+1)                                                            
      CP(2) = -A(2)*D(M+1)                                                      
      CALL TRIVEC (NH,A(2),B(2),C(2),CP(2))                                     
      RETURN                                                                    
C                                                                               
C     N EVEN, M ODD                                                             
C                                                                               
  290 B(2) = 2.*FNNP1-TMSQ-4.                                                   
      C(2) = 12.-FNNP1                                                          
      CP(1) = D(M+1)                                                            
      CP(2) = -A(3)*D(M+1)                                                      
      CALL TRIVEC (NH-1,A(3),B(3),C(3),CP(2))                                   
      RETURN                                                                    
C                                                                               
C     N ODD                                                                     
C                                                                               
  300 NH = (N+1)/2                                                              
      TI = 0.                                                                   
      DO 310 I=1,NH                                                             
         TI = TI+2.                                                             
         A(I) = (TI-2.)*(TI-3.)-FNNP1                                           
         B(I) = 2.*(FNNP1-(TI-1.)**2-TMSQ)                                      
         C(I) = TI*(TI+1.)-FNNP1                                                
  310 CONTINUE                                                                  
      IF (MMOD) 320,320,330                                                     
C                                                                               
C     N ODD, M EVEN                                                             
C                                                                               
  320 B(1) = FNNP1-TMSQ-TMSQ-2.                                                 
      C(1) = 6.-FNNP1                                                           
      CP(1) = D(M+1)                                                            
      CP(2) = -A(2)*D(M+1)                                                      
      CALL TRIVEC (NH-1,A(2),B(2),C(2),CP(2))                                   
      RETURN                                                                    
C                                                                               
C     N ODD, M ODD                                                              
C                                                                               
  330 B(1) = 3.*FNNP1-TMSQ-TMSQ-2.                                              
      C(1) = 6.-FNNP1                                                           
      CP(1) = D(M+1)                                                            
      CP(2) = -A(2)*D(M+1)                                                      
      CALL TRIVEC (NH-1,A(2),B(2),C(2),CP(2))                                   
  340 RETURN                                                                    
      END                                                                       
      FUNCTION QP (L)                                                           
      QP = 1.                                                                   
      FL = L                                                                    
   10 IF (FL)  30, 30, 20                                                       
   20 QP = FL*QP/(FL+1.)                                                        
      FL = FL-2.                                                                
      GO TO  10                                                                 
   30 RETURN                                                                    
      END                                                                       
      SUBROUTINE LFMA1 (N,L,I,P,PZ,P1,A,B)                                      
      DIMENSION       P(1)       ,PZ(1)      ,P1(1)      ,A(1)       ,          
     1                B(1)                                                      
      DATA NLAST /-1/, LLAST /-1/                                               
      IF (N)  10, 10, 20                                                        
   10 P(1) = 1./SQRT(2.)                                                        
      GO TO 100                                                                 
   20 IF (N-NLAST)  40, 30, 40                                                  
   30 IF (L-LLAST)  40, 60, 40                                                  
   40 NLAST = N                                                                 
      LLAST = L                                                                 
      CALL LFK (N,0,P,DUMW)                                                     
      CALL LFP (N,0,L,P,PZ)                                                     
      CALL LFK (N,1,P,DUMW)                                                     
      CALL LFP (N,1,L,P,P1)                                                     
      NM1 = N-1                                                                 
      NP1 = N+1                                                                 
      PI = 3.14159265358979                                                     
      DT = PI/FLOAT(L-1)                                                        
      FNPM = N                                                                  
      FNMM = N+1                                                                
      FMPM = 0.                                                                 
      DO  50 M=1,N                                                              
         FNPM = FNPM+1.                                                         
         FNMM = FNMM-1.                                                         
         FMPM = FMPM-2.                                                         
         B(M) = SQRT(FNMM*FNPM)                                                 
   50 CONTINUE                                                                  
   60 THETA = FLOAT(I-1)*DT                                                     
      COST = COS(THETA)                                                         
      SINT = SIN(THETA)                                                         
      FMPM = 0.                                                                 
      DO  70 M=1,N                                                              
         FMPM = FMPM-2.                                                         
         A(M) = B(M)*SINT                                                       
         P(M+1) = FMPM*COST                                                     
   70 CONTINUE                                                                  
      IF (ABS(PZ(I))-ABS(P1(I)))  90, 90, 80                                    
   80 R = -A(1)*PZ(I)                                                           
      CALL TRIH (N,A,P(2),A(2),R)                                               
      P(1) = PZ(I)                                                              
      GO TO 100                                                                 
   90 R = -A(2)*P1(I)                                                           
      CALL TRIH (NM1,A(2),P(3),A(3),R)                                          
      P(1) = PZ(I)                                                              
      P(2) = P1(I)                                                              
  100 RETURN                                                                    
      END                                                                       
      SUBROUTINE LFMB1 (L,N,I,P,PMZ,DPMZ,P1,PM1,DPM1,P2)                        
      DOUBLE PRECISIONDPMZ       ,DPM1       ,DPMZ1      ,DPMZ2      ,          
     1                DPM11      ,DPM12      ,DTH        ,DPI                   
      DIMENSION       P1(1)      ,P2(1)      ,P(1)       ,PMZ(1)     ,          
     1                PM1(1)     ,DPMZ(1)    ,DPM1(1)                           
      DATA ILAST/-1/, LLAST/-1/                                                 
      IF (I-ILAST)  20, 10, 20                                                  
   10 IF (L-LLAST)  20, 50, 20                                                  
   20 ILAST = I                                                                 
      LLAST = L                                                                 
      LM1 = L-1                                                                 
      LM2 = L-2                                                                 
      DPI = 4.*DATAN(1.D0)                                                      
      DTH = (I-1)*DPI/LM1                                                       
      THETA = DTH                                                               
      COST = COS(THETA)                                                         
      SINT = SIN(THETA)                                                         
      CALL DLFK (LM1,0,P1,DUMW)                                                 
      CALL DLFPT (LM1,0,DTH,P1,DPMZ1)                                           
      CALL DLFK (LM1,1,P1,DUMW)                                                 
      CALL DLFPT (LM1,1,DTH,P1,DPMZ2)                                           
      CALL DLFK (LM2,0,P1,DUMW)                                                 
      CALL DLFPT (LM2,0,DTH,P1,DPM11)                                           
      CALL DLFK (LM2,1,P1,DUMW)                                                 
      CALL DLFPT (LM2,1,DTH,P1,DPM12)                                           
      DPMZ(1) = DPMZ1                                                           
      DPMZ(2) = DPM11                                                           
      CALL DRECN (L,0,DTH,DPMZ,PM1)                                             
      DO  30 NH=1,L                                                             
         P(NH) = DPMZ(NH)                                                       
   30 CONTINUE                                                                  
      DPM1(1) = DPMZ2                                                           
      DPM1(2) = DPM12                                                           
      CALL DRECN (L,1,DTH,DPM1,PMZ)                                             
      DO  40 NH=1,L                                                             
         PMZ(NH) = P(NH)                                                        
         PM1(NH) = DPM1(NH)                                                     
   40 CONTINUE                                                                  
   50 IF (N) 130, 60, 70                                                        
   60 P(1) = PMZ(1)                                                             
      GO TO 130                                                                 
   70 IF (N-1) 130, 80, 90                                                      
   80 P(1) = PMZ(2)                                                             
      P(2) = PM1(2)                                                             
      GO TO 130                                                                 
   90 FNPM = N                                                                  
      FNMM = N+1                                                                
      FMPM = 0.                                                                 
      DO 100 M=1,N                                                              
         FNPM = FNPM+1.                                                         
         FNMM = FNMM-1.                                                         
         FMPM = FMPM-2.                                                         
         P1(M) = SINT*SQRT(FNMM*FNPM)                                           
         P2(M) = FMPM*COST                                                      
  100 CONTINUE                                                                  
      IF (ABS(PMZ(N+1))-ABS(PM1(N+1))) 120,120,110                              
  110 P(2) = -P1(1)*PMZ(N+1)                                                    
      CALL TRIVEC (N,P1,P2,P1(2),P(2))                                          
      P(1) = PMZ(N+1)                                                           
      GO TO 130                                                                 
  120 P(3) = -P1(2)*PM1(N+1)                                                    
      CALL TRIVEC (N-1,P1(2),P2(2),P1(3),P(3))                                  
      P(1) = PMZ(N+1)                                                           
      P(2) = PM1(N+1)                                                           
  130 RETURN                                                                    
      END                                                                       
      SUBROUTINE DRECN (L,M,THETA,B,A)                                          
      DIMENSION       A(1)       ,B(1)                                          
      DOUBLE PRECISIONA          ,B          ,COST       ,R1         ,          
     1                R2         ,FNMM       ,FNPN       ,FNPM       ,          
     2                THETA      ,BHOLD                                         
      COST = DCOS(THETA)                                                        
      LMM = L-M                                                                 
      R1 = B(1)                                                                 
      R2 = B(2)                                                                 
      K = 0                                                                     
      N = L-1                                                                   
      FNMM = N-M+1                                                              
      FNPN = N+N+1                                                              
      FNPM = N+M+1                                                              
   10 N = N-1                                                                   
      IF (N-M)  30, 20, 20                                                      
   20 FNMM = FNMM-1.                                                            
      FNPN = FNPN-2.                                                            
      FNPM = FNPM-1.                                                            
      K = K+1                                                                   
      A(K) = DSQRT(FNMM*FNPM/(FNPN*(FNPN+2.)))                                  
      B(K+1) = -COST                                                            
      GO TO  10                                                                 
   30 IF (DABS(R1) .LT. DABS(R2)) GO TO  40                                     
      B(1) = R1                                                                 
      R1 = -A(1)*R1                                                             
      CALL DTRIH (LMM-1,A,B(2),A(2),R1)                                         
      GO TO  50                                                                 
   40 B(1) = R1                                                                 
      B(2) = R2                                                                 
      R2 = -A(2)*R2                                                             
      CALL DTRIH (LMM-2,A(2),B(3),A(3),R2)                                      
   50 NDO = (L+1)/2                                                             
      DO  60 N=1,NDO                                                            
         N1 = L-N                                                               
         BHOLD = B(N1+1)                                                        
         B(N1+1) = B(N)                                                         
         B(N) = BHOLD                                                           
   60 CONTINUE                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE LFNA1 (M,L,I,P,PZ,P1,A,B)                                      
      DIMENSION       P(1)       ,PZ(1)      ,P1(1)      ,A(1)       ,          
     1                B(1)                                                      
      DATA MLAST /-1/, LLAST /-1/                                               
      DO  10 N=1,L                                                              
         P(N) = 0.                                                              
   10 CONTINUE                                                                  
      IF (L-LLAST)  30, 20, 30                                                  
   20 IF (M-MLAST)  30, 70, 30                                                  
   30 MLAST = M                                                                 
      LLAST = L                                                                 
      LM1 = L-1                                                                 
      LM2 = L-2                                                                 
      LMM = L-M                                                                 
      PI = 3.14159265358979                                                     
      CALL LFK (LM2,M,B,PZ)                                                     
      CALL LFP (LM2,M,L,B,P)                                                    
      CALL LFK (LM1,M,B,PZ)                                                     
      CALL LFP (LM1,M,L,B,PZ)                                                   
      DO  40 I=1,L                                                              
         P1(I) = P(I)                                                           
   40 CONTINUE                                                                  
      K = 0                                                                     
      N = L-1                                                                   
      FNMM = N-M+1                                                              
      FNPN = N+N+1                                                              
      FNPM = N+M+1                                                              
   50 N = N-1                                                                   
      IF (N-M)  70, 60, 60                                                      
   60 FNMM = FNMM-1.                                                            
      FNPN = FNPN-2.                                                            
      FNPM = FNPM-1.                                                            
      K = K+1                                                                   
      B(K) = SQRT(FNMM*FNPM/(FNPN*(FNPN+2.)))                                   
      GO TO  50                                                                 
   70 P(L) = PZ(I)                                                              
      IF (LMM-2) 140, 80, 90                                                    
   80 P(L-1) = P1(I)                                                            
      GO TO 140                                                                 
   90 THETA = FLOAT(I-1)*PI/FLOAT(LM1)                                          
      COST = COS(THETA)                                                         
      DO 100 K=1,LMM                                                            
         P(K) = -COST                                                           
         A(K) = B(K)                                                            
  100 CONTINUE                                                                  
      IF (ABS(PZ(I)) .LT. ABS(P1(I))) GO TO 110                                 
      P(1) = PZ(I)                                                              
      R = -A(1)*P(1)                                                            
      CALL TRIH (LMM-1,A,P(2),A(2),R)                                           
      GO TO 120                                                                 
  110 P(1) = PZ(I)                                                              
      P(2) = P1(I)                                                              
      R = -A(2)*P(2)                                                            
      CALL TRIH (LMM-2,A(2),P(3),A(3),R)                                        
  120 NDO = (L+1)/2                                                             
      DO 130 N=1,NDO                                                            
         N1 = L-N                                                               
         PHOLD = P(N1+1)                                                        
         P(N1+1) = P(N)                                                         
         P(N) = PHOLD                                                           
  130 CONTINUE                                                                  
  140 RETURN                                                                    
      END                                                                       
      SUBROUTINE LFNB1 (M,L,I,P,PZ,P1,A,B)                                      
      DIMENSION       P(1)       ,PZ(1)      ,P1(1)      ,A(1)       ,          
     1                B(1)                                                      
      DATA ILAST /-1/, LLAST /-1/                                               
      DO  10 N=1,L                                                              
         P(N) = 0.                                                              
   10 CONTINUE                                                                  
      IF (L-LLAST)  30, 20, 30                                                  
   20 IF (I-ILAST)  30, 50, 30                                                  
   30 ILAST = I                                                                 
      LLAST = L                                                                 
      LM1 = L-1                                                                 
      LM2 = L-2                                                                 
      PI = 3.14159265358979                                                     
      THETA = FLOAT(I-1)*PI/FLOAT(LM1)                                          
      COST = COS(THETA)                                                         
      LS = 0                                                                    
      DO  40 N=LM2,LM1                                                          
         LS = LS+1                                                              
         CALL LFK (N,0,A,DUMW)                                                  
         CALL LFPT (N,0,THETA,A,PZ(LS))                                         
         CALL LFK (N,1,A,DUMW)                                                  
         CALL LFPT (N,1,THETA,A,P1(LS))                                         
   40 CONTINUE                                                                  
      PHOLD1 = PZ(1)                                                            
      PHOLD2 = P1(2)                                                            
      PZ(1) = PZ(2)                                                             
      P1(2) = P1(1)                                                             
      P1(1) = PHOLD1                                                            
      PZ(2) = PHOLD2                                                            
      CALL RECM (L,LM2,I,P1,A,B)                                                
      CALL RECM (L,LM1,I,PZ,A,B)                                                
   50 K = 0                                                                     
      LMM = L-M                                                                 
      IF (LMM-1)  60, 60, 70                                                    
   60 P(L) = PZ(M+1)                                                            
      GO TO 140                                                                 
   70 N = L-1                                                                   
      FNMM = N-M+1                                                              
      FNPN = N+N+1                                                              
      FNPM = N+M+1                                                              
   80 N = N-1                                                                   
      IF (N-M) 100, 90, 90                                                      
   90 FNMM = FNMM-1.                                                            
      FNPN = FNPN-2.                                                            
      FNPM = FNPM-1.                                                            
      K = K+1                                                                   
      A(K) = SQRT(FNMM*FNPM/(FNPN*(FNPN+2.)))                                   
      B(K) = -COST                                                              
      GO TO  80                                                                 
  100 IF (ABS(PZ(M+1)) .LT. ABS(P1(M+1))) GO TO 110                             
      P(1) = PZ(M+1)                                                            
      P(2) = -A(1)*P(1)                                                         
      CALL TRIVEC (LMM-1,A,B,A(2),P(2))                                         
      GO TO 120                                                                 
  110 P(1) = PZ(M+1)                                                            
      P(2) = P1(M+1)                                                            
      P(3) = -A(2)*P(2)                                                         
      CALL TRIVEC (LMM-2,A(2),B(2),A(3),P(3))                                   
  120 NDO = (L+1)/2                                                             
      DO 130 N=1,NDO                                                            
         N1 = L-N                                                               
         PHOLD = P(N1+1)                                                        
         P(N1+1) = P(N)                                                         
         P(N) = PHOLD                                                           
  130 CONTINUE                                                                  
  140 RETURN                                                                    
      END                                                                       
      SUBROUTINE RECM (L,N,I,P,A,B)                                             
      DIMENSION       P(1)       ,A(1)       ,B(1)                              
      PI = 3.14159265358979                                                     
      DT = PI/FLOAT(L-1)                                                        
      FNPM = N                                                                  
      FNMM = N+1                                                                
      FMPM = 0.                                                                 
      THETA = FLOAT(I-1)*DT                                                     
      COST = COS(THETA)                                                         
      SINT = SIN(THETA)                                                         
      DO  10 M=1,N                                                              
         FNPM = FNPM+1.                                                         
         FNMM = FNMM-1.                                                         
         FMPM = FMPM-2.                                                         
         A(M) = SINT*SQRT(FNMM*FNPM)                                            
         B(M) = FMPM*COST                                                       
   10 CONTINUE                                                                  
      IF (ABS(P(1))-ABS(P(2)))  30, 30, 20                                      
   20 P(2) = -A(1)*P(1)                                                         
      CALL TRIVEC (N,A,B,A(2),P(2))                                             
      GO TO  40                                                                 
   30 P(3) = -A(2)*P(2)                                                         
      NM1 = N-1                                                                 
      CALL TRIVEC (NM1,A(2),B(2),A(3),P(3))                                     
   40 RETURN                                                                    
      END                                                                       
      SUBROUTINE BEL1 (L,N,I,P3,PMZ,DPMZ,P1,PM1,DPM1,P2)                        
      DOUBLE PRECISIONDPMZ       ,DPM1       ,DPMZ1      ,DPMZ2      ,          
     1                DPM11      ,DPM12      ,DTH        ,DPI                   
      DIMENSION       P1(1)      ,P2(1)      ,P3(1)      ,PMZ(1)     ,          
     1                PM1(1)     ,DPMZ(1)    ,DPM1(1)                           
      IF (N) 100, 30, 10                                                        
   10 IF (N-1) 100, 20, 60                                                      
   20 P2(1) = P3(1)                                                             
      P3(1) = PMZ(2)                                                            
      P3(2) = PM1(2)                                                            
      GO TO 100                                                                 
   30 LM1 = L-1                                                                 
      LM2 = L-2                                                                 
      DPI = 4.*DATAN(1.D0)                                                      
      DTH = FLOAT(I-1)*DPI/FLOAT(LM1)                                           
      THETA = DTH                                                               
      COST = COS(THETA)                                                         
      SINT = SIN(THETA)                                                         
      CALL DLFK (LM1,0,P1,DUMW)                                                 
      CALL DLFPT (LM1,0,DTH,P1,DPMZ1)                                           
      CALL DLFK (LM1,1,P1,DUMW)                                                 
      CALL DLFPT (LM1,1,DTH,P1,DPMZ2)                                           
      CALL DLFK (LM2,0,P1,DUMW)                                                 
      CALL DLFPT (LM2,0,DTH,P1,DPM11)                                           
      CALL DLFK (LM2,1,P1,DUMW)                                                 
      CALL DLFPT (LM2,1,DTH,P1,DPM12)                                           
      DPMZ(1) = DPMZ1                                                           
      DPMZ(2) = DPM11                                                           
      CALL DRECN (L,0,DTH,DPMZ,PM1)                                             
      DO  40 NH=1,L                                                             
         P3(NH) = DPMZ(NH)                                                      
   40 CONTINUE                                                                  
      DPM1(1) = DPMZ2                                                           
      DPM1(2) = DPM12                                                           
      CALL DRECN (L,1,DTH,DPM1,PMZ)                                             
      DO  50 NH=1,L                                                             
         PMZ(NH) = P3(NH)                                                       
         PM1(NH) = DPM1(NH)                                                     
   50 CONTINUE                                                                  
      P3(1) = PMZ(1)                                                            
      GO TO 100                                                                 
   60 DO  70 M=1,N                                                              
         P1(M) = P2(M)                                                          
         P2(M) = P3(M)                                                          
   70 CONTINUE                                                                  
      FN = N                                                                    
      TFN = FN+FN                                                               
      CN = (TFN+1.)/(TFN-3.)                                                    
      DN = (TFN+1.)/(TFN-1.)                                                    
      EN = DN                                                                   
      P3(N+1) = SQRT((FN+.5)*QP(N+N-1))*SINT**N                                 
      P3(1) = PMZ(N+1)                                                          
      P3(2) = PM1(N+1)                                                          
      M = 1                                                                     
      FM = M                                                                    
      FNPM = FN+FM                                                              
      FNMM = FN-FM                                                              
   80 M = M+1                                                                   
      IF (N-M) 100,100, 90                                                      
   90 FNPM = FNPM+1.                                                            
      FNMM = FNMM-1.                                                            
      TEMP = (FNPM-1.)/(FNPM*(FNPM-2.))                                         
      CC = SQRT(CN*TEMP*(FNPM-3.))                                              
      DD = SQRT(DN*TEMP*(FNMM+1.))                                              
      EE = SQRT(EN*FNMM/FNPM)                                                   
      P3(M+1) = CC*P1(M-1)-COST*(DD*P2(M-1)-EE*P2(M+1))                         
      GO TO  80                                                                 
  100 RETURN                                                                    
      END                                                                       
      SUBROUTINE TRIH (N,A,B,C,R)                                               
      DIMENSION       A(1)       ,B(1)       ,C(1)                              
      IF (N) 120,120, 10                                                        
   10 IF (N-2)  20, 30, 40                                                      
   20 B(1) = R/B(1)                                                             
      RETURN                                                                    
   30 QIH = A(2)                                                                
      BIH = B(2)                                                                
      GO TO  70                                                                 
   40 QIH = A(N)                                                                
      BIH = B(N)                                                                
      DO  60 IDO=3,N                                                            
         I = N-IDO+2                                                            
         IF (ABS(BIH) .LT. ABS(C(I))) GO TO  50                                 
         RATIO = C(I)/BIH                                                       
         C(I) = 0.                                                              
         B(I+1) = QIH/BIH                                                       
         BIH = B(I)-RATIO*QIH                                                   
         QIH = A(I)                                                             
         GO TO  60                                                              
   50    B(I+1) = B(I)/C(I)                                                     
         C(I) = A(I)/C(I)                                                       
         BIH1 = QIH-BIH*B(I+1)                                                  
         QIH = -BIH*C(I)                                                        
         BIH = BIH1                                                             
   60 CONTINUE                                                                  
   70 IF (ABS(BIH) .LT. ABS(C(1))) GO TO  80                                    
      Q2 = QIH/BIH                                                              
      BIH = B(1)-C(1)/BIH*QIH                                                   
      B(1) = R/BIH                                                              
      B(2) = -Q2*B(1)                                                           
      GO TO  90                                                                 
   80 RATIO = BIH/C(1)                                                          
      BIH = QIH-RATIO*B(1)                                                      
      RIH = -RATIO*R                                                            
      B1 = RIH/BIH                                                              
      B(2) = (R-B(1)*B1)/C(1)                                                   
      B(1) = B1                                                                 
   90 IF (N-3) 120,100,100                                                      
  100 DO 110 I=3,N                                                              
         B(I) = -B(I)*B(I-1)-C(I-1)*B(I-2)                                      
  110 CONTINUE                                                                  
  120 RETURN                                                                    
      END                                                                       
      SUBROUTINE TRIVEC (N,A,B,C,E)                                             
      DIMENSION       A(1)       ,B(1)       ,C(1)       ,E(1)                  
      R = E(1)                                                                  
      IF (N) 120,120, 10                                                        
   10 IF (N-2)  20, 30, 40                                                      
   20 E(1) = R/B(1)                                                             
      RETURN                                                                    
   30 QIH = A(2)                                                                
      BIH = B(2)                                                                
      GO TO  70                                                                 
   40 QIH = A(N)                                                                
      BIH = B(N)                                                                
      DO  60 IDO=3,N                                                            
         I = N-IDO+2                                                            
         IF (ABS(BIH) .LT. ABS(C(I))) GO TO  50                                 
         RATIO = C(I)/BIH                                                       
         C(I) = 0.                                                              
         B(I+1) = QIH/BIH                                                       
C                                                                               
C     P(I+1)=0.                                                                 
C     Q(I+1)=QIH/BIH                                                            
C                                                                               
         BIH = B(I)-RATIO*QIH                                                   
         QIH = A(I)                                                             
         GO TO  60                                                              
   50    B(I+1) = B(I)/C(I)                                                     
         C(I) = A(I)/C(I)                                                       
         BIH1 = QIH-BIH*B(I+1)                                                  
         QIH = -BIH*C(I)                                                        
         BIH = BIH1                                                             
C                                                                               
C     P(I+1)=A(I)/C(I)                                                          
C     Q(I+1)=B(I)/C(I)                                                          
C                                                                               
   60 CONTINUE                                                                  
   70 IF (ABS(BIH) .LT. ABS(C(1))) GO TO  80                                    
      Q2 = QIH/BIH                                                              
C                                                                               
C     Q(2)=QIH/BIH                                                              
C                                                                               
      BIH = B(1)-C(1)/BIH*QIH                                                   
      E(1) = R/BIH                                                              
      E(2) = -Q2*E(1)                                                           
C                                                                               
C     E(2)=-Q(2)*E(1)                                                           
C                                                                               
      GO TO  90                                                                 
   80 RATIO = BIH/C(1)                                                          
      BIH = QIH-RATIO*B(1)                                                      
      RIH = -RATIO*R                                                            
      E(1) = RIH/BIH                                                            
      E(2) = (R-B(1)*E(1))/C(1)                                                 
   90 IF (N-3) 120,100,100                                                      
  100 DO 110 I=3,N                                                              
         E(I) = -B(I)*E(I-1)-C(I-1)*E(I-2)                                      
C                                                                               
C     E(I)=-Q(I)*E(I-1)-P(I)*E(I-2)                                             
C                                                                               
  110 CONTINUE                                                                  
  120 RETURN                                                                    
      END                                                                       
      SUBROUTINE DLFK (N,M,CP,W)                                                
      DOUBLE PRECISIONCP         ,W                                             
      DIMENSION       CP(1)      ,W(1)                                          
C                                                                               
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR         
C                                                                               
      CALL Q8QST4 (4HXLIB,6HALFPAC,4HDLFK,10HVERSION  1)                        
      IDW = N/2+1                                                               
      IW1 = N+2                                                                 
      IW2 = IW1+IDW                                                             
      IW3 = IW2+IDW                                                             
      CALL DLFK1 (N,M,CP,W,W(IW1),W(IW2),W(IW3))                                
      RETURN                                                                    
      END                                                                       
      SUBROUTINE DLFK1 (N,M,CP,D,A,B,C)                                         
      DOUBLE PRECISIONCP         ,A          ,B          ,C          ,          
     1                D          ,FN         ,FK         ,FNPK       ,          
     2                FNMK       ,ENM        ,TK         ,TMSQ       ,          
     3                FNNP1      ,FNNP12     ,FNNP16     ,HFNM2      ,          
     4                TN         ,FM         ,FNMM       ,FNPM       ,          
     5                ALPHA      ,TI         ,DQP                               
      DIMENSION       CP(1)      ,D(1)       ,A(1)       ,B(1)       ,          
     1                C(1)                                                      
      DATA NLAST /-1/                                                           
      IF (M-1)  10, 10,130                                                      
C                                                                               
C     COEFFICIENTS FOR M=0                                                      
C                                                                               
   10 NMOD = MOD(N,2)                                                           
      FN = FLOAT(N)                                                             
      IF (NMOD)  20, 20, 30                                                     
   20 K = -1                                                                    
      CP(1) = 2.*DSQRT(FN+.5)*DQP(N-1)**2                                       
      GO TO  40                                                                 
   30 K = 0                                                                     
      CP(1) = 2.*DSQRT(FN+.5)*(FN+1.)/FN*DQP(N)**2                              
   40 KS = 1                                                                    
      FK = K                                                                    
      FNPK = FN+FK                                                              
      FNMK = FN-FK                                                              
   50 K = K+2                                                                   
      IF (N-K)  70, 70, 60                                                      
   60 KS = KS+1                                                                 
      FNPK = FNPK+2.                                                            
      FNMK = FNMK-2.                                                            
      CP(KS) = FNPK*(FNMK+1.)/FNMK*CP(KS-1)/(FNPK+1.)                           
      GO TO  50                                                                 
   70 IF (M) 340,340, 80                                                        
C                                                                               
C     COEFFICIENTS FOR M=1                                                      
C                                                                               
   80 ENM = 1./DSQRT(FN*(FN+1.))                                                
      IF (NMOD)  90, 90,110                                                     
   90 KSM = KS-1                                                                
      TK = 0.                                                                   
      DO 100 K=1,KSM                                                            
         TK = TK+2.                                                             
         CP(K) = ENM*TK*CP(K+1)                                                 
  100 CONTINUE                                                                  
      GO TO 340                                                                 
  110 TK = -1.                                                                  
      DO 120 K=1,KS                                                             
         TK = TK+2.                                                             
         CP(K) = ENM*TK*CP(K)                                                   
  120 CONTINUE                                                                  
      GO TO 340                                                                 
  130 NP1 = N+1                                                                 
      MSQ = M*M                                                                 
      TMSQ = MSQ+MSQ                                                            
      MMOD = MOD(M,2)                                                           
      IF (N-NLAST) 140,250,140                                                  
  140 NLAST = N                                                                 
      NNP1 = N*(N+1)                                                            
      NMOD = MOD(N,2)                                                           
      FN = FLOAT(N)                                                             
      FNNP1 = FN*(FN+1.)                                                        
      FNNP12 = FNNP1-2.                                                         
      FNNP16 = FNNP1-6.                                                         
      HFNM2 = FNNP12/2.                                                         
      TN = FLOAT(N+N)                                                           
      IF (NMOD) 200,150,200                                                     
C                                                                               
C     STARTING COEFFICIENTS FOR N EVEN                                          
C                                                                               
  150 D(1) = 2.*DSQRT(FN+.5)*DQP(N-1)**2                                        
      IP = 0                                                                    
      MH = 0                                                                    
      FM = 0.                                                                   
      FNMM = FN-FM                                                              
      FNPM = FN+FM                                                              
  160 MH = MH+1                                                                 
      IP = 1-IP                                                                 
      IF (N-MH) 250,170,170                                                     
  170 FM = FM+1.                                                                
      FNMM = FNMM-1.                                                            
      FNPM = FNPM+1.                                                            
      ALPHA = DSQRT(FNPM*(FNMM+1.))                                             
      IF (IP) 180,180,190                                                       
  180 D(MH+1) = HFNM2/(FM-1.)*D(MH)/ALPHA                                       
      GO TO 160                                                                 
  190 D(MH+1) = ALPHA/HFNM2*FM*D(MH)                                            
      GO TO 160                                                                 
C                                                                               
C     STARTING COEFFICIENTS FOR N ODD                                           
C                                                                               
  200 D(1) = 2.*DSQRT(FN+.5)*(FN+1.)/FN*DQP(N)**2                               
      IP = 0                                                                    
      MH = 0                                                                    
      FM = 0.                                                                   
      FNMM = FN-FM                                                              
      FNPM = FN+FM                                                              
  210 MH = MH+1                                                                 
      IP = 1-IP                                                                 
      IF (N-MH) 250,220,220                                                     
  220 FM = FM+1.                                                                
      FNMM = FNMM-1.                                                            
      FNPM = FNPM+1.                                                            
      ALPHA = DSQRT(FNPM*(FNMM+1.))                                             
      IF (IP) 230,230,240                                                       
  230 D(MH+1) = ALPHA/(FM-1.)*D(MH)                                             
      GO TO 210                                                                 
  240 D(MH+1) = FM/ALPHA*D(MH)                                                  
      GO TO 210                                                                 
  250 IF (NMOD) 260,260,300                                                     
C                                                                               
C     N IS EVEN                                                                 
C                                                                               
  260 NH = N/2                                                                  
      NHDO = NH+1                                                               
      TI = -2.                                                                  
      DO 270 IP1=1,NHDO                                                         
         TI = TI+2.                                                             
         A(IP1) = (TI-1.)*(TI-2.)-FNNP1                                         
         B(IP1) = 2.*(FNNP1-TI*TI-TMSQ)                                         
         C(IP1) = (TI+1.)*(TI+2.)-FNNP1                                         
  270 CONTINUE                                                                  
      IF (MMOD) 280,280,290                                                     
C                                                                               
C     N EVEN, M EVEN                                                            
C                                                                               
  280 B(1) = FNNP1-TMSQ                                                         
      C(1) = 2.-FNNP1                                                           
      CP(1) = D(M+1)                                                            
      CP(2) = -A(2)*D(M+1)                                                      
      CALL DTRI (NH,A(2),B(2),C(2),CP(2))                                       
      RETURN                                                                    
C                                                                               
C     N EVEN, M ODD                                                             
C                                                                               
  290 B(2) = 2.*FNNP1-TMSQ-4.                                                   
      C(2) = 12.-FNNP1                                                          
      CP(1) = D(M+1)                                                            
      CP(2) = -A(3)*D(M+1)                                                      
      CALL DTRI (NH-1,A(3),B(3),C(3),CP(2))                                     
      RETURN                                                                    
C                                                                               
C     N ODD                                                                     
C                                                                               
  300 NH = (N+1)/2                                                              
      TI = 0.                                                                   
      DO 310 I=1,NH                                                             
         TI = TI+2.                                                             
         A(I) = (TI-2.)*(TI-3.)-FNNP1                                           
         B(I) = 2.*(FNNP1-(TI-1.)**2-TMSQ)                                      
         C(I) = TI*(TI+1.)-FNNP1                                                
  310 CONTINUE                                                                  
      IF (MMOD) 320,320,330                                                     
C                                                                               
C     N ODD, M EVEN                                                             
C                                                                               
  320 B(1) = FNNP1-TMSQ-TMSQ-2.                                                 
      C(1) = 6.-FNNP1                                                           
      CP(1) = D(M+1)                                                            
      CP(2) = -A(2)*D(M+1)                                                      
      CALL DTRI (NH-1,A(2),B(2),C(2),CP(2))                                     
      RETURN                                                                    
C                                                                               
C     N ODD, M ODD                                                              
C                                                                               
  330 B(1) = 3.*FNNP1-TMSQ-TMSQ-2.                                              
      C(1) = 6.-FNNP1                                                           
      CP(1) = D(M+1)                                                            
      CP(2) = -A(2)*D(M+1)                                                      
      CALL DTRI (NH-1,A(2),B(2),C(2),CP(2))                                     
  340 RETURN                                                                    
      END                                                                       
      FUNCTION DQP(L)                                                           
      DOUBLE PRECISIONDQP        ,FL                                            
      DQP = 1.                                                                  
      FL = L                                                                    
   10 IF (FL)  30, 30, 20                                                       
   20 DQP = FL*DQP/(FL+1.)                                                      
      FL = FL-2.                                                                
      GO TO  10                                                                 
   30 RETURN                                                                    
      END                                                                       
      SUBROUTINE DLFP (N,M,L,CP,P)                                              
      DOUBLE PRECISIONCP         ,P          ,DT         ,TH         ,          
     1                CDT        ,SDT        ,CT         ,ST         ,          
     2                SUM        ,CTH        ,PI                                
      DIMENSION       CP(1)      ,P(1)                                          
C                                                                               
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR         
C                                                                               
      CALL Q8QST4 (4HXLIB,6HALFPAC,4HDLFP,10HVERSION  1)                        
      IF (N)  10, 10, 40                                                        
   10 IF (M)  20, 20, 40                                                        
   20 DO  30 I=1,L                                                              
         P(I) = DSQRT(5.D-1)                                                    
   30 CONTINUE                                                                  
      GO TO 240                                                                 
   40 LS2 = (L+1)/2                                                             
      LM1 = L-1                                                                 
      NP1 = N+1                                                                 
      PI = 4.*DATAN(1.D0)                                                       
      DT = PI/FLOAT(LM1)                                                        
      NMOD = MOD(N,2)                                                           
      MMOD = MOD(M,2)                                                           
      IF (NMOD)  50, 50,120                                                     
   50 IF (MMOD)  60, 60, 90                                                     
   60 KDO = N/2+1                                                               
      DO  80 I=1,LS2                                                            
         TH = FLOAT(I-1)*DT                                                     
         CDT = DCOS(TH+TH)                                                      
         SDT = DSIN(TH+TH)                                                      
      IF(I.EQ.1) SDT=0.                                                         
         CT = 1.                                                                
         ST = 0.                                                                
         SUM = .5*CP(1)                                                         
         DO  70 KP1=2,KDO                                                       
C                                                                               
C     K=KP1-1                                                                   
C                                                                               
            CTH = CDT*CT-SDT*ST                                                 
            ST = SDT*CT+CDT*ST                                                  
            CT = CTH                                                            
            SUM = SUM+CP(KP1)*CT                                                
C                                                                               
C     SUM=SUM+CP(KP1)*COS((K+K)*TH)                                             
C                                                                               
   70    CONTINUE                                                               
         P(I) = SUM                                                             
   80 CONTINUE                                                                  
      GO TO 190                                                                 
   90 KDO = N/2                                                                 
      DO 110 I=1,LS2                                                            
         TH = FLOAT(I-1)*DT                                                     
         CDT = DCOS(TH+TH)                                                      
         SDT = DSIN(TH+TH)                                                      
      IF(I.EQ.1) SDT=0.                                                         
         CT = 1.                                                                
         ST = 0.                                                                
         SUM = 0.                                                               
         DO 100 K=1,KDO                                                         
            CTH = CDT*CT-SDT*ST                                                 
            ST = SDT*CT+CDT*ST                                                  
            CT = CTH                                                            
            SUM = SUM+CP(K)*ST                                                  
C                                                                               
C     SUM=SUM+CP(K)*SIN((K+K)*TH)                                               
C                                                                               
  100    CONTINUE                                                               
         P(I) = SUM                                                             
  110 CONTINUE                                                                  
      GO TO 190                                                                 
  120 KDO = (N+1)/2                                                             
      IF (MMOD) 130,130,160                                                     
  130 DO 150 I=1,LS2                                                            
         TH = FLOAT(I-1)*DT                                                     
         CDT = DCOS(TH+TH)                                                      
         SDT = DSIN(TH+TH)                                                      
      IF(I.EQ.1) SDT=0.                                                         
         CT = DCOS(TH)                                                          
         ST = -DSIN(TH)                                                         
      IF(I.EQ.1) ST=0.                                                          
         SUM = 0.                                                               
         DO 140 K=1,KDO                                                         
            CTH = CDT*CT-SDT*ST                                                 
            ST = SDT*CT+CDT*ST                                                  
            CT = CTH                                                            
            SUM = SUM+CP(K)*CT                                                  
C                                                                               
C     SUM=SUM+CP(K)*COS((K+K-1)*TH)                                             
C                                                                               
  140    CONTINUE                                                               
         P(I) = SUM                                                             
  150 CONTINUE                                                                  
      GO TO 190                                                                 
  160 DO 180 I=1,LS2                                                            
         TH = FLOAT(I-1)*DT                                                     
         CDT = DCOS(TH+TH)                                                      
         SDT = DSIN(TH+TH)                                                      
      IF(I.EQ.1) SDT=0.                                                         
         CT = DCOS(TH)                                                          
         ST = -DSIN(TH)                                                         
      IF(I.EQ.1) ST=0.                                                          
         SUM = 0.                                                               
         DO 170 K=1,KDO                                                         
            CTH = CDT*CT-SDT*ST                                                 
            ST = SDT*CT+CDT*ST                                                  
            CT = CTH                                                            
            SUM = SUM+CP(K)*ST                                                  
C                                                                               
C     SUM=SUM+CP(K)*SIN((K+K-1)*TH)                                             
C                                                                               
  170    CONTINUE                                                               
         P(I) = SUM                                                             
  180 CONTINUE                                                                  
  190 IF (MOD(N-M,2)) 220,200,220                                               
  200 DO 210 I=1,LS2                                                            
         LMI = L-I                                                              
         P(LMI+1) = P(I)                                                        
  210 CONTINUE                                                                  
      GO TO 240                                                                 
  220 DO 230 I=1,LS2                                                            
         LMI = L-I                                                              
         P(LMI+1) = -P(I)                                                       
  230 CONTINUE                                                                  
  240 RETURN                                                                    
      END                                                                       
      SUBROUTINE DLFPT (N,M,TH,CP,P)                                            
      DOUBLE PRECISIONCP         ,P          ,CDT        ,SDT        ,          
     1                CT         ,ST         ,SUM        ,CTH        ,TH        
      DIMENSION       CP(1)                                                     
C                                                                               
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR         
C                                                                               
      CALL Q8QST4 (4HXLIB,6HALFPAC,5HDLFPT,10HVERSION  1)                       
      IF (N)  10, 10, 30                                                        
   10 IF (M)  20, 20, 30                                                        
   20 P = DSQRT(5.D-1)                                                          
      GO TO 140                                                                 
   30 NP1 = N+1                                                                 
      NMOD = MOD(N,2)                                                           
      MMOD = MOD(M,2)                                                           
      IF (NMOD)  40, 40, 90                                                     
   40 IF (MMOD)  50, 50, 70                                                     
   50 KDO = N/2+1                                                               
      CDT = DCOS(TH+TH)                                                         
      SDT = DSIN(TH+TH)                                                         
      IF(TH.EQ.0.) SDT=0.                                                       
      CT = 1.                                                                   
      ST = 0.                                                                   
      SUM = .5*CP(1)                                                            
      DO  60 KP1=2,KDO                                                          
         CTH = CDT*CT-SDT*ST                                                    
         ST = SDT*CT+CDT*ST                                                     
         CT = CTH                                                               
         SUM = SUM+CP(KP1)*CT                                                   
   60 CONTINUE                                                                  
      P = SUM                                                                   
      GO TO 140                                                                 
   70 KDO = N/2                                                                 
      CDT = DCOS(TH+TH)                                                         
      SDT = DSIN(TH+TH)                                                         
      IF(TH.EQ.0.) SDT=0.                                                       
      CT = 1.                                                                   
      ST = 0.                                                                   
      SUM = 0.                                                                  
      DO  80 K=1,KDO                                                            
         CTH = CDT*CT-SDT*ST                                                    
         ST = SDT*CT+CDT*ST                                                     
         CT = CTH                                                               
         SUM = SUM+CP(K)*ST                                                     
   80 CONTINUE                                                                  
      P = SUM                                                                   
      GO TO 140                                                                 
   90 KDO = (N+1)/2                                                             
      IF (MMOD) 100,100,120                                                     
  100 CDT = DCOS(TH+TH)                                                         
      SDT = DSIN(TH+TH)                                                         
      IF(TH.EQ.0.) SDT=0.                                                       
      CT = DCOS(TH)                                                             
      ST = -DSIN(TH)                                                            
      IF(TH.EQ.0.) ST=0.                                                        
      SUM = 0.                                                                  
      DO 110 K=1,KDO                                                            
         CTH = CDT*CT-SDT*ST                                                    
         ST = SDT*CT+CDT*ST                                                     
         CT = CTH                                                               
         SUM = SUM+CP(K)*CT                                                     
  110 CONTINUE                                                                  
      P = SUM                                                                   
      GO TO 140                                                                 
  120 CDT = DCOS(TH+TH)                                                         
      SDT = DSIN(TH+TH)                                                         
      IF(TH.EQ.0.) SDT=0.                                                       
      CT = DCOS(TH)                                                             
      ST = -DSIN(TH)                                                            
      IF(TH.EQ.0.) ST=0.                                                        
      SUM = 0.                                                                  
      DO 130 K=1,KDO                                                            
         CTH = CDT*CT-SDT*ST                                                    
         ST = SDT*CT+CDT*ST                                                     
         CT = CTH                                                               
         SUM = SUM+CP(K)*ST                                                     
  130 CONTINUE                                                                  
      P = SUM                                                                   
  140 RETURN                                                                    
      END                                                                       
      SUBROUTINE DLFMA (N,L,I,P,W)                                              
      DOUBLE PRECISIONP          ,W                                             
      DIMENSION       P(1)       ,W(1)                                          
C                                                                               
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR         
C                                                                               
      CALL Q8QST4 (4HXLIB,6HALFPAC,5HDLFMA,10HVERSION  1)                       
      IW1 = L+1                                                                 
      IW2 = IW1+L                                                               
      IW3 = IW2+L                                                               
      CALL DLFMA1 (N,L,I,P,W,W(IW1),W(IW2),W(IW3))                              
      RETURN                                                                    
      END                                                                       
      SUBROUTINE DLFMA1 (N,L,I,P,PZ,P1,A,B)                                     
      DOUBLE PRECISIONP          ,PZ         ,P1       ,A            ,          
     1                B          ,PI         ,DT         ,FNPM       ,          
     2                FNMM       ,FMPM       ,THETA      ,COST       ,          
     3                SINT       ,R                                             
      DIMENSION       P(1)       ,PZ(1)      ,P1(1)      ,A(1)       ,          
     1                B(1)                                                      
      DATA NLAST /-1/, LLAST /-1/                                               
      IF (N)  10, 10, 20                                                        
   10 P(1) = 1./DSQRT(2.D0)                                                     
      GO TO 100                                                                 
   20 IF (N-NLAST)  40, 30, 40                                                  
   30 IF (L-LLAST)  40, 60, 40                                                  
   40 NLAST = N                                                                 
      LLAST = L                                                                 
      CALL DLFK (N,0,P,DUMW)                                                    
      CALL DLFP (N,0,L,P,PZ)                                                    
      CALL DLFK (N,1,P,DUMW)                                                    
      CALL DLFP (N,1,L,P,P1)                                                    
      NM1 = N-1                                                                 
      NP1 = N+1                                                                 
      PI = 4.*DATAN(1.D0)                                                       
      DT = PI/FLOAT(L-1)                                                        
      FNPM = N                                                                  
      FNMM = N+1                                                                
      FMPM = 0.                                                                 
      DO  50 M=1,N                                                              
         FNPM = FNPM+1.                                                         
         FNMM = FNMM-1.                                                         
         FMPM = FMPM-2.                                                         
         B(M) = DSQRT(FNMM*FNPM)                                                
   50 CONTINUE                                                                  
   60 THETA = FLOAT(I-1)*DT                                                     
      COST = DCOS(THETA)                                                        
      SINT = DSIN(THETA)                                                        
      IF(I.EQ.1) SDT=0.                                                         
      FMPM = 0.                                                                 
      DO  70 M=1,N                                                              
         FMPM = FMPM-2.                                                         
         A(M) = B(M)*SINT                                                       
         P(M+1) = FMPM*COST                                                     
   70 CONTINUE                                                                  
      IF (DABS(PZ(I))-DABS(P1(I)))  90, 90, 80                                  
   80 R = -A(1)*PZ(I)                                                           
      CALL DTRIH (N,A,P(2),A(2),R)                                              
      P(1) = PZ(I)                                                              
      GO TO 100                                                                 
   90 R = -A(2)*P1(I)                                                           
      CALL DTRIH (NM1,A(2),P(3),A(3),R)                                         
      P(1) = PZ(I)                                                              
      P(2) = P1(I)                                                              
  100 RETURN                                                                    
      END                                                                       
      SUBROUTINE DLFMB (N,L,I,P,W)                                              
      DOUBLE PRECISIONP          ,W                                             
      DIMENSION       P(1)       ,W(1)                                          
C                                                                               
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR         
C                                                                               
      CALL Q8QST4 (4HXLIB,6HALFPAC,5HDLFMB,10HVERSION  1)                       
      IW1 = L+1                                                                 
      IW2 = IW1+L                                                               
      IW3 = IW2+L                                                               
      CALL DLFMB1 (L,N,I,P,W,W,W(IW1),W(IW2),W(IW2),W(IW3))                     
      RETURN                                                                    
      END                                                                       
      SUBROUTINE DLFMB1 (L,N,I,P,PMZ,DPMZ,P1,PM1,DPM1,P2)                       
      DOUBLE PRECISIONDPMZ       ,DPM1       ,DPMZ1      ,DPMZ2      ,          
     1                DPM11      ,DPM12      ,DTH        ,DPI        ,          
     2                THETA      ,COST       ,SINT       ,P          ,          
     3                P1         ,P2         ,PMZ        ,PM1        ,          
     4                FNPM       ,FNMM       ,FMPM                              
      DIMENSION       P1(1)      ,P2(1)      ,P(1)       ,PMZ(1)     ,          
     1                PM1(1)     ,DPMZ(1)    ,DPM1(1)                           
      DATA ILAST/-1/, LLAST/-1/                                                 
      IF (I-ILAST)  20, 10, 20                                                  
   10 IF (L-LLAST)  20, 50, 20                                                  
   20 ILAST = I                                                                 
      LLAST = L                                                                 
      LM1 = L-1                                                                 
      LM2 = L-2                                                                 
      DPI = 4.*DATAN(1.D0)                                                      
      DTH = FLOAT(I-1)*DPI/FLOAT(LM1)                                           
      THETA = DTH                                                               
      COST = DCOS(THETA)                                                        
      SINT = DSIN(THETA)                                                        
      IF(I.EQ.1) SDT=0.                                                         
      CALL DLFK (LM1,0,P1,DUMW)                                                 
      CALL DLFPT (LM1,0,DTH,P1,DPMZ1)                                           
      CALL DLFK (LM1,1,P1,DUMW)                                                 
      CALL DLFPT (LM1,1,DTH,P1,DPMZ2)                                           
      CALL DLFK (LM2,0,P1,DUMW)                                                 
      CALL DLFPT (LM2,0,DTH,P1,DPM11)                                           
      CALL DLFK (LM2,1,P1,DUMW)                                                 
      CALL DLFPT (LM2,1,DTH,P1,DPM12)                                           
      DPMZ(1) = DPMZ1                                                           
      DPMZ(2) = DPM11                                                           
      CALL DRECN (L,0,DTH,DPMZ,PM1)                                             
      DO  30 NH=1,L                                                             
         P(NH) = DPMZ(NH)                                                       
   30 CONTINUE                                                                  
      DPM1(1) = DPMZ2                                                           
      DPM1(2) = DPM12                                                           
      CALL DRECN (L,1,DTH,DPM1,PMZ)                                             
      DO  40 NH=1,L                                                             
         PMZ(NH) = P(NH)                                                        
         PM1(NH) = DPM1(NH)                                                     
   40 CONTINUE                                                                  
   50 IF (N) 130, 60, 70                                                        
   60 P(1) = PMZ(1)                                                             
      GO TO 130                                                                 
   70 IF (N-1) 130, 80, 90                                                      
   80 P(1) = PMZ(2)                                                             
      P(2) = PM1(2)                                                             
      GO TO 130                                                                 
   90 FNPM = N                                                                  
      FNMM = N+1                                                                
      FMPM = 0.                                                                 
      DO 100 M=1,N                                                              
         FNPM = FNPM+1.                                                         
         FNMM = FNMM-1.                                                         
         FMPM = FMPM-2.                                                         
         P1(M) = SINT*DSQRT(FNMM*FNPM)                                          
         P2(M) = FMPM*COST                                                      
  100 CONTINUE                                                                  
      IF (DABS(PMZ(N+1))-DABS(PM1(N+1))) 120,120,110                            
  110 P(2) = -P1(1)*PMZ(N+1)                                                    
      CALL DTRI (N,P1,P2,P1(2),P(2))                                            
      P(1) = PMZ(N+1)                                                           
      GO TO 130                                                                 
  120 P(3) = -P1(2)*PM1(N+1)                                                    
      CALL DTRI (N-1,P1(2),P2(2),P1(3),P(3))                                    
      P(1) = PMZ(N+1)                                                           
      P(2) = PM1(N+1)                                                           
  130 RETURN                                                                    
      END                                                                       
      SUBROUTINE DLFNA (M,L,I,P,W)                                              
C                                                                               
      DOUBLE PRECISIONP          ,W                                             
      DIMENSION       P(1)       ,W(1)                                          
C                                                                               
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR         
C                                                                               
      CALL Q8QST4 (4HXLIB,6HALFPAC,5HDLFNA,10HVERSION  1)                       
      IW1 = L+1                                                                 
      IW2 = IW1+L                                                               
      IW3 = IW2+L                                                               
      CALL DLFNA1 (M,L,I,P,W,W(IW1),W(IW2),W(IW3))                              
      RETURN                                                                    
      END                                                                       
      SUBROUTINE DLFNA1 (M,L,I,P,PZ,P1,A,B)                                     
      DOUBLE PRECISIONP          ,PZ         ,P1         ,A          ,          
     1                B          ,FNMM       ,FNPN       ,FNPM       ,          
     2                THETA      ,COST       ,R          ,PI         ,          
     3                PHOLD                                                     
      DIMENSION       P(1)       ,PZ(1)      ,P1(1)      ,A(1)       ,          
     1                B(1)                                                      
      DATA MLAST /-1/, LLAST /-1/                                               
      DO  10 N=1,L                                                              
         P(N) = 0.                                                              
   10 CONTINUE                                                                  
      IF (L-LLAST)  30, 20, 30                                                  
   20 IF (M-MLAST)  30, 70, 30                                                  
   30 MLAST = M                                                                 
      LLAST = L                                                                 
      LM1 = L-1                                                                 
      LM2 = L-2                                                                 
      LMM = L-M                                                                 
      PI = 4.*DATAN(1.D0)                                                       
      CALL DLFK (LM2,M,B,PZ)                                                    
      CALL DLFP (LM2,M,L,B,P)                                                   
      CALL DLFK (LM1,M,B,PZ)                                                    
      CALL DLFP (LM1,M,L,B,PZ)                                                  
      DO  40 I=1,L                                                              
         P1(I) = P(I)                                                           
   40 CONTINUE                                                                  
      K = 0                                                                     
      N = L-1                                                                   
      FNMM = N-M+1                                                              
      FNPN = N+N+1                                                              
      FNPM = N+M+1                                                              
   50 N = N-1                                                                   
      IF (N-M)  70, 60, 60                                                      
   60 FNMM = FNMM-1.                                                            
      FNPN = FNPN-2.                                                            
      FNPM = FNPM-1.                                                            
      K = K+1                                                                   
      B(K) = DSQRT(FNMM*FNPM/(FNPN*(FNPN+2.)))                                  
      GO TO  50                                                                 
   70 P(L) = PZ(I)                                                              
      IF (LMM-2) 140, 80, 90                                                    
   80 P(L-1) = P1(I)                                                            
      GO TO 140                                                                 
   90 THETA = FLOAT(I-1)*PI/FLOAT(LM1)                                          
      COST = DCOS(THETA)                                                        
      DO 100 K=1,LMM                                                            
         P(K) = -COST                                                           
         A(K) = B(K)                                                            
  100 CONTINUE                                                                  
      IF (DABS(PZ(I)) .LT. DABS(P1(I))) GO TO 110                               
      P(1) = PZ(I)                                                              
      R = -A(1)*P(1)                                                            
      CALL DTRIH (LMM-1,A,P(2),A(2),R)                                          
      GO TO 120                                                                 
  110 P(1) = PZ(I)                                                              
      P(2) = P1(I)                                                              
      R = -A(2)*P(2)                                                            
      CALL DTRIH (LMM-2,A(2),P(3),A(3),R)                                       
  120 NDO = (L+1)/2                                                             
      DO 130 N=1,NDO                                                            
         N1 = L-N                                                               
         PHOLD = P(N1+1)                                                        
         P(N1+1) = P(N)                                                         
         P(N) = PHOLD                                                           
  130 CONTINUE                                                                  
  140 RETURN                                                                    
      END                                                                       
      SUBROUTINE DLFNB (M,L,I,P,W)                                              
      DOUBLE PRECISIONP          ,W                                             
      DIMENSION       P(1)       ,W(1)                                          
C                                                                               
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR         
C                                                                               
      CALL Q8QST4 (4HXLIB,6HALFPAC,5HDLFNB,10HVERSION  1)                       
      IW1 = L+1                                                                 
      IW2 = IW1+L                                                               
      IW3 = IW2+L                                                               
      CALL DLFNB1 (M,L,I,P,W,W(IW1),W(IW2),W(IW3))                              
      RETURN                                                                    
      END                                                                       
      SUBROUTINE DLFNB1 (M,L,I,P,PZ,P1,A,B)                                     
      DOUBLE PRECISIONPZ         ,P1         ,A          ,B                     
     1                P          ,PI         ,THETA      ,COST       ,          
     2                PHOLD1     ,PHOLD2     ,FNMM       ,FNPM       ,          
     3                FNPN       ,PHOLD                                         
      DIMENSION       P(1)       ,PZ(1)      ,P1(1)      ,A(1)       ,          
     1                B(1)                                                      
      DATA ILAST /-1/, LLAST /-1/                                               
      DO  10 N=1,L                                                              
         P(N) = 0.                                                              
   10 CONTINUE                                                                  
      IF (L-LLAST)  30, 20, 30                                                  
   20 IF (I-ILAST)  30, 50, 30                                                  
   30 ILAST = I                                                                 
      LLAST = L                                                                 
      LM1 = L-1                                                                 
      LM2 = L-2                                                                 
      PI = 4.*DATAN(1.D0)                                                       
      THETA = FLOAT(I-1)*PI/FLOAT(LM1)                                          
      COST = DCOS(THETA)                                                        
      LS = 0                                                                    
      DO  40 N=LM2,LM1                                                          
         LS = LS+1                                                              
         CALL DLFK (N,0,A,DUMW)                                                 
         CALL DLFPT (N,0,THETA,A,PZ(LS))                                        
         CALL DLFK (N,1,A,DUMW)                                                 
         CALL DLFPT (N,1,THETA,A,P1(LS))                                        
   40 CONTINUE                                                                  
      PHOLD1 = PZ(1)                                                            
      PHOLD2 = P1(2)                                                            
      PZ(1) = PZ(2)                                                             
      P1(2) = P1(1)                                                             
      P1(1) = PHOLD1                                                            
      PZ(2) = PHOLD2                                                            
      CALL DRECM (L,LM2,I,P1,A,B)                                               
      CALL DRECM (L,LM1,I,PZ,A,B)                                               
   50 K = 0                                                                     
      LMM = L-M                                                                 
      IF (LMM-1)  60, 60, 70                                                    
   60 P(L) = PZ(M+1)                                                            
      GO TO 160                                                                 
   70 N = L-1                                                                   
      FNMM = N-M+1                                                              
      FNPN = N+N+1                                                              
      FNPM = N+M+1                                                              
   80 N = N-1                                                                   
      IF (N-M) 100, 90, 90                                                      
   90 FNMM = FNMM-1.                                                            
      FNPN = FNPN-2.                                                            
      FNPM = FNPM-1.                                                            
      K = K+1                                                                   
      A(K) = DSQRT(FNMM*FNPM/(FNPN*(FNPN+2.)))                                  
      B(K) = -COST                                                              
      GO TO  80                                                                 
  100 IF (DABS(PZ(M+1)) .LT. DABS(P1(M+1))) GO TO 120                           
      P(1) = PZ(M+1)                                                            
      P(2) = -A(1)*P(1)                                                         
      IF (LMM-1) 140,140,110                                                    
  110 CALL DTRI (LMM-1,A,B,A(2),P(2))                                           
      GO TO 140                                                                 
  120 P(1) = PZ(M+1)                                                            
      P(2) = P1(M+1)                                                            
      P(3) = -A(2)*P(2)                                                         
      IF (LMM-2) 140,140,130                                                    
  130 CALL DTRI (LMM-2,A(2),B(2),A(3),P(3))                                     
  140 NDO = (L+1)/2                                                             
      DO 150 N=1,NDO                                                            
         N1 = L-N                                                               
         PHOLD = P(N1+1)                                                        
         P(N1+1) = P(N)                                                         
         P(N) = PHOLD                                                           
  150 CONTINUE                                                                  
  160 RETURN                                                                    
      END                                                                       
      SUBROUTINE DRECM (L,N,I,P,A,B)                                            
      DOUBLE PRECISIONP          ,A          ,B          ,PI         ,          
     1                DT         ,FNPM       ,FNMM       ,FMPM       ,          
     2                THETA      ,COST       ,SINT                              
      DIMENSION       P(1)       ,A(1)       ,B(1)                              
      PI = 4.*DATAN(1.D0)                                                       
      DT = PI/(L-1)                                                             
      FNPM = N                                                                  
      FNMM = N+1                                                                
      FMPM = 0.                                                                 
      THETA = FLOAT(I-1)*DT                                                     
      COST = DCOS(THETA)                                                        
      SINT = DSIN(THETA)                                                        
      IF(I.EQ.1) SDT=0.                                                         
      DO  10 M=1,N                                                              
         FNPM = FNPM+1.                                                         
         FNMM = FNMM-1.                                                         
         FMPM = FMPM-2.                                                         
         A(M) = SINT*DSQRT(FNMM*FNPM)                                           
         B(M) = FMPM*COST                                                       
   10 CONTINUE                                                                  
      IF (DABS(P(1))-DABS(P(2)))  30, 30, 20                                    
   20 P(2) = -A(1)*P(1)                                                         
      CALL DTRI (N,A,B,A(2),P(2))                                               
      GO TO  40                                                                 
   30 P(3) = -A(2)*P(2)                                                         
      NM1 = N-1                                                                 
      CALL DTRI (NM1,A(2),B(2),A(3),P(3))                                       
   40 RETURN                                                                    
      END                                                                       
      SUBROUTINE DTRI (N,A,B,C,E)                                               
      DOUBLE PRECISIONA          ,B          ,C          ,E          ,          
     1                R          ,QIH        ,BIH        ,RATIO      ,          
     2                BIH1       ,Q2         ,RIH                               
      DIMENSION       A(1)       ,B(1)       ,C(1)       ,E(1)                  
      R = E(1)                                                                  
      IF (N) 120,120, 10                                                        
   10 IF (N-2)  20, 30, 40                                                      
   20 E(1) = R/B(1)                                                             
      RETURN                                                                    
   30 QIH = A(2)                                                                
      BIH = B(2)                                                                
      GO TO  70                                                                 
   40 QIH = A(N)                                                                
      BIH = B(N)                                                                
      DO  60 IDO=3,N                                                            
         I = N-IDO+2                                                            
         IF (DABS(BIH) .LT. DABS(C(I))) GO TO  50                               
         RATIO = C(I)/BIH                                                       
         C(I) = 0.                                                              
         B(I+1) = QIH/BIH                                                       
C                                                                               
C     P(I+1)=0.                                                                 
C     Q(I+1)=QIH/BIH                                                            
C                                                                               
         BIH = B(I)-RATIO*QIH                                                   
         QIH = A(I)                                                             
         GO TO  60                                                              
   50    B(I+1) = B(I)/C(I)                                                     
         C(I) = A(I)/C(I)                                                       
         BIH1 = QIH-BIH*B(I+1)                                                  
         QIH = -BIH*C(I)                                                        
         BIH = BIH1                                                             
C                                                                               
C     P(I+1)=A(I)/C(I)                                                          
C     Q(I+1)=B(I)/C(I)                                                          
C                                                                               
   60 CONTINUE                                                                  
   70 IF (DABS(BIH) .LT. DABS(C(1))) GO TO  80                                  
      Q2 = QIH/BIH                                                              
C                                                                               
C     Q(2)=QIH/BIH                                                              
C                                                                               
      BIH = B(1)-C(1)/BIH*QIH                                                   
      E(1) = R/BIH                                                              
      E(2) = -Q2*E(1)                                                           
C                                                                               
C     E(2)=-Q(2)*E(1)                                                           
C                                                                               
      GO TO  90                                                                 
   80 RATIO = BIH/C(1)                                                          
      BIH = QIH-RATIO*B(1)                                                      
      RIH = -RATIO*R                                                            
      E(1) = RIH/BIH                                                            
      E(2) = (R-B(1)*E(1))/C(1)                                                 
   90 IF (N-3) 120,100,100                                                      
  100 DO 110 I=3,N                                                              
         E(I) = -B(I)*E(I-1)-C(I-1)*E(I-2)                                      
C                                                                               
C     E(I)=-Q(I)*E(I-1)-P(I)*E(I-2)                                             
C                                                                               
  110 CONTINUE                                                                  
  120 RETURN                                                                    
      END                                                                       
      SUBROUTINE DTRIH (N,A,B,C,R)                                              
      DOUBLE PRECISIONA          ,B          ,C          ,B1         ,          
     1                R          ,QIH        ,BIH        ,RATIO      ,          
     2                BIH1       ,Q2         ,RIH                               
      DIMENSION       A(1)       ,B(1)       ,C(1)                              
      IF (N) 120,120, 10                                                        
   10 IF (N-2)  20, 30, 40                                                      
   20 B(1) = R/B(1)                                                             
      RETURN                                                                    
   30 QIH = A(2)                                                                
      BIH = B(2)                                                                
      GO TO  70                                                                 
   40 QIH = A(N)                                                                
      BIH = B(N)                                                                
      DO  60 IDO=3,N                                                            
         I = N-IDO+2                                                            
         IF (DABS(BIH) .LT. DABS(C(I))) GO TO  50                               
         RATIO = C(I)/BIH                                                       
         C(I) = 0.                                                              
         B(I+1) = QIH/BIH                                                       
         BIH = B(I)-RATIO*QIH                                                   
         QIH = A(I)                                                             
         GO TO  60                                                              
   50    B(I+1) = B(I)/C(I)                                                     
         C(I) = A(I)/C(I)                                                       
         BIH1 = QIH-BIH*B(I+1)                                                  
         QIH = -BIH*C(I)                                                        
         BIH = BIH1                                                             
   60 CONTINUE                                                                  
   70 IF (DABS(BIH) .LT. DABS(C(1))) GO TO  80                                  
      Q2 = QIH/BIH                                                              
      BIH = B(1)-C(1)/BIH*QIH                                                   
      B(1) = R/BIH                                                              
      B(2) = -Q2*B(1)                                                           
      GO TO  90                                                                 
   80 RATIO = BIH/C(1)                                                          
      BIH = QIH-RATIO*B(1)                                                      
      RIH = -RATIO*R                                                            
      B1 = RIH/BIH                                                              
      B(2) = (R-B(1)*B1)/C(1)                                                   
      B(1) = B1                                                                 
   90 IF (N-3) 120,100,100                                                      
  100 DO 110 I=3,N                                                              
         B(I) = -B(I)*B(I-1)-C(I-1)*B(I-2)                                      
  110 CONTINUE                                                                  
  120 RETURN                                                                    
      END                                                                       

c /\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\
c \/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/
c;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
c Start of new code additions - BR
c;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	subroutine fixit(iarry, xarry)
	dimension iarry(72,36), xarry(73,37)
c Convert the integer array 'iarry' into a float and tack on an extra bit
c at the end.
	do j=1,36 
	  jj= 36-(j-1)
	  do i= 1,72
	    xarry(i,jj)= float(iarry(i,j))
	  enddo
	enddo

c Add the extra value.
	xarry(73,37)= xarry(72,36)
	return
	end
c;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
