C**************************************************************************** 
      SUBROUTINE RAYPLT
C**************************************************************************** 
C     MAIN PLOTTING PROGRAM; INITIALIZES, READS INPUT, PLOTS
C     PROJECTIONS OF RAYS ON A VERTICAL OR HORIZONTAL PLANE.
C        ABS(PLT)=1. PLOTS PROJECTION OF RAYPATH ON VERTICAL PLANE
C                    RECTANGULAR EXPANSION BY FACTOR 'PFACTR'
C               =2. PLOTS PROJECTION OF RAYPATH ON GROUND
C               =3. VERICAL PROJECTION USING RADIAL EXPANSION BY FACTOR
C                    'PFACTR'
C     COMMON DECK "FILEC" INSERTED HERE
 
      implicit double precision (a-h,o-z)
      CHARACTER ID*80
      COMMON /FILEC/NPLTDP
      COMMON /PLT/ XL,XR,YB,YT,PRESET
      COMMON/PLT/RMIN,RMAX,ALPHA,APLT
C
C     COMMON DECK "CONST" INSERTED HERE
      COMMON/PCONST/CREF,RGAS,GAMMA
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10
C     COMMON DECK "FLAG" INSERTED HERE
      LOGICAL NEWWR,NEWWP,NEWTRC,PENET
      COMMON /FLG/ NTYP,NEWWR,NEWWP,NEWTRC,PENET,LINES,IHOP,HPUNCH
      COMMON/FLGP/NSET
C     COMMON DECK "RKAM" INSERTED HERE
      REAL*8 KR,KTH,KPH
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)
      COMMON/WW/MAXW,W(1000)
      common/ident/id
      REAL*8 MAXSTP,MAXERR,INTYP,LLAT,LLON
      EQUIVALENCE (EARTHR,W(1)),(RAY,W(2)),(XMTRH,W(3)),(TLAT,W(4)),
     1 (TLON,W(5)),(OW,W(6)),(FBEG,W(7)),(FEND,W(8)),(FSTEP,W(9)),
     2 (AZ1,W(10)),(AZBEG,W(11)),(AZEND,W(12)),(AZSTEP,W(13)),
     3 (BETA,W(14)),(ELBEG,W(15)),(ELEND,W(16)),(ELSTEP,W(17)),
     8   (RUNSUP,W(18)),(RCVRH,W(20)),
     4 (ONLY,W(21)),(HOP,W(22)),(MAXSTP,W(23)),(PLAT,W(24)),(PLON,W(25))
     5,(HMAX,W(26)),(RAYFNC,W(29)),(EXTINC,W(33)),
     6 (HMIN,W(27)),(RGMAX,W(28)),
     8 (INTYP,W(41)),(MAXERR,W(42)),(ERATIO,W(43)),
     6 (STEP1,W(44)),(STPMAX,W(45)),(STPMIN,W(46)),(FACTR,W(47)),
     7 (SKIP,W(71)),(RAYSET,W(72)),(PRTSRP,W(74)),(HITLET,W(75))
     9 ,(BINRAY,W(76)),(PAGLN,W(77)),(PLT,W(81)),(PFACTR,W(82)),
     1 (LLAT,W(83)),(LLON,W(84)),(RLAT,W(85)),(RLON,W(86))
     2,(TIC,W(87)),(HB,W(88)),(HT,W(89)),(TICV,W(96))
      REAL*8 MMODEL,MFORM,MID
C
C     WIND        100-124
      EQUIVALENCE (W(100),UMODEL),(W(101),UFORM),(W(102),UID)
C
C     DELTA WIND        125-149
      EQUIVALENCE (W(125),DUMODEL),(W(126),DUFORM),(W(127),DUID)
C
C     SOUND SPEED 150-174
      EQUIVALENCE (W(150),CMODEL),(W(151),CFORM),(W(152),CID)
      EQUIVALENCE (W(153),REFC)
C
C     DELTA SOUND SPEED 175-199
      EQUIVALENCE (W(175),DCMODEL),(W(176),DCFORM),(W(177),DCID)
C
C     TEMPERATURE    200-224
      EQUIVALENCE (W(200),TMODEL),(W(201),TFORM),(W(202),TID)
C
C     DELTA TEMPERATURE    225-249
      EQUIVALENCE (W(225),DTMODEL),(W(226),DTFORM),(W(227),DTID)
C
C     MOLECULAR   250-274
      EQUIVALENCE (W(250),MMODEL),(W(251),MFORM),(W(252),MID)
C
C     RECEIVER HEIGHT   275-299
      EQUIVALENCE (W(275),RMODEL),(W(276),RFORM),(W(277),RID)
C
C     TOPOGRAPHY   300-324
      EQUIVALENCE (W(300),GMODEL),(W(301),GFORM),(W(302),GID)
C
C     DELTA TOPOGRAPHY   325-349
      EQUIVALENCE (W(325),GUMODEL),(W(326),GUFORM),(W(327),GUID)
C
C     ATMOSPHERIC SURFACE TOPOGRAPHY   350-374
      EQUIVALENCE (W(350),SMODEL),(W(351),SFORM),(W(352),SID)
C
C     ATMOSPHERIC SURFACE TOPOGRAPHY   375-399
      EQUIVALENCE (W(375),SUMODEL),(W(376),SUFORM),(W(377),SUID)
C     PLOT ENHANCEMENTS CONTROL PARAMETERS
C
      EQUIVALENCE (W(490),XFQMDL),(W(491),YFQMDL)
C     ABSORPTION        500-524
      EQUIVALENCE (W(500),AMODEL),(W(501),AFORM),(W(502),AID)
C
C     DELTA ABSORPTION        525-549
      EQUIVALENCE (W(525),DAMODEL),(W(526),DAFORM),(W(527),DAID)
C
C     PRESSURE        550-574
      EQUIVALENCE (W(550),PMODEL),(W(551),PFORM),(W(552),PID)
C
C     DELTA PRESSURE        575-599
      EQUIVALENCE (W(575),DPMODEL),(W(576),DPFORM),(W(577),DPID)
C
      REAL*8 LTIC
C%      INTEGER IWDMP(411)
C%      EQUIVALENCE(IWDMP,ID)
C
      EXTERNAL PLOT,PLTLB,PLTLBN,PLTHLB
C
      DATA NWDSRK/4/
