C File: multisub_1b-template.for

C MULTISUB Version 1b May 2005.
C This is a template to construct your own program

	PROGRAM MULTISUB

C Reference: 
C J.Garcia de la Torre, H.E. Perez Sanchez, A. Ortega,J.G. Hernandez, 
C F.G.Diaz and M.C. Lopez Martinez, "Calculation of the solution 
C properties of flexible macromolecules: methods and applications", 
C Eur. Biophys. J. 32, 477-486 (2003)


C This program is intended to be an example of a generator of
C files for program HYDROSUB, with the purpose of calculating
C solution properties of flexible structures composed of ellipsoidal
c and cylindrical subunits, and strings of beads.


C In this example, the SPECIFIC model consists of two ellipsoids 
C joined by a flexible chain, with small sidechains, See the above 
C reference for details

c We distinguish the sections of the program that are GENERAL (should
c be present in  any MULTISUB-like program, and those which are 
c SPECIFIC of the present example


c GENERAL SECTION 1 +++++++++++++++++++++++++++++++++++++++++++++++++
c
c MAIN INPUT FILE (hydrodynamic data, etc) with fixed 
c    name: multisub-dat.txt, and OUTPUT file with fixed
c    name: multisub-res.txt (Caution: program execution will 
c    produce a multisub-res.txt file that would overwrite a file
c    with the same name, if there is one!!!
c    
c CONTENTS OF MAIN INPUT FILE multisub-dat.txt
c   
c TITLE : title to appear in output; up to 30 characters
c FILENAME5 :generic name -without extension - for output files; just 
c    five characters
c NSIG :number of values of sigma. As described in the HYDROSUB 
c    User's Guide, if you give NSIG = -1, then NSIG, SIGMIN and
c    SIGMAX will be taken automatically
c (If NSIG is not -1) SIGMIN :lowest sigma, amstrongs
c (If NSIG is not -1) SIGMAX :highest sigma, amstrongs
c T : temperature, Kelvin
c ETA: solvent viscosity, poise
c RM:  molecular weight (see Hints and Notes in HYDROSUB User Guide)
c VBAR : partial specific volume (cm3/g)
c RHO : solution density (g/cm3)
c IDIF (flag =1 for full report on 6x6 trans+rot diff.)

c NCONF (number of conformations to be generated)

c All these data (except NCONF) are required by HYDROSUB. Please
c   consult the HYDROSUB User Guide for details


c***********************
c DECLARATIONS

	IMPLICIT NONE

	INTEGER NINP,NOUT,NHYD,NFIC
	INTEGER NSIG,IDIF,NCONF,I
	CHARACTER TITLE*20, FILENAME5*5,NAMFIC*12,NUM*3
	CHARACTER TXT*4,AZ*1,BZ*1,CZ*1
	INTEGER I1,I2,I3,J1,J2,J3,J
	REAL SIGMIN,SIGMAX,TEMP,ETA,RM,VBAR,RHO
	PARAMETER(NINP=8,NOUT=9,NHYD=2,NFIC=7)

      INTEGER NSUB,NESF,NSUBMAX,NESFMAX,ITYPE
	REAL X,Y,Z,E,XC,YC,ZC,SEMILONG,SEMISHORT,T,F
	PARAMETER(NSUBMAX=99,NESFMAX=999)
      DIMENSION X(NESFMAX),Y(NESFMAX),Z(NESFMAX),E(NESFMAX)
	DIMENSION XC(NSUBMAX),YC(NSUBMAX),ZC(NSUBMAX),T(NSUBMAX),F(NSUBMAX)
	DIMENSION ITYPE(NSUBMAX),SEMILONG(NSUBMAX),SEMISHORT(NSUBMAX)

	INTEGER NUNCOOR,NUNCOLOR
	PARAMETER(NUNCOOR=21,NUNCOLOR=22)

c END OF GENERAL SECTION 1+++++++++++++++++++++++++++++++++++++++++++

c SPECIFIC SECTION 2 ------------------------------------------------

c ... place here the declarations for variables used
c ... in your specific pieces of code


c END OF SPECIFIC SECTION 2 -----------------------------------------

	
c GENERAL SECTION 3 +++++++++++++++++++++++++++++++++++++++++++++++++
c Open and read main input file

C**********************
c PROGRAM BODY

	OPEN(UNIT=NINP,FILE='multisub-dat.txt',TYPE='OLD') 

	READ(NINP,91) TITLE
91	FORMAT(A20)
	READ(NINP,92) FILENAME5
92	FORMAT(A5)
	READ(NINP,*) NSIG
	IF(NSIG.NE.-1) READ(NINP,*) SIGMIN
	IF(NSIG.NE.-1) READ(NINP,*) SIGMAX
	READ(NINP,*) TEMP
	READ(NINP,*) ETA
	READ(NINP,*) RM
	READ(NINP,*) VBAR
	READ(NINP,*) RHO
	READ(NINP,*) IDIF
	READ(NINP,*) NCONF

	CLOSE(NINP)

c Open and write header and primary data in main (log) output file
	OPEN(UNIT=NOUT,FILE='multisub-log.txt',TYPE='UNKNOWN') 

	WRITE(*,100)
	WRITE(NOUT,100)
100   FORMAT(/
     *5x,'------------------------------------------------------------'/
     *5x,'                  MULTISUB Version 1.b                      '/
     *5x,'J. Garcia de la Torre, H.E. Perez Sanchez, A. Ortega,       '/
     *5x,'J.G. Hernandez, F.G.Diaz and M.C. Lopez Martinez,           '/      
     *5x,'"Calculation of the solution properties of flexible         '/
     *5x,'macromolecules: methods and applications", Eur. Biophys. J. '/
     *5x,'32, 477-486 (2003)                                          '/
     *5x,'                 HYDROSUB Version 5.b-2                     '/
     *5x,'J. Garcia de la Torre and B. Carrasco, "Hydrodynamic        '/
     *5x,'properties of rigid macromolecules composed of ellipsoidal  '/
     *5x,'and cylindrical subunits, Biopolymers 63, 163-167 (2002).   '/                                          '/
     *5x,'------------------------------------------------------------'/
     &)

	WRITE(*,101) TITLE
101	FORMAT('* Title to appear in output: ',A30)
	WRITE(*,102) FILENAME5
102	FORMAT('* Name -without extension- for output files: ',A5)
	WRITE(*,103) NSIG
103   FORMAT('* Number of values of sigma: ',I2)
	IF(NSIG.NE.-1) WRITE(*,104) SIGMIN 
104   FORMAT('* Lowest sigma (Amstrongs): ',f5.2)
	IF(NSIG.NE.-1) WRITE(*,105) SIGMAX 
105   FORMAT('* Highest sigma (Amstrongs): ',f5.2)
	WRITE(*,106) TEMP 
106	FORMAT('* Temperature (Kelvin): ',f6.1)
	WRITE(*,107) ETA 
107   FORMAT('* Solvent viscosity (poise); ',f8.5)
	WRITE(*,108) RM 
108	FORMAT('* Molecular weight: ',f9.0)
      WRITE(*,109) VBAR 
109	FORMAT('* Partial specific volume (cm3/g): ',f6.3)
      WRITE(*,110) VBAR 
110	FORMAT('* Solution density (g/cm3): ',f6.3)
	WRITE(*,111) IDIF
111   FORMAT('* IDIF (=1 for full report on 6x6 trans+rot diff.): ',I1)
	WRITE(*,112) NCONF
112	FORMAT('* Number of conformations to be generated: ',I6)

	WRITE(NOUT,101) TITLE
	WRITE(NOUT,102) FILENAME5
	WRITE(NOUT,103) NSIG
	IF(NSIG.NE.-1) WRITE(NOUT,104) SIGMIN 
	IF(NSIG.NE.-1) WRITE(NOUT,105) SIGMAX 
	WRITE(NOUT,106) TEMP 
	WRITE(NOUT,107) ETA 
	WRITE(NOUT,108) RM 
      WRITE(NOUT,109) VBAR 
      WRITE(NOUT,110) VBAR 
	WRITE(NOUT,111) IDIF
	WRITE(NOUT,112) NCONF

C Open and write general hydrosub.dat file


	OPEN(UNIT=NHYD,FILE='hydrosub.dat',TYPE='UNKNOWN')

	WRITE(*,151) 
	WRITE(NOUT,151) 
151	FORMAT('* Write file hydrosub.dat')
	DO I=1,NCONF
	  TXT='.txt'
	  I1=I/100
	  I2=(I-I1*100)/10
	  I3=I-100*I1-10*I2
	  J1=I1+48
	  J2=I2+48
	  J3=I3+48
	  AZ=CHAR(J1)
	  BZ=CHAR(J2)
	  CZ=CHAR(J3)
	  NAMFIC=FILENAME5//AZ//BZ//CZ//TXT
	  NUM=AZ//BZ//CZ

	  WRITE(NHYD,50)FILENAME5,NUM
	  WRITE(NHYD,50)FILENAME5,NUM
50	  FORMAT(A5,A3)
	  WRITE(NHYD,60)NAMFIC
60	  FORMAT(A12)
	  WRITE(NHYD,300)NSIG
	  WRITE(NHYD,310)SIGMIN
	  WRITE(NHYD,310)SIGMAX	
	  WRITE(NHYD,320)TEMP
        WRITE(NHYD,330)ETA
	  WRITE(NHYD,340)RM
	  WRITE(NHYD,350)VBAR
	  WRITE(NHYD,330)RHO
	  WRITE(NHYD,300)IDIF
300	  FORMAT(I3,',')
310	  FORMAT(F5.2,',')
320	  FORMAT(F6.2,',')
330	  FORMAT(F4.2,',')
340	  FORMAT(F9.2,',')
350	  FORMAT(F6.3,',')

	ENDDO

	WRITE(NHYD,200)
200	FORMAT('*')

	CLOSE(NHYD)


c END OF GENERAL SECTION 3+++++++++++++++++++++++++++++++++++++++

c SPECIFIC SECTION 4 --------------------------------------------

c ... Here you will provide the main data of your model, which
c ... will be needed for the generation of each conformation

C END OF SPECIFIC SECTION 4 ------------------------------------

C GENERAL SECTION 5 ++++++++++++++++++++++++++++++++++++++++++++

c DO loop over the number, NCONF of conformations to be generated
	
	WRITE(NOUT,131) NCONF
	WRITE(*,131) NCONF
131	FORMAT('* Writing structural files for',i4,' conformations')

      DO I=1,NCONF

	TXT='.txt'

	  I1=I/100
	  I2=(I-I1*100)/10
	  I3=I-100*I1-10*I2
	  J1=I1+48
	  J2=I2+48
	  J3=I3+48
	  AZ=CHAR(J1)
	  BZ=CHAR(J2)
	  CZ=CHAR(J3)
	  NAMFIC=FILENAME5//AZ//BZ//CZ//TXT

	  WRITE(NOUT,132) I,NAMFIC
  	  WRITE(*,132) I,NAMFIC
132	  FORMAT('   * Conformation:',i4,'   File:',A12)

	  OPEN(UNIT=NFIC,FILE=NAMFIC,TYPE='UNKNOWN')

c END OF GENERAL SECTION 5 +++++++++++++++++++++++++++++++++++++++

c SPECIFIC SECTION 6  --------------------------------------------

c ... Here each specific conformation is constructed
c ... Values of ITYPE, SEMILONG, SEMISHORT, XC, YC, ZC, T, F,
C ... NESF, X, Y, Z , E will be asigned for each shell subunit and 
c ... and for each sphere


C END OF SPECIFIC SECTION 6  -----------------------------------------

C GENERAL SECTION 7  +++++++++++++++++++++++++++++++++++++++++++++++++

C... Write data for this conformation to structural file

	  WRITE(NFIC,230) NSUB
230	  FORMAT(I4)
	  DO J=1,NSUB
	    WRITE(NFIC,240)ITYPE(J),XC(J),YC(J),ZC(J),T(J),F(J),
     *      SEMILONG(J),SEMISHORT(J)
	  ENDDO
240	  FORMAT(I4,7F7.1)

	  WRITE(NFIC,210)NESF
210	  FORMAT(I3)
	  DO J=1,NESF
		WRITE(NFIC,220)X(J),Y(J),Z(J),E(J)
220	  FORMAT(4(F9.3,2X))
	  ENDDO

	  CLOSE(NFIC)
	ENDDO

	WRITE(*,139)
	WRITE(NOUT,139)
139	FORMAT('* End of calculations. All files generated')

	END

C END OF SPECIFIC SECTION 7  -----------------------------------------

C GENERAL SECTION 8  +++++++++++++++++++++++++++++++++++++++++++++++++

C****************************************************
C Subroutine GENERCHAIN

C This subroutine generates conformations of a freely jointed chain
c with small side groups

c The chain consists of:
c a) NSKEL skeletal beads of diameter b, joined by NSKEL-1
c skeletal bonds of length b. The bond angles are random except for
c the excluded volume restriction, which makes them greater than 
c 60 degrees
c b) NSKEL-2 beads, representing side groups, joined to all the 
c skeletal beads except the terminal ones. The side bonds are directed 
c along the external bisector of the bond angles.
c
C The total number of beads is NT=2*NSKEL-2
c
c Excluded-volume conditions are introduced so that the distance between
c any two beads must be larger than the diameter b
c  
c INPUT data
c-----------
C ISEED Integer random seed
c NSKEL Number of skeletal beads
c BESF bead diameter and bond length
c NESFMAX Number to which X,Y and Z are dimensioned in the calling program
c OUTPUT results
c---------------
c X, Y, Z : Cartesian coordinates of the NT beads
c---------------------------------------------------------------

	SUBROUTINE GENERCHAIN(ISEED,NSKEL,NESFMAX,BESF,X,Y,Z)

	IMPLICIT NONE

	INTEGER ISEED,NSKEL,NESFMAX,NTOT
      REAL X,Y,Z,BESF,PI,RAN2
      DIMENSION X(NESFMAX),Y(NESFMAX),Z(NESFMAX)
      REAL CT,ST,FI
	INTEGER J,JC
	REAL BJM1X,BJX,BSUMX,BJM1Y,BJY,BSUMY,BJM1Z,BJZ,BSUMZ,RMOD
      REAL DX,DY,DZ,DIS2

	PARAMETER(PI=3.14159)

