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

First Commit

Browse Source
master
Brett Miller 6 years ago
parent 9ad3b9ca9f
commit 2608e8689e
22 changed files with 6691 additions and 0 deletions
Show Stats Download Patch File Download Diff File

1251
source/3dcondits.for
View File

File diff suppressed because it is too large Load Diff

70
source/3dgeneral.for
Unescape Escape View File

@ -0,0 +1,70 @@
c---------------------------------------------------------------calc_ab
SUBROUTINE calc_ab(a, b, E, N)
c---------------------------------------------------------------------c
c purpose: c
c To rotate and translate coordinates from E-N to grid. c
c---------------------------------------------------------------------c
INCLUDE '3dgeom.cb'
REAL E, N, a, b
REAL nt, et
nt= N-goN
et= E-goE
a = (nt*gaS+et*gaC)/dimXY
b = (nt*gaC-et*gaS)/dimXY
RETURN
END
c---------------------------------------------------------------calc_EN
subroutine calc_EN(E,N,x,y)
c---------------------------------------------------------------------c
c Purpose: c
c To rotate and translate coordinates from grid to E-N. c
c---------------------------------------------------------------------c
INCLUDE '3dgeom.cb'
real E,N,x,y
c
N=goN+x*dimXY*gaS+y*dimXY*gaC
E=goE+x*dimXY*gaC-y*dimXY*gaS
return
end
c----------------------------------------------------------------calc_c
SUBROUTINE calc_c(c, Z)
c---------------------------------------------------------------------c
c purpose: c
c To translate z-coordinates from E-N to grid. c
c---------------------------------------------------------------------c
INCLUDE '3dgeom.cb'
REAL c, Z
c=-Z/dimZ
RETURN
END
c--------------------------------------------------------------calc_XYZ
subroutine calc_XYZ(nop,cord,itype,pel,pnl,pzl,pil,pe,pn,pz)
c---------------------------------------------------------------------c
c purpose: c
c To interplate the global coordinates. c
c---------------------------------------------------------------------c
common /etype/inode(4)
common /shape/shap(20),shpx(20),shpy(20),shpz(20)
real pe,pn,pz,pel,pnl,pzl,cord(3,1),shap
integer*2 nop(20,1)
integer pil,it,nen,i,itype(1)
c
it=itype(pil)
nen=inode(it)
pe=0.0
pn=0.0
pz=0.0
call xn3(it,nen,pel,pnl,pzl)
do i=1,nen
pe=pe+shap(i)*cord(1,nop(i,pil))
pn=pn+shap(i)*cord(2,nop(i,pil))
pz=pz+shap(i)*cord(3,nop(i,pil))
enddo
return
end

13
source/3dgeom.cb
Unescape Escape View File

@ -0,0 +1,13 @@
c GEOMETRY COMMON BLOCK.
c Used as part of the random walk model.
c BMM 17/5/93.
REAL goN, goE, gaS, gaC
INTEGER numX, numY, numZ
REAL dimXY, dimZ
INTEGER resXY, resZ
COMMON /GEOM/ bathym, goN, goE, gaS, gaC, numX, numY, numZ,
* dimXY, dimZ, resXY, resZ

2
source/3dmesh.cb
Unescape Escape View File

@ -0,0 +1,2 @@
common /mesh/np,ne,npm,nes,n1_2d
integer np,ne,npm,nes,n1_2d

1014
source/3doutput-block.for
View File

File diff suppressed because it is too large Load Diff

10
source/3dpoints.cb
Unescape Escape View File

@ -0,0 +1,10 @@
c POINTS COMMON BLOCK.
c Used as part of the random walk model.
c BMM 17/5/93.
INTEGER numviews
INTEGER viewpoint(3)
COMMON /POINT/ numviews, viewpoint

14
source/3dpolls.cb
Unescape Escape View File

@ -0,0 +1,14 @@
c POLLUTANTS COMMON BLOCK
c Used as part of the random walk model.
INTEGER poll_type
INTEGER num_pollwins
COMMON /POLLS/ poll_type, num_pollwins
c Descriptions
c ============
c poll_type :
c 1 - Plume particles (no die off)
c 2 - Floatable particles
c

790
source/3drwalk.for
Unescape Escape View File

@ -0,0 +1,790 @@
c----------------------------------------------------------------RWALK_3D
program RWALK_3D
c-----------------------------------------------------------------------c
c Purpose: c
c RWALK_3D - Random Walk Model c
c-----------------------------------------------------------------------c
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
integer t,NUMPARTS,max_numps,v_time
logical flipflag,restout,existing
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 /dbsys/numa,next,idir,ipp(3)
common /etype/inode(4)
common /randm/isd, isd1
common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10,
* nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19
common mtot,npp,ia(ntotal)
integer isd, isd1
integer ADAY,AYEAR,AMONTH,ATIME,EGDS,EGDF
data inode/10,15,20,13/
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)') '-----------------------------------------------'
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----------------------------------------------------------------
mtot=ntotal
idir=mtot
call setcst
c
c-------open required files
c
c 21 : output stream for RWD graphics file
c 22 : output stream for RWC graphics file
c 23 : output stream for RWZ graphics file
c 9 : geometry file
c 2 : RMA10 velocity file
c 3 : input user parameters file
c 4 : input pollutants file
c 7 : input open-boundary file
c 8 : input hard-boundary file
c 10 : input restart file
c 11 : output restart file
c 19 : Floatable Bin Dump file
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)
partl=fname
dimopt='10'
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)
else
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)
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)
endif
fname='usrparam.rw'
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)
fname='no'
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
openbd=.true.
else
openbd=.false.
endif
fname='no'
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.
else
restart=.false.
endif
fname='no'
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)
fltnam=fname
if(fname(1:1).ne.'n'.and.fname(1:1).ne.'N') then
fbindump=.TRUE.
else
fbindump=.false.
endif
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
else
do itemp=1,50
pbindump(itemp)=.false.
enddo
endif
c print*,plmnam,'pbindump=',pbindump(1)
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)
if(existing) then
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.
rr=-1.0
rf=100.0
y0f=y0f*100.0/30.0
y1f=y1f*100.0/30.0
else
settles=.false.
endif
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)
if(existing) then
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.
else
tracking=.false.
endif
goto 20
10 len0=lnblnk(fname)
write(6,'(a)') '*** File: '''//fname(1:len0)//''' not found ***'
stop
c
c-------read data
c
20 write(*,'(/a/)') 'Loading the user parameters ...'
call getusrparams
print *,'plumeonly,floatableonly=',plumeonly,floatableonly
call echoallparameters
write(*,'(/a/)') 'Loading the geometry of FE mesh ...'
rewind 9
read(9) np,ne
rewind 9
call defini('nop ',nnop ,10,ne+4)
call defini('ilst',nilst,1,ne/2+4)
call defini('imat',nimat,1,ne/2+4)
call definr('cord',ncord,3,7*np) !leave more space for 2D case
call definr('ao ',nao ,2,7*np+12)
call defini('nsur',nnsur,1,np*7/2+12)
call definr('widt',nwidt,1,np*7+12)
c BMM extra line for diffusivities
call definr('difs',ndifs,2,ne+4)
print *,'np=',np,npm,ne
call getgeometry(ia(ncord),ia(nnop),ia(nnsur),ia(nimat),
& ia(nwidt),ia(nao),ia(nilst),ia(ndifs))
print *,'np=',np,npm,ne
call refinr('cord',ncord,3,np+800) !leave some space for polygon
call refinr('ao ',nao ,2,np+4)
call refini('nsur',nnsur,1,np/2+4)
call refinr('widt',nwidt,1,np+4)
call refinr('difs',ndifs,2,ne+4)
c
c-------calculate the day of year DAYOFY
c
if(dieoffc) then
READ(start_stamp,'(I4,2I2,I4)') ISY, ISM, ISD, IST
RMIN=-outputTS/60
CALL TINC2(IST,ISD,ISM,ISY,ATIME,ADAY,AMONTH,AYEAR,RMIN,IER,0)
call gdate(1,AMONTH,ADAY,AYEAR,IDWK,IDYR,IDMON,EGDF,IER,0)
call gdate(1,1,1,AYEAR,IDWK,IDYR,IDMON,EGDS,IER,0)
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
call defini('nply',nnply,2,1000)
call defini('nump',nnump,1,10000)
call getareaparams(ia(nnply),ia(nnump),ia(ncord),np,no_poly,
& nlengthp)
print *,'no_poly,nlengthp=',no_poly,nlengthp
call refini('nump',nnump,1,nlengthp)
call refini('nply',nnply,2,no_poly)
call defini('no_p',nno_p,1,no_poly)
call definr('ppm ',nppm ,1,no_poly)
c
c-------create an array storing number of particles in specified area with
c time tag
c
mlp=tot_dead_age/3600+2
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)
if(existing) then
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)
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
write(*,'(/a/)') 'Loading the pollutant windows ...'
read(4,*) poll_type,num_pollwins
call definr('pwin',npwin,4,num_pollwins)
call definr('masA',nmasA,1,num_pollwins)
call getpollwins(ia(npwin),ia(nmasA))
c
c-------divide FE elements into blocks
c
nx=2
ny=5
nz=1
nxyz=nx*ny*nz
call defini('ityp',nityp,1,ne)
call definr('ebox',nebox,6,ne)
call definr('egrp',negrp,6,nxyz)
call defini('nbeg',nnbeg,1,nxyz)
call defini('nend',nnend,1,nxyz)
call defini('nelg',nnelg,3,ne)
call definr('egps',negps,6,2000)
call defini('nbb ',nnbb ,1,2000)
call defini('nee ',nnee ,1,2000)
call defini('nblk',nnblk,2,2000)
call defini('nelb',nnelb,1,1)
call blocks(nx,ny,nz,np,ne,ia(nnop),ia(ncord),nxyz,ia(negrp),
* ia(nnbeg),ia(nnend),ia(nnelg),ia(nnbb),ia(nnee),ia(nnelb),
* ia(nnblk),ia(negps),ia(nityp),ia(nebox),m,num)
call refini('nelb',nnelb,1,num)
call delete('nelg')
call delete('nbeg')
call delete('nend')
call refinr('egrp',negrp,6,nxyz)
call refinr('egps',negps,6,m)
call refini('nbb ',nnbb ,1,m)
call refini('nee ',nnee ,1,m)
call refini('nblk',nnblk,2,m)
call locate('nelb',nnelb,mlt,num)
c
c-------read open boundaries
c
if(openbd) then
read(7,*) no_open
call definr('opbd',nopbd,4,no_open)
call getopenbound(no_open,ia(nopbd))
endif
call bcname
c
c-------setup boundary conditions
c
c
c-------main program loop
c
write(*,*) 'Start the walking ...'
no_settles=1
NUMPARTS=0
max_numpsA=0
max_numpsB=0
write(*,'(/a/)') 'Loading velocity field ...'
call definr('Avel',nAvel,3,np)
call definr('Bvel',nBvel,3,np)
call definr('vel ',nvel ,3,np)
call definr('dudt',ndudt,3,np)
call definr('dhdt',ndhdt,1,np)
call definr('wsel',nwsel,1,np)
call definr('wsl1',nwsl1,1,np)
call defini('ic ',nic ,1,np/2+4)
call definr('Apol',nApol,4,num_pollwins)
call definr('Bpol',nBpol,4,num_pollwins)
call definr('poll',npoll,4,num_pollwins)
call definr('pdot',npdot,1,num_pollwins)
call definr('masp',nmasp,1,num_pollwins)
call clearr(ia(npdot),1,num_pollwins)
call clearr(ia(ndhdt),1,np)
c
c-------initialise dudt(3,np)
c
if(.not.Euler) then
call clearr(ia(ndudt),3,np)
endif
c
c-------skip the first pollutant release for hot start
c
if(.not.psteadystate.and.restart) then
call getpollconditions(ia(nApol),ia(nmasp),ia(nmasA),
* NUMPARTS,max_numpsA,keyp)
if(keyp.eq.0) goto 50
endif
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),
* NUMPARTS,max_numpsA,keyp)
if(keyf.eq.0.or.keyp.eq.0) goto 50
if(.not.restart) 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)
if(keyf.eq.0) goto 50
else
nBvel=nAvel
endif
if(.not.psteadystate) then
call getpollconditions(ia(nBpol),ia(nmasp),ia(nmasA),
* NUMPARTS,max_numpsB,keyp)
if(keyp.eq.0) goto 50
else
nBpol=nApol
max_numpsB=max_numpsA
endif
endif
flipflag=.FALSE.
max_numps=max(max_numpsA,max_numpsB)
print *,'no_grids=',no_grids
do m=1,no_grids
if(no_grids.eq.1) then
cext=' '
else
write(cext,'(i2.2)') m
endif
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),
* 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),
* EXIST=existing)
if(existing) then
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),
* access='direct',status='new',form='unformatted',recl=1)
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),
* status='old')
close(23,status='delete')
endif
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)
gaS=gass(m)
gaC=gacc(m)
numX=numxx(m)
numY=numyy(m)
numZ=numzz(m)
dimXY=dimxyy(m)
dimZ=dimzz(m)
resXY=resxyy(m)
resZ=reszz(m)
call addoutputheader(21)
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),
* EXIST=existing)
if(existing) then
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),
1 STATUS='unknown',FORM='unformatted')
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)
WRITE(20,REC=1) r
r=FLOAT(numY*resXY)
WRITE(20,REC=2) r
r=FLOAT(numZ*resZ)
WRITE(20,REC=3) r
r=float(4+maxoutputs)
WRITE(20,REC=3) int(r)
WRITE(20,REC=4) 0.0
WRITE(20,REC=4+maxoutputs) 0.0
close(20)
ENDIF
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),
* status='old')
close(19,status='delete')
endif
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)
WRITE(19,REC=1) r
r=FLOAT(numY*resXY)
WRITE(19,REC=2) r
r=float(4+maxoutputs)
WRITE(19,REC=3) r
WRITE(19,REC=4) 0.0
WRITE(19,REC=4+maxoutputs) 0.0
close(19)
ENDIF
enddo
if(settles) write(30,rec=1) no_settles,mass_pp
if(restart) then
rewind 10
read(10) NUMPARTS
endif
call defini('prtI',nprtI,1,NUMPARTS+1)
call definr('prtE',nprtE,1,NUMPARTS+1)
call definr('prtN',nprtN,1,NUMPARTS+1)
call definr('prtZ',nprtZ,1,NUMPARTS+1)
call definr('prtM',nprtM,1,NUMPARTS+1)
call definr('prtA',nprtA,1,NUMPARTS+1)
call definr('prtV',nprtV,1,NUMPARTS+1)
c
c-------read data for restart
c
if(restart) then
call inputRS(NUMPARTS,max_numps,ia(nprtI),ia(nprtE),
* ia(nprtN),ia(nprtZ),ia(nprtM),ia(nprtA),ia(nprtV),ia(nBvel),
* ia(nwsl1),ia(nBpol),ia(nmasA))
flipflag=.true.
endif
v_time=0 !sum of dtnew until newcondTS
t_p=0 !sum of time until termination
do i=1,maxoutputs
t=0
c
c----------estimate the maximum number of particles released
c
max_num=max_numps+NUMPARTS
call refini('prtI',nprtI,1,max_num)
call refinr('prtE',nprtE,1,max_num)
call refinr('prtN',nprtN,1,max_num)
call refinr('prtZ',nprtZ,1,max_num)
call refinr('prtM',nprtM,1,max_num)
call refinr('prtA',nprtA,1,max_num)
call refinr('prtV',nprtV,1,max_num)
c
c----------unsteady case
c
if(.not.hsteadystate.or.(.not.psteadystate)) then
30 if(flipflag) then
call walkem(ia(ncord),ia(nnop),ia(nityp),ia(nebox),
* ia(nBvel),ia(nAvel),ia(nvel),ia(ndudt),ia(nao),nxyz,
* ia(negrp),ia(negps),ia(nnblk),ia(nnelb),ia(nnbb),
* ia(nnee),ia(nBpol),ia(nApol),ia(npoll),ia(npdot),
* ia(nmasp),ia(npwin),ia(nprtE),ia(nprtN),ia(nprtZ),
* ia(nprtI),ia(nprtM),ia(nprtA),ia(nprtV),ia(nopbd),
* no_open,NUMPARTS,t,t_p,v_time,ia(ndifs))
if(v_time.eq.newcondTS) then
v_time=0
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
call getpollconditions(ia(nBpol),ia(nmasp),
* ia(nmasA),NUMPARTS,max_numpsB,keyp)
endif
flipflag=.FALSE.
if(max_numps.lt.max_numpsB) then
max_num=max_num + max_numpsB - max_numps
max_numps=max_numpsB
call refini('prtI',nprtI,1,max_num)
call refinr('prtE',nprtE,1,max_num)
call refinr('prtN',nprtN,1,max_num)
call refinr('prtZ',nprtZ,1,max_num)
call refinr('prtM',nprtM,1,max_num)
call refinr('prtA',nprtA,1,max_num)
call refinr('prtV',nprtV,1,max_num)
endif
endif
else
call walkem(ia(ncord),ia(nnop),ia(nityp),ia(nebox),
* ia(nAvel),ia(nBvel),ia(nvel),ia(ndudt),ia(nao),nxyz,
* ia(negrp),ia(negps),ia(nnblk),ia(nnelb),ia(nnbb),
* ia(nnee),ia(nApol),ia(nBpol),ia(npoll),ia(npdot),
* ia(nmasp),ia(npwin),ia(nprtE),ia(nprtN),ia(nprtZ),
* ia(nprtI),ia(nprtM),ia(nprtA),ia(nprtV),ia(nopbd),
* no_open,NUMPARTS,t,t_p,v_time,ia(ndifs))
if(v_time.eq.newcondTS) then
v_time=0
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
call getpollconditions(ia(nApol),ia(nmasp),
* ia(nmasA),NUMPARTS,max_numpsA,keyp)
endif
flipflag=.TRUE.
if(max_numps.lt.max_numpsA) then
max_num=max_num + max_numpsA - max_numps
max_numps=max_numpsA
call refini('prtI',nprtI,1,max_num)
call refinr('prtE',nprtE,1,max_num)
call refinr('prtN',nprtN,1,max_num)
call refinr('prtZ',nprtZ,1,max_num)
call refinr('prtM',nprtM,1,max_num)
call refinr('prtA',nprtA,1,max_num)
call refinr('prtV',nprtV,1,max_num)
endif
endif
endif
if(t.lt.outputTS) then
if(.not.hsteadystate.and.keyf.eq.0) goto 35
if(.not.psteadystate.and.keyp.eq.0) goto 35
goto 30
endif
c
c----------steady-state case
c
else
call walkem(ia(ncord),ia(nnop),ia(nityp),ia(nebox),
* ia(nAvel),ia(nAvel),ia(nAvel),ia(ndudt),ia(nao),nxyz,
* ia(negrp),ia(negps),ia(nnblk),ia(nnelb),ia(nnbb),
* ia(nnee),ia(nApol),ia(nApol),ia(nApol),ia(npdot),
* ia(nmasp),ia(npwin),ia(nprtE),ia(nprtN),ia(nprtZ),
* ia(nprtI),ia(nprtM),ia(nprtA),ia(nprtV),ia(nopbd),
* no_open,NUMPARTS,t,t_p,v_time,ia(ndifs))
endif
c
c----------output graphics
c
35 call refini('prtI',nprtI,1,NUMPARTS)
call refinr('prtE',nprtE,1,NUMPARTS)
call refinr('prtN',nprtN,1,NUMPARTS)
call refinr('prtZ',nprtZ,1,NUMPARTS)
call refinr('prtM',nprtM,1,NUMPARTS)
call refinr('prtA',nprtA,1,NUMPARTS)
call refinr('prtV',nprtV,1,NUMPARTS)
if(.not.settles.and.t.eq.outputTS) then
do m=1,no_grids
if(no_grids.eq.1) then
cext=' '
else
write(cext,'(i2.2)') m
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),
* access='direct',status='old',form='unformatted',recl=1)
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)
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)
endif
goE=goee(m)
goN=gonn(m)
gaS=gass(m)
gaC=gacc(m)
numX=numxx(m)
numY=numyy(m)
numZ=numzz(m)
dimXY=dimxyy(m)
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)
close(21)
close(22)
close(23)
if(fbindump) close(19)
if(pbindump(m)) close(20)
enddo
endif
if(.not.hsteadystate.and.keyf.eq.0) goto 40
if(.not.psteadystate.and.keyp.eq.0) goto 40
enddo
40 if(settles) write(30,rec=1) no_settles,mass_pp
call bcname
c
c-------output restart file
c
if(restout) then
if(hsteadystate.and.psteadystate) then
call outputRS(NUMPARTS,max_numps,ia(nprtI),ia(nprtE),
* ia(nprtN),ia(nprtZ),ia(nprtM),ia(nprtA),ia(nprtV),ia(nAvel),
* ia(nwsl1),ia(nApol),ia(nmasA))
elseif(flipflag) then
call outputRS(NUMPARTS,max_numps,ia(nprtI),ia(nprtE),
* ia(nprtN),ia(nprtZ),ia(nprtM),ia(nprtA),ia(nprtV),ia(nAvel),
* ia(nwsl1),ia(nApol),ia(nmasA))
else
call outputRS(NUMPARTS,max_numps,ia(nprtI),ia(nprtE),
* ia(nprtN),ia(nprtZ),ia(nprtM),ia(nprtA),ia(nprtV),ia(nBvel),
* ia(nwsl1),ia(nBpol),ia(nmasA))
endif
endif
50 close(2)
close(3)
close(4)
close(9)
if(restart) close(10)
if(restout) close(11)
if(settles) close(30)
stop
end
c----------------------------------------------------------------outputRS
subroutine outputRS(NUMPARTS,max_numps,partI,partE,partN,partZ,
* partM,partA,partV,vel,wsel1,poll,masspA)
c-----------------------------------------------------------------------c
c Purpose: c
c To output neccessary data for restart. c
c-----------------------------------------------------------------------c
include '3dpolls.cb'
include '3dmesh.cb'
integer NUMPARTS,max_numps,i,j
integer partI(1)
real partE(1),partN(1),partZ(1),partM(1),partA(1),partV(1)
real poll(4,1),vel(3,1),wsel1(1),masspA(1)
c
rewind 11
write(11) NUMPARTS,max_numps,(partI(i),i=1,NUMPARTS),(partE(i),
1 i=1,NUMPARTS),(partN(i),i=1,NUMPARTS),(partZ(i),i=1,
2 NUMPARTS),(partM(i),i=1,NUMPARTS),(partA(i),i=1,
3 NUMPARTS),(partV(i),i=1,NUMPARTS),((vel(j,i),j=1,3),
4 wsel1(i),i=1,np),((poll(j,i),j=1,4),i=1,num_pollwins),
5 (masspA(i),i=1,num_pollwins)
return
end
c-----------------------------------------------------------------inputRS
subroutine inputRS(NUMPARTS,max_numps,partI,partE,partN,partZ,
* partM,partA,partV,vel,wsel1,poll,masspA)
c-----------------------------------------------------------------------c
c Purpose: c
c To input neccessary data for restart. c
c-----------------------------------------------------------------------c
include '3dpolls.cb'
include '3dmesh.cb'
integer NUMPARTS,max_numps,i,j
integer partI(1)
real partE(1),partN(1),partZ(1),partM(1),partA(1),partV(1)
real poll(4,1),vel(3,1),wsel1(1),masspA(1)
c
rewind 10
read(10) NUMPARTS,max_numps,(partI(i),i=1,NUMPARTS),(partE(i),
1 i=1,NUMPARTS),(partN(i),i=1,NUMPARTS),(partZ(i),i=1,
2 NUMPARTS),(partM(i),i=1,NUMPARTS),(partA(i),i=1,
3 NUMPARTS),(partV(i),i=1,NUMPARTS),((vel(j,i),j=1,3),
4 wsel1(i),i=1,np),((poll(j,i),j=1,4),i=1,num_pollwins),
5 (masspA(i),i=1,num_pollwins)
return
end