C
      IF (.NOT.NEWWR) GO TO 8
      APLT=dABS(PLT)
C
C          NEW W ARRAY -- REINITIALIZE
      NEWWR=.FALSE.
      PRESET=1.d0
C     IF NO ACTIVE PLOTTING,WE ARE STILL DUMPING DATA
      IF(PLT.EQ.0.d0) GO TO 5
C     INITIALIZE ANNOTATION MODEL
C  CONVERT COORDINATES OF VERTICAL PLANE FROM GEOGRAPHIC TO GEOMAGNETIC
      SW=dSIN (PLAT)
      CW=dSIN (PID2-PLAT)
      SLM=dSIN (LLAT)
      CLM=dSIN (PID2-LLAT)
      SRM=dSIN (RLAT)
      CRM=dSIN (PID2-RLAT)
      CDPHI=dSIN (PID2-(LLON-PLON))
      PHL=dATAN2(dSIN (LLON-PLON)*CLM,CDPHI*SW*CLM-CW*SLM)
      CTHL=CDPHI*CW*CLM+SW*SLM
      STHL=dSIN (dACOS (CTHL))
      CDPHI=dSIN (PID2-(RLON-PLON))
      PHR=dATAN2(dSIN (RLON-PLON)*CRM,CDPHI*SW*CRM-CW*SRM)
      CTHR=CDPHI*CW*CRM+SW*SRM
      STHR=dSIN (dACOS (CTHR))
      CLR=CTHL*CTHR+STHL*STHR*dSIN (PID2-(PHL-PHR))
      ALPHA=.5d0*dACOS (CLR)
      SLR=dSQRT (1.d0-CLR**2)
C
      IF (APLT.EQ.2.d0) GO TO 3
C
C     VERTICAL PROJECTIONS ONLY
      T=HB
      HB=DMIN1(T,HT)
      HT=DMAX1(T,HT)
      RMIN=EARTHR+HB
      YT1=YT
      RMAX=EARTHR+HT
C
      IF(APLT.NE.4.d0) GO TO 100
C
C     SCALE FOR CARTESIAN PROJECTION(=4)
      XR=ALPHA
      XR1=XR
      XL=-XR
      XL1=XL
      YB=RMIN
      YT=RMAX
      GO TO 5
C
C     PROJECTIONS 1 AND 3
100   XR1=RMIN*dSIN (ALPHA)
      XL1=-XR1
      YB=RMIN*dSIN (PID2-ALPHA)
      IF(APLT.EQ.3.d0) RMAX=RMIN+(RMAX-RMIN)*PFACTR
      XR=DMAX1(RMAX*DSIN(ALPHA),(RMAX-YB)/2.D0)
      XL=-XR
      YT1=RMAX*(YB/RMIN)
      YT=2.d0*XR
      IF(APLT.EQ.1.d0) YT=YT/PFACTR
      YT=YT+YB
      GO TO 5
C
C     HORIZONTAL PROJECTION(=2)
3     ALPH1=dATAN2(STHR*dSIN (PHR-PHL),(CTHR-CTHL*CLR)/STHL)
      XL=0.d0
      XL1=0.d0
      XR=EARTHR*2.d0*ALPHA
      XR1=XR
C USE 90% OF X-RANGE FOR Y-RANGE
      RMAX=0.5d0*(0.9d0*XR)/PFACTR
      YT=RMAX
      RMIN=-RMAX
      YB=RMIN
C
5     IF(NPLTDP.LE.0) GO TO 8
      WRITE(NPLTDP) 1,NWDSRK,1,411
      WRITE(NPLTDP) id
C
8     NEW=0
      IF (NEWTRC) NEW=1
      NEWTRC=.FALSE.
C
      IF(NPLTDP.LE.0) GO TO 88
      IF(NEW.NE.1) GO TO 84
      WRITE(NPLTDP) 1,NWDSRK,17,25
      WRITE(NPLTDP) id
C
84    WRITE(NPLTDP) 3-NEW,R,TH,PH
C
88    IF(PLT.EQ.0.d0) RETURN
C
      STH=dSIN (TH)
      CTH=dSIN (PID2-TH)
      CR=CTHR*CTH+STHR*STH*dSIN (PID2-(PHR-PH))
      CL=CTHL*CTH+STHL*STH*dSIN (PID2-(PHL-PH))
      CEA=dATAN2(CR-CL*CLR,CL*SLR)
C
      IF(APLT.NE.4.d0) GO TO 150
      CALL PLOT(CEA-ALPHA,R,NEW)
      RETURN
C
150   IF(APLT.EQ.2.d0) GO TO 10
      RX=R
      IF(APLT.EQ.3.d0) RX=RMIN+(R-RMIN)*PFACTR
      CALL PLOT(CEA-ALPHA,RX,NEW)
      RETURN
C
10    SL=DSQRT(DMAX1(0.D0,1.D0-CL**2))
      TMP1=STH*dSIN (PH-PHL)
      TMP2=(CTH-CTHL*CL)/STHL
      ALPH2=0.D0
      IF (TMP1.NE.0.d0.OR.TMP2.NE.0.d0) ALPH2=DATAN2(TMP1,TMP2)
      CALL PLOT (EARTHR*CEA,EARTHR*dASIN(SL*dSIN (ALPH1-ALPH2)),NEW)
      RETURN
C
C          DRAW AXES AND CALL FOR LABELING AND TERMINATION OF THIS PLOT
C**************************************************************************** 
      ENTRY ENDPLT
C**************************************************************************** 
C
C     IF NEWWR IS STILL TRUE, NO PLOTS WHERE PRODUCED
      IF(NEWWR) RETURN
      NEWWR=.TRUE.
C
C     SIGNAL END OF PLOT
      IF(NPLTDP.GT.0) WRITE(NPLTDP) 4,(NWDSRK,I=2,NWDSRK)
C
      IF(PLT.EQ.0.d0) RETURN
C
      TICKX=1.d-2*(YT-YB)
      DTIC=TIC*EARTHR
      CALL SETXY(APLT,-ALPHA,RMIN,ALPHA,RMAX)
C
      IF(APLT.EQ.4.d0) GO TO 200
      IF(APLT.EQ.2.d0) GO TO 25
C
C     CURVLINEAR PROJECTIONS(=1,3)
      CALL ARCTIC(-ALPHA,ALPHA,RMIN,TICKX,PLTHLB)
      CALL ARCTIC(-ALPHA,ALPHA,RMAX,-TICKX,PLTHLB)
      GO TO 300