30	CONTINUE

C First two beads

	X(1)=0.
	Y(1)=0.
	Z(1)=0.
	X(2)=0.
	Y(2)=0.
	Z(2)=BESF

	j=1
	j=2

	DO J=3,NSKEL

c Bead J
	  
20	  CONTINUE

	  CT=-1.+2*RAN2(ISEED)
	  ST=SQRT(1-CT**2)
	  FI=2*PI*RAN2(ISEED)
	  X(J)=X(J-1)+BESF*ST*COS(FI)
        Y(J)=Y(J-1)+BESF*ST*SIN(FI)
	  Z(J)=Z(J-1)+BESF*CT
					
c Check excluded volume for bead J

	  DO 25 JC=1,J-2
		DX=X(JC)-X(J)
		IF(ABS(DX).GT.BESF) GOTO 25
		DY=Y(JC)-Y(J)
		IF(ABS(DY).GT.BESF) GOTO 25
		DZ=Z(JC)-Z(J)
		IF(ABS(DZ).GT.BESF) GOTO 25
		DIS2=DX**2+DY**2+DZ**2
		IF(DIS2.LT.BESF*BESF) GOTO 20
25	  CONTINUE


	ENDDO

	NTOT=NSKEL
	DO J=2,NSKEL-1

C Side bead for skeletal bead J

	  BJM1X=X(J)-X(J-1)
	  BJM1Y=Y(J)-Y(J-1)
	  BJM1Z=Z(J)-Z(J-1)
	  BJX=X(J+1)-X(J)
	  BJY=Y(J+1)-Y(J)
	  BJZ=Z(J+1)-Z(J)

	  BSUMX=BJM1X-BJX
	  BSUMY=BJM1Y-BJY
	  BSUMZ=BJM1Z-BJZ
	  RMOD=SQRT(BSUMX**2+BSUMY**2+BSUMZ**2)
	  BSUMX=BSUMX*BESF/RMOD
	  BSUMY=BSUMY*BESF/RMOD
	  BSUMZ=BSUMZ*BESF/RMOD

	  NTOT=NTOT+1
	  X(NTOT)=X(J)+BSUMX
	  Y(NTOT)=Y(J)+BSUMY
	  Z(NTOT)=Z(J)+BSUMZ
	