35
source/3duparms.cb
Unescape Escape View File

@ -0,0 +1,35 @@
c USER PARAMETERS COMMON BLOCK
c Used as part of the random walk model
INTEGER newcondTS
INTEGER outputTS
INTEGER dt,dtnew
INTEGER tot_dead_age
INTEGER maxoutputs
INTEGER lnum_ps
INTEGER no_grids
REAL hdiff(1:10), vdiff(1:10)
INTEGER ediff(1:10)
REAL mass_pp,c_init
real xscale,yscale,clevel,dlevel
real t90_s,t90_d,tdpth,tofday,dayofy
character*12 start_stamp
character*5 dimopt
logical hsteadystate,psteadystate,restart,fbindump
logical pbindump(50),settles,openbd
logical p_count,dilution,concentration,dieoff,dieoffc
logical plumeonly,floatableonly
logical tracking,euler
logical compact
COMMON /UPARMS/newcondTS,outputTS,dt,dtnew,tot_dead_age,
1 maxoutputs,hdiff,vdiff,ediff,mass_pp,c_init,lnum_ps,no_grids,
2 start_stamp,hsteadystate,psteadystate,restart,fbindump,
3 pbindump,settles,openbd,p_count,dilution,concentration,
4 plumeonly,floatableonly,tracking,euler,compact,dimopt
common /settle/no_settles
common /dimen2/xscale,yscale
common /conserv/t90_s,t90_d,tdpth,tofday,dayofy,dieoff,dieoffc
common /cdlvl/clevel(5),dlevel(5)

772
source/3dwalking.for
Unescape Escape View File