C
C     CARTESIAN PROJECTION PUT IN BOTTOM BOUNDARY
200   CALL TIKLINE(-ALPHA,RMIN,ALPHA,RMIN,TIC,-TICKX,PLTHLB)
C     PUT IN TOP BOUNDARY
      CALL TIKLINE(-ALPHA,RMAX,ALPHA,RMAX,TIC,TICKX,PLTHLB)
C
300   TIKL=2.d-2*ALPHA
      LTIC=TICV
      IF(APLT.EQ.3.d0) LTIC=TICV*PFACTR
C     PUT IN LEFT BOUNDARY
      CALL TIKLINE(-ALPHA,RMIN,-ALPHA,RMAX,LTIC,TIKL,PLTLB)
C     PUT IN RIGHT BOUNDARY
      CALL TIKLINE(ALPHA,RMIN,ALPHA,RMAX,LTIC,-TIKL,PLTLBN)
C
      GO TO 50
C   DRAW TICKS, BOX FOR HORIZONTAL PLOT
25    CALL DRAWTKS(DTIC,TICV,XL,XR,YB,YT,PLOT)
C
50    IF(APLT.EQ.2.d0) THEN
        CALL PLTANOT(ID,OW/PIT2,XFQMDL,YFQMDL,XL,YB,XR,YT,DEGS,
     2     LLAT,LLON,RLAT,RLON,HB,HT,APLT,DTIC,DTIC/PFACTR,PLOT)
      ELSE
        CALL PLTANOT(ID,OW/PIT2,XFQMDL,YFQMDL,-ALPHA,RMIN,ALPHA,RMAX,
     2     DEGS,LLAT,LLON,RLAT,RLON,HB,HT,APLT,DTIC,TICV,PLOT)
      ENDIF
C
      CALL LABPLT
      CALL PLTEND
      RETURN
      END
 
 
 
 
C**************************************************************************** 
      SUBROUTINE OCNHD
C**************************************************************************** 
C     PRINTS PAGE HEADINGS
C     IF W(72) IS NEGATIVE, ONLY ONE HEADER OUTPUT INCLUDED IN RAYSET
 
      implicit double precision (a-h,o-z)
 
      CHARACTER pcc*(*),pc*1,BLANKS*100,DIVIDR*132,banner(10)*80
      CHARACTER NUMSTG*80,STMP*80,typ*20,id*80,strg5*200
      LOGICAL NOPUNCH
      INTEGER STRIM
C
C     TWO ENTRY POINTS ARE PROVIDED. ONE FOR THE FIRST PAGE HEADER
C     OF THE COMPUTATIONAL PRINTOUT AND FOR THE RAYSET FILE. THE
C     SECOND ENTRY IS FOR ALL SUBSEQUENT PAGES OF THE COMPUTATIONAL
C     PRINTOUT.
C
C     COMMON DECK "LL" INSERTED HERE
      real*8 MODL
      COMMON/LL/MODL(4),APH,APHPT,APHPR,APHPTH,APHPPH
      common/ident/id
C     COMMON DECK "PP" INSERTED HERE
      REAL*8 MODP
      COMMON/PP/MODP(4),P,PPT,PPR,PPTH,PPPH
C
C     COMMON DECK "SS" INSERTED HERE
      real*8 MODSURF
      COMMON/SS/ MODSURF(4)
      COMMON/SS/U,PUR,PURR,PURTH,PURPH
      COMMON/SS/PUTH,PUPH,PUTHTH,PUPHPH,PUTHPH,USELECT,UTIME
C     COMMON DECK "CONST" INSERTED HERE
      COMMON/PCONST/CREF,RGAS,GAMMA
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10
      COMMON/WW/MAXW,W(1000)
      REAL*8 MAXSTP,MAXERR,INTYP,LLAT,LLON
      EQUIVALENCE (EARTHR,W(1)),(RAY,W(2)),(XMTRH,W(3)),(TLAT,W(4)),
     1 (TLON,W(5)),(OW,W(6)),(FBEG,W(7)),(FEND,W(8)),(FSTEP,W(9)),
     2 (AZ1,W(10)),(AZBEG,W(11)),(AZEND,W(12)),(AZSTEP,W(13)),
     3 (BETA,W(14)),(ELBEG,W(15)),(ELEND,W(16)),(ELSTEP,W(17)),
     8   (RUNSUP,W(18)),(RCVRH,W(20)),
     4 (ONLY,W(21)),(HOP,W(22)),(MAXSTP,W(23)),(PLAT,W(24)),(PLON,W(25))
     5,(HMAX,W(26)),(RAYFNC,W(29)),(EXTINC,W(33)),
     6 (HMIN,W(27)),(RGMAX,W(28)),
     8 (INTYP,W(41)),(MAXERR,W(42)),(ERATIO,W(43)),
     6 (STEP1,W(44)),(STPMAX,W(45)),(STPMIN,W(46)),(FACTR,W(47)),
     7 (SKIP,W(71)),(RAYSET,W(72)),(PRTSRP,W(74)),(HITLET,W(75))
     9 ,(BINRAY,W(76)),(PAGLN,W(77)),(PLT,W(81)),(PFACTR,W(82)),
     1 (LLAT,W(83)),(LLON,W(84)),(RLAT,W(85)),(RLON,W(86))
     2,(TIC,W(87)),(HB,W(88)),(HT,W(89)),(TICV,W(96))
C     COMMON DECK "GG" INSERTED HERE
      real*8 MODG
      COMMON/GG/MODG(4)
      COMMON/GG/G,PGR,PGRR,PGRTH,PGRPH
      COMMON/GG/PGTH,PGPH,PGTHTH,PGPHPH,PGTHPH,GSELECT,GTIME
C     COMMON DECK "RR" INSERTED HERE
      real*8 MODREC
      COMMON/RR/ MODREC(4)
      COMMON/RR/F,PFR,PFRR,PFRTH,PFRPH
      COMMON/RR/PFTH,PFPH,PFTHTH,PFPHPH,PFTHPH,FSELECT,FTIME
C     COMMON DECK "B9" INSERTED HERE
      INTEGER GMX,GNTBL,GITBL,GFRMTBL,IDSG(20)
      COMMON/B9/GMX,GNTBL(10),GITBL(10),GFRMTBL(10),GGP(113)
      EQUIVALENCE (GGP,IDSG),(ANG,GGP(11))