C Check excluded volume for side bead

	  DO  35 JC=1,NTOT-1
	  IF(JC.NE.J) THEN
		DX=X(JC)-X(NTOT)
		IF(ABS(DX).GT.BESF) GOTO 35
		DY=Y(JC)-Y(NTOT)
		IF(ABS(DY).GT.BESF) GOTO 35
		DZ=Z(JC)-Z(NTOT)
		IF(ABS(DZ).GT.BESF) GOTO 35
		DIS2=DX**2+DY**2+DZ**2
		IF(DIS2.LT.BESF*BESF) GOTO 30
        ENDIF
35	  CONTINUE
	
	ENDDO

	RETURN
	END

c**************

c SUBROUTINE CHECK_ESF1_ESF2

C Check overlapping between two spheres
c Sphere 1 Center at (XC1, YC1, ZC1), radius R1
c Sphere 2 Center at (XC2, YC2, ZC2), radius R2
c Condition: Center-to-center distance must be 
c greater than sum of radii
c IOVERLAP=1 in case of bead overlappin


	SUBROUTINE CHECK_ESF1_ESF2(R1,XC1,YC1,ZC1,
     * R2,XC2,YC2,ZC2,IOVERLAP)

      IOVERLAP=0
	DIS2=(XC1-XC2)**2+(YC1-YC2)**2+(ZC1-ZC2)**2
	DMAX=R1+R2
	IF(DIS2.LT.DMAX**2) IOVERLAP=1

	RETURN
	END