@ -0,0 +1,772 @@
c---------------------------------------------------------------walkem
subroutine walkem(cord,nop,itype,ebox,vel1,vel2,vel,dudt,ao,
* nxyz,egrp,egps,nblk,nelb,nbb,nee,poll1,poll2,
* poll,PDOT,massp,pollwin,partE,partN,partZ,
* partI,partM,partA,partV,opbd,no_open,NUMPARTS,
* t,t_p,v_time,diffs)
c---------------------------------------------------------------------c
c purpose: c
c To perform random walking. c
c---------------------------------------------------------------------c
include '3duparms.cb'
include '3dpolls.cb'
include '3dgeom.cb'
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
common /shape/shap(20),shpx(20),shpy(20),shpz(20)
common /etype/inode(4)
common /randm/isd, isd1
LOGICAL dead, walkjustone
INTEGER isd, isd1, p, pr, m, w, i, j, key, it, nen, num
INTEGER oldparts, newparts, deadparts, NUMPARTS, nxyz, pil
REAL r, ran3, tZ, tR, q, c, pel, pnl, pzl, pa, pv
REAL pe,pn,pz,xi,eta,zeta,shap
INTEGER t,t_p,v_time,no_open
INTEGER*2 nop(20,1)
INTEGER partI(1),nblk(2,1),nbb(1),nee(1),nelb(1),itype(1)
REAL partE(1),partN(1),partZ(1),partM(1),partA(1),vel(3,1)
REAL cord(3,1),egrp(6,1),egps(6,1),poll(4,1),pollwin(4,1)
REAL ebox(6,1),vel1(3,1),vel2(3,1),poll1(4,1),poll2(4,1)
REAL partV(1),opbd(4,1),dudt(3,1),PDOT(1),massp(1)
REAL TOFDAY1,DAYOFY1
REAL*8 ao(1)
REAL diffs(2,1)
oldparts=NUMPARTS
newparts=0
deadparts=0
num_steps=outputTS/dt
if(.not.hsteadystate.or.(.not.psteadystate)) then
num_steps=min(num_steps,newcondTS/dt)
endif
ibd=0
do i=1,num_steps
c
c----------calculate dieoff
c
if(dieoffc) then
TOFDAY=TOFDAY+dt/3600.0
if(TOFDAY.ge.24.0) then
TOFDAY=TOFDAY-24.0
DAYOFY=DAYOFY+1.0
if(DAYOFY.gt.366.0) DAYOFY=DAYOFY-366.0
endif
DELT=float(dt)
TOFDAY1=TOFDAY
DAYOFY1=DAYOFY
call DIEOFFS(TOFDAY1,DAYOFY1,T90_S)
cdrc
call DIEOFFS(0.0,DAYOFY1,T90_D)
c print *,'T90_S=',T90_S
crdc
c print *,'T90_D=',T90_D
T90_S=3600.0*T90_S
crdc
T90_D=3600.0*T90_D
endif
c
c----------Temporal interpolation of velocities
c
if(.not.hsteadystate.or.(.not.psteadystate)) then
call interp_vels(vel1,vel2,vel,dudt,poll1,
* poll2,poll,v_time)
endif
c
c----------Old particles
c
m=1
do while (m.le.NUMPARTS)
c
c---------------The particle is "active" and in the system.
c
pel=partE(m)
pnl=partN(m)
pzl=partZ(m)
pil=partI(m)
pa=partA(m)
pv=partV(m)
pm=partM(m)
it=itype(pil)
nen=inode(it)
call xn3(it,nen,pel,pnl,pzl)
pe=0.0
pn=0.0
pz=0.0
do j=1,nen
pe=pe+shap(j)*cord(1,nop(j,pil))
pn=pn+shap(j)*cord(2,nop(j,pil))
pz=pz+shap(j)*cord(3,nop(j,pil))
enddo
dead=walkjustone(pel,pnl,pzl,pil,pa,pv,pm,cord,itype,
* ebox,nop,vel,dudt,ao,egrp,egps,nblk,nelb,nbb,
* nee,nxyz,pe,pn,pz,no_open,opbd,dt,diffs,ibd)
partE(m)=pel
partN(m)=pnl
partZ(m)=pzl
partI(m)=pil
partA(m)=pa
partV(m)=pv
partM(m)=pm
IF (dead.and.(.not.tracking)) THEN
c
c-----------------------particle is dead or out of system
c
NUMPARTS=NUMPARTS-1
do j=m,NUMPARTS
partI(j)=partI(j+1)
partE(j)=partE(j+1)
partN(j)=partN(j+1)
partZ(j)=partZ(j+1)
partM(j)=partM(j+1)
partA(j)=partA(j+1)
partV(j)=partV(j+1)
enddo
deadparts=deadparts+1
m=m-1
END IF
m=m+1
enddo
c
c----------New Particles.
c
NUMPARTS0=NUMPARTS
do ns=1,int(dt/dtnew)
c
c------------Temporal interpolation of velocities
c
if(ns.ne.1) then
if(.not.hsteadystate.or.(.not.psteadystate)) then
call interp_vels(vel1,vel2,vel,dudt,poll1,
* poll2,poll,v_time)
endif
endif
ccc num_s=num_steps-ns+1
c
c------------Old particles related to time increment dtnew
c
m=1+NUMPARTS0
do while (m.le.NUMPARTS)
c
c---------------The particle is "active" and in the system.
c
pel=partE(m)
pnl=partN(m)
pzl=partZ(m)
pil=partI(m)
pa=partA(m)
pv=partV(m)
pm=partM(m)
it=itype(pil)
nen=inode(it)
call xn3(it,nen,pel,pnl,pzl)
pe=0.0
pn=0.0
pz=0.0
do j=1,nen
pe=pe+shap(j)*cord(1,nop(j,pil))
pn=pn+shap(j)*cord(2,nop(j,pil))
pz=pz+shap(j)*cord(3,nop(j,pil))
enddo
dead=walkjustone(pel,pnl,pzl,pil,pa,pv,pm,cord,itype,
* ebox,nop,vel,dudt,ao,egrp,egps,nblk,nelb,nbb,
* nee,nxyz,pe,pn,pz,no_open,opbd,dtnew,diffs,ibd)
partE(m)=pel
partN(m)=pnl
partZ(m)=pzl
partI(m)=pil
partA(m)=pa
partV(m)=pv
partM(m)=pm
IF (dead.and.(.not.tracking)) THEN
c
c-----------------------particle is dead or out of system
c
NUMPARTS=NUMPARTS-1
do j=m,NUMPARTS
partI(j)=partI(j+1)
partE(j)=partE(j+1)
partN(j)=partN(j+1)
partZ(j)=partZ(j+1)
partM(j)=partM(j+1)
partA(j)=partA(j+1)
partV(j)=partV(j+1)
enddo
deadparts=deadparts+1
m=m-1
END IF
m=m+1
enddo
do w=1,num_pollwins
c
c-------------Establish the number of particles to be released.
c
q=poll(1,w)
c=poll(2,w)
PDOT(w)=PDOT(w)+(c*q*dtnew)/massp(w)
P=INT(PDOT(w))
PDOT(w)=PDOT(w)-P
c p=100000
c if(i.gt.1.or.ns.gt.1) p=0
do pr=1,P
NUMPARTS=NUMPARTS+1
m=NUMPARTS
r=ran3(isd)
pe=pollwin(1,w)+r*(pollwin(3,w)-pollwin(1,w))
pn=pollwin(2,w)+r*(pollwin(4,w)-pollwin(2,w))
partA(m)=0.
pm=massp(w)
tZ=-poll(3,w)
tR=poll(4,w)
r=2.0*(ran3(isd)-0.5)
pz=tZ+r*tR
pz_keep=pz
c
c---------------Check if the particle is in the water column.
c
IF (pz.gt.0.0) pz=0.0
c i.e. if it is above the water surface then put it on the surface
num=-1
call trinvs3(nop,cord,itype,ebox,egrp,nblk,egps,nelb,
* nbb,nee,nxyz,nen,it,pe,pn,pz,xi,eta,zeta,num,key)
num_keep1=num
if(key.eq.1) then
c i.e. the pe,pn,pz coordinate was inside of element num
partE(m)=xi
partN(m)=eta
partZ(m)=zeta
partI(m)=num
else
c try and see if the plan form (pe,pn) is inside an element
num=-1
call trinvs3(nop,cord,itype,ebox,egrp,nblk,egps,
* nelb,nbb,nee,nxyz,nen,it,pe,pn,0.0,xi,eta,
* zeta,num,key)
num_keep2=num
if(key.eq.1) then
c the (pe, pn) point was in element num, but we knew that (pe,pn,pz) wasn't so the
c old pz value must have been below the bed.
pz=0.0
call xn3(it,nen,xi,eta,zeta)
do j=1,nen
pz=pz+shap(j)*ao(nop(j,num))
enddo
pz=0.99*pz
c set the pz value to be at 99% of the depth and all is well.
else
c the (pe,pn) pair was not found to be in an element.
c This may happen if the (pe, pn) pair falls directly onto a element boundary.
write(*,'(a)') ' *** Window for discharging '//
* 'plumes out of the modelling domain (1)***'
write(*,'(i5,3f15.5)') 0,pe,pn,pz
stop
endif
num=-1
call trinvs3(nop,cord,itype,ebox,egrp,nblk,egps,
* nelb,nbb,nee,nxyz,nen,it,pe,pn,pz,xi,eta,
* zeta,num,key)
if(key.eq.1) then
partE(m)=xi
partN(m)=eta
partZ(m)=zeta
partI(m)=num
else
write(*,'(a)') ' *** Window for discharging '//
* 'plumes out of the modelling domain (2)***'
write(*,'(i5,3f15.5)') 1,pe,pn,pz
write(*,*) 'orig pz, num_keep1, num_keep2',pz_keep
* ,num_keep1, num_keep2
stop
endif
endif
pel=partE(m)
pnl=partN(m)
pzl=partZ(m)
pil=partI(m)
pa=partA(m)
it=itype(pil)
nen=inode(it)
call xn3(it,nen,pel,pnl,pzl)
call rise_fall_vel(pv)
c following line for plume only
if(plumeonly) pv=0.0
c following few lines for floatables only
if(floatableonly) then
if(pv.le.0.0) then
NUMPARTS = NUMPARTS - 1
goto 70
endif
endif
c------------------------------------------
if(settles.and.pv.ge.0.0) then
NUMPARTS = NUMPARTS - 1
goto 70
endif
dead=walkjustone(pel,pnl,pzl,pil,pa,pv,pm,cord,itype,
* ebox,nop,vel,dudt,ao,egrp,egps,nblk,nelb,nbb,
* nee,nxyz,pe,pn,pz,no_open,opbd,dtnew,diffs,ibd)
newparts=newparts+1
partE(m)=pel
partN(m)=pnl
partZ(m)=pzl
partI(m)=pil
partA(m)=pa
partV(m)=pv
partm(m)=pm
IF (dead.and.(.not.tracking)) THEN
NUMPARTS=NUMPARTS-1
deadparts=deadparts+1
END IF
70 enddo
enddo
v_time=v_time+dtnew
enddo
c---------------------
c
c----------output tracking path
c
if(tracking.and.t+dt.eq.outputTS) then
ctmp inquire(35,NEXTREC=ntrack)
ctmp write(35,rec=ntrack) NUMPARTS,t_p
write(35,*) NUMPARTS,t_p
ntrack=ntrack+1
do m=1,NUMPARTS
pel=partE(m)
pnl=partN(m)
pzl=partZ(m)
pil=partI(m)
pa=partA(m)
pv=partV(m)
pm=partM(m)
it=itype(pil)
nen=inode(it)
call xn3(it,nen,pel,pnl,pzl)
pe=0.0
pn=0.0
pz=0.0
do j=1,nen
pe=pe+shap(j)*cord(1,nop(j,pil))
pn=pn+shap(j)*cord(2,nop(j,pil))
pz=pz+shap(j)*cord(3,nop(j,pil))
enddo
call calc_ab(aE,bN,pe,pn)
ctmp write(35,rec=ntrack) aE*resXY,bN*resXY
write(35,*) pe,pn,pz,pm
ntrack=ntrack+1
enddo
endif
t=t+dt
t_p=t_p+dt
enddo
print *,'ibd=',ibd
if(.not.hsteadystate.or.(.not.psteadystate)) then
WRITE(*,200) min(outputTS,newcondTS), oldparts,
* newparts,deadparts
else
WRITE(*,200) outputTS, oldparts,
* newparts,deadparts
endif
200 FORMAT('During the last',I5,' seconds:- [ Old:',I8,
1 ' ; New:',I7,' ; Dead:',I7,' ]')
RETURN
END
c-----------------------------------------------------------walkjustone
FUNCTION walkjustone(pel,pnl,pzl,pil,pa,pv,pm,cord,itype,ebox,
* nop,vel,dudt,ao,egrp,egps,nblk,nelb,nbb,nee,nxyz,pe,
* pn,pz,no_open,opbd,dtt,diffs,ibd)
c---------------------------------------------------------------------c
c purpose: c
c To walk one step. c
c---------------------------------------------------------------------c
include '3duparms.cb'
common /shape/shap(20),shpx(20),shpy(20),shpz(20)
common /etype/inode(4)
common /randm/isd, isd1
LOGICAL walkjustone
REAL cord(3,1),vel(3,1),egrp(6,1),egps(6,1),ebox(6,1)
REAL opbd(4,1),dudt(3,1)
REAL pe,pn,pz,pa,pel,pnl,pzl,xi,eta,zeta
REAL peold,pnold,pzold,peold1,pnold1,pzold1
REAL*8 ao(1)
INTEGER*2 nop(20,1)
INTEGER nblk(2,1),nelb(1),nbb(1),nee(1),itype(1)
INTEGER num,pil,key,nen,nxyz,it,no_open
REAL ran3, normdis
LOGICAL boundaries
INTEGER isd, isd1
REAL len,ang
REAL vE,vN,vZ
REAL pv
real*8 a(3,3),b(3),ai(3,3),du(3,3),a1,a2,a3,dudt0,dvdt0,dwdt0
real*8 detj,ass
integer dtt
real diffs(2,1)
walkjustone=.FALSE.
ccc do i=1,num_steps
pa=pa+dtt
if(pa.gt.tot_dead_age.or.pm/mass_pp.le.0.001) then
c if(pa.gt.tot_dead_age) then
walkjustone=.TRUE.
return
endif
c
c----------gradually die-off
c
if(dieoff.and.pv.eq.0.0) then
c print*,tdpth
if(abs(pz).le.tdpth) then
c print*,'surface dieoff=',10.0**(-dtt/t90_s)
c print*,pm
pm=pm*10.0**(-dtt/t90_s)
c print*,pm
else
c print*,'deep dieoff=',10.0**(-dtt/t90_d)
c print*,pm
pm=pm*10.0**(-dtt/t90_d)
c print*,pm
endif
endif
c
c----------settleable particles on bed -- return
c
if(pv.eq.-9999.0) return
c
c----------interpolate velocity from FE mesh
c
vE=0.0
vN=0.0
vZ=0.0
it=itype(pil)
nen=inode(it)
ccc call xn3(it,nen,pel,pnl,pzl)
do j=1,nen
vE=vE+shap(j)*vel(1,nop(j,pil))
vN=vN+shap(j)*vel(2,nop(j,pil))
vZ=vZ+shap(j)*vel(3,nop(j,pil))
c print *,'pil,vel=',pil,j,vel(3,nop(j,pil))
enddo
c print *,'vE,vN,vZ1=',vE,vN,vZ
c
c----------predict the trajectory of the particle
c
if(.not.euler) then
call xn3x(it,nen,pel,pnl,pzl)
c
c----------calculate the Jacobian matrix
c
do j=1,3
a(j,1)=0.0
a(j,2)=0.0
a(j,3)=0.0
du(1,j)=0.0
du(2,j)=0.0
du(3,j)=0.0
do k=1,nen
nele=nop(k,pil)
a(j,1)=a(j,1)+shpx(k)*cord(j,nele)
a(j,2)=a(j,2)+shpy(k)*cord(j,nele)
a(j,3)=a(j,3)+shpz(k)*cord(j,nele)
du(1,j)=du(1,j)+shpx(k)*vel(j,nele)
du(2,j)=du(2,j)+shpy(k)*vel(j,nele)
du(3,j)=du(3,j)+shpz(k)*vel(j,nele)
enddo
enddo
a1=a(2,2)*a(3,3)-a(2,3)*a(3,2)
a2=a(3,2)*a(1,3)-a(3,3)*a(1,2)
a3=a(1,2)*a(2,3)-a(1,3)*a(2,2)
detj=a(1,1)*a1+a(2,1)*a2+a(3,1)*a3
c if(detj.lt.1.0e-10) then
c print *,'detj,pil,nen,pel,pnl,pzl=',detj,pil,nen,pel,pnl,pzl
c stop
c endif
c
c----------calculate the inverse Jacobi matrix
c
ai(1,1)=(a(2,2)*a(3,3)-a(2,3)*a(3,2))/detj
ai(1,2)=(a(3,1)*a(2,3)-a(2,1)*a(3,3))/detj
ai(1,3)=(a(2,1)*a(3,2)-a(3,1)*a(2,2))/detj
ai(2,1)=(a(3,2)*a(1,3)-a(1,2)*a(3,3))/detj
ai(2,2)=(a(1,1)*a(3,3)-a(1,3)*a(3,1))/detj
ai(2,3)=(a(3,1)*a(1,2)-a(1,1)*a(3,2))/detj
ai(3,1)=(a(1,2)*a(2,3)-a(2,2)*a(1,3))/detj
ai(3,2)=(a(2,1)*a(1,3)-a(1,1)*a(2,3))/detj
ai(3,3)=(a(1,1)*a(2,2)-a(1,2)*a(2,1))/detj
do m=1,3
do j=1,3
a(m,j)=0.0
do k=1,3
a(m,j)=a(m,j)+dtt*ai(j,k)*du(k,m)
enddo
c a(m,j)=dtt*a(m,j)
enddo
c print *,'a=',(a(m,j),j=1,3)
enddo
do k=1,3
a(k,k)=a(k,k)-2.0
enddo
dudt0=0.0
dvdt0=0.0
dwdt0=0.0
do k=1,nen
dudt0=dudt0+shap(k)*dudt(1,nop(k,pil))
dvdt0=dvdt0+shap(k)*dudt(2,nop(k,pil))
dwdt0=dwdt0+shap(k)*dudt(3,nop(k,pil))
enddo
b(1)=-2.0*vE-dudt0*dtt
b(2)=-2.0*vN-dvdt0*dtt
b(3)=-2.0*vZ-dwdt0*dtt
c do i=1,3
c ass=max(dabs(a(i,1)),dabs(a(i,2)),dabs(a(i,3)))
c do j=1,3
c a(i,j)=a(i,j)/ass
c enddo
c b(i)=b(i)/ass
c 25 enddo
a1=a(2,2)*a(3,3)-a(2,3)*a(3,2)
a2=a(3,2)*a(1,3)-a(3,3)*a(1,2)
a3=a(1,2)*a(2,3)-a(1,3)*a(2,2)
detj=a(1,1)*a1+a(2,1)*a2+a(3,1)*a3
c detj=detj*dtt
vE=(a1*b(1)+a2*b(2)+a3*b(3))/detj
a1=a(3,1)*a(2,3)-a(3,3)*a(2,1)
a2=a(1,1)*a(3,3)-a(1,3)*a(3,1)
a3=a(2,1)*a(1,3)-a(2,3)*a(1,1)
vN=(a1*b(1)+a2*b(2)+a3*b(3))/detj
a1=a(2,1)*a(3,2)-a(2,2)*a(3,1)
a2=a(3,1)*a(1,2)-a(3,2)*a(1,1)
a3=a(1,1)*a(2,2)-a(1,2)*a(2,1)
vZ=(a1*b(1)+a2*b(2)+a3*b(3))/detj
c print *,'vE,vN,vZ2=',vE,vN,vZ,dtt
c stop
endif
c-------------------------------
peold=pe
pnold=pn
pzold=pz
len=sqrt(6*diffs(1,pil)*dtt)*ran3(isd)
ang=6.2831853*ran3(isd)
angz=6.2831853*ran3(isd)
c if(vdiff.ne.0.0) then
c angz=6.2831853*ran3(isd)
c else
c angz=0.0
c endif
peold1 = pe + vE*dtt
pnold1 = pn + vN*dtt
pe = peold1 + len*cos(ang)*cos(angz)
pn = pnold1 + len*sin(ang)*cos(angz)
if(pv.eq.0.0) then
pzold1 = pz + vZ*dtt
pz = pzold1 + ran3(isd)*sqrt(6*diffs(2,pil)*dtt)*sin(angz)
else
pzold1 = pz + (vZ+pv)*dtt
pz = pzold1
endif
IF(boundaries(pe,pn,pz,peold,pnold,pzold,no_open,opbd)) THEN
walkjustone=.TRUE.
RETURN
ENDIF
c -->
IF (pz.gt.0.0) pz=0.0
num=pil
call trinvs3(nop,cord,itype,ebox,egrp,nblk,egps,nelb,nbb,
* nee,nxyz,nen,it,pe,pn,pz,xi,eta,zeta,num,key)
if(key.eq.1) then
pel=xi
pnl=eta
pzl=zeta
pil=num
else
c
c-------------neglect the random displacements
c
num=pil
if(pnold1.gt.0.0) pnold1=0.0
call trinvs3(nop,cord,itype,ebox,egrp,nblk,egps,nelb,nbb,
* nee,nxyz,nen,it,peold1,pnold1,pzold1,xi,eta,zeta,num,key)
if(key.eq.1) then
pe=peold1
pn=pnold1
pz=pzold1
pel=xi
pnl=eta
pzl=zeta
pil=num
c
c-------------neutral bouyant and floatable particales
c
else if(pv.ge.0.0) then
pe=peold
pn=pnold
pz=pzold
ibd=ibd+1
c
c-------------settleable particles
c
else
num=pil
call trinvs3(nop,cord,itype,ebox,egrp,nblk,egps,nelb,
* nbb,nee,nxyz,nen,it,pe,pn,0.0,xi,eta,zeta,num,key)
if(key.eq.1) then
if(settles) then
no_settles=no_settles+1
write(30,rec=no_settles) pe,pn
walkjustone=.true.
else
pz=0.0
call xn3(it,nen,xi,eta,zeta)
do j=1,nen
pz=pz+shap(j)*ao(nop(j,num))
enddo
pv=-9999.0
endif
else
pe=peold
pn=pnold
pz=pzold
endif
endif
endif
ccc enddo
RETURN
END
c------------------------------------------------------------boundaries
FUNCTION boundaries(pe,pn,pz,peold,pnold,pzold,no_open,opbd)
c---------------------------------------------------------------------c
c purpose: c
c To check whether the particle within the caculation region. c
c---------------------------------------------------------------------c
include '3duparms.cb'
include '3dgeom.cb'
REAL pe, pn, pz, peold, pnold, pzold
REAL a, b
INTEGER no_open
LOGICAL boundaries
real opbd(4,1),s,r
real*8 deta
c This is very rudimentary. It must be modified later.
c A particle can leave the domain in two ways:-
c i) Leave the (x,y) range.
c ii) Enter a part of the (x,y) range that is "land".
c Any particle leaving by mode (i) are considered to be "gone"
c Particles leaving by (ii) are returned to their old position.
c
c To particle is considered to entered a land area if:-
c a) It is beneath the bathymetry.
c b) It is in a grid cell that doesn't have all four
c corners within the ocean area.
c
c The function returns true only if the particle is "gone".
boundaries=.FALSE.
c CALL calc_ab(a, b, pe, pn)
c Check if particle has left the array region. Case i
c IF ((a.le.0.0).or.(a.ge.float(numX)).or.
c 1 (b.le.0.0).or.(b.ge.float(numY))) THEN
c boundaries=.TRUE.
c RETURN
c ENDIF
c-------check if particle has left from specified open boundaries
if(openbd) then
do i=1,no_open
deta=(pnold-pn)*(opbd(3,i)-opbd(1,i))
* -(peold-pe)*(opbd(4,i)-opbd(2,i))
if(dabs(deta).lt.1.0e-8) goto 10
r=((pnold-opbd(2,i))*(opbd(3,i)-opbd(1,i))
* -(peold-opbd(1,i))*(opbd(4,i)-opbd(2,i)))/deta
if(r.lt.-1.0e-5.or.r.gt.1.0+1.0e-5) goto 10
s=((pnold-pn)*(peold-opbd(1,i))
* -(peold-pe)*(pnold-opbd(2,i)))/deta
if(s.le.1.0+1.0e-5.and.s.ge.-1.0e-5) then
boundaries=.TRUE.
return
endif
10 enddo
endif
c Check for particles below the bathymetry. (Case (iia))
RETURN
END
c-----------------------------------------------------------interp_vels
subroutine interp_vels(vel1,vel2,vel,dudt,poll1,poll2,poll,
* v_time)
c---------------------------------------------------------------------c
c purpose: c
c To perform temporal interpolation of velocities. c
c---------------------------------------------------------------------c
INCLUDE '3duparms.cb'
INCLUDE '3dpolls.cb'
INCLUDE '3dmesh.cb'
dimension vel1(3,1),vel2(3,1),vel(3,1),poll1(4,1),poll2(4,1),
* poll(4,1),dudt(3,1)
REAL r
INTEGER i,j,v_time
c
r=float(v_time)/float(newcondTS)
do i=1,np
do j=1,3
vel(j,i)=vel1(j,i)+r*(vel2(j,i)-vel1(j,i))
enddo
enddo
do i=1,num_pollwins
do j=1,4
poll(j,i)=poll1(j,i)+r*(poll2(j,i)-poll1(j,i))
enddo
enddo
c
c-------calculate dudt
c
if(.not.euler) then
do i=1,np
do j=1,3
dudt(j,i)=(vel2(j,i)-vel1(j,i))/float(newcondTS)
enddo
enddo
endif
return
end

25
source/askstr.for
Unescape Escape View File

@ -0,0 +1,25 @@
SUBROUTINE ASKSTR(PROMPT, STRING)
C
C Prompt standard input for a string, and return that string.
C If just <return> is pressed, the passed string is not changed.
CHARACTER*(*) PROMPT,STRING
C
INTEGER STRLEN
CHARACTER*80 TMPSTR
C
STRLEN = LEN(STRING)
DO WHILE(STRING(STRLEN:STRLEN) .EQ. ' ')
IF (STRLEN .EQ. 1) GOTO 150
STRLEN = STRLEN - 1
END DO
150 CONTINUE
cdrc TYPE 10,PROMPT,STRING(1:STRLEN)
PRINT 10,PROMPT,STRING(1:STRLEN)
10 FORMAT(1X,A,' [',A,']: ',$)
READ(*,'(A)') TMPSTR
IF (TMPSTR .NE. ' ') THEN
STRING = TMPSTR
END IF
C
RETURN
END

245
source/blocks.for
Unescape Escape View File

