c----------------------------------------------------------getusrparams SUBROUTINE getusrparams c---------------------------------------------------------------------c c purpose: c c To read user parameters. c c---------------------------------------------------------------------c call gettitle call getgridparams call getgeneralparams call getanalysisoptn call getoutputflags call getconservparams c BMM 6/8/96 call getgridoutputs return END c-----------------------------------------------------echoallparameters SUBROUTINE echoallparameters c---------------------------------------------------------------------c c purpose: c c To display all user parameters. c c---------------------------------------------------------------------c INCLUDE '3duparms.cb' INCLUDE '3dgeom.cb' WRITE(*,*) WRITE(*,*) '--------------------------------------------' WRITE(*,*) 'User parameters are:-' WRITE(*,*) WRITE(*,'(2A15)') 'goE','goN' WRITE(*,20) goE, goN WRITE(*,*) WRITE(*,'(2A15)') 'gaS','gaC' WRITE(*,25) gaS, gaC WRITE(*,*) WRITE(*,'(3A15)') 'numX', 'numY', 'numZ' WRITE(*,30) numX, numY, numZ WRITE(*,*) WRITE(*,'(2A15)') 'dimXY', 'dimZ' WRITE(*,40) dimXY, dimZ WRITE(*,*) WRITE(*,'(2A15)') 'resXY', 'resZ' WRITE(*,42) resXY, resZ WRITE(*,*) WRITE(*,'(5A15)') 'newcondTS','outputTS','dt1','dt2', 1 'maxoutputs' WRITE(*,50) newcondTS, outputTS, dt,dtnew, maxoutputs WRITE(*,*) WRITE(*,'(3A15)') 'tot_dead_age', 'mass_pp', 'c_init' WRITE(*,60) tot_dead_age, mass_pp, c_init WRITE(*,*) WRITE(*,'(A)') 'Diffs: Element types not listed use the first' WRITE(*,'(3A15)') 'Element type', 'HDIFF','VDIFF' i=1 DO WHILE (EDIFF(i).ne.-1) WRITE(*,'(I15,2F15.7)') EDIFF(i),HDIFF(i),VDIFF(i) i=i+1 ENDDO WRITE(*,*) IF (dieoff) THEN WRITE(*,'(A)') 'dieoff = TRUE' ELSE WRITE(*,'(A)') 'dieoff = FALSE' ENDIF IF (dieoffc) THEN WRITE(*,'(A)') 'dieoffc = TRUE' ELSE WRITE(*,'(A)') 'dieoffc = FALSE' ENDIF WRITE(*,'(3A15)') 'T90_S','T90_D','TDPTH' WRITE(*,'(3F15.3)') t90_s,t90_d,tdpth WRITE(*,*) WRITE(*,*) '--------------------------------------------' WRITE(*,*) c c-------check whether the time steps meet the required relations c if(mod(newcondTS,dt).ne.0.or.mod(outputTS,dt).ne.0) goto 70 if(mod(newcondTS,outputTS).ne.0.and. * mod(outputTS,newcondTS).ne.0) goto 70 if(mod(dt,dtnew).ne.0) goto 70 20 FORMAT(2F15.3) 25 FORMAT(2F15.9) 30 FORMAT(3I15) 40 FORMAT(2F15.3) 42 FORMAT(2I15) 50 FORMAT(5I15) 60 FORMAT(I15,2E15.3) RETURN 70 write(*,*) '** Time Steps of newcondTS, outputTS and dt'// * ' not as required **' stop END c-----------------------------------------------------getpollconditions SUBROUTINE getpollconditions(poll,massp,masspA,NUMPARTS, * max_numps,key) c---------------------------------------------------------------------c c purpose: c c To input a new set of pollution conditions. c c---------------------------------------------------------------------c INCLUDE '3duparms.cb' INCLUDE '3dpolls.cb' REAL poll(4,1) c poll(1,i) = Vol of discharge at window i c poll(2,i) = Conc of discharge at window i c poll(3,i) = Depth of discharge at window i c poll(4,i) = Radius of discharge at window i INTEGER i,d,w,max_numps,NUMPARTS,key REAL massp(1),masspA(1),masspB,numd c *** Read in the new pollutant releases. key=0 max_numps=NUMPARTS DO 60 w=1, num_pollwins READ(4,*,ERR=200,end=205) d,(poll(i,d),i=1,4) IF (d.gt.num_pollwins.or.d.lt.1) THEN PRINT *,'ERROR - Pollutant window' STOP ENDIF masspB=mass_pp numd=poll(1,d)*poll(2,d)*min(outputTS,newcondTS) * /mass_pp if(numd.gt.float(lnum_ps)) then masspB=poll(1,d)*poll(2,d)* * min(outputTS,newcondTS)/lnum_ps endif massp(d)=max(masspB,masspA(d)) numd=poll(1,d)*poll(2,d)*min(outputTS,newcondTS) * /massp(d)+25 max_numps=max_numps+int(numd) masspA(d)=masspB 60 CONTINUE key=1 65 FORMAT(I2, F10.5, E10.1, 2F10.3) RETURN 200 WRITE(*,*) 'Error in pollutants!' STOP 205 WRITE(*,*) 'Run out of pollutant data' return END c-----------------------------------------------------getflowconditions SUBROUTINE getflowconditions(nsurf,vel,cord,ao,wsel,wsel1,ic, & dhdt,nop,ielist,key) c---------------------------------------------------------------------c c purpose: c c To input velocity field. c c This reads one time step as called from walkem c c If it is a 2D file then it makes an artificial 3d one c c---------------------------------------------------------------------c INCLUDE '3duparms.cb' INCLUDE '3dmesh.cb' INCLUDE '3dgeom.cb' INCLUDE '3dpolls.cb' common /shape/shap(20),shpx(20),shpy(20),shpz(20) REAL vel(3,1) c vel(1,i) = Velocity to E at node i c vel(2,i) = Velocity to N at node i c vel(3,i) = Velocity to Z at node i INTEGER*2 nsurf(1),nop(20,1),ic(1),ielist(1) INTEGER i,k,j,no3dn,ndg,no3de,key REAL vel1,vel2,depth,vel4,vel5,vel6,tet,dfct,wsel(1),wsel1(1) real dhdt(1),prsm(3,6),brck(3,8),xi(8),eta(8),zeta(8) real cord(3,1),a(2,2) real*8 ao(1),detj data prsm/0.,-1.,-1.,1.,0.,-1.,0.,1.,-1., & 0.,-1.,1.,1.,0.,1.,0.,1.,1./ data brck/-1.,-1.,-1.,1.,-1.,-1.,1.,1.,-1.,-1.,1.,-1., & -1.,-1.,1.,1.,-1.,1.,1.,1.,1.,-1.,1.,1./ c *** Read in the velocity field key=0 if(dimopt(1:2).eq.'10') then cBM2019 updating reads to OUTBNRMA format if(compact) then c read(2,err=210,end=215) TET,no3dn,NDG,no3de, c 1 ((VEL(K,J),K=1,2),wsel(j),VEL(3,j),J=1,No3dn) WRITE(*,*) 'No compact!' STOP else c read(2,err=210,end=215) TET,no3dn,NDG,no3de,IYRR, c 1 ((VEL(K,J),K=1,2),wsel(j),vel4,vel5,vel6, c 2 VEL(3,j),J=1,No3dn),(DFCT,J=1,No3dE) read(2,err=210,end=215) TET,no3dn,NDG,no3de,IYRR, 1 ((VEL(K,J),K=1,2),depth,vel4,vel5,vel6, 2 VEL(3,j),wsel(j),J=1,no3dn), 3 (DFCT,J=1,no3de), 4 (vel7,j=1,no3dn) c WRITE(IRMAFM) TETT,NP,NDG,NE,IYRR, c 1 ((VSING(K,J),K=1,6),VVEL(J),WSLL(J),J=1,NP),(DFCT(J),J=1,NE),(VSING(7,J),J=1,NP) C Output RMA results file contains c 1 time in hours (Julian) c 2 number of nodes c 3 obsolete counter of degrees of freedom (set to 5) NDF=6 c 3a number of elements c 4 year c 5 VSING subscript(1) x-vel by node c 6 VSING subscript(2) y-vel by node c 7 VSING subscript(3) depth by node c 8 VSING subscript(4) salinity by node c 9 VSING subscript(5) temperature by node c 10 VSING subscript(6) sus-sed by node c 11 VVEL w-vel by node c 12 WSLL water surface elevation c 13 DFCT stratification multiplier by element c 14 VSING subscript(7) time derivative of depth by node print*,'Velocity file read at time',TET endif else c read(2,err=210,end=215) TET,no3dn,((vel(k,j),k=1,2), c 1 wsel(j),j=1,no3dn),(vel4,j=1,no3dn) WRITE(*,*) 'No RMA2!' STOP endif c c-------For 2D case - START c-------artificially create a 3D velocity field c-------The 2D field could be from RMA2 or a depth averaged RMA10 if(no3dn.ne.np) then do i=1,np ic(i)=0 enddo do ii=1,n1_2d i=ielist(ii) if(nop(20,i).ne.0) then nen=20 else nen=15 endif do k=1,nen vel(3,nop(k,i))=0.0 enddo enddo do i=no3dn+1,np wsel(i)=wsel(nsurf(i)) enddo if(restart) then do i=1,np dhdt(i)=(wsel(i)-wsel1(i))/float(newcondTS) wsel0=wsel(i) wsel(i)=wsel1(i) wsel1(i)=wsel0 enddo endif do i=no3dn+1,np if(nsurf(i).eq.nsurf(i-1)) then cord(3,i)=-wsel(nsurf(i))/2.0 else cord(3,i)=-wsel(nsurf(i)) endif do k=1,3 vel(k,i)=vel(k,nsurf(i)) enddo enddo c c----------calculate vertical velocity at corner nodes c do ii=1,n1_2d i=ielist(ii) if(nop(6,i).eq.0) goto 10 if(nop(20,i).ne.0) then nen=20 iad=12 it=3 do k=1,8 xi(k)=brck(1,k) eta(k)=brck(2,k) zeta(k)=brck(3,k) enddo else nen=15 iad=9 it=2 do k=1,6 xi(k)=prsm(1,k) eta(k)=prsm(2,k) zeta(k)=prsm(3,k) enddo endif k=0 do n=1,2*nen/5,2 k=k+1 call xn3x(it,nen,xi(k),eta(k),zeta(k)) c c----------------calculate the Jacobian matrix c do m=1,2 a(1,m)=0.0 a(2,m)=0.0 dadx=0.0 dady=0.0 dhdx=0.0 dhdy=0.0 do j=1,nen nod=nop(j,i) a(1,m)=a(1,m)+shpx(j)*cord(m,nod) a(2,m)=a(2,m)+shpy(j)*cord(m,nod) dadx=dadx+shpx(j)*ao(nod) dady=dady+shpy(j)*ao(nod) dhdx=dhdx+shpx(j)*(wsel(nod)+ao(nod)) dhdy=dhdy+shpy(j)*(wsel(nod)+ao(nod)) enddo enddo detj=a(1,1)*a(2,2)-a(1,2)*a(2,1) vel(3,nop(n,i))=vel(3,nop(n,i))+ & vel(1,nop(n,i))*(dadx*a(2,2)-dady*a(1,2))/detj+ & vel(2,nop(n,i))*(dady*a(1,1)-dadx*a(2,1))/detj vel(3,nop(n+iad,i))=vel(3,nop(n+iad,i))+ & vel(1,nop(n+iad,i))*(dhdx*a(2,2)-dhdy*a(1,2))/detj+ & vel(2,nop(n+iad,i))*(dhdy*a(1,1)-dhdx*a(2,1))/detj+ & dhdt(nop(n+iad,i)) ic(nop(n,i))=ic(nop(n,i))+1 ic(nop(n+iad,i))=ic(nop(n+iad,i))+1 enddo 10 enddo do i=1,np if(ic(i).gt.1) then vel(3,i)=vel(3,i)/float(ic(i)) endif enddo c c----------interpolate the mid-side nodal values c do i=1,ne if(nop(6,i).eq.0) goto 20 if(nop(20,i).ne.0) then nen=8 iad=12 else nen=6 iad=9 endif do n=2,nen,2 vel(3,nop(n,i))=(vel(3,nop(n-1,i))+ & vel(3,nop(mod(n+1,nen),i)))/2.0 vel(3,nop(n+iad,i))=(vel(3,nop(n+iad-1,i))+ & vel(3,nop(mod(n+1,nen)+iad,i)))/2.0 enddo do n=1,nen/2 vel(3,nop(nen+n,i))=(vel(3,nop(2*n-1,i))+ & vel(3,nop(2*n-1+iad,i)))/2.0 enddo 20 enddo c c----------calculate dhdt(i) for next time step c cBM2019 WRITE(*,*) 'Shouldnt be here' STOP c if(dimopt(1:2).eq.'10') then if(compact) then read(2,err=30,end=30) TET,no3dn,NDG,no3de, 1 (vel1,vel2,wsel1(j), 2 VEL1,J=1,No3dn) else cDRC 04/06/99 IYRR added for RMA10v6.5 read(2,err=30,end=30) TET,no3dn,NDG,no3de,IYRR, 1 (vel1,vel2,wsel1(j),vel4,vel5,vel6, 2 VEL1,J=1,No3dn),(DFCT,J=1,No3dE) endif else read(2,err=30,end=30) TET,no3dn,(vel1,vel2,wsel1(j), 1 j=1,no3dn),(vel1,j=1,no3dn) endif backspace 2 do i=no3dn+1,np wsel1(i)=wsel1(nsurf(i)) enddo do i=1,np dhdt(i)=(wsel1(i)-wsel(i))/float(newcondTS) enddo 30 continue endif c-------For 2D case - END 40 key=1 RETURN 210 WRITE(*,*) ' *** Error in velocities! *** ' STOP 215 WRITE(*,*) ' *** Run out of velocity data *** ' return END c-----------------------------------------------------------getgeometry subroutine getgeometry(cord,nop,nsurf,imat,width,ao, * ielist,diffs) c---------------------------------------------------------------------c c purpose: c c To read geometric data of FE mesh. c c---------------------------------------------------------------------c include '3duparms.cb' include '3dmesh.cb' integer j,k integer nfix,nfix1,nref integer*2 imat(1),ncorn,iem,ndep,netyp integer*4 nfixh !BM2019 changed from integer*2 real width(1),paxis,alfa,th cBM2019 c the 3dg file now has cord as real*8 c and nop as real*4 integer*2 nop(20,1),nsurf(1) integer nop4(30000,20),nsurf4(30000) integer*2 ielist(1) ! this one OK as it is created internally real*8 cord8(3,30000) ! BM2019 hack because of real*8 real cord(3,1) real spec(3) real*8 ao(1) real diffs(2,1) c if(dimopt(1:2).eq.'10') then read(9,err=10) NP,NE,NPM,NES,((CORD8(k,j),SPEC(k),K=1,3), 1 ALFA,NFIX,NFIX1,AO(J),NSURF4(J),J=1,NP), 2 (NDEP,NREF,J=1,NPM), 3 ((NOP4(k,j),K=1,20),NCORN,IMAT(J),NETYP,TH, 4 NFIXH,J=1,NE),(WIDTH(J),J=1,NP) do j=1,np do k=1,3 cord(k,j)=cord8(k,j) enddo c write(58,'(I5,3F20.3)') j,(cord(k,j),k=1,3) enddo do j=1,ne do k=1,20 nop(k,j)=nop4(k,j) enddo c write(59,'(i5,20I8)') j,(nop(k,j),k=1,20) enddo do j=1,np nsurf(j)=nsurf4(j) enddo else cBM2019 write(*,*) 'Stop: Geometry: No 2D' stop read(9,err=10) NP,NE,((CORD(K,J),K=1,2),ALFA,cord(3,j), 1 J=1,NP),((NOP(K,J),K=1,8),IMAT(j),PAXIS,IEM,J=1,NE) do i=1,np ao(i)=cord(3,i) cord(1,i)=cord(1,i)*xscale cord(2,i)=cord(2,i)*yscale enddo endif print *,'np,npm=',np,npm print *,'ne,nes=',ne,nes do i=1,NE if(nop(20,i).ne.0) then nen=20 goto 11 else if(nop(15,i).ne.0.and.nop(16,i).eq.0) then nen=15 goto 11 else if(nop(13,i).ne.0.and.nop(14,i).eq.0) then nen=13 goto 11 else if(nop(10,i).ne.0.and.nop(11,i).eq.0) then nen=10 goto 11 else if(nop(8,i).ne.0.and.nop(9,i).eq.0) then nen=8 goto 11 else if(nop(6,i).ne.0.and.nop(7,i).eq.0) then nen=6 goto 11 else if(nop(3,i).ne.0.and.nop(4,i).eq.0) then nen=3 goto 11 endif c print *,i,imat(i),nen c print *,i,(nop(k,i),k=1,20) c 11 print *,i,imat(i),nen c print *,i,(nop(k,i),k=1,nen) 11 enddo n1_2d=0 if(np.eq.npm.or.dimopt(1:1).eq.'2') then c c----------generate nodes below surface c cbm2019 write(*,*) 'Stop: no 2d!' stop np3=np do n=1,np cord(3,n)=0.0 do i=1,ne if(nop(6,i).eq.0) goto 15 !1D or 1D to 2D junction elements if(nop(7,i).eq.0) then nen=6 else nen=8 endif do k=1,nen-1,2 if(n.eq.nop(k,i)) goto 20 !corner nodes enddo do k=2,nen,2 if(n.eq.nop(k,i)) goto 30 !mid-side nodes enddo 15 enddo 20 nn=2 nsurf(n)=n goto 40 30 nn=1 nsurf(n)=-n 40 do k=1,nn np3=np3+1 cord(1,np3)=cord(1,n) cord(2,np3)=cord(2,n) cord(3,np3)=ao(n)/float(k) ao(np3)=ao(n) nsurf(np3)=n enddo enddo c c----------form 3D element c do i=1,ne if(nop(6,i).eq.0) goto 60 if(nop(7,i).eq.0) then nen=6 iad=9 else nen=8 iad=12 endif do k=1,nen nop(k+iad,i)=nop(k,i) enddo do k=1,nen do n=np+1,np3 if(nsurf(n).eq.nop(k+iad,i)) goto 50 enddo write(*,*) ' ** Error in code **' stop 50 nop(k,i)=n if(mod(k,2).eq.1) then ik=nen+(k+1)/2 nop(ik,i)=n+1 endif enddo n1_2d=n1_2d+1 ielist(n1_2d)=i c print *,'nop=',i,(nop(k,i),k=1,nen+iad) 60 enddo npm=np np=np3 else c c----------check whether there are any 1D & 2D elements in the 3D mesh c np3=np do i=1,ne if(imat(i).lt.900) then if(nop(10,i).ne.0) goto 65 nen=0 if(nop(8,i).ne.0.and.nop(9,i).eq.0) then nen=8 iad=12 else if(nop(6,i).ne.0.and.nop(7,i).eq.0) then nen=6 iad=9 else if(nop(3,i).ne.0.and.nop(4,i).eq.0) then nen=3 c Modify 22/1/96 ver 2.2 c Works for meshes that don't have all element numbers. else if(nop(1,i).eq.0.and.nop(2,i).eq.0) then nen=0 else nen=5 endif if(nen.eq.0) goto 65 n1_2d=n1_2d+1 ielist(n1_2d)=i if(nen.gt.5) then !2D elements do k=1,nen nop(k+iad,i)=nop(k,i) nsurf(nop(k+iad,i))=nop(k,i) cord(3,nop(k+iad,i))=cord(3,1) ao(nop(k+iad,i))=ao(nop(k,i)) enddo else if(nen.le.5) then !1D or 1D to 2D junction elements print *,'TT',i,nen,nop(1,i),nop(2,i) iad=12 aa=sqrt((cord(2,nop(3,i))-cord(2,nop(1,i)))* & (cord(2,nop(3,i))-cord(2,nop(1,i)))+ & (cord(1,nop(3,i))-cord(1,nop(1,i)))* & (cord(1,nop(3,i))-cord(1,nop(1,i)))) asin=-(cord(2,nop(3,i))-cord(2,nop(1,i)))/aa acos=(cord(1,nop(3,i))-cord(1,nop(1,i)))/aa wd=width(nop(1,i)) x=cord(1,nop(1,i))+asin*wd/2.0 y=cord(2,nop(1,i))+acos*wd/2.0 z=cord(3,1) ao0=ao(nop(1,i)) call existnod(cord,ao,nsurf,nop,np,np3,x,y,z,ao0, & 1,0,i,l) nop(1+iad,i)=l x=cord(1,nop(1,i))-asin*wd/2.0 y=cord(2,nop(1,i))-acos*wd/2.0 z=cord(3,1) ao0=ao(nop(1,i)) call existnod(cord,ao,nsurf,nop,np,np3,x,y,z,ao0, & 1,0,i,l) nop(3+iad,i)=l nop(2+iad,i)=nop(1,i) nsurf(nop(2+iad,i))=nop(1,i) if(nen.eq.5) then !1D to 2D junction elements nop(6+iad,i)=nop(3,i) nsurf(nop(6+iad,i))=nop(3,i) x13=cord(1,nop(4,i))-cord(1,nop(1,i)) y13=cord(2,nop(4,i))-cord(2,nop(1,i)) x23=cord(1,nop(5,i))-cord(1,nop(1,i)) y23=cord(2,nop(5,i))-cord(2,nop(1,i)) area=x13*y23-x23*y13 if(area.gt.0) then nop(5+iad,i)=nop(4,i) nsurf(nop(5+iad,i))=nop(4,i) nop(7+iad,i)=nop(5,i) nsurf(nop(7+iad,i))=nop(5,i) else nop(5+iad,i)=nop(5,i) nsurf(nop(5+iad,i))=nop(5,i) nop(7+iad,i)=nop(4,i) nsurf(nop(7+iad,i))=nop(4,i) endif else !1D elements wd=width(nop(3,i)) x=cord(1,nop(3,i))-asin*wd/2.0 y=cord(2,nop(3,i))-acos*wd/2.0 z=cord(3,1) ao0=ao(nop(3,i)) call existnod(cord,ao,nsurf,nop,np,np3,x,y,z, & ao0,3,0,i,l) nop(5+iad,i)=l x=cord(1,nop(3,i))+asin*wd/2.0 y=cord(2,nop(3,i))+acos*wd/2.0 z=cord(3,1) ao0=ao(nop(3,i)) call existnod(cord,ao,nsurf,nop,np,np3,x,y,z, & ao0,3,0,i,l) nop(7+iad,i)=l nop(6+iad,i)=nop(3,i) nsurf(nop(6+iad,i))=nop(3,i) endif np3=np3+1 nop(4+iad,i)=np3 nsurf(nop(4+iad,i))=nop(2,i) cord(1,np3)=(cord(1,nop(3+iad,i)) & +cord(1,nop(5+iad,i)))/2.0 cord(2,np3)=(cord(2,nop(3+iad,i)) & +cord(2,nop(5+iad,i)))/2.0 cord(3,np3)=(cord(3,nop(3+iad,i)) & +cord(3,nop(5+iad,i)))/2.0 ao(np3)=(ao(nop(3+iad,i))+ao(nop(5+iad,i)))/2.0 np3=np3+1 nop(8+iad,i)=np3 nsurf(nop(8+iad,i))=nop(2,i) cord(1,np3)=(cord(1,nop(1+iad,i)) & +cord(1,nop(7+iad,i)))/2.0 cord(2,np3)=(cord(2,nop(1+iad,i)) & +cord(2,nop(7+iad,i)))/2.0 cord(3,np3)=(cord(3,nop(1+iad,i)) & +cord(3,nop(7+iad,i)))/2.0 ao(np3)=(ao(nop(1+iad,i))+ao(nop(7+iad,i)))/2.0 nen=8 endif c c----------------generate nodes below surface c if(nen.gt.5) then do k=1,nen if(mod(k,2).eq.0) then nn=1 else nn=2 endif do m=1,nn x=cord(1,nop(k+iad,i)) y=cord(2,nop(k+iad,i)) if(m.eq.1) then z=ao(nop(k+iad,i)) else z=(cord(3,nop(k+iad,i)) & +ao(nop(k+iad,i)))/2.0 endif ao0=ao(nop(k+iad,i)) c c-------------------------check whether it exists c call existnod(cord,ao,nsurf,nop,np,np3,x,y,z, & ao0,k,iad,i,l) if(m.eq.1) then nop(k,i)=l else nop(k/2+nen+1,i)=l endif enddo enddo endif endif 65 enddo np=np3 print *,'npm,np=',npm,np endif c print *,'After generation ....' c do i=1,ne c if(nop(20,i).ne.0) then c nen=20 c else if(nop(15,i).ne.0.and.nop(16,i).eq.0) then c nen=15 c else if(nop(13,i).ne.0.and.nop(14,i).eq.0) then c nen=13 c else if(nop(10,i).ne.0.and.nop(11,i).eq.0) then c nen=10 c else if(nop(8,i).ne.0.and.nop(9,i).eq.0) then c nen=8 c else if(nop(6,i).ne.0.and.nop(7,i).eq.0) then c nen=6 c else if(nop(3,i).ne.0.and.nop(4,i).eq.0) then c nen=3 c else c nen=5 c endif c print *,i,imat(i),nen c print *,i,(nop(k,i),k=1,nen) c enddo c stop c c Creating array for diffusivities. k=1 DO WHILE(EDIFF(k).ne.-1) k=k+1 ENDDO DO i=1,ne j=1 DO WHILE ((imat(i).ne.ediff(j)).and.(j.lt.k)) j=j+1 ENDDO if (imat(i).eq.ediff(j)) then diffs(1,i)=hdiff(j) diffs(2,i)=vdiff(j) else diffs(1,i)=hdiff(1) diffs(2,i)=vdiff(1) endif ENDDO return 10 write(*,'(/a/)') ' *** Error in geometry file ***' stop end c---------------------------------------------------------------existnod subroutine existnod(cord,ao,nsurf,nop,np,np3,x,y,z,ao0,k,iad,i, & l) c----------------------------------------------------------------------c c purpose: c c To check whether the created node exists in the extension c c list of nodes. c c----------------------------------------------------------------------c real cord(3,1) real*8 ao(1) integer*2 nsurf(1),nop(20,1) c do l=np+1,np3 if(x.eq.cord(1,l).and.y.eq.cord(2,l). & and.z.eq.cord(3,l)) return enddo np3=np3+1 l=np3 cord(1,np3)=x cord(2,np3)=y cord(3,np3)=z ao(np3)=ao0 nsurf(np3)=nop(k+iad,i) return end c-----------------------------------------------------------getpollwins subroutine getpollwins(pollwin,masspA) c---------------------------------------------------------------------c c purpose: c c To read pollutant windows. c c---------------------------------------------------------------------c INCLUDE '3duparms.cb' include '3dpolls.cb' real pollwin(4,1),masspA(1) c do i=1,num_pollwins read(4,*,err=10) (pollwin(j,i),j=1,4) masspA(i)=mass_pp enddo return 10 write(*,'(/a/)') ' *** Error in pollutant file ***' stop end c-------------------------------------------------------------poll_skip SUBROUTINE poll_skip(poll,key) c---------------------------------------------------------------------c c purpose: c c To skip the first pollutant for hot start in unsteady case. c c---------------------------------------------------------------------c INCLUDE '3duparms.cb' INCLUDE '3dpolls.cb' REAL poll(4,1) c poll(1,i) = Vol of discharge at window i c poll(2,i) = Conc of discharge at window i c poll(3,i) = Depth of discharge at window i c poll(4,i) = Radius of discharge at window i INTEGER d, w, key key=0 c *** Read in the new pollutant releases. DO 60 w=1, num_pollwins READ(4,*,ERR=200,end=205) d,(poll(i,d),i=1,4) IF (d.gt.num_pollwins.or.d.lt.1) THEN PRINT *,'ERROR - Pollutant window not as '// 1 'expected.' STOP ENDIF 60 CONTINUE key=1 65 FORMAT(I2, F10.5, E10.1, 2F10.3) RETURN 200 WRITE(*,*) 'Error in pollutants!' STOP 205 WRITE(*,*) 'Run out of pollutant data' return END c----------------------------------------------------------getopenbound subroutine getopenbound(no_open,opbd) c---------------------------------------------------------------------c c purpose: c c To read open boundaries. c c---------------------------------------------------------------------c c no_open -- number of open boundaries specified. c c opbd(1,i) -- x-component of the first point c c opbd(2,i) -- y-component of the first point c c opbd(3,i) -- x-component of the second point c c opbd(4,i) -- y-component of the second point c c---------------------------------------------------------------------c integer no_open real opbd(4,1) c do i=1,no_open read(7,*) (opbd(k,i),k=1,4) enddo return end c---------------------------------------------------------getareaparams subroutine getareaparams(nply,nump,cord,np,no_poly,ii) c---------------------------------------------------------------------c c purpose: c c To read parameters of areas in which particles are counted. c c---------------------------------------------------------------------c include '3duparms.cb' dimension nump(1),nply(2,1),cord(3,1),xy(2,500) character*1 line common /cline/line(160) common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 c npadd=np ntmp=nin nin=3 rewind nin ii=1 no_poly=0 call find('AREAS',5,key) if(key.eq.0) then p_count=.false. return endif p_count=.true. 10 call free if(line(1).eq.' ') goto 20 call freei('N',nns,1) call freei('R',nump(ii),nns) if(nump(ii).eq.0) then call freer('Y',xy(1,1),2*nns) do k=1,nns npadd=npadd+1 if(npadd-np.gt.500) goto 30 nump(ii+k-1)=npadd cord(1,npadd)=xy(1,k) cord(2,npadd)=xy(2,k) enddo endif no_poly=no_poly+1 nply(1,no_poly)=nns nply(2,no_poly)=ii ii=ii+nns goto 10 20 ii=ii-1 nin=ntmp print *,'nply,nump=',nply(1,1),nply(2,1),(nump(i),i=1,ii) return 30 write(*,'(a)') ' *** Total no of polygon nodes > 500 ***' stop end c--------------------------------------------------------getoutputflags subroutine getoutputflags c---------------------------------------------------------------------c c purpose: c c To read output flags. c c---------------------------------------------------------------------c include '3duparms.cb' character*1 line character*4 flag common /cline/line(160) common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 data dlevel/100.0,200.0,500.0,1000.0,2000.0/ data clevel/100.0,500.0,1000.0,5000.0,20000.0/ c ntmp=nin nin=3 rewind nin call find('OUTPUT',6,key) if(key.eq.0) then dilution=.true. write(*,*) 'Dilution' concentration=.true. write(*,*) 'Concentration' return endif 10 call free if(line(1).eq.' ') goto 20 if(line(1).eq.'P') then call freeh('N',flag,4,1) if(flag(1:4).eq.'CONC') then dilution=.false. write(*,*) 'No Dilution' concentration=.true. write(*,*) 'Concentration' elseif(flag(1:4).eq.'BOTH') then dilution=.true. write(*,*) 'Dilution' concentration=.true. write(*,*) 'Concentration' else dilution=.true. write(*,*) 'Dilution' concentration=.false. write(*,*) 'No Concentration' endif elseif(line(1).eq.'D') then call freer('S',dlevel,5) write(*,'(a,5f12.1)') 'Input Dilution Levels=',dlevel elseif(line(1).eq.'C') then call freer('S',clevel,5) write(*,'(a,5f12.1)') 'Input Concentration Levels=',clevel endif goto 10 20 nin=ntmp return end c------------------------------------------------------getconservparams subroutine getconservparams c---------------------------------------------------------------------c c purpose: c c To read conservative option & parameters. c c---------------------------------------------------------------------c include '3duparms.cb' character*1 line character*4 flag common /cline/line(160) common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 c ntmp=nin nin=3 rewind nin call find('CONSERVATION',12,key) if(key.eq.0) then dieoff=.false. dieoffc=.false. return endif call free call freeh('F',flag,1,1) if(flag(1:1).eq.'Y') then dieoff=.true. dieoffc=.false. else if(flag(1:1).eq.'C') then dieoff=.true. dieoffc=.true. endif call freer('W',t90_s,1) call freer('P',t90_d,1) call freer('H',tdpth,1) if(tdpth.eq.0.0) tdpth=4.0 nin=ntmp return end c---------------------------------------------------------getgridouputs subroutine getgridoutputs c---------------------------------------------------------------------c c purpose: c c To read which grids should have bindump files c c added BMM 6/8/96 c c---------------------------------------------------------------------c include '3duparms.cb' character*1 line character*4 flag common /cline/line(160) common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 c ntmp=nin nin=3 rewind nin call find('GRID OUTPUT',11,key) if(key.eq.0) then return endif do i=1,50 pbindump(i)=.FALSE. enddo 10 call free if(line(1).eq.' ') goto 90 call freei('P',n,1) pbindump(n)=.TRUE. goto 10 90 nin=ntmp return end c--------------------------------------------------------------gettitle subroutine gettitle c---------------------------------------------------------------------c c purpose: c c To read the title of the problem. c c---------------------------------------------------------------------c include '3duparms.cb' CHARACTER*80 dummstr character*1 line common /cline/line(160) common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 c ntmp=nin nin=3 rewind nin call find('TITLE',5,key) if(key.eq.0) return call free call freeh(' ',dummstr,1,80) print *,dummstr(1:lenstr(dummstr)) nin=ntmp return end c---------------------------------------------------------getgridparams subroutine getgridparams c---------------------------------------------------------------------c c purpose: c c To read grid parameters. c c---------------------------------------------------------------------c include '3duparms.cb' INCLUDE '3dgeom.cb' character*1 line integer numxx(50),numyy(50),numzz(50),resxyy(50),reszz(50) real goee(50),gonn(50),gass(50),gacc(50),dimxyy(50),dimzz(50) common /mgird/goee,gonn,gass,gacc,numxx,numyy,numzz,dimxyy, * dimzz,resxyy,reszz common /cline/line(160) common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 c ntmp=nin nin=3 rewind nin no_grids=1 call find('GRID PARAM',10,key) if(key.eq.0) then write(*,*) ' ** Grid parameters are not defined **' stop endif 10 call free if(line(1).eq.' ') goto 90 call freei('R',n,1) if(n.ne.0) then call free else n=1 endif no_grids=max(no_grids,n) call freer('E',goE,1) call freer('N',goN,1) call freer('S',gaS,1) call freer('C',gaC,1) call free call freei('X',numX,1) call freei('Y',numY,1) call freei('Z',numZ,1) call free call freer('Y',dimXY,1) call freer('Z',dimZ,1) call free call freei('Y',resXY,1) call freei('Z',resZ,1) print *,'goE,goN=',goE,goN print *,'gaS,gaC=',gaS,gaC print *,'numX,numY,numZ=',numX,numY,numZ print *,'dimXY,dimZ=',dimXY,dimZ print *,'resXY,resZ=',resXY,resZ goee(n)=goE gonn(n)=goN gass(n)=gaS gacc(n)=gaC numxx(n)=numX numyy(n)=numY numzz(n)=numZ dimxyy(n)=dimXY dimzz(n)=dimZ resxyy(n)=resXY reszz(n)=resZ goto 10 90 nin=ntmp return end c------------------------------------------------------getgeneralparams subroutine getgeneralparams c---------------------------------------------------------------------c c purpose: c c To read general parameters. c c---------------------------------------------------------------------c include '3duparms.cb' character*1 line common /cline/line(160) common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 c ntmp=nin nin=3 rewind nin call find('GENERAL PARAM',13,key) if(key.eq.0) then write(*,*) ' ** General parameters are not defined **' stop endif call free call freei('S',maxoutputs,1) call free call freei('S',newcondTS,1) call free call freei('S',outputTS,1) call free call freei('T',dt,1) if(dt.eq.0) then call freei('1',dt,1) call freei('2',dtnew,1) else dtnew=dt endif call free call freei('E',tot_dead_age,1) cBMM additional lines can be read in for varying diffusivity ver 2.2 do i=1,10 EDIFF(i)=-1 end do call free call freer('Y',HDIFF(1),1) if (HDIFF(1).ge.0.0) then EDIFF(1)=1 call free call freer('Y',VDIFF,1) else j=-1*int(HDIFF(1)) do i=1,j call free call freei('E',EDIFF(i),1) call free call freer('Y',HDIFF(i),1) call free call freer('Y',VDIFF(i),1) enddo endif call free call freer('E',mass_pp,1) call freei('R',lnum_ps,1) if(lnum_ps.eq.0) lnum_ps=2000000000 call free call freer('N',c_init,1) nin=ntmp return end c-------------------------------------------------------getanalysisoptn subroutine getanalysisoptn c---------------------------------------------------------------------c c purpose: c c To read analysis option. c c---------------------------------------------------------------------c include '3duparms.cb' include 'rise_fall.cb' character*1 line character*40 flag common /cline/line(160) common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 c hsteadystate=.true. psteadystate=.true. plumeonly=.false. floatableonly=.false. euler=.true. compact=.false. ntmp=nin nin=3 rewind nin call find('ANALYSIS OPT',12,key) if(key.eq.0) then return endif call free if(line(1).eq.' ') goto 90 call freeh('W',flag,1,8) if(flag(1:8).eq.'UNSTEADY') then hsteadystate=.false. else if(flag(1:7).eq.'COMPACT') then compact=.true. hsteadystate=.false. endif call freeh('E',flag,1,8) if(flag(1:8).eq.'UNSTEADY') then psteadystate=.false. endif call freeh('S',flag,1,5) if(flag(1:5).eq.'PLUME'.and..not.settles) then plumeonly=.true. rr=-1.0 rf=-1.0 elseif(flag(1:5).eq.'FLOAT'.and..not.settles) then floatableonly=.true. rr=100.0 rf=-1.0 y0r=y0r*100.0/5.0 y1r=y1r*100.0/5.0 endif call free if(line(1).eq.' ') goto 90 call freeh('N',flag,1,5) if(flag(1:5).eq.'SECAN') euler=.false. if(.not.plumeonly) then call free if(line(1).eq.' ') goto 90 if(.not.floatableonly) then call freer('F',rf,1) if(rf.eq.0.0) rf=-1.0 call freer('A',y0f,1) call freer('B',y1f,1) endif if(.not.settles) then call freer('R',rr,1) if(rr.eq.0.0) rr=-1.0 call freer('C',y0r,1) call freer('D',y1r,1) endif if(y1r.eq.0.0.or.y1f.eq.0.0) then write(*,*) ' *** b and d cannot be zero. ***' stop endif if(rr.gt.100.0.or.rf.gt.100.0) then write(*,*) ' *** rf and rr having a range of 0 - 100 ***' stop endif endif 90 nin=ntmp return end