C     COMMON DECK "B2" INSERTED HERE
      INTEGER DUMX,DUNTBL,DUITBL,DUFRMTB,IDSDU(20)
      COMMON/B2/DUMX,DUNTBL(10),DUITBL(10),DUFRMTB(10),DUGP(10)
      EQUIVALENCE (DUGP,IDSDU)
C     COMMON DECK "B4" INSERTED HERE
      INTEGER DCMX,DCNTBL,DCITBL,DCFRMTB,IDSDC(20)
      COMMON/B4/DCMX,DCNTBL(10),DCITBL(10),DCFRMTB(10),DCGP(10)
      EQUIVALENCE (DCGP,IDSDC)
C     COMMON DECK "B6" INSERTED HERE
      INTEGER DTMX,DTNTBL,DTITBL,DTFRMTB,IDSDT(20)
      COMMON/B6/DTMX,DTNTBL(10),DTITBL(10),DTFRMTB(10),DTGP(10)
      EQUIVALENCE (DTGP,IDSDT)
C     COMMON DECK "B1" INSERTED HERE
      INTEGER UMX,UNTBL,UITBL,UFRMTBL,IDSU(20)
      COMMON/B1/UMX,UNTBL(10),UITBL(10),UFRMTBL(10),UGP(10)
      EQUIVALENCE (UGP,IDSU)
C     COMMON DECK "B3" INSERTED HERE
      INTEGER CMX,CNTBL,CITBL,CFRMTBL,IDSC(20)
      COMMON/B3/CMX,CNTBL(10),CITBL(10),CFRMTBL(10),CGP(512)
      EQUIVALENCE (CGP,IDSC),(ANC,CGP(11))
C     COMMON DECK "B5" INSERTED HERE
      INTEGER TMX,TNTBL,TITBL,TFRMTBL,IDST(20)
      COMMON/B5/TMX,TNTBL(10),TITBL(10),TFRMTBL(10),TGP(262)
      EQUIVALENCE (TGP,IDST),(ANT,TGP(11))
C     COMMON DECK "B7" INSERTED HERE
      INTEGER MMX,MNTBL,MITBL,MFRMTBL,IDSM(20)
      REAL*8 MGP
      COMMON/B7/MMX,MNTBL(10),MITBL(10),MFRMTBL(10),MGP(10)
      EQUIVALENCE (MGP,IDSM)
C     COMMON DECK "HDR" INSERTED HERE
      CHARACTER INITID*80,DAT*10,TOD*10
      COMMON/HDR/SEC
      COMMON/HDRC/INITID,DAT,TOD
C     COMMON DECK "RINPLEX" INSERTED HERE
      REAL*8 KAY2,KAY2I
      COMPLEX*16 PNP,POLAR,LPOLAR
      LOGICAL SPACE
      CHARACTER DISPM*6,modrin*64
      COMMON/RINPL/DISPM
      COMMON /RIN/ MODRIN,RAYNAME(2,3),TYPE(3),SPACE
      COMMON/RIN/OMEGMIN,OMEGMAX,KAY2,KAY2I
      COMMON/RIN/PNP(10),POLAR,LPOLAR,SGN
C     COMMON DECK "CC" INSERTED HERE
      real*8 MODC
      COMMON/CC/MODC(4),CS,PCST,PCSR,PCSTH,PCSPH
C     COMMON DECK "MM" INSERTED HERE
      REAL*8 M,MODM
      COMMON/MM/MODM(4),M,PMT,PMR,PMTH,PMPH
C     COMMON DECK "FLAG" INSERTED HERE
      LOGICAL NEWWR,NEWWP,NEWTRC,PENET
      COMMON /FLG/ NTYP,NEWWR,NEWWP,NEWTRC,PENET,LINES,IHOP,HPUNCH
      COMMON/FLGP/NSET
C     COMMON DECK "TT" INSERTED HERE
      REAL*8 MODT
      COMMON/TT/MODT(4), T,PTT,PTR,PTTH,PTPH
C     COMMON DECK "UU" INSERTED HERE
      real*8 MODU
      COMMON/UU/MODU(4)
     1  ,V  ,PVT  ,PVR  ,PVTH  ,PVPH
     2  ,VR ,PVRT ,PVRR ,PVRTH ,PVRPH
     3  ,VTH,PVTHT,PVTHR,PVTHTH,PVTHPH
     4  ,VPH,PVPHT,PVPHR,PVPHTH,PVPHPH
C
C     COMMON DECK "B10" INSERTED HERE
      INTEGER DGMX,DGNTBL,DGITBL,DGFRMTB,IDSDG(20)
      COMMON/B10/DGMX,DGNTBL(10),DGITBL(10),DGFRMTB(10),DGGP(10)
      EQUIVALENCE (DGGP,IDSDG)
C     COMMON DECK "B8" INSERTED HERE
      INTEGER RMX,RNTBL,RITBL,RFRMTBL,IDSR(20)
      COMMON/B8/RMX,RNTBL(10),RITBL(10),RFRMTBL(10),RGP(10)
      EQUIVALENCE (RGP,IDSR)
C     COMMON DECK "CB11" INSERTED HERE
      INTEGER SMX,SNTBL,SITBL,SFRMTBL,IDSS(20)
      COMMON/B11/SMX,SNTBL(10),SITBL(10),SFRMTBL(10),SGP(11)
      EQUIVALENCE (SGP,IDSS),(ANS,SGP(11))
C     COMMON DECK "CB12" INSERTED HERE
      INTEGER DSMX,DSNTBL,DSITBL,DSFRMTB,IDSDS(20)
      COMMON/B12/DSMX,DSNTBL(10),DSITBL(10),DSFRMTB(10),DSGP(11)
      EQUIVALENCE (DSGP,IDSDS),(ANDS,DSGP(11))
C     COMMON DECK "CB17" INSERTED HERE
      INTEGER VMX,VNTBL,VITBL,VFRMTBL,IDSV(20)
      COMMON/B17/VMX,VNTBL(10),VITBL(10),VFRMTBL(10),VGP(53)
      EQUIVALENCE (VGP,IDSV),(ANV,VGP(11))
