Old FORTRAN program for estimating worst-case rainout hazards from low-yield tactical nuclear weapons

rainout.f

**************************************************************************           
*                                                                        *
*                              RAINOUT MODEL                             *
*                                                                        *
*  This program computes the maximum dose rate. in Rads per hour,        *
*  produced by a nuclear air burst. The nuclear cloud is washed out by a *
*  simulated rain storm at a specified time t of cloud washout. The      *
*  activity grounded as a result of the radioactive cloud washout,       *
*  A(x,y,t), is used to compute the dose rate. This program will accomo- *
*  date changes in the following Input parameters: particle size dist-   *
*  ribution, weapon yield (for yields ranging from .01 to 1.0 kilotons), *
*  fractionation ratio, fallout particle density (in kilograms per cubic *
*  meter),wind velocity (in meters per second), and wind shear (in       *
*  kilometers per hour per kilometer of atmospheric altitude).           *
*                                                                        *
**************************************************************************

      REAL PI,ALPHA0,BETA0,FV,RHOP,GA,RMEAN(10),YIELD,HC,ALPHA2,ALPHA3
      REAL TEMP,MU,DELTAZ,VPRINT,X0,Y0,AXYT,DDOT(4),DOSE(4)
      REAL V,SIGMAZ(10),DELTAT,ZC(10,0:50),TIME(0:50),UTRD(4)
      REAL TC,SIGZ,SIG0,EXPO,SIGX,TSTAPR,SIGY,NUMBER,SHEARX,SHEARY,VX,VY
      REAL X(4),Y,FX,FY,EXPOFX,EXPOFY
      INTEGER I,J,K,CODE,K1

C     INPUT OF ALL PROBLEM PARAMETERS.

      PARAMETER (PI=3.14159, ALPHA0=LOG(.075), BETA0=LOG(2.0), FV=1.0, 
     &     RHOP=5500.0, GA=9.80, YIELD=1.0, VX=4.0, VY=0.0,        
     &     SHEARX=10.0, SHEARY=1.0, Y=0.0, CODE=0)

      OPEN (3,FILE='dosedata')
      REWIND 3

C     SPECIFY THE MEAN PARTICLE RADII FOR THE TEN EQUAL ACTIVITY GROUPS.

      CALL MEAN(ALPHA0,BETA0,ALPHA2,ALPHA3,RMEAN)