C******************

c SUBROUTINE CHECK_ESF_ELI2

C Check overlapping between:
c    Spherical bead with coords. XX,YY,ZZ
c    Ellipsoid centered at XC2,YC2,ZC2, with symmetry and perpendicular
c       semiaxes A2 and B2, and polar angles of symmetry axis FI2 (phi)
c       and T2 (theta)
c RESULT: ICELIP2 =0 (no EV) or =1 (EV detected)
c
c NOTE: Bead size is ignored. What is tested is wheter the CENTER of
c the beads is within the ellipsoid

	SUBROUTINE CHECK_ESF_ELI2(XX,YY,ZZ,A2,B2,
     + XC2,YC2,ZC2,FI2,T2,ICELIP2)

	IMPLICIT NONE
      REAL A2,B2,XC2,YC2,ZC2,XCR,YCR,ZCR,FI2,T2,XX,YY,ZZ
	INTEGER ICELIP2
	REAL XC,YC,ZC,AAA,BBB,CCC,DENTRO

C Translate center of coords. to center of ellipsoid
	XC=XX-XC2
	YC=YY-YC2
	ZC=ZZ-ZC2

C Rotate axis of coords. so that symmetry axis of ellipsoid is along z
c The rotation angles are those of z-axis of ellipsoid 2 with a change
c of sign
			
      CALL ROTA(XC,YC,ZC,XCR,YCR,ZCR,-T2,-FI2)