C     COMMON DECK "CB18" INSERTED HERE
      INTEGER DVMX,DVNTBL,DVITBL,DVFRMTB,IDSDV(20)
      COMMON/B18/DVMX,DVNTBL(10),DVITBL(10),DVFRMTB(10),DVGP(11)
      EQUIVALENCE (DVGP,IDSDV),(ANDV,DVGP(11))
C     COMMON DECK "CB19" INSERTED HERE
      INTEGER PRMX,PRNTBL,PRITBL,PRFRMTB,IDSPR(20)
      COMMON/B19/PRMX,PRNTBL(10),PRITBL(10),PRFRMTB(10),PRGP(11)
      EQUIVALENCE (PRGP,IDSPR),(ANP,PRGP(11))
C     COMMON DECK "CB20" INSERTED HERE
      INTEGER DPMX,DPNTBL,DPITBL,DPFRMTB,IDSDP(20)
      COMMON/B20/DPMX,DPNTBL(10),DPITBL(10),DPFRMTB(10),DPGP(11)
      EQUIVALENCE (DPGP,IDSDP),(ANDP,DPGP(11))
C
      PARAMETER (nbnrls=10,nblns=10)
C
      DATA BLANKS/' '/
      DATA BANNER/
     1  '<<<<<   H A R P O   >>>>>'
     2 ,'HAMILTONIAN ACOUSTIC RAY-TRACING PROGRAM FOR THE OCEAN'
     3 ,' '
     4 ,'BY'
     5 ,'R. M. JONES,  J. P. RILEY  AND  T. M. GEORGES'
     6 ,'WAVE PROPAGATION LABORATORY'
     7 ,'NATIONAL OCEANIC AND ATMOSPHERIC ADMINISTRATION'
     8 ,'BOULDER, COLORADO 80303'
     9 ,' '
     1 ,'PC VERSION 7/1//90 BY T. M. GEORGES, R. M. JONES AND R. S. LAWR
     2ENCE'/     
      DATA NOPUNCH/.FALSE./
      DATA IBLK/1H /
C
C**************************************************************************** 
      ENTRY HEADER1
C**************************************************************************** 
C
C     COMPUTE EFFECTIVE LINES COUNT BASED ON FIXED PAGE SIZE
      CALL NEWPAG(NPAG,INT(PAGLN),PC)
      CALL SFILL(DIVIDR,LEN(DIVIDR),'_')
      DIVIDR(1:1)=' '
      NTYP=2
      IF(RAY.NE.0.D0) NTYP=2.D0+DSIGN(1.D0,RAY)
