diff --git a/source/3drwalk.for b/source/3drwalk.for index d5953c8..b6a5379 100644 --- a/source/3drwalk.for +++ b/source/3drwalk.for @@ -2,18 +2,36 @@ c----------------------------------------------------------------RWALK_3D program RWALK_3D c-----------------------------------------------------------------------c c Purpose: c -c RWALK_3D - Random Walk Model c +c RWALK_3D - Random Walk Model +c +c History +c Originally written for the Sydney Deepwater Outfall +c Original version written by Brett Miller +c Updated by You Cong Wang circa 1996 (version 2.2) +c Used between 1996 and 2007 on many studies (up to Burwood) +c 2019 Modifications +c Starting with Version 2.5 (last used on Burwood Beach) +c Modifications to read RMA lastest file formats +c +c Planned modifications +c Modify to make use of Fortran 2000 dynamic arrays +c Modify to use class structure variables to improve readability c-----------------------------------------------------------------------c + +c Define each commonblock in external .cb files to ensure they are +c the same in the main block and the subroutines. + include '3duparms.cb' include '3dpolls.cb' include '3dgeom.cb' include '3dmesh.cb' include '3dpoints.cb' include 'rise_fall.cb' + parameter (ntotal=50000000) - character fname*100,buff*100,partl*100,fltnam*100 - character plmnam*100,cext*5 - real r + character fname*100,buff*100,partl*100,fltnam*100 + character plmnam*100,cext*5 + real r integer t,NUMPARTS,max_numps,v_time logical flipflag,restout,existing integer numxx(50),numyy(50),numzz(50),resxyy(50),reszz(50) @@ -32,47 +50,20 @@ c-----------------------------------------------------------------------c data isd/0/, isd1/0/ data y0f/10.6/, y1f/-11.34/, y0r/1.5951/, y1r/-0.97525/ data rf/30.0/, rr/5.0/ -c + write(*,'(a)') '-----------------------------------------------' - write(*,'(a)') '--- ---' - write(*,'(a)') '--- 3D_RWALK Version 2.2 ---' - write(*,'(a)') '--- 06/08/96 ---' - write(*,'(a)') '--- ---' - write(*,'(a)') '--- Y.C. Wang & B. Miller ---' - write(*,'(a)') '--- Water Research Lab. U.N.S.W. ---' - write(*,'(a)') '--- Compiled with ifort 12/7/07 ---' - write(*,'(a)') '--- ---' + write(*,'(a)') '--- 3D_RWALK v2.6 ---' + write(*,'(a)') '--- UNSW Sydney ---' + write(*,'(a)') '--- Water Research Laboratory ---' write(*,'(a)') '-----------------------------------------------' -c------------------------------------------------------------------ -c Use Defined variables to echo where and when this program was compiled -c added to code by AF 12/8/04 -c -c WRITE(*,'(2A)') 'This program was compiled in: ' -c 1 ,PATH -c WRITE(*,'(2A)') ' with filename : ' -c 1 ,__FILE__ -c WRITE(*,'(2A)') 'It was compiled on : ' -c 1 ,__DATE__ -c WRITE(*,'(2A)') ' at time : ' -c 1 ,__TIME__ -c -c IF (INDEX -c 1 (PATH,'/wrlnmr/') -c 2 .eq.0) THEN -c PRINT *,'This is not being run from the QA area' -c ELSE -c PRINT*, 'This is being run from the QA area' -c ENDIF -c -c print *,'Running the program' -c---------------------------------------------------------------- +c Dynamic array setup. Don't change mtot=ntotal idir=mtot call setcst -c -c-------open required files -c + +c Filestreams used + c 21 : output stream for RWD graphics file c 22 : output stream for RWC graphics file c 23 : output stream for RWZ graphics file @@ -89,43 +80,45 @@ c 20 : Plume Bin Dump file c 30 : settling particles file c 35 : particle tracking path file c 36 : particle counting file -c - WRITE(*,*) ' Enter start-stamp (yyyymmddhhmm): ' - READ(*,'(A12)') start_stamp - WRITE(*,'(A18,A12,A1)') ' START STAMP is "',start_stamp,'"' - WRITE(*,*) - fname='920820' - call askstr(' Enter output graphics filename',fname) +c Read in the user paramters + + WRITE(*,*) ' Enter start-stamp (yyyymmddhhmm): ' + READ(*,'(A12)') start_stamp + WRITE(*,'(A18,A12,A1)') ' START STAMP is "',start_stamp,'"' + WRITE(*,*) + + fname='920820' + call askstr(' Enter output graphics filename',fname) partl=fname dimopt='10' - call askstr(' Enter RMA-2 or RMA-10 option',dimopt) + call askstr(' Enter RMA-2 or RMA-10 option',dimopt) if(dimopt(1:2).eq.'10') then - fname='sydney.3dg' - call askstr(' Enter 3d geometry filename',fname) - open(9,file=fname,status='old',form='unformatted',err=10) - fname='sydney.res' - call askstr(' Enter RMA-10 velocity filename',fname) - open(2,file=fname,status='old',form='unformatted',err=10) + fname='sydney.3dg' + call askstr(' Enter 3d geometry filename',fname) + open(9,file=fname,status='old',form='unformatted',err=10) + fname='sydney.res' + call askstr(' Enter RMA-10 velocity filename',fname) + open(2,file=fname,status='old',form='unformatted',err=10) else - fname='sydney.geo' - call askstr(' Enter 2D geometry filename',fname) - open(9,file=fname,status='old',form='unformatted',err=10) + fname='sydney.geo' + call askstr(' Enter 2D geometry filename',fname) + open(9,file=fname,status='old',form='unformatted',err=10) buff='1.0,1.0' - call askstr(' Enter scales in x and y directions',buff) + call askstr(' Enter scales in x and y directions',buff) read(buff,*) xscale,yscale - fname='sydney.res' - call askstr(' Enter RMA-2 velocity filename',fname) - open(2,file=fname,status='old',form='unformatted',err=10) + fname='sydney.res' + call askstr(' Enter RMA-2 velocity filename',fname) + open(2,file=fname,status='old',form='unformatted',err=10) endif fname='usrparam.rw' - call askstr(' Enter user parameters filename',fname) - open(3,file=fname,status='old',form='formatted',err=10) + call askstr(' Enter user parameters filename',fname) + open(3,file=fname,status='old',form='formatted',err=10) fname='polls' - call askstr(' Enter pollutants filename',fname) - open(4,file=fname,status='old',form='formatted',err=10) + call askstr(' Enter pollutants filename',fname) + open(4,file=fname,status='old',form='formatted',err=10) fname='no' - call askstr(' Enter open boundary filename',fname) + call askstr(' Enter open boundary filename',fname) if(fname(1:1).ne.'n'.and.fname(1:1).ne.'N') then open(7,file=fname,status='old',form='formatted',err=10) rewind 7 @@ -134,7 +127,7 @@ c openbd=.false. endif fname='no' - call askstr(' Enter input restart filename',fname) + call askstr(' Enter input restart filename',fname) if(fname(1:1).ne.'n'.and.fname(1:1).ne.'N') then open(10,file=fname,status='unknown',form='unformatted',err=10) restart=.true. @@ -142,45 +135,44 @@ c restart=.false. endif fname='no' - call askstr(' Enter output restart filename',fname) + call askstr(' Enter output restart filename',fname) if(fname(1:1).ne.'n'.and.fname(1:1).ne.'N') then open(11,file=fname,status='unknown',form='unformatted',err=10) restout=.true. else restout=.false. endif - fname='no' - CALL askstr(' Enter Floatables Bin-Dump filename',fname) + fname='no' + CALL askstr(' Enter Floatables Bin-Dump filename',fname) fltnam=fname if(fname(1:1).ne.'n'.and.fname(1:1).ne.'N') then - fbindump=.TRUE. + fbindump=.TRUE. else fbindump=.false. endif - fname='no' - CALL askstr(' Enter Plumes Bin-Dump filename',fname) + fname='no' + CALL askstr(' Enter Plumes Bin-Dump filename',fname) plmnam=fname if(fname(1:1).ne.'n'.and.fname(1:1).ne.'N') then - do itemp=1,50 - pbindump(itemp)=.TRUE. - enddo + do itemp=1,50 + pbindump(itemp)=.TRUE. + enddo else - do itemp=1,50 - pbindump(itemp)=.false. - enddo + do itemp=1,50 + pbindump(itemp)=.false. + enddo endif -c print*,plmnam,'pbindump=',pbindump(1) - fname='no' - CALL askstr(' Enter settleable particles filename',fname) + fname='no' + CALL askstr(' Enter settleable particles filename',fname) if(fname(1:1).ne.'n'.and.fname(1:1).ne.'N') then - inquire(file=fname(1:lenstr(fname)),EXIST=existing) + inquire(file=fname(1:lenstr(fname)),EXIST=existing) if(existing) then - open(30,file=fname(1:lenstr(fname)),status='old') + open(30,file=fname(1:lenstr(fname)),status='old') close(30,status='delete') endif - OPEN(30,FILE=fname, ACCESS='direct', - 1 STATUS='new',FORM='unformatted',RECL=8, ERR=10) - settles=.TRUE. + OPEN(30,FILE=fname, ACCESS='direct', + 1 STATUS='new',FORM='unformatted',RECL=8, ERR=10) + settles=.TRUE. rr=-1.0 rf=100.0 y0f=y0f*100.0/30.0 @@ -188,19 +180,17 @@ c print*,plmnam,'pbindump=',pbindump(1) else settles=.false. endif - fname='no' - CALL askstr(' Enter particle tracking path filename',fname) + fname='no' + CALL askstr(' Enter particle tracking path filename',fname) if(fname(1:1).ne.'n'.and.fname(1:1).ne.'N') then - inquire(file=fname(1:lenstr(fname)),EXIST=existing) + inquire(file=fname(1:lenstr(fname)),EXIST=existing) if(existing) then - open(35,file=fname(1:lenstr(fname)),status='old') + open(35,file=fname(1:lenstr(fname)),status='old') close(35,status='delete') endif -ctmp OPEN(35,FILE=fname, ACCESS='direct', -ctmp 1 STATUS='new',FORM='unformatted',RECL=8, ERR=10) - OPEN(35,FILE=fname, - 1 STATUS='new',FORM='formatted',ERR=10) - tracking=.TRUE. + OPEN(35,FILE=fname, + 1 STATUS='new',FORM='formatted',ERR=10) + tracking=.TRUE. else tracking=.false. endif @@ -217,7 +207,7 @@ c call echoallparameters write(*,'(/a/)') 'Loading the geometry of FE mesh ...' - rewind 9 + rewind 9 read(9) np,ne rewind 9 call defini('nop ',nnop ,10,ne+4) @@ -250,7 +240,6 @@ c DAYOFY=EGDF-EGDS+1 TOFDAY=int(ATIME/100)+mod(ATIME,100)/60.0 call DIEOFFS(0.0,DAYOFY,T90_D) -c print *,'T90_D=',T90_D T90_D=3600.0*T90_D endif print *,'Euler=',euler @@ -272,22 +261,22 @@ c call defini('nopt',nnopt,no_poly,mlp) if(p_count) then - fname=partl(1:lenstr(partl))//'.par' - CALL askstr(' Enter particle counting filename',fname) - inquire(file=fname(1:lenstr(fname)),EXIST=existing) + fname=partl(1:lenstr(partl))//'.par' + CALL askstr(' Enter particle counting filename',fname) + inquire(file=fname(1:lenstr(fname)),EXIST=existing) if(existing) then - open(36,file=fname(1:lenstr(fname)),status='old') + open(36,file=fname(1:lenstr(fname)),status='old') close(36,status='delete') endif - OPEN(36,FILE=fname, - 1 STATUS='new',FORM='formatted', ERR=10) + OPEN(36,FILE=fname, + 1 STATUS='new',FORM='formatted', ERR=10) write(36,'(a)') ' Start_stamp='//start_stamp write(36,'(a,i15)') ' outputTS=',outputTS write(36,'(a,i15)') 'tot_dead_age=',tot_dead_age write(36,'(a,i15)') 'No. of Areas=',no_poly - endif + endif write(*,'(/a/)') 'Loading the pollutant windows ...' read(4,*) poll_type,num_pollwins @@ -378,14 +367,14 @@ c if(keyp.eq.0) goto 50 endif - call getflowconditions(ia(nnsur),ia(nAvel),ia(ncord), + call getflowconditions(ia(nnsur),ia(nAvel),ia(ncord), * ia(nao),ia(nwsel),ia(nwsl1),ia(nic), * ia(ndhdt),ia(nnop),ia(nilst),keyf) - call getpollconditions(ia(nApol),ia(nmasp),ia(nmasA), + call getpollconditions(ia(nApol),ia(nmasp),ia(nmasA), * NUMPARTS,max_numpsA,keyp) if(keyf.eq.0.or.keyp.eq.0) goto 50 if(.not.restart) then - if(.not.hsteadystate) then + if(.not.hsteadystate) then call getflowconditions(ia(nnsur), * ia(nBvel),ia(ncord),ia(nao),ia(nwsel),ia(nwsl1), * ia(nic),ia(ndhdt),ia(nnop),ia(nilst),keyf) @@ -393,7 +382,7 @@ c else nBvel=nAvel endif - if(.not.psteadystate) then + if(.not.psteadystate) then call getpollconditions(ia(nBpol),ia(nmasp),ia(nmasA), * NUMPARTS,max_numpsB,keyp) if(keyp.eq.0) goto 50 @@ -412,32 +401,32 @@ c else write(cext,'(i2.2)') m endif - inquire(file=partl(1:lenstr(partl))//'.rwd'//cext(1:2), + inquire(file=partl(1:lenstr(partl))//'.rwd'//cext(1:2), * EXIST=existing) if(existing) then - open(21,file=partl(1:lenstr(partl))//'.rwd'//cext(1:2), + open(21,file=partl(1:lenstr(partl))//'.rwd'//cext(1:2), * status='old') close(21,status='delete') endif open(21,file=partl(1:lenstr(partl))//'.rwd'//cext(1:2), * access='direct',status='new',form='unformatted',recl=1) - inquire(file=partl(1:lenstr(partl))//'.rwc'//cext(1:2), + inquire(file=partl(1:lenstr(partl))//'.rwc'//cext(1:2), * EXIST=existing) if(existing) then - open(22,file=partl(1:lenstr(partl))//'.rwc'//cext(1:2), + open(22,file=partl(1:lenstr(partl))//'.rwc'//cext(1:2), * status='old') close(22,status='delete') endif - open(22,file=partl(1:lenstr(partl))//'.rwc'//cext(1:2), + open(22,file=partl(1:lenstr(partl))//'.rwc'//cext(1:2), * access='direct',status='new',form='unformatted',recl=1) - inquire(file=partl(1:lenstr(partl))//'.rwz'//cext(1:2), + inquire(file=partl(1:lenstr(partl))//'.rwz'//cext(1:2), * EXIST=existing) if(existing) then - open(23,file=partl(1:lenstr(partl))//'.rwz'//cext(1:2), + open(23,file=partl(1:lenstr(partl))//'.rwz'//cext(1:2), * status='old') close(23,status='delete') endif - open(23,file=partl(1:lenstr(partl))//'.rwz'//cext(1:2), + open(23,file=partl(1:lenstr(partl))//'.rwz'//cext(1:2), * access='direct',status='new',form='unformatted',recl=1) goE=goee(m) goN=gonn(m) @@ -451,27 +440,27 @@ c resXY=resxyy(m) resZ=reszz(m) call addoutputheader(21) - call addoutputheader(22) - call addoutputheader(23) + call addoutputheader(22) + call addoutputheader(23) close(21) close(22) close(23) - IF (pbindump(m)) THEN - inquire(file=plmnam(1:lenstr(plmnam))//cext(1:2), + IF (pbindump(m)) THEN + inquire(file=plmnam(1:lenstr(plmnam))//cext(1:2), * EXIST=existing) if(existing) then - open(20,file=plmnam(1:lenstr(plmnam))//cext(1:2), + open(20,file=plmnam(1:lenstr(plmnam))//cext(1:2), * status='old') close(20,status='delete') endif c BMM 070810 Addition of another file type for storing the bin file data - OPEN(37,FILE=plmnam(1:lenstr(plmnam))//'-block'//cext(1:2), + OPEN(37,FILE=plmnam(1:lenstr(plmnam))//'-block'//cext(1:2), 1 STATUS='unknown',FORM='unformatted') - OPEN(20,FILE=plmnam(1:lenstr(plmnam))//cext(1:2), + OPEN(20,FILE=plmnam(1:lenstr(plmnam))//cext(1:2), cdrc 1 ACCESS='direct',STATUS='new',FORM='unformatted',RECL=1) 1 ACCESS='direct',STATUS='new',FORM='unformatted',RECL=4) r=FLOAT(numX*resXY) @@ -488,15 +477,15 @@ cdrc 1 ACCESS='direct',STATUS='new',FORM='unformatted',RECL=1) ENDIF - IF (fbindump) THEN - inquire(file=fltnam(1:lenstr(fltnam))//cext(1:2), + IF (fbindump) THEN + inquire(file=fltnam(1:lenstr(fltnam))//cext(1:2), * EXIST=existing) if(existing) then - open(19,file=fltnam(1:lenstr(fltnam))//cext(1:2), + open(19,file=fltnam(1:lenstr(fltnam))//cext(1:2), * status='old') close(19,status='delete') endif - OPEN(19,FILE=fltnam(1:lenstr(fltnam))//cext(1:2), + OPEN(19,FILE=fltnam(1:lenstr(fltnam))//cext(1:2), cdrc 1 ACCESS='direct',STATUS='new',FORM='unformatted',RECL=1) 1 ACCESS='direct',STATUS='new',FORM='unformatted',RECL=4) r=FLOAT(numX*resXY) @@ -508,7 +497,7 @@ cdrc 1 ACCESS='direct',STATUS='new',FORM='unformatted',RECL=1) WRITE(19,REC=4) 0.0 WRITE(19,REC=4+maxoutputs) 0.0 close(19) - ENDIF + ENDIF enddo if(settles) write(30,rec=1) no_settles,mass_pp @@ -571,12 +560,12 @@ c * no_open,NUMPARTS,t,t_p,v_time,ia(ndifs)) if(v_time.eq.newcondTS) then v_time=0 - if(.not.hsteadystate) then + if(.not.hsteadystate) then call getflowconditions(ia(nnsur), * ia(nBvel),ia(ncord),ia(nao),ia(nwsel),ia(nwsl1), * ia(nic),ia(ndhdt),ia(nnop),ia(nilst),keyf) endif - if(.not.psteadystate) then + if(.not.psteadystate) then call getpollconditions(ia(nBpol),ia(nmasp), * ia(nmasA),NUMPARTS,max_numpsB,keyp) endif @@ -603,12 +592,12 @@ c * no_open,NUMPARTS,t,t_p,v_time,ia(ndifs)) if(v_time.eq.newcondTS) then v_time=0 - if(.not.hsteadystate) then + if(.not.hsteadystate) then call getflowconditions(ia(nnsur), * ia(nAvel),ia(ncord),ia(nao),ia(nwsel),ia(nwsl1), * ia(nic),ia(ndhdt),ia(nnop),ia(nilst),keyf) endif - if(.not.psteadystate) then + if(.not.psteadystate) then call getpollconditions(ia(nApol),ia(nmasp), * ia(nmasA),NUMPARTS,max_numpsA,keyp) endif @@ -665,20 +654,17 @@ c endif open(21,file=partl(1:lenstr(partl))//'.rwd'//cext(1:2), * access='direct',status='old',form='unformatted',recl=1) - open(22,file=partl(1:lenstr(partl))//'.rwc'//cext(1:2), + open(22,file=partl(1:lenstr(partl))//'.rwc'//cext(1:2), * access='direct',status='old',form='unformatted',recl=1) - open(23,file=partl(1:lenstr(partl))//'.rwz'//cext(1:2), + open(23,file=partl(1:lenstr(partl))//'.rwz'//cext(1:2), * access='direct',status='old',form='unformatted',recl=1) if(fbindump) then - OPEN(19,FILE=fltnam(1:lenstr(fltnam))//cext(1:2), -cdrc 1 ACCESS='direct',STATUS='old',FORM='unformatted',RECL=1) - 1 ACCESS='direct',STATUS='old',FORM='unformatted',RECL=4) + OPEN(19,FILE=fltnam(1:lenstr(fltnam))//cext(1:2), + 1 ACCESS='direct',STATUS='old',FORM='unformatted',RECL=4) endif -c print*,plmnam,'pbindump=',pbindump(m) if(pbindump(m)) then - OPEN(20,FILE=plmnam(1:lenstr(plmnam))//cext(1:2), -crdc 1 ACCESS='direct',STATUS='old',FORM='unformatted',RECL=1) - 1 ACCESS='direct',STATUS='old',FORM='unformatted',RECL=4) + OPEN(20,FILE=plmnam(1:lenstr(plmnam))//cext(1:2), + 1 ACCESS='direct',STATUS='old',FORM='unformatted',RECL=4) endif goE=goee(m) goN=gonn(m) @@ -691,12 +677,11 @@ crdc 1 ACCESS='direct',STATUS='old',FORM='unformatted',RECL=1) dimZ=dimzz(m) resXY=resxyy(m) resZ=reszz(m) + call outputgraphicscodes(m) -c print*,'here' call addpointers(21) -c print*,'here2' - call addpointers(22) - call addpointers(23) + call addpointers(22) + call addpointers(23) close(21) close(22) close(23) @@ -736,7 +721,7 @@ c close(9) if(restart) close(10) if(restout) close(11) - if(settles) close(30) + if(settles) close(30) stop end