      SUBROUTINE BQWFWC                                                 BQWC001 

C********* CALCULATES A HAMILTONIAN H                                   BQWC002 

C********* (= BOOKER QUARTIC FOR VERTICAL INCIDENCE, S=0, C=1)          BQWC003 

C********* AND ITS PARTIAL DERIVATIVES WITH RESPECT TO                  BQWC004 

C********* TIME, R, THETA, PHI, OMEGA, KR, KTHETA, AND KPHI.            BQWC005 

C********* WITH FIELD, WITH COLLISIONS                                  BQWC006 

      COMMON /CONST/ PI,PIT2,PID2,DEGS,RADIAN,K,C,LOGTEN                BQWC007 

      COMMON /RIN/ MODRIN(3),COLL,FIELD,SPACE,KAY2,H,PHPT,PHPR,PHPTH,   BQWC008 

     1             PHPPH,PHPOM,PHPKR,PHPKTH,PHPKPH,KPHPK,POLAR,LPOLAR,  BQWC009 

     2             SGN                                                  BQWC010 

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  BQWC011 

      COMMON /YY/ MODY,Y,PYPR,PYPTH,PYPPH,YR,PYRPR,PYRPT,PYRPP,YTH,PYTPRBQWC012 

     1           ,PYTPT,PYTPP,YPH,PYPPR,PYPPT,PYPPP                     BQWC013 

      COMMON /ZZ/ MODZ,Z,PZPR,PZPTH,PZPPH                               BQWC014 

      COMMON R,TH,PH,KR,KTH,KPH   /WW/ ID(10),W0,W(400)                 BQWC015 

      COMMON /RK/ N,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART       BQWC016 

      COMMON /FLG/ NTYP,NEWWR,NEWWP,PENET,LINES,IHOP,HPUNCH             BQWC017 

      EQUIVALENCE (RAY,W(1)),(F,W(6))                                   BQWC018 

      LOGICAL SPACE                                                     BQWC019 

      REAL KR,KTH,KPH,K2,KDOTY,K4,KDOTY2                                BQWC020 

      COMPLEX KAY2,H,PHPT,PHPR,PHPTH,PHPPH,PHPOM,PHPKR,PHPKTH,PHPKPH,   BQWC021 

     1        POLAR,LPOLAR,I,U,RAD,D,PNPPS,PNPX,PNPY,PNPZ,UX,UX2,D2,    BQWC022 

     2        KPHPK,U2,A,B,ALPHA,BETA,GAMMA,PHPX,PHPY2,PHPK2,PHPU,PHPZ, BQWC023 

     3        N2,PNPR,PNPTH,PNPPH,PNPVR,PNPVTH,PNPVPH,NNP,PNPT          BQWC024 

      DATA (MODRIN=8HBOOKER Q,8HUARTIC, ,8HS=0, C=1),(COLL=1.),         BQWC025 

     1     (FIELD=1.),                                                  BQWC026 

     2     (X=0.),(PXPR=0.),(PXPTH=0.),(PXPPH=0.),(PXPT=0.),            BQWC027 

     3     (Y=0.),(PYPR=0.),(PYPTH=0.),(PYPPH=0.),(YR=0.),(PYRPR=0.),   BQWC028 

     4     (PYRPT=0.),(PYRPP=0.),(YTH=0.),(PYTPR=0.),(PYTPT=0.),        BQWC029 

     5     (PYTPP=0.),(YPH=0.),(PYPPR=0.),(PYPPT=0.),(PYPPP=0.)         BQWC030 

     6    ,(Z=0.),(PZPR=0.),(PZPTH=0.),(PZPPH=0.),                      BQWC031 

     7     (I=(0.,1.)),(ABSLIM=1.E-5),(SGN=1.)                          BQWC032 

      ENTRY RINDEX                                                      BQWC033 

      OM=PIT2*1.E6*F                                                    BQWC034 

      C2=C*C                                                            BQWC035 

      K2=KR*KR+KTH*KTH+KPH*KPH                                          BQWC036 

      OM2=OM*OM                                                         BQWC037 

      CALL ELECTX                                                       BQWC038 

      IF(X.LT..1) GO TO 2                                               BQWC039 

      K4=K2*K2                                                          BQWC040 

      OM4=OM2*OM2                                                       BQWC041 

      C4=C2*C2                                                          BQWC042 

      CALL MAGY                                                         BQWC043 

      Y2=Y*Y                                                            BQWC044 

      KDOTY=KR*YR+KTH*YTH+KPH*YPH                                       BQWC045 

      KDOTY2=KDOTY*KDOTY                                                BQWC046 

      CALL COLFRZ                                                       BQWC047 

      U=CMPLX(1.,-Z)                                                    BQWC048 

      U2=U*U                                                            BQWC049 

      UX=U-X                                                            BQWC050 

      UX2=UX*UX                                                         BQWC051 

      A=UX*U2-U*Y2                                                      BQWC052 

      B=-2.*U*UX2+Y2*(2.*U-X)                                           BQWC053 

      ALPHA=A*C4*K4+X*KDOTY2*C4*K2                                      BQWC054 

      BETA=B*C2*OM2*K2-X*KDOTY2*C2*OM2                                  BQWC055 

      GAMMA=(UX2-Y2)*UX*OM4                                             BQWC056 

      H=ALPHA+BETA+GAMMA                                                BQWC057 

      PHPX=-U2*C4*K4+KDOTY2*C4*K2+(4.*U*UX-Y2)*C2*OM2*K2-KDOTY2*C2*OM2+ BQWC058 

     1 (-3.*UX2+Y2)*OM4                                                 BQWC059 

      PHPY2=-U*C4*K4+(2.*U-X)*C2*OM2*K2-UX*OM4                          BQWC060 

      PHPKY2   =X*C2*(C2*K2-OM2)                                        BQWC061 

      PHPU=(2.*U*UX+U2-Y2)*C4*K4+2.*(Y2-UX2-2.*U*UX)*C2*K2*OM2+(3.*UX2  BQWC062 

     1 -Y2)*OM4                                                         BQWC063 

      PHPZ=-I*PHPU                                                      BQWC064 

      PHPK2=2.*A*C4*K2+X*KDOTY2*C4+B*C2*OM2                             BQWC065 

      PHPT=PHPX*PXPT                                                    BQWC066 

      PHPR =PHPX*PXPR +PHPY2*2.*Y*PYPR +PHPKY2  *2.*KDOTY*              BQWC067 

     1 (KR*PYRPR+KTH*PYTPR+KPH*PYPPR) +PHPZ*PZPR                        BQWC068 

      PHPTH=PHPX*PXPTH+PHPY2*2.*Y*PYPTH+PHPKY2  *2.*KDOTY*              BQWC069 

     1 (KR*PYRPT+KTH*PYTPT+KPH*PYPPT) +PHPZ*PZPTH                       BQWC070 

      PHPPH=PHPX*PXPPH+PHPY2*2.*Y*PYPPH+PHPKY2  *2.*KDOTY*              BQWC071 

     1 (KR*PYRPP+KTH*PYTPP+KPH*PYPPP) +PHPZ*PZPPH                       BQWC072 

      PHPOM=(2.*BETA+4.*GAMMA)/OM                                       BQWC073 

     1 -2.*PHPX*X/OM-2.*PHPY2*Y2/OM-2.*PHPKY2  *KDOTY2/OM -PHPZ*Z/OM    BQWC074 

      PHPKR= 2.*PHPK2*KR +2.*KDOTY*PHPKY2  *YR                          BQWC075 

      PHPKTH=2.*PHPK2*KTH+2.*KDOTY*PHPKY2  *YTH                         BQWC076 

      PHPKPH=2.*PHPK2*KPH+2.*KDOTY*PHPKY2  *YPH                         BQWC077 

      KAY2=K2*(-BETA+SGN*RAY*CSQRT(BETA**2-4.*ALPHA*GAMMA))/(2.*ALPHA)  BQWC078 