@ -0,0 +1,245 @@
c----------------------------------------------------------------blocks
subroutine blocks(nx,ny,nz,np,ne,nop,cord,nxyz,egrp,nbeg,nend,
* nelg,nbb,nee,nelb,nblk,egps,itype,ebox,m,num)
c---------------------------------------------------------------------c
c Purpose: c
c To divide FE mesh into blocks. c
c---------------------------------------------------------------------c
common /etype/inode(4)
integer*2 nop
dimension nop(20,1),cord(3,1),egrp(6,1),nbeg(1),nend(1),
* nelg(1),nee(1),nbb(1),nelb(1),nblk(2,1),egps(6,1),
* itype(1),ebox(6,1)
c
xmax=cord(1,1)
xmin=cord(1,1)
ymax=cord(2,1)
ymin=cord(2,1)
zmax=cord(3,1)
zmin=cord(3,1)
disd=0.0
do i=1,ne
if(nop(9,i).eq.0) goto 10
c
c----------store element type
c
if(nop(20,i).ne.0) then
itype(i)=3
else if(nop(15,i).ne.0) then
itype(i)=2
else if(nop(13,i).ne.0) then
itype(i)=4
else
itype(i)=1
endif
c
c----------find B.L. and T.R. conners of each element
c
ebox(1,i)=cord(1,nop(1,i))
ebox(4,i)=cord(1,nop(1,i))
ebox(2,i)=cord(2,nop(1,i))
ebox(5,i)=cord(2,nop(1,i))
ebox(3,i)=cord(3,nop(1,i))
ebox(6,i)=cord(3,nop(1,i))
do k=2,inode(itype(i))
ebox(1,i)=min(ebox(1,i),cord(1,nop(k,i)))
ebox(4,i)=max(ebox(4,i),cord(1,nop(k,i)))
ebox(2,i)=min(ebox(2,i),cord(2,nop(k,i)))
ebox(5,i)=max(ebox(5,i),cord(2,nop(k,i)))
ebox(3,i)=min(ebox(3,i),cord(3,nop(k,i)))
ebox(6,i)=max(ebox(6,i),cord(3,nop(k,i)))
enddo
ebox(1,i)=ebox(1,i)-1.0e-5*abs(ebox(1,i))
ebox(4,i)=ebox(4,i)+1.0e-5*abs(ebox(4,i))
ebox(2,i)=ebox(2,i)-1.0e-5*abs(ebox(2,i))
ebox(5,i)=ebox(5,i)+1.0e-5*abs(ebox(5,i))
ebox(3,i)=ebox(3,i)-1.0e-5*abs(ebox(3,i))
ebox(6,i)=ebox(6,i)+1.0e-5*abs(ebox(6,i))
c
c----------calculate the maximum diagonal distance
c
disd0=sqrt((ebox(4,i)-ebox(1,i))*(ebox(4,i)-ebox(1,i))+
* (ebox(5,i)-ebox(2,i))*(ebox(5,i)-ebox(2,i)))
disd=max(disd,disd0)
c
c-------find B.L. and T.R. conners of the FE mesh
c
xmax=max(xmax,ebox(4,i))
xmin=min(xmin,ebox(1,i))
ymax=max(ymax,ebox(5,i))
ymin=min(ymin,ebox(2,i))
zmax=max(zmax,ebox(6,i))
zmin=min(zmin,ebox(3,i))
10 enddo
m=1
dx=(xmax-xmin)/float(nx)
if(dx.lt.disd) then
nx=(xmax-xmin)/disd
if(nx.lt.1) nx=1
dx=(xmax-xmin)/float(nx)
endif
dy=(ymax-ymin)/float(ny)
if(dy.lt.disd) then
ny=(ymax-ymin)/disd
if(ny.lt.1) ny=1
dy=(ymax-ymin)/float(ny)
endif
dz=(zmax-zmin)/float(nz)
x=xmin
c
c-------divide into blocks
c
do i=1,nx
y=ymin
do j=1,ny
z=zmin
do k=1,nz
egrp(1,m)=x-dx*0.01
egrp(2,m)=y-dy*0.01
egrp(3,m)=z-dz*0.01
egrp(4,m)=x+dx*1.01
egrp(5,m)=y+dy*1.01
egrp(6,m)=z+dz*1.01
z=z+dz
m=m+1
enddo
y=y+dy
enddo
x=x+dx
enddo
nxyz=nx*ny*nz
nxyz0=nxyz
num=1
do i=1,nxyz
nbeg(i)=num
do j=1,ne
if(nop(9,j).ne.0) then
do k=1,inode(itype(j))
if(cord(1,nop(k,j)).gt.egrp(4,i).or.
* cord(1,nop(k,j)).lt.egrp(1,i)) goto 30
if(cord(2,nop(k,j)).gt.egrp(5,i).or.
* cord(2,nop(k,j)).lt.egrp(2,i)) goto 30
if(cord(3,nop(k,j)).gt.egrp(6,i).or.
* cord(3,nop(k,j)).lt.egrp(3,i)) goto 30
nelg(num)=j
num=num+1
goto 40
30 enddo
endif
40 enddo
nend(i)=num-1
enddo
i=1
45 if(nend(i)-nbeg(i).lt.0) then
nxyz=nxyz-1
do j=i,nxyz
nbeg(j)=nbeg(j+1)
nend(j)=nend(j+1)
do k=1,6
egrp(k,j)=egrp(k,j+1)
enddo
enddo
else
i=i+1
endif
if(i.le.nxyz) goto 45
c
c-------divide into sub-blocks
c
m=1
num=1
do n=1,nxyz
num_elem=nend(n)-nbeg(n)+1
if(num_elem.gt.nxyz+nxyz) then
a_multiplier=sqrt(float(num_elem))/float(nx*ny)
nxs=int(nx*a_multiplier)
if(nxs.lt.1) nxs=1
nys=int(ny*a_multiplier)
if(nys.lt.1) nys=1
dx=(egrp(4,n)-egrp(1,n))/float(nxs)
dy=(egrp(5,n)-egrp(2,n))/float(nys)
if(dx.lt.disd) then
nxs=(egrp(4,n)-egrp(1,n))/disd
if(nxs.lt.1) nxs=1
dx=(egrp(4,n)-egrp(1,n))/float(nxs)
endif
if(dy.lt.disd) then
dy=disd
nys=(egrp(5,n)-egrp(2,n))/disd
if(nys.lt.1) nys=1
dy=(egrp(5,n)-egrp(2,n))/float(nys)
endif
nxyzs=nxs*nys
if(nxyzs.le.1) goto 70
dz=(egrp(6,n)-egrp(3,n))
nblk(1,n)=nxyzs
nblk(2,n)=m
x=egrp(1,n)
do i=1,nxs
y=egrp(2,n)
do j=1,nys
z=egrp(3,n)
egps(1,m)=x-0.01*dx
egps(2,m)=y-0.01*dy
egps(3,m)=z-0.01*dz
egps(4,m)=x+dx*1.01
egps(5,m)=y+dy*1.01
egps(6,m)=z+dz*1.01
y=y+dy
m=m+1
enddo
x=x+dx
enddo
do i=nblk(2,n),nblk(2,n)+nblk(1,n)-1
nbb(i)=num
do j=nbeg(n),nend(n)
if(nop(9,nelg(j)).ne.0) then
do k=1,inode(itype(nelg(j)))
if(cord(1,nop(k,nelg(j))).gt.egps(4,i).or.
* cord(1,nop(k,nelg(j))).lt.egps(1,i)) goto 50
if(cord(2,nop(k,nelg(j))).gt.egps(5,i).or.
* cord(2,nop(k,nelg(j))).lt.egps(2,i)) goto 50
if(cord(3,nop(k,nelg(j))).gt.egps(6,i).or.
* cord(3,nop(k,nelg(j))).lt.egps(3,i)) goto 50
nelb(num)=nelg(j)
num=num+1
goto 60
50 enddo
endif
60 enddo
nee(i)=num-1
enddo
i=nblk(2,n)
65 if(nee(i)-nbb(i).lt.0) then
nblk(1,n)=nblk(1,n)-1
m=m-1
do j=i,nblk(2,n)+nblk(1,n)-1
nbb(j)=nbb(j+1)
nee(j)=nee(j+1)
do k=1,6
egps(k,j)=egps(k,j+1)
enddo
enddo
else
i=i+1
endif
if(i.le.nblk(2,n)+nblk(1,n)-1) goto 65
else
70 nblk(1,n)=1
nblk(2,n)=m
nbb(m)=num
do i=nbeg(n),nend(n)
nelb(num)=nelg(i)
num=num+1
enddo
do j=1,6
egps(j,m)=egrp(j,n)
enddo
nee(m)=num-1
m=m+1
endif
enddo
return
end

43
source/dieoffs.for
Unescape Escape View File

@ -0,0 +1,43 @@
SUBROUTINE DIEOFFS(TOFDAY,DAYOFY,DELT,T90)
C Calculates faecal coliform dieoff T90 based on solar radiation during
C the time interval DELT(sec).
C Input Time of day (hrs) TOFDAY, Day of year (1-366) DAYOFY.
C Output T90 (hrs)
REAL TOFDAY,DAYOFY,LAT,LONG,ELEVAT,DELT,DELTH,STANM
REAL TSOLHR,T90
COMMON/INPUTS/LAT,LONG,ELEVAT,STANM
COMMON/MET/WETBLB,ATMPR,DRYBLB,CLOUD,DAT
LAT=-33.8
LONG=151.3
ELEVAT=0.0 ! Elevation (m)
DELTH=DELT/3600.
STANM=150.0 ! Longitude corresponding to standard time
WETBLB=20.0 ! Dry bulb temp
DRYBLB=17.5 ! Wet bulb temp (depends on humidity)
CLOUD=0.0 ! Cloud cover (0-1)
ATMPR=1015. ! Atm pressure
DAT=0.07 ! Dust attenuation
C Dieoff-Solar Radiation Parameters
C Ref. Caldwell-Connell Engineers (1980) Fig. 7.6
A=3.3
B=-0.7
CALL HEATEX(TOFDAY,DAYOFY,DELT,TSOLHR)
C Convert radiation in DELT to hourly solar radiation in MJ/m2
TSOLHR=TSOLHR/1000./DELTH
C PRINT*,TSOLHR
C Calculate the T90
IF (TSOLHR.LT.0.08) THEN
T90=20.
ELSE
T90=A*(TSOLHR**(B))
ENDIF
RETURN
END

875
source/flib.for
Unescape Escape View File