C Check for overlapping
	AAA=ABS(XCR/A2)
	IF(AAA.GE.1)GOTO 456
	BBB=ABS(YCR/A2)
	IF(BBB.GE.1)GOTO 456
	CCC=ABS(ZCR/B2)
	IF(CCC.GE.1)GOTO 456
	DENTRO=AAA**2+BBB**2+CCC**2
	IF(DENTRO.LT.1)THEN
	  ICELIP2=1	
        RETURN
	ENDIF
456	CONTINUE

	ICELIP2=0

	RETURN
	END



C****************************


c SUBROUTINE CHECK_ELI1_ELI2

C Check overlapping between two ellipsoids such that:
c 
c    Ellipsoid 2, with symmetry and perpendicular semiaxes B1 and A1, is 
c    centered at an arbitrary point, (XC2,YC2,ZC2), and its symmetry
c    axis makes polar FI2 (phi) and T2 (theta) - in radians - with 
c    the Cartesian axes of reference

c    Ellipsoid 1,  centered at XC1=0, YC1=0, and ZC1 variable, with 
c    symmetry and perpendicular semiaxes B1 and A1, and its symmetry
c    axis along Z

c    Check whether some point on the surface of ellipsoid 2 is within
c    the interior of ellipsoid 1

c RESULT: ICELIPS =0 (no overlap) or =1 (overlap detected)

	SUBROUTINE CHECK_ELI1_ELI2(A1,B1,A2,B2,XC1,YC1,ZC1,
     * XC2,YC2,ZC2,FI2,T2,ICELIPS)

	IMPLICIT NONE
	REAL A1,B1,A2,B2,XC1,YC1,ZC1,XC2,YC2,ZC2,FI2,T2,PI
      INTEGER ICELIPS
	PARAMETER(PI=3.14159)
	INTEGER IT2,IFI2
	REAL TE2,F2
	REAL AAA,BBB,CCC,DENTRO
	REAL XXX2,YYY2,ZZZ2,XX2,YY2,ZZ2,X2,Y2,Z2