C                                                                       BQWC079 

      IF(RSTART.EQ.0.) GO TO 1                                          BQWC080 

      SCALE=SQRT((-REAL(BETA)+SGN*RAY*SQRT(REAL(BETA)**2                BQWC081 

     1 -4.*REAL(ALPHA)*REAL(GAMMA)))/(2.*REAL(ALPHA)))                  BQWC082 

      KR =SCALE*KR                                                      BQWC083 

      KTH=SCALE*KTH                                                     BQWC084 

      KPH=SCALE*KPH                                                     BQWC085 

  1   CONTINUE                                                          BQWC086 

C********* THE FOLLOWING 3 CARDS USED FOR RAY TRACING IN COMPLEX SPACE  BQWC087 

C     IF(CABS((-BETA-SGN*RAY*CSQRT(BETA**2-4.*ALPHA*GAMMA))/ALPHA-2.).  BQWC088 

C    1LT.CABS((-BETA+SGN*RAY*CSQRT(BETA**2-4.*ALPHA*GAMMA))/ALPHA-2.)   BQWC089 

C    2 .AND.RSTART.EQ.0.) SGN=-SGN                                      BQWC090 

      KPHPK=4.*ALPHA+2.*BETA                                            BQWC091 

      SPACE=CABS(C2*KAY2/OM2-1.).LT.ABSLIM                              BQWC092 

      POLAR =SQRT(K2)*(U+X*OM2/(C2*KAY2-OM2))/KDOTY*I                   BQWC093 

      LPOLR = SQRT(Y2-KDOTY2/K2)/UX*(1.-C2*KAY2/OM2)*I                  BQWC094 

      RETURN                                                            BQWC095 