C
1600  FORMAT(2(2(A8,F7.1,1X),2X))
C%      CALL PUTKST(3,'1'
C%     1  //DAT//TOD//BLANKS(:100)//'PAGE'//NUMSTG(NPAG,1,'(I5)'))
      call pcodes (pc,typ)
      call putkst(3,typ//pc)
      strg5=blanks(:3)//dat//tod//blanks(:100)//'PAGE'
     1 //numstg(npag,1,'(i4)')
      call putkst(3,strg5)
      CALL PUTKST(3,DIVIDR)
      CALL PUTKBK(3,1)
      DO 15 I=1,NBNRLS
15    CALL PUTKCT(3,BANNER(I))
      CALL PUTKBK(3,1)
      CALL PUTKST(3,DIVIDR)
      CALL PUTKBK(3,NBLNS)
C
C%      CALL PUTKST(3,BLANKS(:57)//'RUN SET NUMBER'
C%     1  //NUMSTG(NSET,1,'(a5)'))
      strg5=blanks(:57)//'RUN SET NUMBER'//numstg(nset,1,'(i5)')
      call putkst (3,strg5)
      CALL PUTKBK(3,1)
C*      WRITE(STMP,'(10A8)') ID
      WRITE(STMP,'(A80)') ID
      CALL PUTKST(3,BLANKS(:57)//'OCEAN MODEL ID: '//STMP(:3))
      CALL PUTKBK(3,1)
      CALL PUTKST(3,
     1 BLANKS(:25)//'OCEAN MODEL DESCRIPTION:  '//STMP(7:))
      CALL PUTKBK(3,1)
C
      WRITE(STMP,'(2A8)') (RAYNAME(I,NTYP),I=1,2)
      CALL PUTKST(3,STMP)
      CALL PUTKBK(3,1)
      CALL PUTKST(3,DIVIDR)
      CALL PUTKBK(3,1)
      CALL PUTKST(3,
     1  BLANKS(:8)//'MODEL                 SUBROUTINE    DATA SET'
     1  //'       DESCRIPTION')
      CALL PUTKST(3,
     1  BLANKS(:8)//'TYPE                     NAME         ID')
      CALL PUTKBK(3,1)
      CALL PUTKST(3,DIVIDR)
      CALL PUTKBK(3,1)
      CALL PUTKST(3,' DISPERSION RELATION           '//DISPM
     1  //BLANKS(:16)//MODRIN )
C*     1  //BLANKS(:16)//NUMSTG(MODRIN,8,'(8A8)'))
      CALL PUTDES(3,'BACKGROUND CURRENT VELOCITY',MODU,IDSU)
      CALL PUTDES(3,'CURRENT VELOCITY PERTURBATION',MODU(3),IDSDU)
      CALL PUTDES(3,'BACKGROUND SOUND SPEED',MODC,IDSC)
      CALL PUTDES(3,'SOUND SPEED PERTURBATION',MODC(3),IDSDC)
      CALL PUTDES(3,'BACKGROUND BOTTOM',MODG,IDSG)
      CALL PUTDES(3,'BOTTOM PERTURBATION',MODG(3),IDSDG)
      CALL PUTDES(3,'ABSORPTION',MODL,IDSV)
      CALL PUTDES(3,'ABSORPTION PERTURBATION',MODL(3),IDSDV)
      CALL PUTDES(3,'RECEIVER SURFACE',MODREC,IDSR)
      CALL PUTDES(3,'OCEAN SURFACE',MODSURF,IDSS)
      CALL PUTDES(3,'OCEAN SURFACE PERTURBATION',MODSURF(3),IDSDS)
      CALL PUTKST(3,DIVIDR)
C
      NOPUNCH=NOPUNCH.AND.RAYSET.LT.0.0
      IF(NOPUNCH) RETURN
      NOPUNCH=RAYSET.LT.0.0
C
C*    WRITE(9,1200) ID,DAT,TOD
      WRITE(9,1201) ID,DAT,TOD
1201  FORMAT(A80,2(A10,2X))
      WRITE(9,1600) MODU,MODC,MODG,MODL,MODREC(1),MODREC(2),MODSURF
C
      IF(IDSU(1) .NE. IBLK) WRITE(9,1205) IDSU
      IF(IDSDU(1) .NE. IBLK) WRITE(9,1205) IDSDU
      IF(IDSC(1) .NE. IBLK) WRITE(9,1205) IDSC
      IF(IDSDC(1) .NE. IBLK) WRITE(9,1205) IDSDC
      IF(IDST(1) .NE. IBLK) WRITE(9,1205) IDST
      IF(IDSDT(1) .NE. IBLK) WRITE(9,1205) IDSDT
      IF(IDSM(1) .NE. IBLK) WRITE(9,1205) IDSM
      IF(IDSG(1) .NE. IBLK) WRITE(9,1205) IDSG
C
1000  format (A1/A80,21X,2A10,' PAGE',I4)
1200  FORMAT(10A8,2(A10,2x))
1205  format (20A4)
C
      RETURN
C**************************************************************************** 
      ENTRY HEADER2
C**************************************************************************** 
C     COMPUTE EFFECTIVE LINES COUNT BASED ON FIXED PAGE SIZE
      CALL NEWPAG(NPAG,INT(PAGLN),PC)
      LINES=LINES+5
C
      WRITE(3,1000) PC,ID,DAT,TOD,NPAG
      WRITE(3,2400) SNGL(AZ1*DEGS),SNGL(TLAT*DEGS),SNGL(OW/PIT2),
     1         SNGL(BETA*DEGS)
     1        ,SNGL(TLON*DEGS),MAXERR
 2400 FORMAT (
     1  /' AZIMUTH ANGLE OF TRANSMISSION   =',F12.6,' DEG'
     2  ,'      TRANSMITTER  LATITUDE  =',F12.6,' DEG'
     3  ,'    FREQUENCY =',F18.6,' HZ'
     4  /' ELEVATION ANGLE OF TRANSMISSION =',F12.6,' DEG'
     5  ,'      TRANSMITTER  LONGITUDE =',F12.6,' DEG'
     6  ,'    SINGLE STEP ERROR =',1PG13.6/)
C
      RETURN
C**************************************************************************** 
      ENTRY PUTDVR(NUNIT)
C**************************************************************************** 
      CALL PUTKST(NUNIT,DIVIDR)
      RETURN
C
C
C**************************************************************************** 
      ENTRY PUTHDR(NUNIT,PCC,NP)
C**************************************************************************** 
C%      CALL PUTKST(NUNIT,
C%     1  PCC//DAT//TOD//BLANKS(:100)//'PAGE'//NUMSTG(NP,1,'(I5)'))
      strg5=pcc//blanks(:2)//dat//tod//blanks(:100)//'PAGE'
     1 //numstg(np,1,'(i4)')
      call putkst(nunit,strg5)
      RETURN
      END
 
 
 
 
C****************************************************************************
      SUBROUTINE PUTDES(NUNIT,DES,MOD,ids)
C*****************************************************************************
C     TITLING AID, OUTPUT MODULE DESCRIPTION AS HOLLERITH NAME
C     AND NUMERIC IDENTIFIER
      implicit double precision (a-h,o-z)
      CHARACTER DES*(*),TITLE*30,SMODL*20,NUMSTG*80,strg5*200
      REAL*8 MOD(2)
      TITLE=DES
      WRITE(SMODL,'(A8,F10.2)')MOD
      CALL PUTKST(NUNIT, ' '//TITLE//SMODL//NUMSTG(ids,20,'(20a4)'))
      return
      END
 
 
 
 
C**************************************************************************** 
      FUNCTION NUMSTG(V,N,FRM)
C****************************************************************************
C     CONVERTS AN integer v(n) to a string numstg*80
C     USING FORMAT frm PROVIDED IN CALL
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER NUMSTG*80,FRM*(*)
      INTEGER V(N)
C
      NUMSTG=' '
      WRITE(NUMSTG,FRM) V
      RETURN
      END
 
 
 
 
 
C**************************************************************************** 
      SUBROUTINE SFILL(STG,LN,C)
C**************************************************************************** 
C     FILLS A STRING WITH N SPECIFIED CHARACTERS
      implicit double precision (a-h,o-z)
      CHARACTER STG*(*),C*1
C
      DO 10 I=1,LN
10    STG(I:I)=C
      RETURN
      END
 
 
 
 
 
C**************************************************************************** 
      INTEGER FUNCTION STRIM(C)
C**************************************************************************** 
C     DETERMINES POSITION OF LAST NONBLANK CHARACTER OF A STRING
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER*(*) C
C
      L=LEN(C)
      DO 100 I=L,1,-1
100   IF(C(I:I) .NE. ' ') GO TO 200
      I=0
200   STRIM=I
      return
      END
 
 
 
 
 
C**************************************************************************** 
      FUNCTION RERR(V,ERT,ELAB,NKV,PREF)
C**************************************************************************** 
C     RETURNS RELATIVE ERROR OF VARIABLE RKVAR(NKV) IF PREF<>0 AND
C     THE ERROR IS GREATER THAN PREVIOUS ERROR 'V'.
C     COMMON DECK "RKAM" INSERTED HERE
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      character*(*) elab
      REAL*8 KR,KTH,KPH
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)
C
      RERR=V
      IF(NKV.LE.0) RETURN
      IF(PREF.EQ.0.d0) RETURN
      V1=RKVARS(NKV)/PREF-1.d0
      IF(dABS(V1).LT.dABS(V)) RETURN
C
      ERT=ELAB
      RERR=V1
      END
 
 
 
 
 
C**************************************************************************** 
      SUBROUTINE RERROR(ROUTIN,STR,VAL)
C**************************************************************************** 
C     REPORTS ERROR CONDITIONS AND STOPS PROGRAM.
C
C%      implicit double precision (a-h,o-z)
	character routin*(*),str*(*)
      PRINT 10,ROUTIN,STR,VAL
10    FORMAT(' ERROR CONDITION IN RAYTRACE ROUTINE <'
     1  ,A6,'> DUE TO ',A10,/,10x,' VALUE=  ',F8.2)
      STOP
      END
 
 
 
 
 
C**************************************************************************** 
      SUBROUTINE STOPIT(A)
C**************************************************************************** 
C     PRINTS CONDITION AND STOPS PROGRAM
C     AFTER CALLING THE SYSTEM POST MORTEM DUMP.
      implicit double precision (a-h,o-z)
      CHARACTER A*(*)
C
      PRINT 100, A
      CALL MORTEM
100   FORMAT('*** STOPIT WITH CONDITION <',A,'>')
      STOP
      END
 
 
 
 
 
C**************************************************************************** 
      SUBROUTINE PUTKST(NUNIT,STRG)
C**************************************************************************** 
C     WRITE LINE OF OUTPUT TO PRINTER UNIT ADDING TO LINE COUNT
      implicit double precision (a-h,o-z)
      CHARACTER STRG*(*)
      INTEGER STRIM
C     COMMON DECK "FLAG" INSERTED HERE
      LOGICAL NEWWR,NEWWP,NEWTRC,PENET
      COMMON /FLG/ NTYP,NEWWR,NEWWP,NEWTRC,PENET,LINES,IHOP,HPUNCH
      COMMON/FLGP/NSET
C
C     PUT OUT A STRING WITH LINE COUNT INCREMENT
      LINES=LINES+1
      LN=MAX0(1,MIN0(STRIM(STRG),132))
      WRITE(NUNIT,'(A)') STRG(1:LN)
      RETURN
      END
C
 
 
 
 
 
C**************************************************************************** 
      SUBROUTINE PUTKCT(NUNIT,STRG)
C**************************************************************************** 
      implicit double precision (a-h,o-z)
C     PUT OUT A CENTERED LINE AND COUNT
      CHARACTER STRG*(*)
      CHARACTER BLANKS*100,STMP*80
      INTEGER STRIM
      CHARACTER pcc*(*),PC*1,typ*20
C
C     COMMON DECK "FLAG" INSERTED HERE
      LOGICAL NEWWR,NEWWP,NEWTRC,PENET
      COMMON /FLG/ NTYP,NEWWR,NEWWP,NEWTRC,PENET,LINES,IHOP,HPUNCH
      COMMON/FLGP/NSET
C
      DATA BLANKS/' '/
C
      NTRM=MAX0(1,STRIM(STRG))
      NBLKS=66-(NTRM+1)/2
      LINES=LINES+1
      LN=MIN0(NTRM,132-NBLKS)
      STMP=STRG
      WRITE(NUNIT,'(A)') BLANKS(:NBLKS)//STMP(:LN)
      RETURN
C
C**************************************************************************** 
      ENTRY PUTKBK(NUNIT,NBKS)
C**************************************************************************** 
C     PUT NBLKS BLANK LINES TO UNIT NUNIT
C
      LINES=LINES+1
      DO 100 I=1,NBKS
100      WRITE(NUNIT,'(1X)')
      RETURN
C
C**************************************************************************** 
      ENTRY NEWPAG(NXP,LIP,PCC)
C**************************************************************************** 
C............NXP= page number (input)
C............LIP= lines per page (input)
C............PC= new page control character (output)
C     THIS ENTRY COMPUTES NEXT PAGE NUMBER AND LINE COUNT TO END
C     OF CURRENT PAGE
C     RETAIN LINES PER PAGE FOR USE BY TSTPAG ENTRY
      LINSPP=LIP
C     THE LOGIC HERE ASSUMES THAT A FORM FEED IS ALWAYS NEEDED.
C     A MORE SOPHISTICATED VERSION WOULD DETERMINE IF A NEW PAGE
C     IS NEEDED AND RETURN PC=' ' IF NOT. THE CALLING PROGRAM COULD
C     USE OR IGNORE.
C%    PC='1'
      call pcodes (pcc,typ)
C     COMPUTE CURRENT PAGE NUMBER
      NPAG=(LINES+LINSPP-1)/LINSPP
C     COMPUTE LINE COUNT AT END OF CURRENT PAGE
      LINES=NPAG*LINSPP
      NXP=NPAG+1
C     COMPUTE LINE COUNT AT END OF NEXT PAGE
      LNSXPG=LINES+LINSPP-1
      RETURN
C
C**************************************************************************** 
      ENTRY TSTPAG(NUNIT)
C**************************************************************************** 
C     OUTPUT HEADER IF HAVE REACHED LAST LINE OF PAGE
      IF(LINES.GE.LNSXPG) THEN
C     COMPUTE CURENT PAGE NUMBER
         NPAG=(LINES+LINSPP-1)/LINSPP
C     AND ITS LAST LINE
         LINES=NPAG*LINSPP
         LNSXPG=LINES+LINSPP-1
         NPAG=NPAG+1
         call pcodes (pc,typ)
         CALL PUTHDR(NUNIT,pc,NPAG)
         WRITE(NUNIT,'(1X)')
         LINES=LINES+1
      ENDIF
      END
 
 
 
 
 
C**************************************************************************** 
      SUBROUTINE OPNREP(IUN,FNAME)
C**************************************************************************** 
      implicit double precision (a-h,o-z)
C     THIS OPEN OPERATION ALLOWS FOR AN EXISTING VERSION OF
C     A FILE. IF THE FILE EXISTS IT IS OVERWRITTEN.
C     THE LOGIC USED IS HOPEFULLY MACHINE INDEPENDENT.
      CHARACTER*(*) FNAME
      CALL ZAPFIL(IUN,FNAME)
      OPEN(IUN,FILE=FNAME,STATUS='NEW')
      REWIND IUN
      RETURN
C**************************************************************************** 
      ENTRY OPNURP(IUN,FNAME)
C**************************************************************************** 
      CALL ZAPFIL(IUN,FNAME)
      OPEN(IUN,FILE=FNAME,STATUS='NEW',FORM='UNFORMATTED')
      REWIND IUN
      END
 
 
 
 
 
C**************************************************************************** 
      SUBROUTINE ZAPFIL(IUN,FNAME)
C**************************************************************************** 
      implicit double precision (a-h,o-z)
C     TEST FOR EXISTENCE OF FILE AND DELETE
      CHARACTER FNAME*(*)
      LOGICAL EX
      INQUIRE(FILE=FNAME,EXIST=EX)
      IF(EX) THEN
         OPEN(UNIT=IUN,FILE=FNAME,STATUS='OLD')
         CLOSE(UNIT=IUN,STATUS='DELETE')
      ENDIF
      END
 
 
 
 
 
C**************************************************************************** 
      SUBROUTINE OVERRD(VAR,TST,DEFT,NFLG,NFVLEQ,NFVLNE)
C**************************************************************************** 
      implicit double precision (a-h,o-z)
C     OVERRIDE SUPPORT ROUTINE
C     TEST 'VAR' AGAINST 'TST', IF EQUAL SET TO DEFAULT 'DEFT'
C     AND ALSO SET INTEGER FLAG 'NFLG' TO VALUE 'NFVLEQ' ELSE 'NFVLNE'
C
      IF(VAR.EQ.TST) THEN
         VAR=DEFT
         NFLG=NFVLEQ
      ELSE
         NFLG=NFVLNE
      ENDIF
      RETURN
      END
 
 
 
 
 
C**************************************************************************** 
      SUBROUTINE SFILTR(C,S,KSET)
C**************************************************************************** 
C     FILTERS EXTRANEOUS CHARACTERS FROM PLOT LABELS
      implicit double precision (a-h,o-z)
      CHARACTER*(*) C,S,KSET
C
      LN=LEN(C)
      J=0
      DO 10 I=1,LN
        IF(INDEX(KSET,C(I:I)).EQ.0) THEN
           J=J+1
           S(J:J)=C(I:I)
        ENDIF
10    CONTINUE
      RETURN
      END
 
 
 
 
 
C**************************************************************************** 
      FUNCTION ALCOSH(X)
C**************************************************************************** 
C     COMPUTE LOG(COSH(X)) AND USE LARGE ARGUMENT APPROXIMATION
C     WHEN POSSIBLE.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DATA ALOG2/.6931471805599453d0/
C
      IF(dABS(X).GT.50.d0) GO TO 10
      EX=dEXP(X)
      ALCOSH=DLOG((EX+1.D0/EX)*.5D0)
      RETURN
10    ALCOSH=DABS(X)-alog2
      END
 
 
 
 
 
C**************************************************************************** 
      SUBROUTINE GAUSEL (C,NRD,NRR,NCC,NSF)
C**************************************************************************** 
C     MATRIX INVERSION BY METHOD OF GAUSSIAN PIVOT
C     SOLVES NCC-NRR SETS OF NRR EQUATIONS
C********* SAME AS SUBROUTINE GAUSSEL WRITTEN BY L. DAVID LEWIS ********
      implicit double precision (a-h,o-z)
      DIMENSION C(NRD,NCC),L(128,2)
C          BITS = 2.**-18
      DATA BITS/3.8146972656d-6/
      NR=NRR
      NC=NCC
      IF(NC.LT.NR.OR.NR.GT.128.OR.NR.LE.0) STOP 5000
C
C          INITIALIZE.
      NSF=0
      NRM=NR-1
      NRP=NR+1
      D=1.
      LSD=1
      DO 1 KR=1,NR
      L(KR,1)=KR
1     L(KR,2)=0
      IF(NR.EQ.1) GO TO 42
C
C          ELIMINATION PHASE.
      DO 41 KP=1,NRM
      KPP=KP+1
      PM=0.d0
      MPN=0
C
C          SEARCH COLUMN KP FROM DIAGONAL DOWN FOR MAX PIVOT.
      DO 2 KR=KP,NR
      LKR=L(KR,1)
      PT=dABS(C(LKR,KP))
      IF(PT.LE.PM) GO TO 2
      PM=PT
      MPN=KR
      LMP=LKR
2     CONTINUE
C
C          IF MAX PIVOT IS ZERO, MATRIX IS SINGULAR.
      IF(MPN.EQ.0) GO TO 9
      NSF=NSF+1
      IF(MPN.EQ.KP) GO TO 3
C
C          NEW ROW NUMBER KP HAS MAX PIVOT.
      LSD=-LSD
      L(KP,2)=L(KP,1)
      L(MPN,1)=L(KP,1)
      L(KP,1)=LMP
C
C          ROW OPERATIONS TO ZERO COLUMN KP BELOW DIAGONAL.
3     MKP=L(KP,1)
      P=C(MKP,KP)
      D=D*P
      DO 41 KR=KPP,NR
      MKR=L(KR,1)
      Q=C(MKR,KP)/P
      IF(Q.EQ.0.d0) GO TO 41
C
C          SUBTRACT Q * PIVOT ROW FROM ROW KR.
      DO 4 LC=KPP,NC
      R=Q*C(MKP,LC)
      C(MKR,LC)=C(MKR,LC)-R
4     IF(dABS(C(MKR,LC)).LT.dABS(R)*BITS) C(MKR,LC)=0.d0
41    CONTINUE
C
C          LOWER RIGHT HAND CORNER.
42    LNR=L(NR,1)
      P=C(LNR,NR)
      IF(P.EQ.0.d0) GO TO 9
      NSF=NSF+1
      D=D*P*LSD
      IF(NR.EQ.NC) GO TO 8
C
C          BACK SOLUTION PHASE.
      DO 61 MC=NRP,NC
      C(LNR,MC)=C(LNR,MC)/P
      IF(NR.EQ.1) GO TO 61
      DO 6 LL=1,NRM
      KR=NR-LL
      MR=L(KR,1)
      KRP=KR+1
      DO 5 MS=KRP,NR
      LMS=L(MS,1)
      R=C(MR,MS)*C(LMS,MC)
      C(MR,MC)=C(MR,MC)-R
5     IF(dABS(C(MR,MC)).LT.dABS(R)*BITS) C(MR,MC)=0.d0
6     C(MR,MC)=C(MR,MC)/C(MR,KR)
61    CONTINUE
C
C          SHUFFLE SOLUTION ROWS BACK TO NATURAL ORDER.
      DO 71 LL=1,NRM
      KR=NR-LL
      MKR=L(KR,2)
      IF(MKR.EQ.0) GO TO 71
      MKP=L(KR,1)
      DO 7 LC=NRP,NC
      Q=C(MKR,LC)
      C(MKR,LC)=C(MKP,LC)
7     C(MKP,LC)=Q
71    CONTINUE
C
C          NORMAL AND SINGULAR RETURNS.  GOOD SOLUTION COULD HAVE D=0.
8     C(1,1)=D
      GO TO 91
9     C(1,1)=0.d0
   91 RETURN
         END
 
 
 
 
 