c... Generate a number of points on the surface of ellipsoid 2

	DO IT2=0,50
	  TE2=(PI/50.)*IT2
	  DO IFI2=0,100
	    F2=((2*PI)/100.)*IFI2

c 1st, point P, if ellipsoid 2 were centered at origin and aligned
		XXX2=A2*SIN(TE2)*COS(F2)
		YYY2=A2*SIN(TE2)*SIN(F2)
		ZZZ2=B2*COS(TE2)
c 2nd, point P, if ellipsoid 2 is centered at origin and rotated
		CALL ROTA(XXX2,YYY2,ZZZ2,XX2,YY2,ZZ2,T2,FI2)
c 3rd, point P, if ellipsoid 2 is translated to its center
	    X2=XX2+XC2
	    Y2=YY2+YC2
	    Z2=ZZ2+ZC2

c Now check whether point is within ellipsoid 1

		AAA=ABS(X2/A1)
		IF(AAA.GE.1)GOTO 457
		BBB=ABS(Y2/A1)
		IF(BBB.GE.1)GOTO 457
		CCC=ABS((Z2-ZC1)/B1)
		IF(CCC.GE.1)GOTO 457
		DENTRO=AAA**2+BBB**2+CCC**2
		IF(DENTRO.LT.1)THEN
		  ICELIPS=1
		  RETURN
		ENDIF
	
457 			CONTINUE 
			ENDDO
		ENDDO


	ICELIPS=0
	RETURN

	END


c************

	SUBROUTINE  ROTA(X,Y,Z,XX,YY,ZZ,DTHETA,DFI)

	PI=ACOS(-1.)
	DR=0.

C	PASAMOS DE CARTESIANAS A POLARES

	R=SQRT(X**2+Y**2+Z**2)
	IF(R.EQ.0)RETURN

	THETA=ACOS(Z/R)

	IF(X.GE.0.AND.Y.GE.0)GOTO 123

	AAA=ABS(Y)/ABS(X)

	FI=ATAN(AAA)

	IF(X.LE.0.AND.Y.GE.0)THEN 
		FI=PI-FI
	ELSE
		IF(X.LE.0.AND.Y.LE.0)THEN 
			FI=FI+PI
			ELSE
				IF(X.GE.0.AND.Y.LE.0)THEN 
					FI=2*PI-FI
				ENDIF
		ENDIF
	ENDIF

123	CONTINUE

	R=R+DR
	THETA=THETA+DTHETA
	FI=FI+DFI

	XX=R*SIN(THETA)*COS(FI)
	YY=R*SIN(THETA)*SIN(FI)
	ZZ=R*COS(THETA)

	RETURN
	END

c********************************************************************
c This is a machine-independent random number generator
c If you prefer to user your compiler's function, then 
c look for RAN(ISEED) in the code above and replace it for
c your system's call, usually RAN(ISEED)

      REAL FUNCTION ran2(idum)

      INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
      REAL AM,EPS,RNMX
      PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,
     *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,
     *NTAB=322,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
      INTEGER idum2,j,k,iv(NTAB),iy
      SAVE iv,iy,idum2
      DATA idum2/123456789/, iv/NTAB*0/, iy/0/
      if (idum.le.0) then
        idum=max(-idum,1)
        idum2=idum
        do 11 j=NTAB+8,1,-1
          k=idum/IQ1
          idum=IA1*(idum-k*IQ1)-k*IR1
          if (idum.lt.0) idum=idum+IM1
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ1
      idum=IA1*(idum-k*IQ1)-k*IR1
      if (idum.lt.0) idum=idum+IM1
      k=idum2/IQ2
      idum2=IA2*(idum2-k*IQ2)-k*IR2
      if (idum2.lt.0) idum2=idum2+IM2
      j=1+iy/NDIV
      iy=iv(j)-idum2
      iv(j)=idum
      if(iy.lt.1)iy=iy+IMM1
      ran2=min(AM*iy,RNMX)
      return
      END

C END OF SPECIFIC SECTION 7  -----------------------------------------

C END OF PROGRAM *****************************************************