C          CALCULATES THE REFRACTIVE INDEX AND ITS GRADIENT USING THE   BQWC096 

C          APPLETON-HARTREE FORMULA WITH FIELD, WITH COLLISIONS         BQWC097 

  2   CONTINUE                                                          BQWC098 

      VR =C/OM*KR                                                       BQWC099 

      VTH=C/OM*KTH                                                      BQWC100 

      VPH=C/OM*KPH                                                      BQWC101 

      CALL MAGY                                                         BQWC102 

      V2=VR**2+VTH**2+VPH**2                                            BQWC103 

      VDOTY=VR*YR+VTH*YTH+VPH*YPH                                       BQWC104 

      YLV=VDOTY/V2                                                      BQWC105 

      YL2=VDOTY**2/V2                                                   BQWC106 

      YT2=Y**2-YL2                                                      BQWC107 

      YT4=YT2*YT2                                                       BQWC108 

      CALL COLFRZ                                                       BQWC109 

      U=CMPLX(1.,-Z)                                                    BQWC110 

      UX=U-X                                                            BQWC111 

      UX2=UX*UX                                                         BQWC112 

      RAD=SGN*RAY*CSQRT(YT4+4.*YL2*UX2)                                 BQWC113 

      D=2.*U*UX-YT2+RAD                                                 BQWC114 

      D2=D*D                                                            BQWC115 

      N2=1.-2.*X*UX/D                                                   BQWC116 

      PNPPS=2.*X*UX*(-1.+(YT2-2.*UX2)/RAD)/D2                           BQWC117 

      PPSPR= YL2/Y*PYPR -(VR*PYRPR+VTH*PYTPR+VPH*PYPPR)*YLV             BQWC118 

      PPSPTH=YL2/Y*PYPTH-(VR*PYRPT+VTH*PYTPT+VPH*PYPPT)*YLV             BQWC119 

      PPSPPH=YL2/Y*PYPPH-(VR*PYRPP+VTH*PYTPP+VPH*PYPPP)*YLV             BQWC120 

      PNPX=-(2.*U*UX2-YT2*(U-2.*X)+(YT4*(U-2.*X)+4.*YL2*UX*UX2)/RAD)/D2 BQWC121 

      PNPY=2.*X*UX*(-YT2+(YT4+2.*YL2*UX2)/RAD)/(D2*Y)                   BQWC122 

      PNPZ=I*X*(-2.*UX2-YT2+YT4/RAD)/D2                                 BQWC123 

      NNP=N2-(2.*X*PNPX+Y*PNPY+Z*PNPZ)                                  BQWC124 

      PNPR =PNPX*PXPR +PNPY*PYPR +PNPZ*PZPR +PNPPS*PPSPR                BQWC125 

      PNPTH=PNPX*PXPTH+PNPY*PYPTH+PNPZ*PZPTH+PNPPS*PPSPTH               BQWC126 

      PNPPH=PNPX*PXPPH+PNPY*PYPPH+PNPZ*PZPPH+PNPPS*PPSPPH               BQWC127 

      PNPVR =PNPPS*(VR *YL2-VDOTY*YR )/V2                               BQWC128 

      PNPVTH=PNPPS*(VTH*YL2-VDOTY*YTH)/V2                               BQWC129 

      PNPVPH=PNPPS*(VPH*YL2-VDOTY*YPH)/V2                               BQWC130 

      PNPT=PNPX*PXPT                                                    BQWC131 

      SPACE=REAL(N2).EQ.1..AND.ABS(AIMAG(N2)).LT.ABSLIM                 BQWC132 

      POLAR=-I*SQRT(V2)*(-YT2+RAD)/(2.*VDOTY*UX)                        BQWC133 

      GAM=(-YT2+RAD)/(2.*UX)                                            BQWC134 

      LPOLAR=I*X*SQRT(YT2)/(UX*(U+GAM))                                 BQWC135 

      KAY2=OM2/C2*N2                                                    BQWC136 

      IF(RSTART.EQ.0.) GO TO 3                                          BQWC137 

      SCALE=SQRT(REAL(KAY2)/K2)                                         BQWC138 

      KR =SCALE*KR                                                      BQWC139 

      KTH=SCALE*KTH                                                     BQWC140 

      KPH=SCALE*KPH                                                     BQWC141 

  3   CONTINUE                                                          BQWC142 

      H=.5*(C2*K2/OM2-N2)                                               BQWC143 

      PHPT =-PNPT                                                       BQWC144 

      PHPR =-PNPR                                                       BQWC145 

      PHPTH=-PNPTH                                                      BQWC146 

      PHPPH=-PNPPH                                                      BQWC147 

      PHPOM=-NNP/OM                                                     BQWC148 

      PHPKR =C2/OM2*KR -C/OM*PNPVR                                      BQWC149 

      PHPKTH=C2/OM2*KTH-C/OM*PNPVTH                                     BQWC150 

      PHPKPH=C2/OM2*KPH-C/OM*PNPVPH                                     BQWC151 

      KPHPK=N2                                                          BQWC152 

      RETURN                                                            BQWC153 

         END                                                            BQWC154 

      SUBROUTINE BQWFNC                                                 BQNC001 