@ -0,0 +1,875 @@
c----------------------------------------------------------------setcst
subroutine setcst
c---------------------------------------------------------------------c
c purpose: c
c To set consats. c
c---------------------------------------------------------------------c
ccc implicit real*8 (a-h,o-z)
common mtot,np,ia(1)
common /dbsys/numa,next,idir,ipp(3)
common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10
* ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19
c
ipp(1)=4
ipp(2)=8
ipp(3)=1
c mtot=100000
numa=0
next=1
c idir=mtot
ntm=6
ntr=5
nin=1
not=2
nsp=3
nfl=4
nt7=7
nt8=8
nt9=9
nt10=10
nt11=11
nt12=12
nt13=13
nt14=14
nt15=15
nt16=16
nt17=17
nt18=18
nt19=19
c
return
end
c-----------------------------------------------------------------upper
subroutine upper
c---------------------------------------------------------------------c
c purpose: c
c This routine converts lower case to upper case c
c---------------------------------------------------------------------c
ccc implicit real*8 (a-h,o-z)
character*1 line,sp,sc
c
common /cline/line(160)
common /iline/ii
c
data sp/' '/,sc/':'/
c
c------check end of line and covert to upper case
c
isp=ichar(sp)
do 80 i=1,ii
if(line(i).eq.sc) goto 90
nn=ichar(line(i))
if(nn.gt.96) line(i)=char(nn-isp)
80 continue
c
c------remove blanks at end of line
c
jj=0
do 50 j=1,ii
if(line(j).ne.sp) goto 70
jj=jj+1
50 continue
70 if(jj.eq.0) goto 85
if(jj.eq.ii) goto 100
kk=jj+1
do 60 j=1,ii-jj
line(j)=line(kk)
kk=kk+1
60 continue
ii=ii-jj-1
85 if(line(ii).ne.sp) goto 95
ii=ii-1
goto 85
90 ii=i
c
95 return
100 ii=1
return
end
c----------------------------------------------------------------bcmane
subroutine bcname
c---------------------------------------------------------------------c
c purpose: c
c To print arry names in 'ia'. c
c---------------------------------------------------------------------c
ccc implicit real*8 (a-h,o-z)
character*1 name(4)
common mtot,np,ia(1)
common /dbsys/numa,next,idir,ipp(3)
common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10
* ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19
c
k=idir
do 10 i=1,numa
do 20 j=1,4
name(j)=char(ia(k+j-1))
20 continue
ccc write(not,1000) i,name,ia(k+4),ia(k+5),ia(k+7),ia(k+8)
write(ntm,1000) i,name,ia(k+4),ia(k+5),ia(k+7),ia(k+8)
1000 format(1x,i7,5x,4a1,',nr=',i7,',nc=',i7,',na=',i12,
* ',size=',i12)
k=k+10
10 continue
return
end
c----------------------------------------------------------------definr
subroutine definr(name,na,nr,nc)
c---------------------------------------------------------------------c
c purpose: c
c To reserve storage for REAL data. c
c---------------------------------------------------------------------c
ccc implicit real*8 (a-h,o-z)
character*1 name(4)
common mtot,np,ia(1)
c
np=1
call defin(name,na,nr,nc)
return
end
c----------------------------------------------------------------defini
subroutine defini(name,na,nr,nc)
c---------------------------------------------------------------------c
c purpose: c
c To reserve storage for INTEGER data. c
c---------------------------------------------------------------------c
ccc implicit real*8 (a-h,o-z)
character*1 name(4)
common mtot,np,ia(1)
c
np=1
call defin(name,na,nr,nc)
return
end
c-----------------------------------------------------------------defin
subroutine defin(name,na,nr,nc)
c---------------------------------------------------------------------c
c purpuse: c
c To define and reserve storage for array. c
c---------------------------------------------------------------------c
ccc implicit real*8 (a-h,o-z)
character*1 name(4)
common mtot,np,ia(1)
common /dbsys/numa,next,idir,ipp(3)
common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10
* ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19
c
c---------------------------------------------------------------------c
c where name = name of array--4 logicals maxmum c
c na = location of array if in blank common c
c nr = number of rows c
c nc = number of columns c
c mtot = end of directory c
c numa = number of arrays in data base c
c next = next avialable storage location c
c idir = start of directory in data base c
c ipp = number of logicals contained in data type c
c lenr = number of logicals in physical record c
c np = type of data c
c = 1 integer data c
c = 2 real data c
c = 3 logical data c
c------directory definition for core or sequential files--------------c
c idir(i,n) = name of array --iname(4,char)
c idir(5,n) = number of rows -- nr c
c idir(6,n) = number of colomns -- nc c
c idir(7,n) = type of data -- np c
c idir(8,n) = incore address -- na c
c = -1 if sequential file on disk c
c = -2 if direct access on disk c
c idir(9,n) = size of array c
c idir(10,n) = 0 if in core storage c
c------directory defintion for direct access files--------------------c
c idir(5,n) = number of integers c
c idir(6,n) = number of real words c
c idir(7,n) = number of logicals c
c idir(8,n) = number of logical records c
c idir(9,n) = logical record number c
c idir(10,n) = lun if on logical unit lun c
c---------------------------------------------------------------------c
c
c------evaluate storage requirements
c
if((np.ne.1).and.(np.ne.3)) np=2
nsize=(nr*nc*ipp(np)-1)/(ipp(1)*2)
nsize=nsize*2+2
na=next
next=next+nsize
c
c------set up new directory
c
numa=numa+1
idir=idir-10
i=idir
c
c------check storage limits
c
if(i.ge.next) goto 100
i=next-i+mtot-1
write(ntm,2000) i,mtot,(name(k),k=1,4)
stop
100 call icon(name,ia(i))
ia(i+4)=nr
ia(i+5)=nc
ia(i+6)=np
ia(i+7)=na
ia(i+8)=nsize
ia(i+9)=0
900 return
c---------------------------------------------------------------------c
2000 format(/,'storage fault in defin',/,'storage required =',i9,
* /,'storage available =',i9,
* /,'name = ',4a)
end
c
c----------------------------------------------------------------locate
subroutine locate(name,na,nr,nc)
c---------------------------------------------------------------------c
c purpose: c
c To locate array NAME in IA. c
c---------------------------------------------------------------------c
ccc implicit real*8 (a-h,o-z)
character*1 name
dimension name(4),iname(4)
common mtot,np,ia(1)
common /dbsys/numa,next,idir,ipp(3)
common /iolist/ntm,ntr,nun,not,nsp,nfl,nt7,nt8,nt9,nt10
* ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19
c
c------locate and return properties on array
c
na=0
call icon(name,iname)
i=ifind(iname,0)
if(i.eq.0) goto 900
c
c------return array properties
c
na=ia(i+7)
nr=ia(i+4)
nc=ia(i+5)
np=ia(i+6)
return
c---------------------------------------------------------------------c
900 write(ntm,1000) name
1000 format(/,' NO SUCH ARRAY ',4a1,' IN IA ',/)
return
end
c
c----------------------------------------------------------------delete
subroutine delete(name)
c---------------------------------------------------------------------c
c purpose: c
c To delete array NAME. c
c---------------------------------------------------------------------c
ccc implicit real*8 (a-h,o-z)
character*1 name
dimension name(4),iname(4)
common mtot,np,ia(1)
common /dbsys/numa,next,idir,ipp(3)
common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10
* ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19
c
c------delete array from storage
c
100 call icon(name,iname)
i=ifind(iname,0)
if(i.eq.0) goto 900
c
c------check on storage location
c
200 nsize=ia(i+8)
c
c------set size of array
c
next=next-nsize
numa=numa-1
na=ia(i+7)
c
c------check if out of core or direct access
c
if(na.gt.0) goto 500
write(ntm,1000) name
goto 800
500 if(na.eq.next) goto 800
c
c------compact storage
c
ii=na+nsize
nnxt=next-1
do 700 j=na,nnxt
ia(j)=ia(ii)
ii=ii+1
700 continue
c
c------compact and update directory
c
800 na=i-idir
idir=idir+10
if(na.eq.0) goto 900
na=na/10
do 860 k=1,na
ii=i+9
do 850 j=1,10
ia(ii)=ia(ii-10)
ii=ii-1
850 continue
if(ia(i+7).le.0) goto 860
if(ia(i+9).eq.0) ia(i+7)=ia(i+7)-nsize
i=i-10
860 continue
c
900 return
c---------------------------------------------------------------------c
1000 format(/,' NAME ',4a1,' IS BEING USED FOR AN OUT OF CORE FILE'/)
end
c
c------------------------------------------------------------------icon
subroutine icon(name,iname)
c---------------------------------------------------------------------c
c purpose: c
c To convert logicals to integer data. c
c---------------------------------------------------------------------c
ccc implicit real*8 (a-h,o-z)
character*1 name(4)
dimension iname(4)
c
do 100 i=1,4
iname(i)=ichar(name(i))
100 continue
c
return
end
c-----------------------------------------------------------------ifind
function ifind(iname,lun)
c---------------------------------------------------------------------c
c purpose: c
c This function is used to find location of array. c
c---------------------------------------------------------------------c
ccc implicit real*8 (a-h,o-z)
dimension iname(4)
common mtot,np,ia(1)
common /dbsys/numa,next,idir,ipp(3)
common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10
* ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19
c
i=idir
do 100 n=1,numa
if(lun.ne.ia(i+9)) goto 300
if(iname(1).ne.ia(i )) goto 300
if(iname(2).ne.ia(i+1)) goto 300
if(iname(3).ne.ia(i+2)) goto 300
if(iname(4).eq.ia(i+3)) goto 200
300 i=i+10
100 continue
i=0
200 ifind=i
c
return
end
c----------------------------------------------------------------refinr
subroutine refinr(name,na,nr,nc)
c---------------------------------------------------------------------c
c purpose: c
c To reserve storage for Real array. c
c---------------------------------------------------------------------c
ccc implicit real*8 (a-h,o-z)
character*1 name(4)
common mtot,np,ia(1)
c
np=1
call refin(name,na,nr,nc)
return
end
c----------------------------------------------------------------refini
subroutine refini(name,na,nr,nc)
c---------------------------------------------------------------------c
c purpose: c
c To reserve storage for Integer array. c
c---------------------------------------------------------------------c
ccc implicit real*8 (a-h,o-z)
character*1 name(4)
common mtot,np,ia(1)
c
np=1
call refin(name,na,nr,nc)
return
end
c-----------------------------------------------------------------refin
subroutine refin(name,na,nr,nc)
c---------------------------------------------------------------------c
c purpose: c
c To refinde and reserve new storage for array. c
c---------------------------------------------------------------------c
ccc implicit real*8 (a-h,o-z)
character*1 name(4)
dimension iname(4)
common mtot,np,ia(1)
common /dbsys/numa,next,idir,ipp(3)
common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10
* ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19
c
c------locate origin array
c
call icon(name,iname)
i=ifind(iname,0)
nr1=ia(i+4)
nc1=ia(i+5)
na=ia(i+7)
nsize1=ia(i+8)
nsize=(nr*nc*ipp(np)-1)/(ipp(1)*2)
nsize=nsize*2+2
nszicr=nsize-nsize1
nrincr=nr-nr1
ncincr=nc-nc1
if(nszicr.eq.0) goto 450
c
c------set size of array
c
nexto=next-1
next=next+nszicr
if(i.gt.next) goto 450
i=next-i+mtot-1
print *,'name,nr,nc=',name,nr,nc
write(ntm,1000) i,mtot,(name(k),k=1,4)
stop
450 ia(i+4)=nr
ia(i+5)=nc
ia(i+8)=nsize
if(nszicr.eq.0) goto 900
c
c------compact storage
c
nextn=next-1
nnew=na+nsize
nold=na+nsize1
if(np.eq.3) goto 500
if(nrincr.eq.0) goto 500
nrincr=nrincr*np
nexto=nexto-nrincr
nnew=na
nold=na
do 10 k=1,nc1
nnew=nnew+nr*np
nold=nnew-nrincr
nexto=nexto+nrincr
if(i.gt.nexto) goto 480
i=nexto-i+mtot
print *,'name,nr,nc=',name,nr,nc
write(ntm,1000) i,mtot,(name(m),m=1,4)
stop
480 if(nrincr.gt.0) goto 510
ij=nnew
do 20 j=nold,nexto
ia(ij)=ia(j)
ij=ij+1
20 continue
goto 10
510 ij1=nexto+nrincr
do 50 j=nold,nexto
ii1=nexto-j+nold
ia(ij1)=ia(ii1)
ij1=ij1-1
50 continue
10 continue
nnew=na+nsize
nexto=nexto+nrincr
nold=nnew-(nextn-nexto)
500 if(i.eq.idir) goto 900
if(nextn.eq.nexto) goto 600
if(nnew.gt.nold) goto 520
ij=nold
do 30 j=nnew,nextn
ia(j)=ia(ij)
ij=ij+1
30 continue
goto 600
520 ij1=nexto
do 60 j=nnew,nextn
ii1=nextn-j+nnew
ia(ii1)=ia(ij1)
ij1=ij1-1
60 continue
c
c------compact and update directory
c
600 ij=(i-idir)/10
if(ij.eq.0) goto 900
i=i-10
do 40 j=1,ij
if(ia(i+7).le.0) goto 40
if(ia(i+9).eq.0) ia(i+7)=ia(i+7)+nszicr
i=i-10
40 continue
c
900 return
c---------------------------------------------------------------------c
1000 format(/,' STORAGE FAULT IN REFIN',/,' STORAGE REQUIRED =',I9,
* /,' STORAGE AVIABLE =',I9,
* /,' NAME = ',4a)
end
c------------------------------------------------------------------find
subroutine find(sep,nostr,nn)
c---------------------------------------------------------------------c
c input parameters: c
c sep--separator c
c nn--flag whether there is separator or not c
c---------------------------------------------------------------------c
c implicit real*8 (a-h,o-z)
character line*160,sep*80
common /cline/line
common /iline/ii
common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10
* ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19
c
c------find separator in input file
c
nn=0
50 read(nin,1000,err=900,end=900) line(1:80)
ii=80
c
c------convert to upper case
c
call upper
if(sep(1:nostr).ne.line(1:nostr)) goto 50
nn=1
900 return
c---------------------------------------------------------------------c
1000 format(a80)
end
c------------------------------------------------------------------free
subroutine free
c---------------------------------------------------------------------c
c purpose: c
c To read a line of free field data. c
c---------------------------------------------------------------------c
c implicit real*8 (a-h,o-z)
character*1 line,sp,bs,sc,c
common /cline/line(160)
common /iline/ ii
common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10
* ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19
data sp/' '/,c/'%'/,sc/':'/,bs/'\\'/
c
c------read line of free field data
c
do 40 i=1,160
line(i)=sp
40 continue
50 i=1
ii=80
60 read(nin,1000,err=100) (line(k),k=i,ii)
c
c------check for additional line
c
70 do 80 k=i,ii
if(line(k).ne.bs) goto 80
i=k
ii=k+79
if(ii.gt.160) ii=160
goto 60
80 continue
call upper
c
c------check for comment
c
if(line(1).ne.c) goto 900
write(ntm,2000) (line(i),i=1,ii)
goto 50
c
c------error in data
c
100 write(ntm,2100)
call freept
c
900 return
c---------------------------------------------------------------------c
1000 format(80a1)
2000 format(1x,80a1)
2100 format(/'* Error in reading input line *'/)
end
c----------------------------------------------------------------preept
subroutine freept
c---------------------------------------------------------------------c
c purpose: c
c To write record to file NOT and NTM. c
c---------------------------------------------------------------------c
c implicit real*8 (a-h,o-z)
character*1 line
c
common /iline/ii
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
c------echo of free field information
c
write(ntm,1000) (line(i),i=1,ii)
return
1000 format(1x,80a1)
end
c-----------------------------------------------------------------freeh
subroutine freeh(ic,idata,nc,num)
c---------------------------------------------------------------------c
c purpose: c
c To extract hollerith data. c
c---------------------------------------------------------------------c
c num -- negative read strings separated by blank space or comma c
c num -- positive read nc number of characters c
c---------------------------------------------------------------------c
c implicit real*8 (a-h,o-z)
character*1 line,ic,blk,comma,io,eqs,colin,idata
dimension idata(nc,num)
common /cline/line(160)
common /iline/ii
data blk/' '/,comma/','/,io/'0'/,eqs/'='/,colin/':'/
c
c------find hollerith string
c
num0=abs(num)
90 i=0
do 210 j=1,num0
do 210 n=1,nc
idata(n,j)=blk
210 continue
if(ic.eq.blk) goto 200
do 100 i=2,ii
if((line(i-1).eq.ic).and.(line(i).eq.eqs)) goto 200
100 continue
return
c
c------extract hollerith string
c
200 do 300 j=1,num0
260 i=i+1
if(i.gt.ii) goto 400
do 290 n=1,nc
if(line(i).eq.eqs) goto 400
if(line(i).eq.colin) goto 300
if(num.gt.0) goto 280
if(line(i).eq.blk) goto 300
if(line(i).eq.comma) goto 300
280 idata(n,j)=line(i)
if(n.eq.nc) goto 290
i=i+1
290 continue
300 continue
c
400 return
end
c-----------------------------------------------------------------freei
subroutine freei(ic,idata,num)
c---------------------------------------------------------------------c
c purpose: c
c To extract integer data. c
c---------------------------------------------------------------------c
c input parameter: c
c ic--flag of data ic=' ' unidentifier c
c ic=else identifier c
c num--number of data c
c---------------------------------------------------------------------c
c output parameter: c
c idata--integer array c
c---------------------------------------------------------------------c
c implicit real*8 (a-h,o-z)
character*1 line,ic,blk,comma,io,eqs,neg,colin,lne,pls
dimension idata(num)
common /cline/line(160)
common /iline/ii
common /ilist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10
data blk/' '/,comma/','/,io/'0'/,eqs/'='/,neg/'-'/,colin/':'/
data pls/'+'/
c
c------find integer string
c
90 i=0
do 210 j=1,num
idata(j)=0
210 continue
if(ic.eq.blk) goto 200
do 100 i=1,ii
if(line(i).eq.ic.and.line(i+1).eq.eqs) goto 200
100 continue
return
c
c------zero integer string
c
200 if(line(i+1).eq.eqs) i=i+1
do 250 j=1,num
isign=1
c
c------skip blanks between integers
c
215 if(line(i+1).ne.blk) goto 220
i=i+1
if(i.gt.ii) goto 900
goto 215
220 i=i+1
if(i.gt.ii) goto 230
c
c------check for sign
c
if(line(i).ne.pls) goto 235
i=i+1
if(i.gt.ii) goto 230
235 lne=line(i)
if(lne.ne.neg) goto 225
isign=-1
goto 220
c
c------extract integer
c
225 if(lne.eq.blk) goto 230
if(lne.eq.comma) goto 230
if(lne.eq.colin) goto 230
nn=ichar(lne)-ichar(io)
if((nn.lt.0).or.(nn.gt.9)) goto 900
idata(j)=10*idata(j)+nn
goto 220
c
c------set sign
c
230 idata(j)=idata(j)*isign
250 continue
900 return
end
c-----------------------------------------------------------------freer
subroutine freer(ic,data,num)
c---------------------------------------------------------------------c
c purpose: c
c To extract real data. c
c---------------------------------------------------------------------c
c implicit real*8 (a-h,o-z)
character*1 line,blk,ic,mul,div,add,sub,eqs,e
dimension data(num)
common /cline/line(160)
common /iline/ii
c
data blk/' '/,mul/'*'/,div/'\/'/,add/'+'/,sub/'-'/,e/'E'/,
* eqs/'='/
c
c------find real string
c
90 i=0
do 260 j=1,num
data(j)=0.0
260 continue
if(ic.eq.blk) goto 250
do 100 i=1,ii
if((line(i).eq.ic).and.(line(i+1).eq.eqs)) goto 250
100 continue
return
c
c------extract real data
c
250 do 300 j=1,num
jj=0
270 if(i.gt.ii) goto 300
call rdata(i,xx,nn)
if(jj.ne.0) goto 275
data(j)=xx
goto 290
c
c------arithmetric statement
c
275 if(jj.eq.1) data(j)=data(j)*xx
if(jj.eq.2) data(j)=data(j)/xx
if(jj.eq.3) data(j)=data(j)+xx
if(jj.eq.4) data(j)=data(j)-xx
if(jj.ne.5) goto 290
c
c------exponential data
c
jj=abs(xx)
if(jj.eq.0) goto 290
do 280 k=1,jj
if(xx.lt.0.0) data(j)=data(j)/10.0
if(xx.gt.0.0) data(j)=data(j)*10.0
280 continue
c
c------set type of statement
c
290 jj=0
if(line(i).eq.mul) jj=1
if(line(i).eq.div) jj=2
if(line(i).eq.add) jj=3
if(line(i).eq.sub) jj=4
if(line(i).eq.e) jj=5
if(line(i+1).eq.eqs) jj=0
if(jj.ne.0) goto 270
if(nn.gt.9) return
300 continue
c---------------------------------------------------------------------c
return
end
c-----------------------------------------------------------------rdata
subroutine rdata(i,xx,nn)
c---------------------------------------------------------------------c
c purpose: c
c To convert string to real floating point number. c
c---------------------------------------------------------------------c
ccc implicit real*8 (a-h,o-z)
character*1 line,blk,comma,io,dot,neg,eqs,pls
common /cline/line(160)
common /iline/ii
data blk/' '/,comma/','/,io/'0'/,dot/'.'/,neg/'-'/,eqs/'='/
data pls/'+'/
c
c------converts string to real floating point number
c
if(line(i+1).eq.eqs) i=i+1
y=0
is=1
xx=0.0
250 if(line(i+1).ne.blk) goto 255
i=i+1
if(i.gt.ii) goto 300
goto 250
255 if(line(i+1).ne.neg) goto 260
is=-1
i=i+1
if(i.gt.ii) goto 300
260 if(line(i+1).ne.pls) goto 265
i=i+1
if(i.gt.ii) goto 300
265 if(line(i+1).ne.blk) goto 270
i=i+1
if(i.gt.ii) goto 300
goto 265
270 i=i+1
if(i.gt.ii) goto 300
if(line(i).eq.blk.and.line(i+1).eq.blk) goto 270
nn=ichar(line(i))-ichar(io)
xn=isign(nn,is)
if(line(i).ne.dot) goto 275
y=1.0
goto 270
275 if(line(i).eq.blk) goto 300
if(line(i).eq.comma) goto 300
if(nn.lt.0.or.nn.gt.9) goto 300
if(y.eq.0) goto 280
y=y/10.0
xn=xn*y
xx=xx+xn
goto 270
280 xx=10.0*xx+xn
goto 270
300 return
end
c
c---------------------------------------------------------------clearr
subroutine clearr(adata,nr,nc)
c--------------------------------------------------------------------c
c purpose: c
c This routine is used to initialize real array ADATA. c
c--------------------------------------------------------------------c
ccc implicit real*8 (a-h,o-z)
dimension adata(nr,nc)
c
do 10 i=1,nr
do 20 j=1,nc
adata(i,j)=0.0
20 continue
10 continue
return
end

209
source/gdate.for
Unescape Escape View File

@ -0,0 +1,209 @@
C
C ANALYST : D.GORHAM
C INDEX NO. : UT 95.00
C TITLE : EXTENDED GREGORIAN DATE SUBROUTINE FROM YEAR 1/1/1
C DATE : 08/03/82
C REV.01 : 19-SEP-83 BY DJG: MODIFIED FOR VAX
C
C THIS SUBROUTINE IS ADAPTED BY DJG FROM THAT IN HP COMMUNICATOR
C 1981 VOL. 5 ISSUE 3 WITH INTEGER*4 TYPE BEING CHANGED TO REAL
C PLUS ASSOCIATED CHANGES.
C
C PARAMETER LIST :
C OPT = 1 TO CONVERT CALENDAR DATE TO EXTENDED GREGORIAN DATE (EGD)
C = 2 TO CONVERT EXTENDED GREGORIAN DATE (EGD) TO CALENDAR DATE
C MON = MONTH (1-12)
C DAY = DAY OF MONTH
C YEAR= CALENDAR YEAR EG. 1982
C DWK = DAY OF WEEK (1=SUNDAY, 7=SATURDAY)
C DYR = DAY OF YEAR
C DMON= DAYS IN MONTH
C EGD = EXTENDED GREGORIAN DATE SINCE 1/1/1
C IER = 0 FOR NO ERRORS
C 1 IF ERRORS DETECTED
C LU = 0 FOR NO ERROR MESSAGE PRINTOUT
C = N FOR DISPLAY OF ERROR MESSAGE ON LOGICAL UNIT 'N'
C
C TITLE TRUE GEORGIAN DATE ROUTINE ISS. 2 801126 R.A.G. (MWR)
C
SUBROUTINE GDATE (OPT,MON,DAY,YEAR,DWK,DYR,DMON,EGD,IER,LU)
C
C GREGORIAN DATE ROUTINE 801126
C
IMPLICIT INTEGER*4 (A-Z)
C
C
C --------------------------------------------------------------
C
C REVISION LIST
C
C --DATE-- ---BY--- --DESCRIPTION--
C
C 19/03/80 R.A.G. -ORIGINAL ISSUE
C 26/11/80 R.A.G. -CONVERTED REAL NUMBERS TO DOUBLE INTEGER AND
C ADDED DAY OF THE WEEK, DAY OF THE YEAR (JULIAN
C DATE), AND NUMBER OF DAYS IN THE MONTH.
C
C ---------------------------------------------------------------
C
C
INTEGER*4 M(12)
C
LOGICAL LYEAR
C
DATA M/31,28,31,30,31,30,31,31,30,31,30,31/
C
C INITIALIZE SOME VARIABLES
C
M(2)=28
DYR=0
IER=0
C
C CHECK FOR OPTION
C 1=CALLER SUPPLIES DAY MONTH YEAR
C >1=CALLER SUPPLIES EXTENDED GREGORIAN DATE
C
IF (OPT.GT.1) GO TO 120
C
C THIS SECTION CONVERTS TO A GREGORIAN DATE
C
C TEST ARGUMENTS FOR VALIDITY
C
IF (MON.LT.1.OR.MON.GT.12.OR.DAY.LT.1.OR.DAY.GT.31.
1 OR.YEAR.LT.1) THEN
cdrc TYPE 5000
PRINT 5000
5000 FORMAT(' **GDATE ERROR')
cdrc TYPE '(A,I6)',' MONTH =',MON
cdrc TYPE '(A,I6)',' DAY =',DAY
cdrc TYPE '(A,I6)',' YEAR =',YEAR
PRINT '(A,I6)',' MONTH =',MON
PRINT '(A,I6)',' DAY =',DAY
PRINT '(A,I6)',' YEAR =',YEAR
IER=1
GO TO 180
END IF
C
IF (LYEAR(YEAR)) M(2)=29
C
IF (DAY.GT.M(MON)) THEN
cdrc TYPE 5000
cdrc TYPE '(A,I6)',' DAY =',DAY
PRINT 5000
PRINT '(A,I6)',' DAY =',DAY
IER=2
GO TO 180
END IF
C
C CALCULATE GREG. DATE TO 1ST OF REQUESTED YEAR
C
Y=YEAR-1
EGD=GDATS(Y)
C
C CALCULATE TO CURRENT GREGORIAN DATE
C
J=MON-1
IF (J.EQ.0) GO TO 110
DO 100 I=1,J
DYR=DYR+M(I)
100 CONTINUE
EGD=EGD+DYR
110 EGD=EGD+DAY
DYR=DYR+DAY
GO TO 170
C
C THIS SECTION CONVERTS FROM A GREGORIAN DATE
C
120 IF (EGD.LT.1) THEN
IER=3
cdrc TYPE 5000
cdrc TYPE '(A,I6)',' EGD =',EGD
PRINT 5000
PRINT '(A,I6)',' EGD =',EGD
GO TO 180
END IF
C
C CALCULATE CURRENT DATE (DD/MM/YYYY)
C
YEAR=(EGD/366)-1
130 YEAR=YEAR+1
EEGD=GDATS(YEAR)
IF (EGD-EEGD-368) 140,140,130
140 YEAR=YEAR+1
DYR=EGD-EEGD
C
IF (LYEAR(YEAR)) M(2)=29
C
DO 150 MON=1,12
EEGD=EEGD+M(MON)
IF (EGD.LE.EEGD) GO TO 160
150 CONTINUE
M(2)=28
GO TO 140
C
C CALCULATE THE REMAINING ARGUMENTS
C
160 DAY=EGD+M(MON)-EEGD
170 DMON=M(MON)
DWK=MOD(EGD,7)+1
RETURN
C
180 RETURN
END
C
C
C TITLE FUNCTION LYEAR(LEAP YEAR) ISS. 1 801126 R.A.G. (MWR)
FUNCTION LYEAR(YEAR)
C
IMPLICIT INTEGER*4 (A-Z)
C
C THIS FUNCTION WILL TEST A GIVEN YEAR AND RETURN A TRUE/FALSE
C INDICATION
C
C --------------------------------------------------------------
C
C REVISION LIST
C
C --DATE-- ---BY--- --DESCRIPTION--
C
C 26/11/80 R.A.G. -ORIGINAL ISSUE
C
C --------------------------------------------------------------
C
LYEAR=0
IF (MOD(YEAR,4).EQ.0.AND.MOD(YEAR,100).NE.0.OR.
& MOD(YEAR,400).EQ.0) LYEAR= -1
RETURN
END
C TITLE FUNCTION GDATS(GDATE) ISS. 1 801126 R.A.G. (MWR)
FUNCTION GDATS(YEAR)
C
C
IMPLICIT INTEGER*4 (A-Z)
C
C
C THIS FUNCTION IS PART OF THE GDATE SUBROUTINE, AND RETURNS A
C GREGORIAN DATE TO THE 1ST OF THE YEAR BASED ON THAT YEAR
C
C -----------------------------------------------------------------
C
C REVISION LIST
C
C --DATE-- ---BY--- --DESCRIPTION
C
C 26/11/80 R.A.G. -ORIGINAL ISSUE
C
C --------------------------------------------------------------------
C
C **NOTE**
C
C GDATS MUST BE DECLARED AS AN INTEGER*4 FUNCTION!!!!!!!
C
CON=365
GDATS=CON*YEAR
GDATS=GDATS+(24*(YEAR/100))
GDATS=GDATS+(YEAR/400)
GDATS=GDATS+(MOD(YEAR,100)/4)
RETURN
END

