brettm
/
3DRWALK
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Code Tidy Up

Browse Source
master
Brett Miller 5 years ago
parent 2608e8689e
commit 1f24a75037
1 changed files with 131 additions and 146 deletions
Show Stats Download Patch File Download Diff File

277
source/3drwalk.for
Unescape Escape View File

@ -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

Write Preview
Loading…
Cancel
Save