C  AXIAL DOSAGE/ SEMI-CONT/ VAR SOURCE/ VAR MET ORG CGW
      SUBROUTINE DDS
      CHARACTER*1 ST,MTC,AA1,AB1
      CHARACTER*2 AGNT,WT,AA2,CST,I2
      CHARACTER*3 MUNT,REL,LOCT,SEAT,SUR,I3,I4,I5
      CHARACTER*12 ADH,RL
      COMMON NQI,QT(6),TWL(6),D(10),DL(10)
      COMMON PR(1),DXT,HT,HML,SXS,SYS,SZS,TIVCH,UT,BR,SF,TMP,ALFA,SY100,
     1BETA,SZ100,Z,RC,V,QS
      COMMON DDM(12),FP,DDM1(8),Z0,LOCT(1),SEAT,MUNT,AGNT,AA1,REL,
     $       MTC,AA2,SUR,WT,AB1,ADH,ADR,AD2
      COMMON IPR(1),ND,IPO,I2MCS,IMA,IPC,IMM,IDD,IHR,NOV,INP,MRL,ID1,
     1ID2,IDEP,IMTCH,IM,IR,IL,IRL,ISM,IVD,K33,K42
      DIMENSION SY100T(6),SZ100T(6),SY100I(6),ALFAT(6),BETAT(6),ST(11),
     1SY100C(2,3),SZ20C(2,3),BETAC(2,3),TWLS(6),AR(10),Y(10),VF(2),CI(2)
     2,S(2),RL(3),DS(51,2)
      DATA ALFAT/1.,1.,1.,.9,.8,.7/
      DATA BETAT/1.4,1.,.9,.85,.8,.75/
      DATA SY100I/9.,6.33,4.8,4.,3.,2./
      DATA SY100T/27.,19.,12.5,8.,6.1,4./
      DATA SZ100T/14.,11.,7.5,4.5,3.5,2.5/
      DATA SY100C/41.19,31.18,66.56,30.98,26.17,29.33/
      DATA SZ20C/3.,1.652,.797,1.934,.705,1.242/
      DATA BETAC/1.344,.755,1.218,.949,1.182,1./
      DATA ST/'A','B','C','D','E','F','N','I','U','S','W'/
      DATA VF,CI,S/.13,.01,.454,.262,2.38,2.24/
      DATA I4,I5/'EDI','EDS'/
      DATA RL/'NO EFFECTS  ','NO DEATHS   ','1% LETHALITY'/
      IF (ISM.GT.1) GO TO 2
      DO 1 I=1,51   
 1    DS(I,2)=(DS(I,1)+DS(I,2))*ISM 
      RETURN
 2    IF (IMTCH) 3,3,7  
 3    DX=DXT
      DO 4 I=1,ND   
      AR(I)=0.
 4    DL(I)=ALOG(D(I))  
      K=1   
      DLDG=1.   
      DLDGS=1.  
      X=0.  
      TWHML=HML+HML 
      IC=ND 
      MXLF=0
      XX=0. 
      XY=0. 
      XZ=0.
      XS=0. 
      DPMX=0.   
      IRTP=5-IR 
      DPL=-87.5 
      DPLS=-87.5
      I2MC=I2MCS
      TQS=0.
      QRMX=0.   
      SDEPL=1.  
      QTTL=0.   
      DO 5 I=1,NQI  
      QTTL=QTTL+QT(I)
      TWLS(I)=TWL(I)*60.
      QR=QT(I)/(TWLS(I)-TQS)
      IF (QRMX.GT.QR) GO TO 5   
      QRMX=QR   
      TQ2=TWLS(I)   
      TQ1=TQS   
 5    TQS=TWLS(I)   
      TOMXS=(TQ2-TQ1)/2.+TQ1
      IF (IDEP) 6,7,7   
 6    CALL HD42 
      IF (QT(1).EQ.0.) RETURN
      GO TO 11  
 7    IMTCH=0   
      IF (IM.EQ.9.OR.IM.EQ.11) GO TO 11 
      IF (IM.LE.6) GO TO 9  
      IF (IM.LE.8) GO TO 10 
      WRITE (*,8)   
 8    FORMAT (' MET CODE NOT DEFINED')  
      RETURN
 9    IF (TOMXS.LE.2.4) SY100=SY100I(IM)
      IF (TOMXS.GT.2.4) SY100=SY100T(IM)
      SZ100=SZ100T(IM)  
      ALFA=ALFAT(IM)
      BETA=BETAT(IM)
      GO TO 11  
 10   IF (UT.LE.2.235) I=1  
      IF (UT.GT.2.235.AND.UT.LE.4.47) I=2   
      IF (UT.GT.4.47) I=3   
      MC=IM-6   
      SY100=SY100C(MC,I)
      SZ100=SZ20C(MC,I)*5**BETAC(MC,I)  
      ALFA=.5   
      BETA=BETAC(MC,I)  
 11   U=UT
      XCH=1.E36
      IF (TIVCH.NE.1.E36) XCH=U*TIVCH*60.+X
      IF (X.LE.0.) X=1.
      IF (SXS.GT.0.) XX=(SXS/.1522)**(1./.9294)-X
      IF (SYS.GT.0.) XY=100.*(SYS/SY100)**(1./ALFA)-X
      IF (SZS.GT.0.) XZ=100.*(SZS/SZ100)**(1./BETA)-X
      WRITE (*,12)
 12   FORMAT (/4X,'Q(MG)',3X,'TS(MIN)',2X,'HTS(M)',3X,'HML(M)',3X,
     1'SXS(M)',3X,'SYS(M)',3X,'SZS(M)')
      CST=ST(IM)
      IF(IM.EQ.11) CST=WT
      WRITE(*,13) QT(1),TWL(1),HT,HML,SXS,SYS,SZS,CST
 13   FORMAT (1X,1PE9.3,6(1X,E8.2),2X,A2)
      IF (NQI.EQ.1) GO TO 15
      DO 14 I=2,NQI
      WRITE(*,13) QT(I),TWL(I)
 14   CONTINUE
 15   IF (IVD.EQ.1) CALL VDPL(QTTL,SY,SZ,D(1),X,DX,DXA,XS,XCH,DP,
     1DPA,SDEPL,1)
      IF(IDEP.GE.1) WRITE(*,18) BR,SF
      IF(I2MC.GT.0) WRITE(*,103)
      IF(ISM.EQ.2) WRITE(*,105)
      IF (IPO.LT.1.OR.IPO.EQ.4) GO TO 25
      WRITE(*,16)
 16   FORMAT (/3X,'ALFA',5X,'SYR',4X,'BETA',5X,'SZR',5X,'XY',6X,'XZ')   
      WRITE(*,17) ALFA,SY100,BETA,SZ100,XY,XZ,UT
 17   FORMAT (1X,7(1PE8.2)) 
 18   FORMAT (/' BRT=',F4.0,'   SKF=',1PE8.2)   
      I2='DP'   
      IF (IDEP.GT.0) I2='ED'
      I3='   '  
      IF (I2MC.GT.0) I3='2MC'   
      IF (IPO.NE.3) GO TO 21
      WRITE(*,19) (D(I),I=1,ND) 
 19   FORMAT(/6X,'DOSAGE CONTOURS',10F5.0)  
      WRITE(*,20) I2,I3
 20   FORMAT(/7X,'X',7X,A2,A1,3X,'CONTOUR HALF-WIDTH')  
      GO TO 25  
 21   IF (IDEP.EQ.0) WRITE(*,24) I2,I3,I3   
      IF (IDEP.GT.0) WRITE(*,24) I2,I3,I3,I4,I5 
 24   FORMAT(/7X,'X',7X,A2,A1,7X,'RF',7X,A3,7X,A3,7X,A3)
      IF (ISM.EQ.2) WRITE(*,100)
 25   IF (X.GT.XCH) X=XCH   
      B=X/U 
      IIND=0
      IF (IRTP.GT.0) GO TO 26
      CALL PLRS (U,TMP,PMM,IL,IM,IR,X,HT,HML,IPC,IRTP)
 26   DXA=DX
      IF (X.GE.(DX*10.)) DX=DX*10.
      DXA=DXA+DX
      SX=.1522*(XX+X)**.9294
      SY=SY100*((X+XY)*.01)**ALFA
      SZ=SZ100*((X+XZ)*.01)**BETA
      IF (IVD.EQ.1.AND.X.EQ.1) CALL VDPL(QTTL,SY,SZ,D(1),X,DX,DXA,XS,
     1XCH,DP,DPA,SDEPL,2)
      DP=0.
      IF (I2MC.EQ.0) GO TO 27
      A=1./(1.414*SX)
      G=.7979*SX/U
      H=.5/(SX*SX)
      TSX=SX/U
      IT=0
      TT0=B+TOMXS
      TT1=TT0-60.
      TT2=TT1+120.
 27   TSZSQ=SZ*SZ*2.
      HPZ=HT+Z
      HMZ=HT-Z
      VT=V*X/U
      FAC=0.
      HML2=1.E36
      ARG=(HMZ-VT)**2/TSZSQ
      IF (ARG.LT.87.) FAC=FAC+EXP(-ARG) 
      ARG=(HPZ-VT)**2/TSZSQ 
      IF (ARG.LT.87.) FAC=FAC+RC/EXP(ARG)   
      ZFAC=0.   
      IF (HML.GT.1.E10.OR.MXLF.EQ.1) GO TO 30   
      DO 28 JJ=1,20
      SMHML=TWHML*JJ
      ARG=(SMHML-HPZ+VT)**2/TSZSQ   
      IF (NOV.EQ.4) WRITE(*,*) ' ARG',ARG   
      IF (ARG.GT.87.) GO TO 29
      ZFAC=ZFAC+RC**(JJ-1)/EXP(ARG) 
      ARG=(SMHML-HMZ+VT)**2/TSZSQ   
      IF (ARG.LT.87.) ZFAC=ZFAC+RC**JJ/EXP(ARG) 
      ARG=(SMHML+HMZ-VT)**2/TSZSQ   
      IF (ARG.LT.87.) ZFAC=ZFAC+RC**JJ/EXP(ARG) 
      ARG=(SMHML+HPZ-VT)**2/TSZSQ
 28   IF (ARG.LT.87.) ZFAC=ZFAC+RC**(JJ+1)/EXP(ARG) 
 29   IF ((FAC+ZFAC).NE.0.) HML2=(2.5066283*SZ)/(FAC+ZFAC)  
      IF (HML.GT.HML2.AND.(HML-HML2).LT.1) MXLF=1
 30   RF=(FAC+ZFAC)/2.
      IF (IIND) 33,33,31
 31   TT0=B-TSX-TSX
      TT1=TT0
      TT3=B+TWLS(NQI)+TSX+TSX   
      TT3S=TT3  
      TT3=TT3-120.  
      DTT=(TT3-TT0)/30. 
      IF (DTT.LT.10.) DTT=10.   
      TT2=TT1+120.
      DOS=0.
      IF (DTT.LT.120.) GO TO 33 
 32   TT2=TT1+DTT   
 33   TWLT=0.
      DLDP=0.   
      DO 63 ITWL=1,NQI  
      TS=TWLS(ITWL)-TWLT
      IF (MXLF) 35,35,34
 34   C=QT(ITWL)/(300.79539*U*TS*SY*HML)
      GO TO 36
 35   C=QT(ITWL)*RF/(376.991*U*TS*SY*SZ)
 36   IF (I2MC.GT.0) GO TO 37   
      DP=DP+C*2.*TS
      GO TO 63
 37   T1=TT1-TWLT   
      IF (T1) 38,39,39
 38   T1=0.
 39   T2=TT2-TWLT   
      IF (T2) 63,63,40  
 40   I=3   
      PT1=T1
      PT2=T2
      P=T2-TS
      R=T1-TS
      IF (P) 42,42,41   
 41   IF (R) 44,43,43   
 42   I=1
      GO TO 45  
 43   I=2   
      GO TO 45  
 44   T2=TS 
      R=0.  
 45   E=X-U*T2
      EE=E*E
      F=X-U*T1
      FF=F*F
      W=(T2-B)*ERFF(E*A)
      CC=(T1-B)*ERFF(F*A)
      ARG1=EE*H
      ARG2=FF*H
      IF (ARG1-82.6) 47,47,46
 46   ARG1=0.
      GO TO 49
 47   ARG1=EXP(-ARG1)
 48   ARG2=EXP(-ARG2)
      GO TO 51
 49   IF (ARG2-82.6) 48,48,50
 50   ARG2=0.
 51   DD=G*(ARG1-ARG2)
      GO TO (52,54,52,54), I
 52   QQ=(T2-T1)*ERFF(X*A)
      XTERM=C*(QQ-W+CC+DD)
      IF (I-3) 62,53,94
 53   PART=XTERM
      I=4
      T2=PT2
      T1=TS
      GO TO 45
 54   ARG3=X-U*P
      ARG4=X-U*R
      HH=(T2-B-TS)*ERFF(ARG3*A)
      PP=(T1-B-TS)*ERFF(ARG4*A)
      ARG3=(ARG3**2)*H
      ARG4=(ARG4**2)*H  
      IF (ARG3-82.6) 56,56,55   
 55   ARG3=0.   
      GO TO 58  
 56   ARG3=EXP(-(ARG3)) 
 57   ARG4=EXP(-(ARG4))
      GO TO 60
 58   IF (ARG4-82.6) 57,57,59   
 59   ARG4=0.   
 60   GG=G*(ARG3-ARG4)  
      XTERM=C*(HH-PP-W+CC-GG+DD)
      IF (I-4) 62,61,94 
 61   T1=PT1
      I=3   
      XTERM=XTERM+PART  
 62   DLDP=DLDP+XTERM   
      TWLT=TWLS(ITWL)   
 63   CONTINUE  
      DPA=DP
      IF (I2MC.EQ.0) GO TO 75   
      IF (IIND) 64,64,65
 64   IIND=1
      TDLDP=DLDP+0.
      GO TO 31  
 65   DP=DP+DLDP
      IF (IT) 66,66,67  
 66   IT=1  
      TT2C=TT2-TT0  
      GO TO 71
 67   TT10=TT1-TT0
      TT20=TT2-TT0  
      TT21=TT2-TT1  
      TT121=TT2C+TT2-TT1
      TT10P=TT2C**.274  
      TT20P=TT121**.274 
      TT21P=TT21**.274  
      DELDC=.2696*(D0*(TT20P-TT10P))
      IF (DELDC-DLDP) 69,68,70  
 68   TT2C=TT121
      GO TO 71  
 69   DPC=DPS+DELDC 
      TT2C=(((DPC*TT20P)+((DLDP-DELDC)*TT21P))/DP)**3.6496
      GO TO 71  
 70   TT2C=((DLDP/TT21)*TT20*TT20P)+((DPS-(DLDP/TT21)*TT10)*TT10P)  
      TT2C=(TT2C/DP)**3.6496
 71   DPS=DP
      D0=DP/(.2696*(TT2C**.274))
      IF (DOS-D0) 72,72,73  
 72   TT2CS=TT2C+0. 
      DOS=D0
 73   TT1=TT2   
      IF (TT2-TT3S) 32,74,74
 74   DLDG=.2696*TT2CS**.274
      IF (DLDG.LT.1.) DLDG=1.
      DPA=DP
      DP=DP/DLDG
      IF (DP.LT.TDLDP) DP=TDLDP
 75   IF (IDEP.LT.1) GO TO 76
      DPA=DPA*VF(IDEP)
      DOSI=DP*VF(IDEP)
      DOSIS=QT(1)*SF*CI(IDEP)*(U/X)**S(IDEP)*1000/BR
      DP=(DOSI+DOSIS)
 76   IF (IVD.EQ.1) CALL VDPL(QTTL,SY,SZ,D(1),X,DX,DXA,XS,
     1XCH,DP,DPA,SDEPL,3)
      IF (K.LT.52) DS(K,1)=DP
      IF (ISM.EQ.2) DP=DS(K,2)+DP   
      IF (DP.LE.0.) GO TO 83
      DPL=ALOG(DP)
      IF (XS.EQ.0..OR.DPL.GE.DPLS) GO TO 83
      IND=1 
      IF (IC-MRL) 84,84,77  
 77   IF (DPMX.LT.D(IC)) GO TO 78   
      IF (DP.GT.D(IC).OR.DPLS.LT.DL(IC)) GO TO 84   
      XL=ALOG(X)
      XLS=ALOG(XS)  
      XINT=EXP(XLS+(XL-XLS)*(DL(IC)-DPLS)/(DPL-DPLS))
      IF (IVD.EQ.1) WRITE(*,'(F5.2)') SDEPL 
      IF (IRL.EQ.0) WRITE(*,101) XINT,D(IC) 
      IF (IRL.EQ.1) WRITE(*,102) XINT,RL(IC)
 78   IF (D(IC).EQ.1.E36) WRITE(*,79)   
 79   FORMAT (' RESP NOT DEF')  
      IC=IC-1   
      IF (IC.GT.1.OR.I2MC.LT.2) GO TO 81
      DPLS=ALOG(EXP(DPLS)*DLDGS)
      DPL=ALOG(DP*DLDG) 
      I2MC=0
      IF (IC.GT.MRL) WRITE(*,80)
 80   FORMAT (' W/O 2-MINUTE CORRECTION')   
 81   IF (IC-MRL) 84,84,77  
 83   IND=0 
      I2MC=I2MCS
      IF (DP.LT.DPMX) GO TO 84
      DPMX=DP   
      XDMX=X
 84   IPOP1=IPO+1   
      GO TO (92,85,87,88,92), IPOP1 
 85   IF (RF.EQ.1..OR.MXLF.EQ.1) GO TO 86   
      WRITE(*,98) X,DP,RF   
      GO TO 92  
 86   WRITE(*,98) X,DP  
      GO TO 92  
 87   IF (IDEP.EQ.0) WRITE(*,98) X,DP,RF,DLDG   
      IF (IDEP.GE.1) WRITE(*,98) X,DP,RF,DLDG,DOSI,DOSIS
      GO TO 92  
 88   IF (DP.LT.D(1)) GO TO 92  
      DO 89 IY=1,ND 
      ARG=DPL-DL(IY)
      IF (ARG.LT.0.) GO TO 90   
      Y(IY)=1.41421*SY*SQRT(ARG)
      IF (X.NE.XCH) AR(IY)=AR(IY)+Y(IY)*DXA
      IF (X.EQ.XCH) AR(IY)=AR(IY)-Y(IY)*(XS+DX-XCH) 
 89   CONTINUE  
      IY=ND+1   
 90   IY=IY-1
      WRITE(*,91) X,DP,(Y(I),I=1,IY)
 91   FORMAT (1X,F10.0,1PE10.3,0P10F5.0)
 92   XS=X  
      DPLS=DPL  
      DLDGS=DLDG
      IF (X.NE.XCH) GO TO 93
      IMTCH=1   
      SXS=SX
      SYS=SY
      SZS=SZ
      RETURN
 93   IF (X.EQ.1..AND.DX.GT.1.) X=0.
      X=X+DX
      K=K+1 
      IF (X.GT.9.E6) RETURN 
      IF (IND.EQ.0.OR.DP.GT.D(MRL+1).OR.IRTP.LT.1) GO TO 25 
      IF (DPMX.LT.D(MRL+1)) WRITE(*,104) XDMX,DPMX  
 94   WRITE(*,*)
      IF (IPO.EQ.3) WRITE(*,96)
 96   FORMAT (4X,'DOSAGE',5X,'AREA')
      IF (IPO.EQ.3) WRITE(*,99) (D(I),AR(I),I=ND,1,-1)
      IF (IVD.EQ.1.AND.IPO.EQ.4) CALL VDPL(QTTL,SY,SZ,D(1),X,DX,DXA,XS,
     1XCH,DP,DPA,SDEPL,4)
      IF (K.GT.51) RETURN
      DO 97 I=K,51
 97   DS(I,1)=0.
      RETURN
 98   FORMAT (1X,F10.0,1P5E10.3)
 99   FORMAT (1X,1P2E10.3)
 100  FORMAT (15X,'SUM')
 101  FORMAT (/1X,F10.0,' (M) IS DISTANCE TO ',E10.3,' (MG-MIN/M^3)')
 102  FORMAT (/1X,F10.0,' (M) IS DISTANCE TO ',A12/)
 103  FORMAT (/' W/2-MINUTE CORRECTION')
 104  FORMAT (/' MINIMUM DOSAGE NOT ATTAINED',//4X,'XDMAX',9X,'DMAX',
     $       /F10.0,6X,E10.3)
 105  FORMAT (/' DOSAGE IS BEING SUMMED')
      END