364
source/heatex.for
Unescape Escape View File

@ -0,0 +1,364 @@
SUBROUTINE HEATEX(TOFDAY,DAYOFY,DELT,TSOLHR)
C
C HEATEX COMPUTES THE NET AMOUNT OF HEAT
C RADIATION FLUX BEING TRANSFERRED ACROSS
C THE AIR-WATER INTERFACE BASED ON AN
C ENERGY BUDGET WHICH CONSIDERS SOLAR
C RADIATION, ATMOSPHERIC RADIATION, BACK
C RADIATION, CONDUCTION, AND EVAPORATION.
C SAVE
INTEGER ITIME,INWDAY,METRIC
REAL CON1,CON2,CON3,CON4,CON5,CON6,DELTSL
REAL SOLCON,TSOLHR,HA,CS,SOLAR,ALPHT
REAL ELEVAT,LAT,LONG
REAL DELTH,DELT,TOFDAY,DAYOFY,TLEFT
REAL DECLIN,REARTH,EQTIME,DECLON
REAL PI,STR,STB,STE,STS,ST1,ST2
COMMON/INPUTS/LAT,LONG,ELEVAT,STANM
COMMON/MET/WETBLB,ATMPR,DRYBLB,CLOUD,DAT
C For metric units, radiation is in kJ/m2*hr
METRIC=1
LOUT=6
CC
CCC
CCC NCASI Commentary, HEATEX Section A. (QUAL2 Step 1-0)
CCC A. Compute and/or define required constants.
CCC
C
CCC NCASI Commentary, HEATEX Section F.
CCC F. Test for beginning of a new day.
CCC
CJFD cheap fix to prevent incrementing time for each material type ...
CJFD jump passed time calcs, this only works because local values are saved
C
C DATA ITIME/0/
C
CIPK MAY94 CHANGED ORDER
C
ITIME=0
IF(ITIME .EQ. 0) THEN
PI=4.0*atan(1.)
c PI=3.141628
CON1=2.0*PI/365.0
CON2=PI/180.0*LAT
CON3=180.0/PI
CON4=23.45*PI/180.0
CON5=PI/12.0
CON6=12.0/PI
DELTSL=(LONG-STANM)/15.0
CIPK MAR94
IF(METRIC .EQ. 0) THEN
SOLCON=438.0
ELEXP=EXP(-ELEVAT/2532.0)
ELSE
SOLCON=4974.4
c SOLCON=4870.8 !Original value incorrect DRC
c ELEXP=EXP(-ELEV*3.2808/2532.0)
ELEXP=EXP(-ELEVAT/771.76)
ENDIF
CIPK MAR94 END CHANGES
ENDIF
TLEFT=0.
TSOLHR=0.
55 CONTINUE
ITIME=ITIME+1
IF(ITIME .EQ. 1) THEN
INWDAY=1
ELSE
INWDAY=0
ENDIF
C INWDAY =1 signals a new day or start of simulation
C
C Get Delta Time in hours
C
DELTH = DELT/3600.0
TOFDAY = TOFDAY+DELTH
C
IF(DELTH .EQ. 0.) THEN
DELTH=1.0
TOFDAY=0.00
TLEFT=24.
ENDIF
C
C write(*,*) TOFDAY, dayofy
c WRITE(IOT,*) 'TOFDAY DAYOFY',TOFDAY,DAYOFY
56 CONTINUE
IF(DELT .EQ. 0.) THEN
TLEFT=TLEFT-DELTH*1.0001
TOFDAY=TOFDAY+DELTH
GO TO 85
ENDIF
IF(TLEFT .GT. 0.0) THEN
TOFDAY=TLEFT
DAYOFY=DAYOFY + 1.0
INWDAY=1
IF (TOFDAY .GT. 24.001) THEN
TLEFT=TOFDAY-24.00
TOFDAY=24.00
DELTH=24.00
ELSE
DELTH=TLEFT
TLEFT=0.0
ENDIF
GO TO 85
ENDIF
IF (TOFDAY .GT. 24.001) THEN
TLEFT=TOFDAY-24.00
TOFDAY=24.00
DELTH=DELTH-TLEFT
INWDAY=0
ENDIF
85 CONTINUE
CDRC write(LOUT,*)'at 85 tleft,tofday,delth,inwday',tleft,tofday,delth
CDRC + ,inwday
cDRC CALL GETMET(TOFDAY)
IF(INWDAY .EQ. 1) THEN
CCC
CCC NCASI Commentary, HEATEX Section B. (QUAL2 Step 2-0)
CCC B. Begin computations for calculating the
CCC net solar radiation term.
CCC
CCC B.1 Test for beginning of a new day.
CCC
CCC
CCC B.1a Calculate seasonal and daily position
CCC of the sun relative to the location
CCC of the basin on the earth's
CCC surface. (QUAL2 Step 2-1)
CCC
REARTH=1.0+0.017*COS(CON1*(186.0-DAYOFY))
DECLIN=CON4*COS(CON1*(172.0-DAYOFY))
RR=REARTH**2
EQTIME=0.000121-0.12319*SIN(CON1*(DAYOFY-1.0)-0.07014)
* -0.16549*SIN(2.0*CON1*(DAYOFY-1.0)+0.3088)
DECLON=ABS(DECLIN)
CC
CC Replace TAN function with SIN/COS.
CC
TANA = SIN(CON2)/COS(CON2)
TANB = SIN(DECLON)/COS(DECLON)
ACS = TANA*TANB
CC
IF (ACS .NE. 0.0) THEN
XX=SQRT(1.0-ACS*ACS)
XX=XX/ACS
cipk oct94 ACS=ATAN(XX)
if(xx .gt. 0.) then
ACS=abs(ATAN(XX))
IF (DECLIN.GT.0.0) ACS=PI-ACS
else
acs=abs(atan(xx))
if (declin.lt.0.0) acs=pi-acs
endif
ELSE
ACS=PI/2.0
ENDIF
CCC
CCC B.1a Calculate the standard time of
CCC sunrise (STR) and sunset (STS).
CCC (QUAL2 Step 2-2)
CCC
C
STR=12.0-CON6*ACS+DELTSL
STS=24.0-STR+2.0*DELTSL
DAYLEN=STS-STR
STB=TOFDAY-DELTH
STE=STB+DELTH
cDRC WRITE(LOUT,*) 'AFTER 2-2 STR,STS,STB,STE',STR,STS,STB,STE
CCC
CCC B.2 Increment the variables that define the
CCC time of the beginning(STB) and the
CCC end (STE) of the time interval.
CCC
ELSE
STB=STB+DELTH
STE=STB+DELTH
ENDIF
CCC
CCC B.3 Test if time to read in local
CCC climatological data. (QUAL2 Step 2-3)
CCC
CJFD IF (TRLCD.NE.1.0) GO TO 82
CCC
CCC B.7 Compute vapor pressures (VPWB and
CCC VPAIR), dew point (DEWPT), AND
CCC dampening effect of clouds (CNS
CCC and CNL). (QUAL2 Step 2-4)
CCC
CIPK MAR94
IF(METRIC .EQ. 0) THEN
VPWB=0.1001*EXP(0.03*WETBLB)-0.0837
VPAIR=VPWB-0.000367*ATMPR*(DRYBLB-WETBLB)
* *(1.0+(WETBLB-32.0)/1571.0)
DEWPT=ALOG((VPAIR+0.0837)/0.1001)/0.03
ELSE
c VPWB=(0.1001*EXP(0.03*(WETBLB(NN)*1.8+32.))-0.0837)/29.53*1000.
c VPAIR=VPWB-0.000367*ATMPR(NN)*(DRYBLB(NN)-WETBLB(NN))*1.8
c * *(1.0+(WETBLB(NN)*1.8)/1571.0)
c DEWPT=((ALOG((VPAIR/1000.*29.53+0.0837)/0.1001)/0.03)-32.0)/1.8
c IF(NN .EQ. 1) WRITE(75,*) 'VPWB',VPWB,VPAIR,DEWPT
vpwb=8.8534*exp(0.054*wetblb)-2.8345
VPAIR=VPWB-0.0006606*ATMPR*(DRYBLB-WETBLB)
* *(1.0+WETBLB/872.78)
c vpwb in millibars
c vpair in millibars
DEWPT=ALOG((VPAIR+2.8345)/8.8534)/0.054
c dewpt in deg C
c IF(NN .EQ. 1) WRITE(75,*) 'VPWB',VPWB,VPAIR,DEWPT
ENDIF
CIPK MAR94 END CHANGES
CS=1.0-0.65*CLOUD**2
IF (CLOUD.GT.0.9) CS=0.50
CNL=CLOUD*10.0+1.0
NL=CNL
82 CONTINUE
IF (STS.LE.STB.OR.STR.GE.STE) GO TO 35
C IF(STR.GT.STB.AND.STR.LT.STE) GO TO 41
C IF (STS.LT.STE.AND.STS.GT.STB) GO TO 42
ST1=STB
ST2=STE
IF(STB .LT. STR) ST1=STR
IF(STE .GT. STS) ST2=STS
CCC
CCC NCASI Commentary, HEATEX Section C. (QUAL2 Step 2-5)
CCC C. Continue with calculations for solar
CCC radiation.
CCC
CCC C.1 Calculate hour angles (TB and TE).
CCC
C TB=STB-12.0-DELTSL+EQTIME
C TE=STE-12.0-DELTSL+EQTIME
C WRITE(LOUT,*) '40 ,DELTSL,EQTIME,STB,STE,TB,TE',
C +DELTSL,EQTIME,STB,STE,TB,TE
C GO TO 43
C 41 TB=STR-12.0-DELTSL+EQTIME
C TE=STE-12.0-DELTSL+EQTIME
C WRITE(LOUT,*) '41 ,DELTSL,EQTIME,STR,STE,TB,TE',
C +DELTSL,EQTIME,STR,STE,TB,TE
C GO TO 43
C 42 TB=STB-12.0-DELTSL+EQTIME
C TE=STS-12.0-DELTSL+EQTIME
C WRITE(LOUT,*) '42 ,DELTSL,EQTIME,STB,STS,TB,TE',
C +DELTSL,EQTIME,STB,STR,TB,TE
C 43 CONTINUE
TB=ST1-12.0-DELTSL+EQTIME
TE=ST2-12.0-DELTSL+EQTIME
CDRC WRITE(LOUT,*) '43 ,DELTSL,EQTIME,ST1,ST1,TB,TE',
CDRC +DELTSL,EQTIME,ST1,ST2,TB,TE
TALT=(TB+TE)/2.0
CDRC WRITE(LOUT,*) 'TB,TE',TB,TE
CCC
CCC C.2 Compute amount of clear sky, solar
CCC radiation(SOLAR), and altitude of
CCC the sun (ALPHT). (QUAL2 Step 2-6)
CCC
SOLAR=SOLCON/RR*(SIN(CON2)*SIN(DECLIN)*(TE-TB)+CON6*COS(CON2)*
* COS(DECLIN)*(SIN(CON5*TE)-SIN(CON5*TB)))
ALPHT=SIN(CON2)*SIN(DECLIN)+COS(CON2)*COS(DECLIN)*COS(CON5*TALT)
IF (ABS(ALPHT).EQ.1.0) GO TO 4
Y=SQRT(1.0-ALPHT*ALPHT)
Y=ALPHT/Y
ALPHT=ATAN(Y)
GO TO 5
4 IF (ALPHT.EQ.-1.0) GO TO 6
ALPHT=PI/2.0
GO TO 5
6 ALPHT=-PI/2.0
5 CONTINUE
CDRC write(LOUT,*) 'alpht',alpht
IF (ALPHT.LT.0.01) GO TO 35
CCC
CCC C.3 Compute absorption and scattering due
CCC to atmospheric conditions. (QUAL2
CCC Step 2-7)
CCC
CIPK MAR94
IF(METRIC .EQ. 0) THEN
PWC=0.00614*EXP(0.0489*DEWPT)
ELSE
c PWC=0.00614*EXP(0.0489*(DEWPT*1.8+32.0))
c IF(NN.EQ. 1) WRITE(75,*) 'PWC',PWC
PWC=0.02936*EXP(0.08802*DEWPT)
c IF(NN.EQ. 1) WRITE(75,*) 'PWC',PWC
ENDIF
CIPK MAR94 END CHANGES
OAM=ELEXP/(SIN(ALPHT)+0.15*(ALPHT*CON3+3.885)**(-1.253))
A1=EXP(-(0.465+0.0408*PWC)*(0.129+0.171*EXP(-0.880*OAM))*OAM)
A2=EXP(-(0.465+0.0408*PWC)*(0.179+0.421*EXP(-0.721*OAM))*OAM)
CCC
CCC C.4 Compute reflectivity coefficient (RS).
CCC (QUAL2 Step 2-8)
CCC
GO TO (30,31,31,31,31,31,32,32,32,32,33), NL
30 AR=1.18
BR=-0.77
GO TO 34
31 AR=2.20
BR=-0.97
GO TO 34
32 AR=0.95
BR=-0.75
GO TO 34
33 AR=0.35
BR=-0.45
34 CONTINUE
RS=AR*(CON3*ALPHT)**BR
CDRC write(LOUT,*) 'rs',rs
CC
CC Add test for RS greater than 1.0.
CC
IF(RS.GE.1.0) GO TO 35
CC
CCC
CCC C.5 Compute atmospheric transmission term (ATC).
CCC
ATC=(A2+0.5*(1.0-A1-DAT))/(1.0-0.5*RS*(1.0-A1+DAT))
CCC
CCC C.6 Compute net solar radiaiont for the time
CCC interval delta t. (QUAL2 Step 2-9)
CCC
CDRC write(LOUT,*) 'solar,atc,cs',solar,atc,cs
TSOLHR = TSOLHR+SOLAR*ATC*CS*(1.0-RS)
GO TO 36
35 TSOLHR = TSOLHR+0.0
36 CONTINUE
CLC=1.0+0.17*CLOUD**2
CCC
CCC NCASI Commentary, HEATEX Section D. (QUAL2 Step 3-0)
CCC D. Compute heat fluxes from other terms.
CCC
CCC D.1 Long wave atmospheric radiation (HA).
CCC
CIPK MAR94
IF(METRIC .EQ. 0) THEN
HA = HA+
+ 0.97*1.73E-09*2.89E-06*(DRYBLB+460.0)**6*CLC*DELTH
ELSE
c
c HA(NN) = HA(NN) +
c + 0.97*1.73E-09*2.89E-06*(DRYBLB(NN)*1.8+32.+460.0)**6
c + *CLC*DELTH*4870.8/438.
c IF(NN.EQ. 1) WRITE(75,*) 'HA',HA(NN)
HA = HA +
+ 0.97*9.37e-06*2.0412E-07*(DRYBLB+273.0)**6
+ *CLC*DELTH
c IF(NN.EQ. 1) WRITE(75,*) 'HA',HA
ENDIF
CIPK MAR94 END CHANGES
C WRITE(LOUT,*) 'tofday,tsolhr,ha',
C + TOFDAY, TSOLHR, HA, STB, STE, STR, STS
IF(TLEFT .GT. 0.) GO TO 56
CJFD cheap fix to prevent incrementing time for each material type ...
CC
DELTH = DELT/3600.
RETURN
END