C********* CALCULATES A HAMILTONIAN H                                   BQNC002 

C********* (= BOOKER QUARTIC FOR VERTICAL INCIDENCE, S=0, C=1)          BQNC003 

C********* AND ITS PARTIAL DERIVATIVES WITH RESPECT TO                  BQNC004 

C********* TIME, R, THETA, PHI, OMEGA, KR, KTHETA, AND KPHI.            BQNC005 

C********* WITH FIELD, NO COLLISIONS                                    BQNC006 

      COMMON /CONST/ PI,PIT2,PID2,DEGS,RADIAN,K,C,LOGTEN                BQNC007 

      COMMON /RIN/ MODRIN(3),COLL,FIELD,SPACE,KAY2,KAY2I,               BQNC008 

     1             H,HI,PHPT,PHPTI,PHPR,PHPRI,PHPTH,PHPTHI,PHPPH,PHPPHI,BQNC009 

     2             PHPOM,PHPOMI,PHPKR,PHPKRI,PHPKTH,PHPKTI,PHPKPH,PHPKPIBQNC010 

     3            ,KPHPK,KPHPKI,POLAR,POLARI,LPOLAR,LPOLRI,SGN          BQNC011 

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  BQNC012 

      COMMON /YY/ MODY,Y,PYPR,PYPTH,PYPPH,YR,PYRPR,PYRPT,PYRPP,YTH,PYTPRBQNC013 

     1           ,PYTPT,PYTPP,YPH,PYPPR,PYPPT,PYPPP                     BQNC014 

      COMMON /ZZ/ MODZ,Z(4)                                             BQNC015 

      COMMON /RK/ N,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART       BQNC016 

      COMMON R,TH,PH,KR,KTH,KPH   /WW/ ID(10),W0,W(400)                 BQNC017 

      EQUIVALENCE (RAY,W(1)),(F,W(6))                                   BQNC018 

      LOGICAL SPACE                                                     BQNC019 

      REAL N2,NNP,LPOLAR,LPOLRI,KR,KTH,KPH,K2,KDOTY,K4,KDOTY2,          BQNC020 

     1     KPHPK,KPHPKI,KAY2,KAY2I                                      BQNC021 

      DATA (MODRIN=8HBOOKER Q,8HUARTIC, ,8HS=0, C=1),(COLL=0.),         BQNC022 

     1     (FIELD=1.),(KAY2I=0.),(HI=0.),(PHPTI=0.),(PHPRI=0.),         BQNC023 

     2     (PHPTHI=0.),(PHPPHI=0.),(PHPOMI=0.),(PHPKRI=0.),(PHPKTI=0.), BQNC024 

     3     (PHPKPI=0.),(KPHPKI=0.),(POLAR=0.),(LPOLAR=0.),              BQNC025 

     4     (X=0.),(PXPR=0.),(PXPTH=0.),(PXPPH=0.),(PXPT=0.),            BQNC026 

     5     (Y=0.),(PYPR=0.),(PYPTH=0.),(PYPPH=0.),(YR=0.),(PYRPR=0.),   BQNC027 

     6     (PYRPT=0.),(PYRPP=0.),(YTH=0.),(PYTPR=0.),(PYTPT=0.),        BQNC028 

     7     (PYTPP=0.),(YPH=0.),(PYPPR=0.),(PYPPT=0.),(PYPPP=0.),        BQNC029 

     8     (MODZ=1H ),(U=1.),(U2=1.)                                    BQNC030 

      ENTRY RINDEX                                                      BQNC031 

      OM=PIT2*1.E6*F                                                    BQNC032 

      C2=C*C                                                            BQNC033 

      K2=KR*KR+KTH*KTH+KPH*KPH                                          BQNC034 

      OM2=OM*OM                                                         BQNC035 

      CALL ELECTX                                                       BQNC036 

      IF(X.LT..1) GO TO 2                                               BQNC037 

      K4=K2*K2                                                          BQNC038 

      OM4=OM2*OM2                                                       BQNC039 

      C4=C2*C2                                                          BQNC040 

      CALL MAGY                                                         BQNC041 

      Y2=Y*Y                                                            BQNC042 

      KDOTY=KR*YR+KTH*YTH+KPH*YPH                                       BQNC043 

      KDOTY2=KDOTY*KDOTY                                                BQNC044 

      UX=U-X                                                            BQNC045 

      UX2=UX*UX                                                         BQNC046 

      A=UX*U2-U*Y2                                                      BQNC047 

      B=-2.*U*UX2+Y2*(2.*U-X)                                           BQNC048 

      ALPHA=A*C4*K4+X*KDOTY2*C4*K2                                      BQNC049 

      BETA=B*C2*OM2*K2-X*KDOTY2*C2*OM2                                  BQNC050 

      GAMMA=(UX2-Y2)*UX*OM4                                             BQNC051 

      H=ALPHA+BETA+GAMMA                                                BQNC052 

      PHPX=-U2*C4*K4+KDOTY2*C4*K2+(4.*U*UX-Y2)*C2*OM2*K2-KDOTY2*C2*OM2+ BQNC053 

     1 (-3.*UX2+Y2)*OM4                                                 BQNC054 

      PHPY2=-U*C4*K4+(2.*U-X)*C2*OM2*K2-UX*OM4                          BQNC055 

      PHPKY2   =X*C2*(C2*K2-OM2)                                        BQNC056 

      PHPK2=2.*A*C4*K2+X*KDOTY2*C4+B*C2*OM2                             BQNC057 

      PHPT=PHPX*PXPT                                                    BQNC058 

      PHPR =PHPX*PXPR +PHPY2*2.*Y*PYPR +PHPKY2  *2.*KDOTY*              BQNC059 

     1 (KR*PYRPR+KTH*PYTPR+KPH*PYPPR)                                   BQNC060 

      PHPTH=PHPX*PXPTH+PHPY2*2.*Y*PYPTH+PHPKY2  *2.*KDOTY*              BQNC061 

     1 (KR*PYRPT+KTH*PYTPT+KPH*PYPPT)                                   BQNC062 

      PHPPH=PHPX*PXPPH+PHPY2*2.*Y*PYPPH+PHPKY2  *2.*KDOTY*              BQNC063 

     1 (KR*PYRPP+KTH*PYTPP+KPH*PYPPP)                                   BQNC064 

      PHPOM=(2.*BETA+4.*GAMMA)/OM                                       BQNC065 

     1 -2.*PHPX*X/OM-2.*PHPY2*Y2/OM-2.*PHPKY2  *KDOTY2/OM               BQNC066 

      PHPKR= 2.*PHPK2*KR +2.*KDOTY*PHPKY2  *YR                          BQNC067 

      PHPKTH=2.*PHPK2*KTH+2.*KDOTY*PHPKY2  *YTH                         BQNC068 

      PHPKPH=2.*PHPK2*KPH+2.*KDOTY*PHPKY2  *YPH                         BQNC069 

      KAY2 = K2   *(-BETA+RAY*SQRT(BETA**2-4.*ALPHA*GAMMA))/(2.*ALPHA)  BQNC070 

      IF(RSTART.EQ.0.) GO TO 1                                          BQNC071 

      SCALE=SQRT(KAY2/K2)                                               BQNC072 

      KR =SCALE*KR                                                      BQNC073 

      KTH=SCALE*KTH                                                     BQNC074 

      KPH=SCALE*KPH                                                     BQNC075 

  1   CONTINUE                                                          BQNC076 

      KPHPK=4.*ALPHA+2.*BETA                                            BQNC077 

      SPACE=KAY2.EQ.OM2/C2                                              BQNC078 

      POLARI=SQRT(K2)*(U+X*OM2/(C2*KAY2-OM2))/KDOTY                     BQNC079 

      LPOLRI= SQRT(Y2-KDOTY2/K2)/UX*(1.-C2*KAY2/OM2)                    BQNC080 

      RETURN                                                            BQNC081 

