home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / FORTRAN / D2SOURCE.ZIP / 37G3.FOR next >
Encoding:
Text File  |  1985-12-20  |  11.6 KB  |  409 lines
  1. C  AXIAL DOSAGE/ SEMI-CONT/ VAR SOURCE/ VAR MET ORG CGW
  2.       SUBROUTINE DDS
  3.       CHARACTER*1 ST,MTC,AA1,AB1
  4.       CHARACTER*2 AGNT,WT,AA2,CST,I2
  5.       CHARACTER*3 MUNT,REL,LOCT,SEAT,SUR,I3,I4,I5
  6.       CHARACTER*12 ADH,RL
  7.       COMMON NQI,QT(6),TWL(6),D(10),DL(10)
  8.       COMMON PR(1),DXT,HT,HML,SXS,SYS,SZS,TIVCH,UT,BR,SF,TMP,ALFA,SY100,
  9.      1BETA,SZ100,Z,RC,V,QS
  10.       COMMON DDM(12),FP,DDM1(8),Z0,LOCT(1),SEAT,MUNT,AGNT,AA1,REL,
  11.      $       MTC,AA2,SUR,WT,AB1,ADH,ADR,AD2
  12.       COMMON IPR(1),ND,IPO,I2MCS,IMA,IPC,IMM,IDD,IHR,NOV,INP,MRL,ID1,
  13.      1ID2,IDEP,IMTCH,IM,IR,IL,IRL,ISM,IVD,K33,K42
  14.       DIMENSION SY100T(6),SZ100T(6),SY100I(6),ALFAT(6),BETAT(6),ST(11),
  15.      1SY100C(2,3),SZ20C(2,3),BETAC(2,3),TWLS(6),AR(10),Y(10),VF(2),CI(2)
  16.      2,S(2),RL(3),DS(51,2)
  17.       DATA ALFAT/1.,1.,1.,.9,.8,.7/
  18.       DATA BETAT/1.4,1.,.9,.85,.8,.75/
  19.       DATA SY100I/9.,6.33,4.8,4.,3.,2./
  20.       DATA SY100T/27.,19.,12.5,8.,6.1,4./
  21.       DATA SZ100T/14.,11.,7.5,4.5,3.5,2.5/
  22.       DATA SY100C/41.19,31.18,66.56,30.98,26.17,29.33/
  23.       DATA SZ20C/3.,1.652,.797,1.934,.705,1.242/
  24.       DATA BETAC/1.344,.755,1.218,.949,1.182,1./
  25.       DATA ST/'A','B','C','D','E','F','N','I','U','S','W'/
  26.       DATA VF,CI,S/.13,.01,.454,.262,2.38,2.24/
  27.       DATA I4,I5/'EDI','EDS'/
  28.       DATA RL/'NO EFFECTS  ','NO DEATHS   ','1% LETHALITY'/
  29.       IF (ISM.GT.1) GO TO 2
  30.       DO 1 I=1,51   
  31.  1    DS(I,2)=(DS(I,1)+DS(I,2))*ISM 
  32.       RETURN
  33.  2    IF (IMTCH) 3,3,7  
  34.  3    DX=DXT
  35.       DO 4 I=1,ND   
  36.       AR(I)=0.
  37.  4    DL(I)=ALOG(D(I))  
  38.       K=1   
  39.       DLDG=1.   
  40.       DLDGS=1.  
  41.       X=0.  
  42.       TWHML=HML+HML 
  43.       IC=ND 
  44.       MXLF=0
  45.       XX=0. 
  46.       XY=0. 
  47.       XZ=0.
  48.       XS=0. 
  49.       DPMX=0.   
  50.       IRTP=5-IR 
  51.       DPL=-87.5 
  52.       DPLS=-87.5
  53.       I2MC=I2MCS
  54.       TQS=0.
  55.       QRMX=0.   
  56.       SDEPL=1.  
  57.       QTTL=0.   
  58.       DO 5 I=1,NQI  
  59.       QTTL=QTTL+QT(I)
  60.       TWLS(I)=TWL(I)*60.
  61.       QR=QT(I)/(TWLS(I)-TQS)
  62.       IF (QRMX.GT.QR) GO TO 5   
  63.       QRMX=QR   
  64.       TQ2=TWLS(I)   
  65.       TQ1=TQS   
  66.  5    TQS=TWLS(I)   
  67.       TOMXS=(TQ2-TQ1)/2.+TQ1
  68.       IF (IDEP) 6,7,7   
  69.  6    CALL HD42 
  70.       IF (QT(1).EQ.0.) RETURN
  71.       GO TO 11  
  72.  7    IMTCH=0   
  73.       IF (IM.EQ.9.OR.IM.EQ.11) GO TO 11 
  74.       IF (IM.LE.6) GO TO 9  
  75.       IF (IM.LE.8) GO TO 10 
  76.       WRITE (*,8)   
  77.  8    FORMAT (' MET CODE NOT DEFINED')  
  78.       RETURN
  79.  9    IF (TOMXS.LE.2.4) SY100=SY100I(IM)
  80.       IF (TOMXS.GT.2.4) SY100=SY100T(IM)
  81.       SZ100=SZ100T(IM)  
  82.       ALFA=ALFAT(IM)
  83.       BETA=BETAT(IM)
  84.       GO TO 11  
  85.  10   IF (UT.LE.2.235) I=1  
  86.       IF (UT.GT.2.235.AND.UT.LE.4.47) I=2   
  87.       IF (UT.GT.4.47) I=3   
  88.       MC=IM-6   
  89.       SY100=SY100C(MC,I)
  90.       SZ100=SZ20C(MC,I)*5**BETAC(MC,I)  
  91.       ALFA=.5   
  92.       BETA=BETAC(MC,I)  
  93.  11   U=UT
  94.       XCH=1.E36
  95.       IF (TIVCH.NE.1.E36) XCH=U*TIVCH*60.+X
  96.       IF (X.LE.0.) X=1.
  97.       IF (SXS.GT.0.) XX=(SXS/.1522)**(1./.9294)-X
  98.       IF (SYS.GT.0.) XY=100.*(SYS/SY100)**(1./ALFA)-X
  99.       IF (SZS.GT.0.) XZ=100.*(SZS/SZ100)**(1./BETA)-X
  100.       WRITE (*,12)
  101.  12   FORMAT (/4X,'Q(MG)',3X,'TS(MIN)',2X,'HTS(M)',3X,'HML(M)',3X,
  102.      1'SXS(M)',3X,'SYS(M)',3X,'SZS(M)')
  103.       CST=ST(IM)
  104.       IF(IM.EQ.11) CST=WT
  105.       WRITE(*,13) QT(1),TWL(1),HT,HML,SXS,SYS,SZS,CST
  106.  13   FORMAT (1X,1PE9.3,6(1X,E8.2),2X,A2)
  107.       IF (NQI.EQ.1) GO TO 15
  108.       DO 14 I=2,NQI
  109.       WRITE(*,13) QT(I),TWL(I)
  110.  14   CONTINUE
  111.  15   IF (IVD.EQ.1) CALL VDPL(QTTL,SY,SZ,D(1),X,DX,DXA,XS,XCH,DP,
  112.      1DPA,SDEPL,1)
  113.       IF(IDEP.GE.1) WRITE(*,18) BR,SF
  114.       IF(I2MC.GT.0) WRITE(*,103)
  115.       IF(ISM.EQ.2) WRITE(*,105)
  116.       IF (IPO.LT.1.OR.IPO.EQ.4) GO TO 25
  117.       WRITE(*,16)
  118.  16   FORMAT (/3X,'ALFA',5X,'SYR',4X,'BETA',5X,'SZR',5X,'XY',6X,'XZ')   
  119.       WRITE(*,17) ALFA,SY100,BETA,SZ100,XY,XZ,UT
  120.  17   FORMAT (1X,7(1PE8.2)) 
  121.  18   FORMAT (/' BRT=',F4.0,'   SKF=',1PE8.2)   
  122.       I2='DP'   
  123.       IF (IDEP.GT.0) I2='ED'
  124.       I3='   '  
  125.       IF (I2MC.GT.0) I3='2MC'   
  126.       IF (IPO.NE.3) GO TO 21
  127.       WRITE(*,19) (D(I),I=1,ND) 
  128.  19   FORMAT(/6X,'DOSAGE CONTOURS',10F5.0)  
  129.       WRITE(*,20) I2,I3
  130.  20   FORMAT(/7X,'X',7X,A2,A1,3X,'CONTOUR HALF-WIDTH')  
  131.       GO TO 25  
  132.  21   IF (IDEP.EQ.0) WRITE(*,24) I2,I3,I3   
  133.       IF (IDEP.GT.0) WRITE(*,24) I2,I3,I3,I4,I5 
  134.  24   FORMAT(/7X,'X',7X,A2,A1,7X,'RF',7X,A3,7X,A3,7X,A3)
  135.       IF (ISM.EQ.2) WRITE(*,100)
  136.  25   IF (X.GT.XCH) X=XCH   
  137.       B=X/U 
  138.       IIND=0
  139.       IF (IRTP.GT.0) GO TO 26
  140.       CALL PLRS (U,TMP,PMM,IL,IM,IR,X,HT,HML,IPC,IRTP)
  141.  26   DXA=DX
  142.       IF (X.GE.(DX*10.)) DX=DX*10.
  143.       DXA=DXA+DX
  144.       SX=.1522*(XX+X)**.9294
  145.       SY=SY100*((X+XY)*.01)**ALFA
  146.       SZ=SZ100*((X+XZ)*.01)**BETA
  147.       IF (IVD.EQ.1.AND.X.EQ.1) CALL VDPL(QTTL,SY,SZ,D(1),X,DX,DXA,XS,
  148.      1XCH,DP,DPA,SDEPL,2)
  149.       DP=0.
  150.       IF (I2MC.EQ.0) GO TO 27
  151.       A=1./(1.414*SX)
  152.       G=.7979*SX/U
  153.       H=.5/(SX*SX)
  154.       TSX=SX/U
  155.       IT=0
  156.       TT0=B+TOMXS
  157.       TT1=TT0-60.
  158.       TT2=TT1+120.
  159.  27   TSZSQ=SZ*SZ*2.
  160.       HPZ=HT+Z
  161.       HMZ=HT-Z
  162.       VT=V*X/U
  163.       FAC=0.
  164.       HML2=1.E36
  165.       ARG=(HMZ-VT)**2/TSZSQ
  166.       IF (ARG.LT.87.) FAC=FAC+EXP(-ARG) 
  167.       ARG=(HPZ-VT)**2/TSZSQ 
  168.       IF (ARG.LT.87.) FAC=FAC+RC/EXP(ARG)   
  169.       ZFAC=0.   
  170.       IF (HML.GT.1.E10.OR.MXLF.EQ.1) GO TO 30   
  171.       DO 28 JJ=1,20
  172.       SMHML=TWHML*JJ
  173.       ARG=(SMHML-HPZ+VT)**2/TSZSQ   
  174.       IF (NOV.EQ.4) WRITE(*,*) ' ARG',ARG   
  175.       IF (ARG.GT.87.) GO TO 29
  176.       ZFAC=ZFAC+RC**(JJ-1)/EXP(ARG) 
  177.       ARG=(SMHML-HMZ+VT)**2/TSZSQ   
  178.       IF (ARG.LT.87.) ZFAC=ZFAC+RC**JJ/EXP(ARG) 
  179.       ARG=(SMHML+HMZ-VT)**2/TSZSQ   
  180.       IF (ARG.LT.87.) ZFAC=ZFAC+RC**JJ/EXP(ARG) 
  181.       ARG=(SMHML+HPZ-VT)**2/TSZSQ
  182.  28   IF (ARG.LT.87.) ZFAC=ZFAC+RC**(JJ+1)/EXP(ARG) 
  183.  29   IF ((FAC+ZFAC).NE.0.) HML2=(2.5066283*SZ)/(FAC+ZFAC)  
  184.       IF (HML.GT.HML2.AND.(HML-HML2).LT.1) MXLF=1
  185.  30   RF=(FAC+ZFAC)/2.
  186.       IF (IIND) 33,33,31
  187.  31   TT0=B-TSX-TSX
  188.       TT1=TT0
  189.       TT3=B+TWLS(NQI)+TSX+TSX   
  190.       TT3S=TT3  
  191.       TT3=TT3-120.  
  192.       DTT=(TT3-TT0)/30. 
  193.       IF (DTT.LT.10.) DTT=10.   
  194.       TT2=TT1+120.
  195.       DOS=0.
  196.       IF (DTT.LT.120.) GO TO 33 
  197.  32   TT2=TT1+DTT   
  198.  33   TWLT=0.
  199.       DLDP=0.   
  200.       DO 63 ITWL=1,NQI  
  201.       TS=TWLS(ITWL)-TWLT
  202.       IF (MXLF) 35,35,34
  203.  34   C=QT(ITWL)/(300.79539*U*TS*SY*HML)
  204.       GO TO 36
  205.  35   C=QT(ITWL)*RF/(376.991*U*TS*SY*SZ)
  206.  36   IF (I2MC.GT.0) GO TO 37   
  207.       DP=DP+C*2.*TS
  208.       GO TO 63
  209.  37   T1=TT1-TWLT   
  210.       IF (T1) 38,39,39
  211.  38   T1=0.
  212.  39   T2=TT2-TWLT   
  213.       IF (T2) 63,63,40  
  214.  40   I=3   
  215.       PT1=T1
  216.       PT2=T2
  217.       P=T2-TS
  218.       R=T1-TS
  219.       IF (P) 42,42,41   
  220.  41   IF (R) 44,43,43   
  221.  42   I=1
  222.       GO TO 45  
  223.  43   I=2   
  224.       GO TO 45  
  225.  44   T2=TS 
  226.       R=0.  
  227.  45   E=X-U*T2
  228.       EE=E*E
  229.       F=X-U*T1
  230.       FF=F*F
  231.       W=(T2-B)*ERFF(E*A)
  232.       CC=(T1-B)*ERFF(F*A)
  233.       ARG1=EE*H
  234.       ARG2=FF*H
  235.       IF (ARG1-82.6) 47,47,46
  236.  46   ARG1=0.
  237.       GO TO 49
  238.  47   ARG1=EXP(-ARG1)
  239.  48   ARG2=EXP(-ARG2)
  240.       GO TO 51
  241.  49   IF (ARG2-82.6) 48,48,50
  242.  50   ARG2=0.
  243.  51   DD=G*(ARG1-ARG2)
  244.       GO TO (52,54,52,54), I
  245.  52   QQ=(T2-T1)*ERFF(X*A)
  246.       XTERM=C*(QQ-W+CC+DD)
  247.       IF (I-3) 62,53,94
  248.  53   PART=XTERM
  249.       I=4
  250.       T2=PT2
  251.       T1=TS
  252.       GO TO 45
  253.  54   ARG3=X-U*P
  254.       ARG4=X-U*R
  255.       HH=(T2-B-TS)*ERFF(ARG3*A)
  256.       PP=(T1-B-TS)*ERFF(ARG4*A)
  257.       ARG3=(ARG3**2)*H
  258.       ARG4=(ARG4**2)*H  
  259.       IF (ARG3-82.6) 56,56,55   
  260.  55   ARG3=0.   
  261.       GO TO 58  
  262.  56   ARG3=EXP(-(ARG3)) 
  263.  57   ARG4=EXP(-(ARG4))
  264.       GO TO 60
  265.  58   IF (ARG4-82.6) 57,57,59   
  266.  59   ARG4=0.   
  267.  60   GG=G*(ARG3-ARG4)  
  268.       XTERM=C*(HH-PP-W+CC-GG+DD)
  269.       IF (I-4) 62,61,94 
  270.  61   T1=PT1
  271.       I=3   
  272.       XTERM=XTERM+PART  
  273.  62   DLDP=DLDP+XTERM   
  274.       TWLT=TWLS(ITWL)   
  275.  63   CONTINUE  
  276.       DPA=DP
  277.       IF (I2MC.EQ.0) GO TO 75   
  278.       IF (IIND) 64,64,65
  279.  64   IIND=1
  280.       TDLDP=DLDP+0.
  281.       GO TO 31  
  282.  65   DP=DP+DLDP
  283.       IF (IT) 66,66,67  
  284.  66   IT=1  
  285.       TT2C=TT2-TT0  
  286.       GO TO 71
  287.  67   TT10=TT1-TT0
  288.       TT20=TT2-TT0  
  289.       TT21=TT2-TT1  
  290.       TT121=TT2C+TT2-TT1
  291.       TT10P=TT2C**.274  
  292.       TT20P=TT121**.274 
  293.       TT21P=TT21**.274  
  294.       DELDC=.2696*(D0*(TT20P-TT10P))
  295.       IF (DELDC-DLDP) 69,68,70  
  296.  68   TT2C=TT121
  297.       GO TO 71  
  298.  69   DPC=DPS+DELDC 
  299.       TT2C=(((DPC*TT20P)+((DLDP-DELDC)*TT21P))/DP)**3.6496
  300.       GO TO 71  
  301.  70   TT2C=((DLDP/TT21)*TT20*TT20P)+((DPS-(DLDP/TT21)*TT10)*TT10P)  
  302.       TT2C=(TT2C/DP)**3.6496
  303.  71   DPS=DP
  304.       D0=DP/(.2696*(TT2C**.274))
  305.       IF (DOS-D0) 72,72,73  
  306.  72   TT2CS=TT2C+0. 
  307.       DOS=D0
  308.  73   TT1=TT2   
  309.       IF (TT2-TT3S) 32,74,74
  310.  74   DLDG=.2696*TT2CS**.274
  311.       IF (DLDG.LT.1.) DLDG=1.
  312.       DPA=DP
  313.       DP=DP/DLDG
  314.       IF (DP.LT.TDLDP) DP=TDLDP
  315.  75   IF (IDEP.LT.1) GO TO 76
  316.       DPA=DPA*VF(IDEP)
  317.       DOSI=DP*VF(IDEP)
  318.       DOSIS=QT(1)*SF*CI(IDEP)*(U/X)**S(IDEP)*1000/BR
  319.       DP=(DOSI+DOSIS)
  320.  76   IF (IVD.EQ.1) CALL VDPL(QTTL,SY,SZ,D(1),X,DX,DXA,XS,
  321.      1XCH,DP,DPA,SDEPL,3)
  322.       IF (K.LT.52) DS(K,1)=DP
  323.       IF (ISM.EQ.2) DP=DS(K,2)+DP   
  324.       IF (DP.LE.0.) GO TO 83
  325.       DPL=ALOG(DP)
  326.       IF (XS.EQ.0..OR.DPL.GE.DPLS) GO TO 83
  327.       IND=1 
  328.       IF (IC-MRL) 84,84,77  
  329.  77   IF (DPMX.LT.D(IC)) GO TO 78   
  330.       IF (DP.GT.D(IC).OR.DPLS.LT.DL(IC)) GO TO 84   
  331.       XL=ALOG(X)
  332.       XLS=ALOG(XS)  
  333.       XINT=EXP(XLS+(XL-XLS)*(DL(IC)-DPLS)/(DPL-DPLS))
  334.       IF (IVD.EQ.1) WRITE(*,'(F5.2)') SDEPL 
  335.       IF (IRL.EQ.0) WRITE(*,101) XINT,D(IC) 
  336.       IF (IRL.EQ.1) WRITE(*,102) XINT,RL(IC)
  337.  78   IF (D(IC).EQ.1.E36) WRITE(*,79)   
  338.  79   FORMAT (' RESP NOT DEF')  
  339.       IC=IC-1   
  340.       IF (IC.GT.1.OR.I2MC.LT.2) GO TO 81
  341.       DPLS=ALOG(EXP(DPLS)*DLDGS)
  342.       DPL=ALOG(DP*DLDG) 
  343.       I2MC=0
  344.       IF (IC.GT.MRL) WRITE(*,80)
  345.  80   FORMAT (' W/O 2-MINUTE CORRECTION')   
  346.  81   IF (IC-MRL) 84,84,77  
  347.  83   IND=0 
  348.       I2MC=I2MCS
  349.       IF (DP.LT.DPMX) GO TO 84
  350.       DPMX=DP   
  351.       XDMX=X
  352.  84   IPOP1=IPO+1   
  353.       GO TO (92,85,87,88,92), IPOP1 
  354.  85   IF (RF.EQ.1..OR.MXLF.EQ.1) GO TO 86   
  355.       WRITE(*,98) X,DP,RF   
  356.       GO TO 92  
  357.  86   WRITE(*,98) X,DP  
  358.       GO TO 92  
  359.  87   IF (IDEP.EQ.0) WRITE(*,98) X,DP,RF,DLDG   
  360.       IF (IDEP.GE.1) WRITE(*,98) X,DP,RF,DLDG,DOSI,DOSIS
  361.       GO TO 92  
  362.  88   IF (DP.LT.D(1)) GO TO 92  
  363.       DO 89 IY=1,ND 
  364.       ARG=DPL-DL(IY)
  365.       IF (ARG.LT.0.) GO TO 90   
  366.       Y(IY)=1.41421*SY*SQRT(ARG)
  367.       IF (X.NE.XCH) AR(IY)=AR(IY)+Y(IY)*DXA
  368.       IF (X.EQ.XCH) AR(IY)=AR(IY)-Y(IY)*(XS+DX-XCH) 
  369.  89   CONTINUE  
  370.       IY=ND+1   
  371.  90   IY=IY-1
  372.       WRITE(*,91) X,DP,(Y(I),I=1,IY)
  373.  91   FORMAT (1X,F10.0,1PE10.3,0P10F5.0)
  374.  92   XS=X  
  375.       DPLS=DPL  
  376.       DLDGS=DLDG
  377.       IF (X.NE.XCH) GO TO 93
  378.       IMTCH=1   
  379.       SXS=SX
  380.       SYS=SY
  381.       SZS=SZ
  382.       RETURN
  383.  93   IF (X.EQ.1..AND.DX.GT.1.) X=0.
  384.       X=X+DX
  385.       K=K+1 
  386.       IF (X.GT.9.E6) RETURN 
  387.       IF (IND.EQ.0.OR.DP.GT.D(MRL+1).OR.IRTP.LT.1) GO TO 25 
  388.       IF (DPMX.LT.D(MRL+1)) WRITE(*,104) XDMX,DPMX  
  389.  94   WRITE(*,*)
  390.       IF (IPO.EQ.3) WRITE(*,96)
  391.  96   FORMAT (4X,'DOSAGE',5X,'AREA')
  392.       IF (IPO.EQ.3) WRITE(*,99) (D(I),AR(I),I=ND,1,-1)
  393.       IF (IVD.EQ.1.AND.IPO.EQ.4) CALL VDPL(QTTL,SY,SZ,D(1),X,DX,DXA,XS,
  394.      1XCH,DP,DPA,SDEPL,4)
  395.       IF (K.GT.51) RETURN
  396.       DO 97 I=K,51
  397.  97   DS(I,1)=0.
  398.       RETURN
  399.  98   FORMAT (1X,F10.0,1P5E10.3)
  400.  99   FORMAT (1X,1P2E10.3)
  401.  100  FORMAT (15X,'SUM')
  402.  101  FORMAT (/1X,F10.0,' (M) IS DISTANCE TO ',E10.3,' (MG-MIN/M^3)')
  403.  102  FORMAT (/1X,F10.0,' (M) IS DISTANCE TO ',A12/)
  404.  103  FORMAT (/' W/2-MINUTE CORRECTION')
  405.  104  FORMAT (/' MINIMUM DOSAGE NOT ATTAINED',//4X,'XDMAX',9X,'DMAX',
  406.      $       /F10.0,6X,E10.3)
  407.  105  FORMAT (/' DOSAGE IS BEING SUMMED')
  408.       END
  409.