18
source/lenstr.for
Unescape Escape View File

@ -0,0 +1,18 @@
integer function lenstr(string)
character*(*) string
integer val,min,max
min=ichar('!')
max=ichar('z')
lenstr = len(string)
val=ichar(string(lenstr:lenstr))
do while(((val.lt.min).or.(val.gt.max)).and.(lenstr.gt.0))
lenstr = lenstr - 1
val=ichar(string(lenstr:lenstr))
end do
return
end

133
source/random.for
Unescape Escape View File

@ -0,0 +1,133 @@
c FILE : random.f
c ----------------------------------------------------------------------
c ------ Uniform Distribution Random Number Generator ------------------
c ----------------------------------------------------------------------
c From Numerical Recipies.
FUNCTION ran3(idum)
INTEGER*2 i,ii,k,inext,inextp,iff
INTEGER*4 idum,mbig,mseed,mz,mj,ma(55),mk
REAL*4 fac, ran3
parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9)
save inext,inextp,ma,iff
data iff /0/
C
if(idum.lt.0.or.iff.eq.0)then
iff=1
mj=mseed-iabs(idum)
mj=mod(mj,mbig)
ma(55)=mj
mk=1
do i=1,54
ii=mod(21*i,55)
ma(ii)=mk
mk=mj-mk
if(mk.lt.mz) mk=mk+mbig
mj=ma(ii)
end do
do k=1,4
do i=1,55
ma(i)=ma(i)-ma(1+mod(i+30,55))
if(ma(i).lt.mz)ma(i)=ma(i)+mbig
end do
end do
inext=0
inextp=31
idum=1
endif
inext=inext+1
if(inext.eq.56)inext=1
inextp=inextp+1
if(inextp.eq.56)inextp=1
mj=ma(inext)-ma(inextp)
if(mj.lt.mz)mj=mj+mbig
ma(inext)=mj
ran3=real(mj)*fac
return
end
FUNCTION ran4(idum)
INTEGER*2 i,ii,k,inext,inextp,iff
INTEGER*4 idum,mbig,mseed,mz,mj,ma(55),mk
REAL*4 fac, ran4
parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9)
save inext,inextp,ma,iff
data iff /0/
C
if(idum.lt.0.or.iff.eq.0)then
iff=1
mj=mseed-iabs(idum)
mj=mod(mj,mbig)
ma(55)=mj
mk=1
do i=1,54
ii=mod(21*i,55)
ma(ii)=mk
mk=mj-mk
if(mk.lt.mz) mk=mk+mbig
mj=ma(ii)
end do
do k=1,4
do i=1,55
ma(i)=ma(i)-ma(1+mod(i+30,55))
if(ma(i).lt.mz)ma(i)=ma(i)+mbig
end do
end do
inext=0
inextp=31
idum=1
endif
inext=inext+1
if(inext.eq.56)inext=1
inextp=inextp+1
if(inextp.eq.56)inextp=1
mj=ma(inext)-ma(inextp)
if(mj.lt.mz)mj=mj+mbig
ma(inext)=mj
ran4=real(mj)*fac
return
end
c ---------------------------------------------------------------------
c -------- Normally Distributed Random number generator ---------------
c ---------------------------------------------------------------------
c From numerical recipies
c Mean=0, Stdev=1
FUNCTION normdis(idum)
INTEGER iset, idum
REAL v1, v2, r, fac, normdis, ran4, gset
SAVE gset, iset
DATA iset/0/
c idum=0
IF (iset.eq.0) THEN
1 v1=2.*RAN4(idum)-1.
v2=2.*RAN4(idum)-1.
r=v1**2+v2**2
IF (r.ge.1.) GOTO 1
fac=SQRT(-2.*LOG(r)/r)
gset=v1*fac
normdis=v2*fac
iset=1
ELSE
normdis=gset
iset=0
END IF
RETURN
END
c ---------------------------------------------------------------
FUNCTION rbit()
INTEGER rbit, isd
REAL r, ran3
r=ran3(isd)
IF (r.le.0.5) THEN
rbit=1
ELSE
rbit=-1
ENDIF
RETURN
END

2
source/rise_fall.cb
Unescape Escape View File

@ -0,0 +1,2 @@
real y0r,y1r,y0f,y1f,rr,rf
common /risefall/y0r,y1r,y0f,y1f,rr,rf

44
source/rise_fall_vel.for
Unescape Escape View File

@ -0,0 +1,44 @@
c--------------------------------------------------------rise_fall_vel
subroutine rise_fall_vel(vv)
c--------------------------------------------------------------------c
c vv > 0 : rise and vv < 0 : fall vv = 0.0 : nuetral c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c This program is to assign a velocity to one of the c
c sewage particles which have various sizes and densities and c
c paricle size and density are uniformly distributed among c
c the particles c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c yo and y1 are coefficients of speed regression equation
c which has form 10.0**((m-y0)/y1)
c particles which fall have possitive velocity
c particles which rise have negantive velocity
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
include 'rise_fall.cb'
REAL rnumb, vv, ran3, fspe, rspe
common /randm/isd,isd1
INTEGER isd,isd1
c
rnumb=100.0*ran3(isd)
vv=0.0
c
if(rnumb.le.rf) vv=fspe(rnumb)
if(rnumb.ge.(100.0-rr)) vv=rspe(rnumb)
return
end
function fspe(rnumb)
include 'rise_fall.cb'
real fspe, rnumb
fspe=-10.0**((rnumb-y0f)/y1f)/3600.0
return
end
function rspe(rnumb)
include 'rise_fall.cb'
real rspe, rnumb
rspe=10.0**(((100.0-rnumb)-y0r)/y1r)/3600.0
return
end

49
source/tinc2.for
Unescape Escape View File

@ -0,0 +1,49 @@
C
SUBROUTINE TINC2(IHRMNS,IDAYS,MONS,IYEARS,IHRMNF,IDAYF,MONF,
& IYEARF,RMIN,IER,LU)
C
C AUTHOR : D.GORHAM
C DATE : 4-MAR-83
C TITLE : GIVEN START TIME/DATE AND NO. OF MINUTE INTERVALS
C FINISH TIME/DATE
C REV.01 : 22-SEP-83 BY DJG: MODIFIED FOR VAX
C REV.02 : 21-FEB-84 BY DJG: MODIFIED SUCH THAT IF 'RMIN' IS
C NEGATIVE CORRECT TIME IS ALWAYS RETURNED
C
INTEGER*4 EGDS,EGDF
CALL GDATE(1,MONS,IDAYS,IYEARS,IDWK,IDYR,IDMON,EGDS,IER,LU)
IF(IER.EQ.0) GO TO 10
IF(LU.EQ.0) GO TO 10
WRITE(LU,1000)
1000 FORMAT(' **ERROR CALLING START DATE/TIME IN ROUTINE TINC2')
C
10 IHRS=IHRMNS/100
MINS=IHRMNS-IHRS*100
TMINS=RMIN+MINS+IHRS*60
ITMINS=TMINS
NDAYS=ITMINS/1440
IREM=ITMINS-NDAYS*1440
NHRS=IREM/60
NMINS=IREM-NHRS*60
C
IF(NMINS.LT.0) THEN
NHRS=NHRS-1
NMINS=60+NMINS
END IF
C
IF(NHRS.LT.0) THEN
NDAYS=NDAYS-1
NHRS=24+NHRS
END IF
C
EGDF=EGDS+NDAYS
IHRMNF=NHRS*100+NMINS
print *,rmin,nhrs,nmins
C
CALL GDATE(2,MONF,IDAYF,IYEARF,IDWK,IDYR,IDMON,EGDF,IER,LU)!M.W.S.
IF(IER.EQ.0) GO TO 40
C
WRITE(LU,1100)
1100 FORMAT(' **ERROR CALLING FINISH DATE/TIME IN ROUTINE TINC2')
40 RETURN
END

713
source/trinvs.for
Unescape Escape View File