C     COMPUTE THE ALTITUDE OF THE STABILIZED CLOUD CENTER IN METERS.

      HC=(1609./5.28)*(8.30476298+.63632921*LOG(YIELD)-.39763749*LOG(YIE
     &LD)*LOG(YIELD)-.01364081*LOG(YIELD)**3+.00560846*LOG(YIELD)**4)
      TC=(12.*(HC/60.)*(5.28/1609.)-2.5*(HC/60.)*(5.28/1609.)*(5.28/1609
     &     .))*(1.0-.5*EXP((HC/25.)* (5.28/1609.)* (5.28/1609.)))
      SIGZ=.18*HC
      EXPO=.7+(1./3.)*LOG(YIELD/1000.)-(3.25/(4.+(LOG(YIELD/1000.)+5.4)
     &     *(LOG(YIELD/1000.)+5.4)))
      SIG0=1609.*EXP(EXPO)

      PRINT 1, YIELD
 1    FORMAT (/,"YIELD = ",F6.4," KILOTON AIRBURST")
      PRINT 2, HC
 2    FORMAT (/,"STABILIZED CLOUD CENTER HEIGHT = ",F5.0," METERS")
      PRINT 25, TC
 25   FORMAT (/,"TC - ",F5.2," HOURS")
      PRINT 26, SIG0
 26   FORMAT (/,"SIGMAO = ",F5.0," METERS")
      PRINT 27, SIGZ
 27   FORMAT (/,"SIGMAZ =",F4.0,"METERS")
      PRINT 28, VX
 28   FORMAT (/,"DOWNWIND VELOCITY = ",F4.1," METERS PER SECOND")
      PRINT 3, FV
 3    FORMAT (/,"VOLUMETRIC FRACTIONATION RATIO = ",F4.2)
      PRINT 35, RHOP
 35   FORMAT (/,"PARTICLE MASS DENSITY = ",F5.0," KILOGRAMS PER CUBIC ME
     &TER")
      PRINT 4, ALPHA0,BETA0
 4    FORMAT (/,"ALPHA0 = LN(.075) = ",F6.3,"     BETA0 = LN(2.0) = ",
     &F5.3)
      PRINT 5, ALPHA2,ALPHA3
 5    FORMAT (/,"ALPHA2 = LN(.196) = " F6.3,"    ALPHA3 = LN(.317) = " 
     &     ,F6.3,//)

C     COMPUTE THE ALTITUDE FOR THE CENTER OF EACH PARTICLE SIZE GROUP.
      DO 7 I=1,10
         ZC(I,0)=HC
         SIGMAZ(I)=.18*ZC(I,0)
 7    CONTINUE

C     COMPUTE THE FALL TIMES FOR EACH SIZE GROUP CENTER USING STOKE'S LAW
C     FALL MECHANICS.
      TIME(0)=0.0
      DELTAT=1.0
      DO 13 J=1,36
         TIME(J)=TIME(J-1)+DELTAT
         DO 12 I=1,10
            DELTAT=DELTAT*3600.0
            TEMP=288.15-.0065*ZC(I,J-1)
            MU=(1.458E-6*TEMP**1.5)/(TEMP+110.4)
            V=(2.0/9.0)*((RMEAN(I)*1E-6)**2*RHOP*GA)/MU
            DELTAZ=V*DELTAT
            ZC(I,J)=ZC(I,J-1)-DELTAZ
            IF(ZC(I,J).LT.0.0)THEN
               ZC(I,J)=0.0
               GO TO 105
            END IF
 105        DELTAT=DELTAT/3600.0
            VPRINT=V*100.0
 12      CONTINUE
 13   CONTINUE

C     COMPUTE THE VALUE OF THE DOSE RATE IN R/hr.
      IF(CODE.EQ.0)THEN
         PRINT*,' '
         PRINT*,'DOWNWIND     SIGMAX       TIME OF                    DO
     &SE TO'
         PRINT*,'DISTANCE    DISTANCE      WASH OUT     DOSE RATE     IN
     &FINITY'
         PRINT*,'  (km)        (km)         (hours)     (Rads/hr)      (
     &Rads)'
         PRINT*,' '
         GO TO 200
      ELSE IF (CODE.EQ.1)THEN
         PRINT*,' '
         PRINT*,'DOWNWIND     SIGMAX       TIME OF                    DOSE
     & TO'
         PRINT*,'DISTANCE    DISTANCE      WASH OUT     DOSE RATE     24 H
     &OURS'
         PRINT*,'  (km)        (km)         (hours)     (Rads/hr)      (Ra
     &ds)'
         PRINT*,' '
         GO TO 200
      ELSE
         PRINT*,' '
         PRINT*,'DOWNWIND     SIGMAX       TIME OF                    DOSE
     & TO'
         PRINT*,'DISTANCE    DISTANCE      WASH OUT     DOSE RATE     24 H
     &OURS'
         PRINT*,'  (km)        (km)         (hours)     (Rads/hr)      (Ra
     &ds)'
         PRINT*,' '
         GO TO 200
      END IF

 200  DO 17 J=1,36
         IF(TIME(J).LE.3.0)THEN
            TSTAR=TIME(J)
            GO TO 135
         END IF
         TSTAR=3.0
 135     SIGX=SIG0*SQRT(1.0+8.0*TSTAR/TC)
         NUMBER=SIG0*SIG0*(1.0+8.0*TSTAR/TC)+(SIGZ*SHEARY*TIME(J))*
     &        (SIGZ*SHEARY*TIME(J))
         SIGY=SQRT(NUMBER)
         X0=VX*TIME(J)*3600.0
         Y0=VY*TIME(J)*3600.0
         DO 14 K=1,4
            X(K)=X0+(K-1)*SIGX
            EXPOFX=-.5*((X(K)-X0)/SIGX)*((X(K)-X0)/SIGX)
            EXPOFY=-.5*((Y-Y0)/SIGY)*((Y-Y0)/SIGY)
            FX=(1.0/(SQRT(2.0*PI)*SIGX))*EXP(EXPOFX)
            FY=(1.0/(SQRT(2.0*PI)*SIGY))*EXP(EXPOFY)
            AXYT=(530.E6*YIELD*TIME(J)**(-1.2))*FX*FY
            DDOT(K)=14.4*AXYT
            UTRD(K)=DDOT(K)*TIME(J)**1.2
            IF(CODE.EQ.0)THEN
               DOSE(K)=5.0*UTRD(K)*TIME(J)**(-.2)
               GO TO 14
            ELSE IF(CODE.EQ.1)THEN
               DOSE(K)=5.0*UTRD(K)*(TIME(J)**(-.2)-(TIME(J)+24.0)**
     &              (-.2))
               GO TO 14
            ELSE
               DOSE(K)=5.0*UTRD(K)*(TIME(J)**(-.2)-(TIME(J)+4.0)**(-.2))
               GO TO 14
            END IF
 14      CONTINUE
         SIGX=SIGX/1000.0
         WRITE (3,136,IOSTAT=K1,ERR=17) TIME(J),DDOT(1),DDOT(2),DDOT(3),
     &        DDOT(4),DOSE(1),DOSE(2),DOSE(3),DOSE(4),SIGX
 136     FORMAT (F4.1,1X,E11.6,1X,E11.6,1X,E11.6,1X,E11.6,1X,E11.6,1X,
     &        E11.6,1X,E11.6,1X,E11.6,1X,F8.2)
         PRINT 15, (X(K)/1000.0,SIGX,TIME(J),DDOT(K),DOSE(K),K=1,4)
 15      FORMAT (F8.3,4X,F8.3,9X,F4.1,7X,F7.2,5X,F8.2)
 17   CONTINUE
      ENDFILE (3)
      CLOSE (3)
      END
      
      
**************************************************************************            
*                                                                        *
*                    SUBROUTINE FOR MEAN PARTICLE RADII                  *
*                                                                        *
**************************************************************************

      SUBROUTINE MEAN(ALPHA0,BETA0,ALPHA2,ALPHA3,RMEAN)
      REAL ALPHA0,BETA0,ALPHA2,ALPHA3,Z(10),RMEAN(10)
      INTEGER I
      DATA (Z(I),I=1,10)/-1.64,-1.03,-.672,-.384,-.126,.126,.384,.672, 
     &     1.03,1.64/
      ALPHA2=ALPHA0+2.0*BETA0*BETA0
      ALPHA3=ALPHA0+3.0*BETA0*BETA0
      DO 10 I=1,10
         RMEAN(I)=EXP(Z(I)*BETA0+ALPHA3)
 10   CONTINUE
      RETURN
      END