C          CALCULATES THE REFRACTIVE INDEX AND ITS GRADIENT USING THE   BQNC082 

C          APPLETON-HARTREE FORMULA WITH FIELD, NO COLLISIONS           BQNC083 

  2   CONTINUE                                                          BQNC084 

      VR =C/OM*KR                                                       BQNC085 

      VTH=C/OM*KTH                                                      BQNC086 

      VPH=C/OM*KPH                                                      BQNC087 

      CALL MAGY                                                         BQNC088 

      V2=VR**2+VTH**2+VPH**2                                            BQNC089 

      VDOTY=VR*YR+VTH*YTH+VPH*YPH                                       BQNC090 

      YLV=VDOTY/V2                                                      BQNC091 

      YL2=VDOTY**2/V2                                                   BQNC092 

      YT2=Y**2-YL2                                                      BQNC093 

      YT4=YT2*YT2                                                       BQNC094 

      UX=U-X                                                            BQNC095 

      UX2=UX*UX                                                         BQNC096 

      RAD=RAY*SQRT(YT4+4.*YL2*UX2)                                      BQNC097 

      D=2.*UX-YT2+RAD                                                   BQNC098 

      D2=D*D                                                            BQNC099 

      N2=1.-2.*X*UX/D                                                   BQNC100 

      PNPPS=2.*X*UX*(-1.+(YT2-2.*UX2)/RAD)/D2                           BQNC101 

      PPSPR= YL2/Y*PYPR -(VR*PYRPR+VTH*PYTPR+VPH*PYPPR)*YLV             BQNC102 

      PPSPTH=YL2/Y*PYPTH-(VR*PYRPT+VTH*PYTPT+VPH*PYPPT)*YLV             BQNC103 

      PPSPPH=YL2/Y*PYPPH-(VR*PYRPP+VTH*PYTPP+VPH*PYPPP)*YLV             BQNC104 

      PNPX=-(2.*UX2-YT2*(U-2.*X)+(YT4*(U-2.*X)+4.*YL2*UX*UX2)/RAD)/D2   BQNC105 

      PNPY=2.*X*UX*(-YT2+(YT4+2.*YL2*UX2)/RAD)/(D2*Y)                   BQNC106 

      NNP=N2-(2.*X*PNPX+Y*PNPY)                                         BQNC107 

      PNPR =PNPX*PXPR +PNPY*PYPR +PNPPS*PPSPR                           BQNC108 

      PNPTH=PNPX*PXPTH+PNPY*PYPTH+PNPPS*PPSPTH                          BQNC109 

      PNPPH=PNPX*PXPPH+PNPY*PYPPH+PNPPS*PPSPPH                          BQNC110 

      PNPVR =PNPPS*(VR *YL2-VDOTY*YR )/V2                               BQNC111 

      PNPVTH=PNPPS*(VTH*YL2-VDOTY*YTH)/V2                               BQNC112 

      PNPVPH=PNPPS*(VPH*YL2-VDOTY*YPH)/V2                               BQNC113 

      PNPT=PNPX*PXPT                                                    BQNC114 

      SPACE=N2.EQ.1.                                                    BQNC115 

      POLARI=SQRT(V2)*(YT2-RAD)/(2.*VDOTY*UX)                           BQNC116 

      GAM=(-YT2+RAD)/(2.*UX)                                            BQNC117 

      LPOLRI=X*SQRT(YT2)/(UX*(U+GAM))                                   BQNC118 

      KAY2=OM2/C2*N2                                                    BQNC119 

      IF(RSTART.EQ.0.) GO TO 3                                          BQNC120 

      SCALE=SQRT(KAY2/K2)                                               BQNC121 

      KR =SCALE*KR                                                      BQNC122 

      KTH=SCALE*KTH                                                     BQNC123 

      KPH=SCALE*KPH                                                     BQNC124 

  3   CONTINUE                                                          BQNC125 

      H=.5*(C2*K2/OM2-N2)                                               BQNC126 

      PHPT =-PNPT                                                       BQNC127 

      PHPR =-PNPR                                                       BQNC128 

      PHPTH=-PNPTH                                                      BQNC129 

      PHPPH=-PNPPH                                                      BQNC130 

      PHPOM=-NNP/OM                                                     BQNC131 

      PHPKR =C2/OM2*KR -C/OM*PNPVR                                      BQNC132 

      PHPKTH=C2/OM2*KTH-C/OM*PNPVTH                                     BQNC133 

      PHPKPH=C2/OM2*KPH-C/OM*PNPVPH                                     BQNC134 

      KPHPK=N2                                                          BQNC135 

      RETURN                                                            BQNC136 

         END                                                            BQNC137-