@ -0,0 +1,713 @@
c----------------------------------------------------------------trinvs3
subroutine trinvs3(nel,cno,itype,ebox,egrp,nblk,egps,nelb,nbb,
* nee,nxyz,nen,it,x,y,z,xi,eta,zeta,num,nn)
c----------------------------------------------------------------------c
c purpose: c
c To locate the element the coordinate belongs to and to c
c calculate the shape functions. c
c----------------------------------------------------------------------c
ccc implicit real*8 (a-h,o-z)
common /etype/inode(4)
common /shape/shap(20),shpx(20),shpy(20),shpz(20)
integer*2 nel(20,1)
dimension pos(3),cno(3,1),nbb(1),nee(1),itype(1),ebox(6,1),
* nelb(1),egps(6,1),nblk(2,1),egrp(6,1)
real shap,xi,eta,zeta
c
pos(1)=x
pos(2)=y
pos(3)=z
nn=0
if(num.gt.0) then
c
c----------box test
c
if(pos(1).lt.ebox(1,num).or.pos(1).gt.ebox(4,num)) goto 1
if(pos(2).lt.ebox(2,num).or.pos(2).gt.ebox(5,num)) goto 1
if(pos(3).lt.ebox(3,num).or.pos(3).gt.ebox(6,num)) goto 1
it=itype(num)
nen=inode(it)
call invs_el3(nel,cno,pos,xi,eta,zeta,nen,it,num,key,ii)
if(key.eq.1) goto 30
endif
c call trinvs1(nel,cno,egrp,nblk,egps,nelb,nbb,nee,nxyz,nen,
c * it,x,y,z,xi,eta,zeta,num,nn)
c if(nn.eq.1) goto 30
c goto 35
1 do n=1,nxyz
if(pos(1).gt.egrp(4,n).or.pos(1).lt.egrp(1,n)) goto 10
if(pos(2).gt.egrp(5,n).or.pos(2).lt.egrp(2,n)) goto 10
if(pos(3).gt.egrp(6,n).or.pos(3).lt.egrp(3,n)) goto 10
if(nblk(1,n).gt.1) then
do j=nblk(2,n),nblk(2,n)+nblk(1,n)-1
if(pos(1).gt.egps(4,j).or.pos(1).lt.egps(1,j)) goto 5
if(pos(2).gt.egps(5,j).or.pos(2).lt.egps(2,j)) goto 5
if(pos(3).gt.egps(6,j).or.pos(3).lt.egps(3,j)) goto 5
ng=j
goto 20
5 enddo
else
ng=nblk(2,n)
goto 20
endif
10 enddo
c
c------out of the finite element mesh
c
goto 900
20 do 35 i=nbb(ng),nee(ng)
num=nelb(i)
c
c----------box test
c
if(pos(1).lt.ebox(1,num).or.pos(1).gt.ebox(4,num)) goto 35
if(pos(2).lt.ebox(2,num).or.pos(2).gt.ebox(5,num)) goto 35
if(pos(3).lt.ebox(3,num).or.pos(3).gt.ebox(6,num)) goto 35
it=itype(num)
nen=inode(it)
call invs_el3(nel,cno,pos,xi,eta,zeta,nen,it,num,key,ii)
if(key.eq.1) goto 30
35 continue
c write(*,'(/a/)') 'The input coordinate is out of the FE mesh.'
goto 900
30 nn=1
900 return
1000 format(2i10)
1100 format(i10,3f10.0)
1200 format(21i5)
1300 format(i6,6e12.4,i6)
end
c--------------------------------------------------------------invs_el3
subroutine invs_el3(nel,cno,pos,xi,eta,zeta,nen,it,num,key,ii)
c---------------------------------------------------------------------c
c Purpose: c
c To perform inverse mapping for element i. c
c---------------------------------------------------------------------c
common /shape/shap(20),shpx(20),shpy(20),shpz(20)
integer*2 nel(20,1)
dimension nele(20),cnod(3,20),pos(3),cno(3,1)
real shap,xi,eta,zeta
real a(3,3),b(3),a1,a2,a3
real*8 detj,ass
c
do 40 j=1,nen
nele(j)=nel(j,num)
cnod(1,j)=cno(1,nele(j))
cnod(2,j)=cno(2,nele(j))
cnod(3,j)=cno(3,nele(j))
40 continue
c
c------Newton-Raphson iteration method
c
ii=0
xi0=0.0
eta0=0.0
zeta0=0.0
c
c------calculate the values of shape funcations and their derivatives
c
1 ii=ii+1
call xn3x(it,nen,xi0,eta0,zeta0)
c
c------calculate the Jacobian matrix
c
do 20 j=1,3
a(j,1)=0.0
a(j,2)=0.0
a(j,3)=0.0
b(j)=-pos(j)
do 20 i=1,nen
a(j,1)=a(j,1)+shpx(i)*cnod(j,i)
a(j,2)=a(j,2)+shpy(i)*cnod(j,i)
a(j,3)=a(j,3)+shpz(i)*cnod(j,i)
b(j)=b(j)+shap(i)*cnod(j,i)
20 continue
if(abs(b(1)).lt.1.0e-5.and.abs(b(2)).lt.1.0e-5.and.
& abs(b(3)).lt.1.0e-5) goto 30
do i=1,3
ass=max(abs(a(i,1)),abs(a(i,2)),abs(a(i,3)))
if(dabs(ass).lt.1.0e-10) goto 25
do j=1,3
a(i,j)=a(i,j)/ass
enddo
b(i)=b(i)/ass
25 enddo
a1=a(2,2)*a(3,3)-a(2,3)*a(3,2)
a2=a(3,2)*a(1,3)-a(3,3)*a(1,2)
a3=a(1,2)*a(2,3)-a(1,3)*a(2,2)
detj=a(1,1)*a1+a(2,1)*a2+a(3,1)*a3
if(dabs(detj).lt.1.0e-10) goto 900
xi=xi0-(a1*b(1)+a2*b(2)+a3*b(3))/detj
a1=a(3,1)*a(2,3)-a(3,3)*a(2,1)
a2=a(1,1)*a(3,3)-a(1,3)*a(3,1)
a3=a(2,1)*a(1,3)-a(2,3)*a(1,1)
eta=eta0-(a1*b(1)+a2*b(2)+a3*b(3))/detj
a1=a(2,1)*a(3,2)-a(2,2)*a(3,1)
a2=a(3,1)*a(1,2)-a(3,2)*a(1,1)
a3=a(1,1)*a(2,2)-a(1,2)*a(2,1)
zeta=zeta0-(a1*b(1)+a2*b(2)+a3*b(3))/detj
if(abs(eta-eta0).le.0.5e-4.and.abs(xi-xi0).le.0.5e-4.and.
* abs(zeta-zeta0).le.0.5e-4) goto 30
if(ii.gt.11) goto 30
xi0=xi
eta0=eta
zeta0=zeta
goto 1
c 900 write(6,*) 'Jacobian Determinant = 0',num,nen,it,xi0,eta0,zeta0,
c * detj
c stop
900 goto 35
c goto 35
c
c------determine whether it is inside the element
c
30 goto (50,60,70,80), it
50 if(xi.lt.-5.0e-3.or.eta.lt.-5.0e-3.or.zeta.lt.-5.0e-3)
* goto 35
if((xi+eta+zeta).gt.1.005) goto 35
goto 90
60 if(xi.lt.-5.0e-3.or.xi.gt.1.005) goto 35
if(eta.lt.-1.005.or.eta.gt.1.005) goto 35
if(zeta.lt.-1.005.or.zeta.gt.1.005) goto 35
if((xi+abs(eta)).gt.1.005) goto 35
goto 90
70 if(xi.lt.-1.005.or.xi.gt.1.005) goto 35
if(eta.lt.-1.005.or.eta.gt.1.005) goto 35
if(zeta.lt.-1.005.or.zeta.gt.1.005) goto 35
goto 90
80 if(xi.lt.-5.0e-3.or.xi.gt.1.005) goto 35
if(eta.lt.-1.005.or.eta.gt.1.005) goto 35
if(zeta.lt.-1.005.or.zeta.gt.1.005) goto 35
if((xi+abs(eta)).gt.1.005) goto 35
if((xi+abs(zeta)).gt.1.005) goto 35
90 key=1
return
35 key=0
return
end
c-------------------------------------------------------------------xn3
subroutine xn3(it,nen,x,y,z)
c---------------------------------------------------------------------c
c
c
c subroutine to evaluate shape function
c
c IT is element type 1 10 point tetrahedron
c type 2 15 point prism
c type 3 20 point parallelipiped
c type 4 13 point rectangular base pyramid
c nen is number of shape functions
c X, Y, Z are cordinates of point to be evaluated in local coord
c
c---------------------------------------------------------------------c
ccc implicit real*8 (a-h,o-z)
save
common /shape/shap(20),shpx(20),shpy(20),shpz(20)
dimension irf(10,2)
dimension a(5),b(5),c(5),xt(3),yt(3),jx(3),kx(3),xm(5)
1 ,sn(2),shpp(15,4),ilokup(13)
c
data ilokup/5,6,1,2,3,4,9,7,14,15,10,11,13/
data irf/1,1,2,2,3,3,1,2,3,4
1 ,1,2,2,3,3,1,4,4,4,4/
data xt/0.0,1.0,0.0/,yt/-1.0,0.0,1.0/,jx/2,3,1/,kx/3,1,2/
data shpp/0.0,0.0,1.0,12*0.0,
+ 0.25,-1.0,3.0,-1.0,0.25,4*0.0,0.25,-1.0,0.0,-1.0,0.25,0.0,
+ 0.25,-1.0,0.0,1.0,-0.25,4*0.0,0.25,-1.0,0.0,1.0,-0.25,0.0,
+ 0.25,-1.0,0.0,-1.0,0.25,4*0.0,-0.25,1.0,0.0,1.0,-0.25,0.0/
data ncall/0/
c
if(it .eq. 4) go to 80
if( it - 2 ) 500,80,300
c-
c-----shape functions for right prism and pyramid.....
c-
80 do j0=1,nen
if(it.eq.2) then
i=j0
else
i=ilokup(j0)
endif
ncall = ncall + 1
if( ncall .gt. 1 ) go to 125
c-
c-----calculate invarient triangular functions.....
c-
n = 0
do 100 j = 1,5,2
n = n + 1
jj = jx(n)
kk = kx(n)
a(j) = xt(jj)*yt(kk) - xt(kk)*yt(jj)
b(j) = yt(jj) - yt(kk)
c(j) = xt(kk) - xt(jj)
100 continue
c-
c-.....shape function calculations.....
c-
125 do 130 k = 1, 5, 2
xm(k) = 0.5*(a(k)+b(k)*x+c(k)*y)
130 continue
c-
c-.....set indexes and normalize coordinates.....
c-
if( i .le. 6 ) then
zooo = -z
else
zooo = z
endif
if(it .eq. 2) go to 132
IF(X .EQ. 1.0) GO TO 250
c if(x .eq. 1.0) x=0.999
if(abs(1.0-x) .le. 1.0e-6) then
c print *,'x=',x
goto 250
endif
zooo=zooo/(1.-x)
c if(x .eq. 1.) zooo=sign(x,zooo)
if(i .eq. 3) go to 200
132 if( i .gt. 6 .and. i .lt. 10 ) go to 170
l = i
if( i .gt. 9 ) l = i - 9
if( mod(l,2) .eq. 0 ) go to 150
c-
c-.....corner nodes.....
c-
shap(j0)=xm(l)*((2.0*xm(l)-1.0)*(1.0+zooo)-(1.0-zooo**2))/2.0
if(it .eq. 4) shap(j0)=shap(j0)+xm(l)/2.0*x*(1.-zooo**2)
goto 240
c-
c-.....mid side node.....
c-
150 n1 = l - 1
n2 = mod(l+1,6)
shap(j0) = 2.0*xm(n1)*xm(n2)*(1.0+zooo)
goto 240
c-
c-.....mid side rectangle.....
c-
170 n1 = 1
if( i .eq. 8 ) n1 = 3
if( i .eq. 9 ) n1 = 5
shap(j0) = xm(n1)*(1.0-zooo**2)
if(it .eq. 4) shap(j0)=shap(j0)*(1.-x)
180 goto 240
c-
c......special case for no 3 shape function of pyramid
c-
200 l=i
shap(j0)=xm(l)*(2.*xm(l)-1.0)
goto 240
c
c-----Special case when x exactly equals 1.0
c
250 continue
shap(j0)=shpp(i,1)
240 enddo
return
c-
c-.....shape functions for cube.....
c-
300 zm1=0.125*(1.0-z)
zp1=0.125*(1.0+z)
c-
c......corner nodes
c-
shap(1)=(1.-x)*(1.-y)*zm1*(-x-y-z-2.)
shap(3)=(1.+x)*(1.-y)*zm1*(x-y-z-2.)
shap(5)=(1.+x)*(1.+y)*zm1*(x+y-z-2.)
shap(7)=(1.-x)*(1.+y)*zm1*(-x+y-z-2.)
shap(13)=(1.-x)*(1.-y)*zp1*(-x-y+z-2.)
shap(15)=(1.+x)*(1.-y)*zp1*(x-y+z-2.)
shap(17)=(1.+x)*(1.+y)*zp1*(x+y+z-2.)
shap(19)=(1.-x)*(1.+y)*zp1*(-x+y+z-2.)
c-
c......mid-side nodes
c-
x21=0.25*(1.0-x*x)
y21=0.25*(1.0-y*y)
z21=0.25*(1.0-z*z)
shap(2)=x21*(1.-y)*(1.-z)
shap(4)=(1.+x)*y21*(1.-z)
shap(6)=x21*(1.+y)*(1.-z)
shap(8)=(1.-x)*y21*(1.-z)
shap(9)=(1.-x)*(1.-y)*z21
shap(10)=(1.+x)*(1.-y)*z21
shap(11)=(1.+x)*(1.+y)*z21
shap(12)=(1.-x)*(1.+y)*z21
shap(14)=x21*(1.-y)*(1.+z)
shap(16)=(1.+x)*y21*(1.+z)
shap(18)=x21*(1.+y)*(1.+z)
shap(20)=(1.-x)*y21*(1.+z)
return
500 continue
c-
c......section for tetrahedron
c-
do i=1,nen
i1=irf(i,1)
i2=irf(i,2)
ia=i1
do 550 k=1,2
go to (505,515,525,535),ia
505 sn(k)=1.-x-y-z
go to 550
515 sn(k)=z
go to 550
525 sn(k)=x
go to 550
535 sn(k)=y
550 ia=i2
c-
c......test for corner nodes
c-
if(i1 .ne. i2) then
c-
c......for midsides evaluate function etc
c-
shap(i)=4.*sn(1)*sn(2)
else
c-
c......corner node evaluation
c-
shap(i)=sn(1)*(2.*sn(1)-1.)
endif
enddo
c-
c......final step
c-
return
end
c------------------------------------------------------------------xn3x
subroutine xn3x(it,nen,x,y,z)
c---------------------------------------------------------------------c
c
c
c subroutine to evaluate shape function and its derivatives
c for three dimensions
c
c IT is element type 1 10 point tetrahedron
c type 2 15 point prism
c type 3 20 point parallelipiped
c type 4 13 point rectangular base pyramid
c I is shape function number
c X, Y, Z are cordinates of point to be evaluated in local coord
c
c---------------------------------------------------------------------c
ccc implicit real*8 (a-h,o-z)
save
common /shape/shap(20),shpx(20),shpy(20),shpz(20)
dimension irf(10,2)
dimension a(5),b(5),c(5),xt(3),yt(3),jx(3),kx(3),xm(5)
1 ,sn(2),snx(2),sny(2),snz(2),shpp(15,4),ilokup(13)
c
data ilokup/5,6,1,2,3,4,9,7,14,15,10,11,13/
data irf/1,1,2,2,3,3,1,2,3,4
1 ,1,2,2,3,3,1,4,4,4,4/
data xt/0.0,1.0,0.0/,yt/-1.0,0.0,1.0/,jx/2,3,1/,kx/3,1,2/
data shpp/0.0,0.0,1.0,12*0.0,
+ 0.25,-1.0,3.0,-1.0,0.25,4*0.0,0.25,-1.0,0.0,-1.0,0.25,0.0,
+ 0.25,-1.0,0.0,1.0,-0.25,4*0.0,0.25,-1.0,0.0,1.0,-0.25,0.0,
+ 0.25,-1.0,0.0,-1.0,0.25,4*0.0,-0.25,1.0,0.0,1.0,-0.25,0.0/
data ncall/0/
c
if(it .eq. 4) go to 80
if( it - 2 ) 500,80,300
c-
c-----shape functions for right prism and pyramid.....
c-
80 do j0=1,nen
if(it.eq.2) then
i=j0
else
i=ilokup(j0)
endif
ncall = ncall + 1
if( ncall .gt. 1 ) go to 125
c-
c-----calculate invarient triangular functions.....
c-
n = 0
do 100 j = 1,5,2
n = n + 1
jj = jx(n)
kk = kx(n)
a(j) = xt(jj)*yt(kk) - xt(kk)*yt(jj)
b(j) = yt(jj) - yt(kk)
c(j) = xt(kk) - xt(jj)
100 continue
c-
c-.....shape function calculations.....
c-
125 do 130 k = 1, 5, 2
xm(k) = 0.5*(a(k)+b(k)*x+c(k)*y)
130 continue
c-
c-.....set indexes and normalize coordinates.....
c-
if( i .le. 6 ) then
zooo = -z
else
zooo = z
endif
if(it .eq. 2) go to 132
IF(X .EQ. 1.0) GO TO 250
c if(x .eq. 1.0) x=0.999
if(abs(1.0-x) .le. 1.0e-6) then
c print *,'x=',x
goto 250
endif
zooo=zooo/(1.-x)
c if(x .eq. 1.) zooo=sign(x,zooo)
if(i .eq. 3) go to 200
132 if( i .gt. 6 .and. i .lt. 10 ) go to 170
l = i
if( i .gt. 9 ) l = i - 9
if( mod(l,2) .eq. 0 ) go to 150
c-
c-.....corner nodes.....
c-
shap(j0)=xm(l)*((2.0*xm(l)-1.0)*(1.0+zooo)-(1.0-zooo**2))/2.0
if(it .eq. 4) shap(j0)=shap(j0)+xm(l)/2.0*x*(1.-zooo**2)
shpx(j0)=0.5*b(l)*((1.+zooo)*(2.*xm(l)-0.5)-0.5*(1.0-zooo**2))
if(it .eq. 4) then
shpx(j0)=shpx(j0)+xm(l)/2.*((2.*xm(l)-1.)*zooo
* /(1.-x)+zooo**2+1.)+b(l)/4.*(1.-zooo**2)*x
endif
shpy(j0)=0.5*c(l)*((1.+zooo)*(2.*xm(l)-0.5)-0.5*(1.0-zooo**2))
if(it .eq. 4) shpy(j0)=shpy(j0)+c(l)/4.*(1-zooo**2)*x
shpz(j0) = xm(l)*(xm(l)+zooo-0.5)
if(it .eq. 4) then
shpz(j0)=(shpz(j0)-zooo*xm(l))/(1.-x)+xm(l)*zooo
endif
if(i .le. 6) shpz(j0) = -shpz(j0)
goto 240
c-
c-.....mid side node.....
c-
150 n1 = l - 1
n2 = mod(l+1,6)
shap(j0) = 2.0*xm(n1)*xm(n2)*(1.0+zooo)
shpx(j0) = (1.0+zooo)*(xm(n1)*b(n2) + xm(n2)*b(n1))
if(it .eq. 4) then
shpx(j0)=shpx(j0)+2.*xm(n1)*xm(n2)*zooo/(1.-x)
endif
shpy(j0) = (1.0+zooo)*( xm(n1)*c(n2)+xm(n2)*c(n1))
shpz(j0) = 2.0*xm(n1)*xm(n2)
if(i .le. 6) shpz(j0) = -shpz(j0)
if(it .eq. 2) goto 240
shpz(j0)=shpz(j0)/(1.-x)
goto 240
c-
c-.....mid side rectangle.....
c-
170 n1 = 1
if( i .eq. 8 ) n1 = 3
if( i .eq. 9 ) n1 = 5
shap(j0) = xm(n1)*(1.0-zooo**2)
if(it .eq. 4) shap(j0)=shap(j0)*(1.-x)
shpx(j0) = 0.5*(1.0-zooo**2)*b(n1)
if(it .eq. 4) shpx(j0)=shpx(j0)*(1.-x)-xm(n1)*(zooo**2+1.)
shpy(j0) = 0.5*(1.0-zooo**2)*c(n1)
if(it .eq. 4) shpy(j0)=shpy(j0)*(1.-x)
shpz(j0) = -2.0*xm(n1)*zooo
180 goto 240
c-
c......special case for no 3 shape function of pyramid
c-
200 l=i
shap(j0)=xm(l)*(2.*xm(l)-1.0)
shpx(j0)=b(l)/2.*(4.*xm(l)-1.0)
shpy(j0)=c(l)/2.*(4.*xm(l)-1.0)
shpz(j0)=0.0
goto 240
c
c--------Special case when x exactly equals 1.0
c
250 continue
shap(j0)=shpp(i,1)
shpx(j0)=shpp(i,2)
shpy(j0)=shpp(i,3)
shpy(j0)=shpp(i,4)
240 enddo
return
c-
c-.....shape functions for cube.....
c-
300 zm1=0.125*(1.0-z)
zp1=0.125*(1.0+z)
ym1=0.125*(1.0-y)
yp1=0.125*(1.0+y)
c-
c......corner nodes
c-
shap(1)=(1.-x)*(1.-y)*zm1*(-x-y-z-2.)
shpx(1)=(1.-y)*zm1*(2.*x+y+z+1.)
shpy(1)=(1.-x)*zm1*(x+2.*y+z+1.)
shpz(1)=(1.-x)*ym1*(x+y+2.*z+1.)
shap(3)=(1.+x)*(1.-y)*zm1*(x-y-z-2.)
shpx(3)=(1.-y)*zm1*(2.*x-y-z-1.)
shpy(3)=(1.+x)*zm1*(-x+2.*y+z+1.)
shpz(3)=(1.+x)*ym1*(-x+y+2.*z+1.)
shap(5)=(1.+x)*(1.+y)*zm1*(x+y-z-2.)
shpx(5)=(1.+y)*zm1*(2.*x+y-z-1.)
shpy(5)=(1.+x)*zm1*(x+2.*y-z-1.)
shpz(5)=(1.+x)*yp1*(-x-y+2.*z+1.)
shap(7)=(1.-x)*(1.+y)*zm1*(-x+y-z-2.)
shpx(7)=(1.+y)*zm1*(2.*x-y+z+1.)
shpy(7)=(1.-x)*zm1*(-x+2.*y-z-1.)
shpz(7)=(1.-x)*yp1*(x-y+2.*z+1.)
shap(13)=(1.-x)*(1.-y)*zp1*(-x-y+z-2.)
shpx(13)=(1.-y)*zp1*(2.*x+y-z+1.)
shpy(13)=(1.-x)*zp1*(x+2.*y-z+1.)
shpz(13)=(1.-x)*ym1*(-x-y+2.*z-1.)
shap(15)=(1.+x)*(1.-y)*zp1*(x-y+z-2.)
shpx(15)=(1.-y)*zp1*(2.*x-y+z-1.)
shpy(15)=(1.+x)*zp1*(-x+2.*y-z+1.)
shpz(15)=(1.+x)*ym1*(x-y+2.*z-1.)
shap(17)=(1.+x)*(1.+y)*zp1*(x+y+z-2.)
shpx(17)=(1.+y)*zp1*(2.*x+y+z-1.)
shpy(17)=(1.+x)*zp1*(x+2.*y+z-1.)
shpz(17)=(1.+x)*yp1*(x+y+2.*z-1.)
shap(19)=(1.-x)*(1.+y)*zp1*(-x+y+z-2.)
shpx(19)=(1.+y)*zp1*(2.*x-y-z+1.)
shpy(19)=(1.-x)*zp1*(-x+2.*y+z-1.)
shpz(19)=(1.-x)*yp1*(-x+y+2.*z-1.)
c-
c......mid-side nodes
c-
x21=0.25*(1.0-x*x)
y21=0.25*(1.0-y*y)
z21=0.25*(1.0-z*z)
x2=0.5*x
y2=0.5*y
z2=0.5*z
shap(2)=x21*(1.-y)*(1.-z)
shpx(2)=-x2*(1.-y)*(1.-z)
shpy(2)=-x21*(1.-z)
shpz(2)=-x21*(1.-y)
shap(4)=(1.+x)*y21*(1.-z)
shpx(4)=y21*(1.-z)
shpy(4)=-y2*(1.+x)*(1.-z)
shpz(4)=-y21*(1.+x)
shap(6)=x21*(1.+y)*(1.-z)
shpx(6)=-x2*(1.+y)*(1.-z)
shpy(6)=x21*(1.-z)
shpz(6)=-x21*(1.+y)
shap(8)=(1.-x)*y21*(1.-z)
shpx(8)=-y21*(1.-z)
shpy(8)=-y2*(1.-x)*(1.-z)
shpz(8)=-y21*(1.-x)
shap(9)=(1.-x)*(1.-y)*z21
shpx(9)=-(1.-y)*z21
shpy(9)=-(1.-x)*z21
shpz(9)=-z2*(1.-x)*(1.-y)
shap(10)=(1.+x)*(1.-y)*z21
shpx(10)=(1.-y)*z21
shpy(10)=-(1.+x)*z21
shpz(10)=-z2*(1.+x)*(1.-y)
shap(11)=(1.+x)*(1.+y)*z21
shpx(11)=(1.+y)*z21
shpy(11)=(1.+x)*z21
shpz(11)=-z2*(1.+x)*(1.+y)
shap(12)=(1.-x)*(1.+y)*z21
shpx(12)=-(1.+y)*z21
shpy(12)=(1.-x)*z21
shpz(12)=-z2*(1.-x)*(1.+y)
shap(14)=x21*(1.-y)*(1.+z)
shpx(14)=-x2*(1.-y)*(1.+z)
shpy(14)=-x21*(1.+z)
shpz(14)=x21*(1.-y)
shap(16)=(1.+x)*y21*(1.+z)
shpx(16)=y21*(1.+z)
shpy(16)=-y2*(1.+x)*(1.+z)
shpz(16)=y21*(1.+x)
shap(18)=x21*(1.+y)*(1.+z)
shpx(18)=-x2*(1.+y)*(1.+z)
shpy(18)=x21*(1.+z)
shpz(18)=x21*(1.+y)
shap(20)=(1.-x)*y21*(1.+z)
shpx(20)=-y21*(1.+z)
shpy(20)=-y2*(1.-x)*(1.+z)
shpz(20)=y21*(1.-x)
return
500 continue
c-
c......section for tetrahedron
c-
do i=1,nen
i1=irf(i,1)
i2=irf(i,2)
ia=i1
do 550 k=1,2
go to (505,515,525,535),ia
505 sn(k)=1.-x-y-z
snx(k)=-1.
sny(k)=-1.
snz(k)=-1.
go to 550
515 sn(k)=z
snx(k)=0.
sny(k)=0.
snz(k)=1.
go to 550
525 sn(k)=x
snx(k)=1.
sny(k)=0.
snz(k)=0.
go to 550
535 sn(k)=y
snx(k)=0.
sny(k)=1.
snz(k)=0.
550 ia=i2
c-
c......test for corner nodes
c-
if(i1 .ne. i2) then
c-
c......for midsides evaluate function etc
c-
shap(i)=4.*sn(1)*sn(2)
shpx(i)=4.*(snx(1)*sn(2)+sn(1)*snx(2))
shpy(i)=4.*(sny(1)*sn(2)+sn(1)*sny(2))
shpz(i)=4.*(snz(1)*sn(2)+sn(1)*snz(2))
else
c-
c......corner node evaluation
c-
shap(i)=sn(1)*(2.*sn(1)-1.)
shpx(i)=snx(1)*(4.*sn(1)-1.)
shpy(i)=sny(1)*(4.*sn(1)-1.)
shpz(i)=snz(1)*(4.*sn(1)-1.)
endif
enddo
c-
c......final step
c-
return
end
Write Preview
Loading…
Cancel
Save