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

Compare commits

base: brettm:eaea2bb3e935e75b92881c3d07446d606ca9e8da
Branches Tags
brettm:master
..
compare: brettm:1f24a75037ecd97df6f66b7b074a102c7f0e85ec
Branches Tags
brettm:master

No commits in common. 'eaea2bb3e935e75b92881c3d07446d606ca9e8da' and '1f24a75037ecd97df6f66b7b074a102c7f0e85ec' have entirely different histories.
eaea2bb3e9 ... 1f24a75037

3 changed files with 796 additions and 950 deletions
Show Stats

133
source/3dcondits.for
Unescape Escape View File

@ -45,14 +45,18 @@ c---------------------------------------------------------------------c
WRITE(*,'(2A15)') 'resXY', 'resZ' WRITE(*,'(2A15)') 'resXY', 'resZ'
WRITE(*,42) resXY, resZ WRITE(*,42) resXY, resZ
WRITE(*,*) WRITE(*,*)
WRITE(*,'(5A15)') 'newcondTS','outputTS','dt1','dt2', WRITE(*,'(5A15)') 'newcondTS','outputTS','dt1','dt2','maxoutputs'
1 'maxoutputs'
WRITE(*,50) newcondTS, outputTS, dt,dtnew, maxoutputs WRITE(*,50) newcondTS, outputTS, dt,dtnew, maxoutputs
WRITE(*,*) WRITE(*,*)
c WRITE(*,'(5A15)') 'tot_dead_age', 'HDIFF', 'VDIFF', 'mass_pp'
c * ,'c_init'
c WRITE(*,60) tot_dead_age, HDIFF, VDIFF, mass_pp,c_init
c BMM modifications made 8/1/96 to print out extra diffusivities ver 2.2
WRITE(*,'(3A15)') 'tot_dead_age', 'mass_pp', 'c_init' WRITE(*,'(3A15)') 'tot_dead_age', 'mass_pp', 'c_init'
WRITE(*,60) tot_dead_age, mass_pp, c_init WRITE(*,60) tot_dead_age, mass_pp, c_init
WRITE(*,*) WRITE(*,*)
WRITE(*,'(A)') 'Diffs: Element types not listed use the first' WRITE(*,'(A)') 'Diffusivities: Element types not listed use the'
* //' first'
WRITE(*,'(3A15)') 'Element type', 'HDIFF','VDIFF' WRITE(*,'(3A15)') 'Element type', 'HDIFF','VDIFF'
i=1 i=1
DO WHILE (EDIFF(i).ne.-1) DO WHILE (EDIFF(i).ne.-1)
@ -125,7 +129,8 @@ c *** Read in the new pollutant releases.
DO 60 w=1, num_pollwins DO 60 w=1, num_pollwins
READ(4,*,ERR=200,end=205) d,(poll(i,d),i=1,4) READ(4,*,ERR=200,end=205) d,(poll(i,d),i=1,4)
IF (d.gt.num_pollwins.or.d.lt.1) THEN IF (d.gt.num_pollwins.or.d.lt.1) THEN
PRINT *,'ERROR - Pollutant window' PRINT *,'ERROR - Pollutant window not as '//
1 'expected.'
STOP STOP
ENDIF ENDIF
masspB=mass_pp masspB=mass_pp
@ -157,8 +162,6 @@ c-----------------------------------------------------getflowconditions
c---------------------------------------------------------------------c c---------------------------------------------------------------------c
c purpose: c c purpose: c
c To input velocity field. c c To input velocity field. c
c This reads one time step as called from walkem c
c If it is a 2D file then it makes an artificial 3d one c
c---------------------------------------------------------------------c c---------------------------------------------------------------------c
INCLUDE '3duparms.cb' INCLUDE '3duparms.cb'
INCLUDE '3dmesh.cb' INCLUDE '3dmesh.cb'
@ -173,7 +176,7 @@ c vel(3,i) = Velocity to Z at node i
INTEGER*2 nsurf(1),nop(20,1),ic(1),ielist(1) INTEGER*2 nsurf(1),nop(20,1),ic(1),ielist(1)
INTEGER i,k,j,no3dn,ndg,no3de,key INTEGER i,k,j,no3dn,ndg,no3de,key
REAL vel1,vel2,depth,vel4,vel5,vel6,tet,dfct,wsel(1),wsel1(1) REAL vel1,vel2,vel4,vel5,vel6,tet,dfct,wsel(1),wsel1(1)
real dhdt(1),prsm(3,6),brck(3,8),xi(8),eta(8),zeta(8) real dhdt(1),prsm(3,6),brck(3,8),xi(8),eta(8),zeta(8)
real cord(3,1),a(2,2) real cord(3,1),a(2,2)
real*8 ao(1),detj real*8 ao(1),detj
@ -187,57 +190,23 @@ c *** Read in the velocity field
key=0 key=0
if(dimopt(1:2).eq.'10') then if(dimopt(1:2).eq.'10') then
cBM2019 updating reads to OUTBNRMA format
if(compact) then if(compact) then
c read(2,err=210,end=215) TET,no3dn,NDG,no3de, read(2,err=210,end=215) TET,no3dn,NDG,no3de,
c 1 ((VEL(K,J),K=1,2),wsel(j),VEL(3,j),J=1,No3dn) 1 ((VEL(K,J),K=1,2),wsel(j),VEL(3,j),J=1,No3dn)
WRITE(*,*) 'No compact!'
STOP
else else
cDRC 04/06/99 IYRR added for RMA10v6.5
c read(2,err=210,end=215) TET,no3dn,NDG,no3de,IYRR,
c 1 ((VEL(K,J),K=1,2),wsel(j),vel4,vel5,vel6,
c 2 VEL(3,j),J=1,No3dn),(DFCT,J=1,No3dE)
read(2,err=210,end=215) TET,no3dn,NDG,no3de,IYRR, read(2,err=210,end=215) TET,no3dn,NDG,no3de,IYRR,
1 ((VEL(K,J),K=1,2),depth,vel4,vel5,vel6, 1 ((VEL(K,J),K=1,2),wsel(j),vel4,vel5,vel6,
2 VEL(3,j),wsel(j),J=1,no3dn), 2 VEL(3,j),J=1,No3dn),(DFCT,J=1,No3dE)
3 (DFCT,J=1,no3de),
4 (vel7,j=1,no3dn)
c WRITE(IRMAFM) TETT,NP,NDG,NE,IYRR,
c 1 ((VSING(K,J),K=1,6),VVEL(J),WSLL(J),J=1,NP),(DFCT(J),J=1,NE),(VSING(7,J),J=1,NP)
C Output RMA results file contains
c 1 time in hours (Julian)
c 2 number of nodes
c 3 obsolete counter of degrees of freedom (set to 5) NDF=6
c 3a number of elements
c 4 year
c 5 VSING subscript(1) x-vel by node
c 6 VSING subscript(2) y-vel by node
c 7 VSING subscript(3) depth by node
c 8 VSING subscript(4) salinity by node
c 9 VSING subscript(5) temperature by node
c 10 VSING subscript(6) sus-sed by node
c 11 VVEL w-vel by node
c 12 WSLL water surface elevation
c 13 DFCT stratification multiplier by element
c 14 VSING subscript(7) time derivative of depth by node
print*,'Velocity file read at time',TET print*,'Velocity file read at time',TET
endif endif
else else
c read(2,err=210,end=215) TET,no3dn,((vel(k,j),k=1,2), read(2,err=210,end=215) TET,no3dn,((vel(k,j),k=1,2),
c 1 wsel(j),j=1,no3dn),(vel4,j=1,no3dn) 1 wsel(j),j=1,no3dn),(vel4,j=1,no3dn)
WRITE(*,*) 'No RMA2!'
STOP
endif endif
c c
c-------For 2D case - START c-------for 2D case
c-------artificially create a 3D velocity field c
c-------The 2D field could be from RMA2 or a depth averaged RMA10
if(no3dn.ne.np) then if(no3dn.ne.np) then
do i=1,np do i=1,np
ic(i)=0 ic(i)=0
@ -366,11 +335,7 @@ c
c c
c----------calculate dhdt(i) for next time step c----------calculate dhdt(i) for next time step
c c
c goto 40
cBM2019
WRITE(*,*) 'Shouldnt be here'
STOP
c
if(dimopt(1:2).eq.'10') then if(dimopt(1:2).eq.'10') then
if(compact) then if(compact) then
read(2,err=30,end=30) TET,no3dn,NDG,no3de, read(2,err=30,end=30) TET,no3dn,NDG,no3de,
@ -397,9 +362,12 @@ cDRC 04/06/99 IYRR added for RMA10v6.5
30 continue 30 continue
endif endif
c-------For 2D case - END
40 key=1 40 key=1
c do i=1,No3dn
c do j=1,3
c vel(j,i)=0.0
c enddo
c enddo
RETURN RETURN
210 WRITE(*,*) ' *** Error in velocities! *** ' 210 WRITE(*,*) ' *** Error in velocities! *** '
@ -419,55 +387,21 @@ c---------------------------------------------------------------------c
include '3dmesh.cb' include '3dmesh.cb'
integer j,k integer j,k
integer nfix,nfix1,nref integer nfix,nfix1,nref
integer*2 imat(1),ncorn,iem,ndep,netyp integer*2 imat(1),ncorn,nfixh,iem,ndep,netyp
integer*4 nfixh !BM2019 changed from integer*2
real width(1),paxis,alfa,th real width(1),paxis,alfa,th
integer*2 nop(20,1),nsurf(1),ielist(1)
cBM2019 real cord(3,1),spec(3)
c the 3dg file now has cord as real*8
c and nop as real*4
integer*2 nop(20,1),nsurf(1)
integer nop4(30000,20),nsurf4(30000)
integer*2 ielist(1) ! this one OK as it is created internally
real*8 cord8(3,30000) ! BM2019 hack because of real*8
real cord(3,1)
real spec(3)
real*8 ao(1) real*8 ao(1)
real diffs(2,1) real diffs(2,1)
c c
if(dimopt(1:2).eq.'10') then if(dimopt(1:2).eq.'10') then
read(9,err=10) NP,NE,NPM,NES,((CORD(k,j),SPEC(k),K=1,3),
read(9,err=10) NP,NE,NPM,NES,((CORD8(k,j),SPEC(k),K=1,3), 1 ALFA,NFIX,NFIX1,AO(J),NSURF(J),J=1,NP),
1 ALFA,NFIX,NFIX1,AO(J),NSURF4(J),J=1,NP),
2 (NDEP,NREF,J=1,NPM), 2 (NDEP,NREF,J=1,NPM),
3 ((NOP4(k,j),K=1,20),NCORN,IMAT(J),NETYP,TH, cdrc16/6/99 3 ((NOP(k,j),K=1,20),NCORN,IMAT(J),IMAT(J),TH,
3 ((NOP(k,j),K=1,20),NCORN,IMAT(J),NETYP,TH,
4 NFIXH,J=1,NE),(WIDTH(J),J=1,NP) 4 NFIXH,J=1,NE),(WIDTH(J),J=1,NP)
do j=1,np
do k=1,3
cord(k,j)=cord8(k,j)
enddo
c write(58,'(I5,3F20.3)') j,(cord(k,j),k=1,3)
enddo
do j=1,ne
do k=1,20
nop(k,j)=nop4(k,j)
enddo
c write(59,'(i5,20I8)') j,(nop(k,j),k=1,20)
enddo
do j=1,np
nsurf(j)=nsurf4(j)
enddo
else else
cBM2019
write(*,*) 'Stop: Geometry: No 2D'
stop
read(9,err=10) NP,NE,((CORD(K,J),K=1,2),ALFA,cord(3,j), read(9,err=10) NP,NE,((CORD(K,J),K=1,2),ALFA,cord(3,j),
1 J=1,NP),((NOP(K,J),K=1,8),IMAT(j),PAXIS,IEM,J=1,NE) 1 J=1,NP),((NOP(K,J),K=1,8),IMAT(j),PAXIS,IEM,J=1,NE)
do i=1,np do i=1,np
@ -512,9 +446,6 @@ c print *,i,(nop(k,i),k=1,nen)
c c
c----------generate nodes below surface c----------generate nodes below surface
c c
cbm2019
write(*,*) 'Stop: no 2d!'
stop
np3=np np3=np
do n=1,np do n=1,np
cord(3,n)=0.0 cord(3,n)=0.0

351
source/3doutput-block.for
Unescape Escape View File

@ -6,8 +6,6 @@ c - FILE: output.f
c - c -
c - University of NSW - Water Research Laboratory. c - University of NSW - Water Research Laboratory.
c - Brett Miller. c - Brett Miller.
c
c 2019 reworking to simply output a large array in text format
c - c -
c -------------------------------------------------------------- c --------------------------------------------------------------
@ -212,12 +210,6 @@ c ------------------------------------------------------------------
INCLUDE '3dgeom.cb' INCLUDE '3dgeom.cb'
INCLUDE '3dpoints.cb' INCLUDE '3dpoints.cb'
cBMM2019
character gridfname*100
integer gridnumber
common /bmhack/gridfname,gridnumber
character fname*100
INTEGER*2 nop(20,1) INTEGER*2 nop(20,1)
INTEGER itype(1) INTEGER itype(1)
INTEGER partI(1) INTEGER partI(1)
@ -247,19 +239,40 @@ cBMM2019
BYTE bot_pv(MAXXRES,MAXYRES) BYTE bot_pv(MAXXRES,MAXYRES)
BYTE dilution_scale BYTE dilution_scale
cBM2019 n1=50
REAL, DIMENSION(:,:,:),allocatable :: grid c print*,'Call 2'
call in4bytes(numviews,n1,21)
allocate (grid(0:MAXXRES,0:MAXYRES,0:MAXZRES)) numviews=numviews+1
if(p_count.and.ngrid.eq.1) then
do x=1,numX do i=1,no_poly
do y=1,numY no_p(i)=0
do z=1,numZ ppm(i)=0.0
grid(x,y,z)=0.0 do k=1,tot_dead_age/3600+1
no_p_time(i,k)=0
enddo
enddo enddo
endif
c *************************
c STEP 1 - The lookup table
c *************************
do x=1,MAXXRES
do y=1, MAXYRES
toplayer(x,y)=0.0
botlayer(x,y)=0.0
end do end do
end do end do
DO 5 m=1, NUMPARTS
nn0=0
DO 10 y=1, numY*resXY
ylookup(y)=-1
10 CONTINUE
table_next=1
c No need to initialise the table as the list will only point to
c used portions.
DO 20 m=1, NUMPARTS
c IF (partE(m).ne.-1.) THEN
pel=partE(m) pel=partE(m)
pnl=partN(m) pnl=partN(m)
pzl=partZ(m) pzl=partZ(m)
@ -269,198 +282,125 @@ cBM2019
pm=partM(m) pm=partM(m)
call calc_XYZ(nop,cord,itype,pel,pnl,pzl,pil,pe,pn,pz) call calc_XYZ(nop,cord,itype,pel,pnl,pzl,pil,pe,pn,pz)
CALL calc_ab(a, b, pe, pn) CALL calc_ab(a, b, pe, pn)
if(a.lt.0.0.or.a.ge.float(numX)) goto 20 !added by Y.C. WANG
if(b.lt.0.0.or.b.ge.float(numY)) goto 20 !added by Y.C. WANG
c
c---------------counts particles in the specified areas
c
if(p_count.and.ngrid.eq.1) call particles_c(cord,nply,
* nump,no_p,no_p_time,ppm,no_poly,pe,pn,pa,pm)
CALL calc_c(c, pz) CALL calc_c(c, pz)
if(a.lt.0.0.or.a.ge.float(numX)) goto 5 IF (pvv.ne.0.0) THEN
if(b.lt.0.0.or.b.ge.float(numY)) goto 5 IF (partV(m).eq.-9999.0) THEN
if(c.lt.0.0.or.c.ge.float(numY)) goto 5
x=INT(a*resXY)+1 x=INT(a*resXY)+1
y=INT(b*resXY)+1 y=INT(b*resXY)+1
z=INT(c*resXY)+1 botlayer(x,y)=botlayer(x,y)+partM(m)
grid(x,y,z)=grid(x,y,z)+pm END IF
5 enddo IF (pz.gt.-0.1.and.pvv.gt.0.0) THEN
x=INT(a*resXY)+1
y=INT(b*resXY)+1
WRITE(fname,'(A,I5.5,A)') gridfname(1:lenstr(gridfname)), toplayer(x,y)=toplayer(x,y)+partM(m)
1 gridnumber,'.grid' nn0=nn0+1
OPEN(UNIT=93,file=fname) END IF
do x=1,numX ELSE
do y=1,numY IF (c.ge.0.0.and.c.lt.float(numZ)) THEN
do z=1,numZ x=INT(a*resXY)+1
call calc_EN(a,b,float(x),float(y)) y=INT(b*resXY)+1
write(93,'(2F9.0,F6.2,F15.0)') a,b,z*dimZ, z=INT(c*resZ)+1
1 grid(x,y,z) table_x(table_next)=x
enddo table_z(table_next)=z
table_M(table_next)=partM(m)
table_point(table_next)=ylookup(y)
ylookup(y)=table_next
table_next=table_next+1
END IF
END IF
c END IF
20 CONTINUE
c
c-------print the proportions of particles in the specified areas
c
if(p_count.and.ngrid.eq.1) then
write(36,'(/a/)') ' Area No of Particles'//
& ' % mass'
do i=1,no_poly
if(NUMPARTS.ne.0) then
write(36,'(2i10,f23.2,e12.3)') i,no_p(i),
& 100.0*no_p(i)/NUMPARTS,ppm(i)
else
write(36,'(2i10,f23.2,e12.3)') i,no_p(i),0.0,ppm(i)
endif
write(36,'(/a,i10/)') ' Total No of Particles=',NUMPARTS
write(36,'(a/)') ' Area Hour No of Particles'
do k=1,tot_dead_age/3600
write(36,'(2i5,i10)') i,k,no_p_time(i,k)
enddo enddo
enddo enddo
write(*,*)
endif
c ************************
c STEP 2 - Run the slices
c ************************
left=1
centre=2
right=3
CALL clearslice(slice,left, MAXXRES, MAXZRES)
CALL clearslice(slice,centre, MAXXRES, MAXZRES)
CALL clearslice(slice,right, MAXXRES, MAXZRES)
y=1
point=ylookup(y)
DO 200 WHILE (point.ne.-1)
x=table_x(point)
z=table_z(point)
mass=table_M(point)
CALL addyslice(x,y,z,mass,slice,left,centre,right,
* MAXXRES,MAXZRES)
point=table_point(point)
200 END DO
DO 800 y=2,numY*resXY
point=ylookup(y)
DO 300 WHILE (point.ne.-1)
x=table_x(point)
z=table_z(point)
mass=table_M(point)
CALL addyslice(x,y,z,mass,slice,left,centre,
* right,MAXXRES,MAXZRES)
point=table_point(point)
300 CONTINUE
CALL establishview(y-1, pv, el, slice, left, MAXXRES,
* MAXYRES,MAXZRES,ngrid)
dummy=left
left=centre
centre=right
right=dummy
CALL clearslice(slice,right, MAXXRES, MAXZRES)
800 CONTINUE
y=numY*resXY
CALL establishview(y, pv, el, slice, left, MAXXRES,
* MAXYRES, MAXZRES,ngrid)
if(dilution) then
CALL encodegraphics(21, 1, pv, el, MAXXRES, MAXYRES, MAXZRES)
endif
if(concentration) then
CALL encodegraphics(22, 2, pv, el, MAXXRES, MAXYRES, MAXZRES)
endif
DO x=1, numX*resXY
DO y=1, numY*resXY
top_pv(x,y)=dilution_scale(toplayer(x,y))
bot_pv(x,y)=dilution_scale(botlayer(x,y))
END DO
END DO
CALL encodedoublegraphics(23, top_pv, bot_pv, MAXXRES, MAXYRES)
if(fbindump) call fbindumproutine(toplayer, MAXXRES, MAXYRES)
CLOSE(93)
gridnumber=gridnumber+1
deallocate(grid)
ccc
ccc n1=50
cccc print*,'Call 2'
ccc call in4bytes(numviews,n1,21)
ccc numviews=numviews+1
ccc if(p_count.and.ngrid.eq.1) then
ccc do i=1,no_poly
ccc no_p(i)=0
ccc ppm(i)=0.0
ccc do k=1,tot_dead_age/3600+1
ccc no_p_time(i,k)=0
ccc enddo
ccc enddo
ccc endif
cccc *************************
cccc STEP 1 - The lookup table
cccc *************************
ccc
ccc do x=1,MAXXRES
ccc do y=1, MAXYRES
ccc toplayer(x,y)=0.0
ccc botlayer(x,y)=0.0
ccc end do
ccc end do
ccc
ccc nn0=0
ccc
ccc DO 10 y=1, numY*resXY
ccc ylookup(y)=-1
ccc10 CONTINUE
ccc table_next=1
cccc No need to initialise the table as the list will only point to
cccc used portions.
ccc DO 20 m=1, NUMPARTS
cccc IF (partE(m).ne.-1.) THEN
ccc pel=partE(m)
ccc pnl=partN(m)
ccc pzl=partZ(m)
ccc pil=partI(m)
ccc pa=partA(m)
ccc pvv=partV(m)
ccc pm=partM(m)
ccc call calc_XYZ(nop,cord,itype,pel,pnl,pzl,pil,pe,pn,pz)
ccc CALL calc_ab(a, b, pe, pn)
ccc if(a.lt.0.0.or.a.ge.float(numX)) goto 20 !added by Y.C. WANG
ccc if(b.lt.0.0.or.b.ge.float(numY)) goto 20 !added by Y.C. WANG
cccc
cccc---------------counts particles in the specified areas
cccc
ccc if(p_count.and.ngrid.eq.1) call particles_c(cord,nply,
ccc * nump,no_p,no_p_time,ppm,no_poly,pe,pn,pa,pm)
ccc
ccc CALL calc_c(c, pz)
ccc IF (pvv.ne.0.0) THEN
ccc IF (partV(m).eq.-9999.0) THEN
ccc x=INT(a*resXY)+1
ccc y=INT(b*resXY)+1
ccc botlayer(x,y)=botlayer(x,y)+partM(m)
ccc END IF
ccc IF (pz.gt.-0.1.and.pvv.gt.0.0) THEN
ccc x=INT(a*resXY)+1
ccc y=INT(b*resXY)+1
ccc toplayer(x,y)=toplayer(x,y)+partM(m)
ccc nn0=nn0+1
ccc END IF
ccc ELSE
ccc IF (c.ge.0.0.and.c.lt.float(numZ)) THEN
ccc x=INT(a*resXY)+1
ccc y=INT(b*resXY)+1
ccc z=INT(c*resZ)+1
ccc table_x(table_next)=x
ccc table_z(table_next)=z
ccc table_M(table_next)=partM(m)
ccc table_point(table_next)=ylookup(y)
ccc ylookup(y)=table_next
ccc table_next=table_next+1
ccc END IF
ccc END IF
cccc END IF
ccc20 CONTINUE
ccc
cccc
cccc-------print the proportions of particles in the specified areas
cccc
ccc if(p_count.and.ngrid.eq.1) then
ccc write(36,'(/a/)') ' Area No of Particles'//
ccc & ' % mass'
ccc do i=1,no_poly
ccc if(NUMPARTS.ne.0) then
ccc write(36,'(2i10,f23.2,e12.3)') i,no_p(i),
ccc & 100.0*no_p(i)/NUMPARTS,ppm(i)
ccc else
ccc write(36,'(2i10,f23.2,e12.3)') i,no_p(i),0.0,ppm(i)
ccc endif
ccc write(36,'(/a,i10/)') ' Total No of Particles=',NUMPARTS
ccc write(36,'(a/)') ' Area Hour No of Particles'
ccc do k=1,tot_dead_age/3600
ccc write(36,'(2i5,i10)') i,k,no_p_time(i,k)
ccc enddo
ccc enddo
ccc write(*,*)
ccc endif
ccc
ccc
cccc ************************
cccc STEP 2 - Run the slices
cccc ************************
ccc
ccc left=1
ccc centre=2
ccc right=3
ccc CALL clearslice(slice,left, MAXXRES, MAXZRES)
ccc CALL clearslice(slice,centre, MAXXRES, MAXZRES)
ccc CALL clearslice(slice,right, MAXXRES, MAXZRES)
ccc y=1
ccc point=ylookup(y)
ccc DO 200 WHILE (point.ne.-1)
ccc x=table_x(point)
ccc z=table_z(point)
ccc mass=table_M(point)
ccc CALL addyslice(x,y,z,mass,slice,left,centre,right,
ccc * MAXXRES,MAXZRES)
ccc point=table_point(point)
ccc200 END DO
ccc DO 800 y=2,numY*resXY
ccc point=ylookup(y)
ccc DO 300 WHILE (point.ne.-1)
ccc x=table_x(point)
ccc z=table_z(point)
ccc mass=table_M(point)
ccc CALL addyslice(x,y,z,mass,slice,left,centre,
ccc * right,MAXXRES,MAXZRES)
ccc point=table_point(point)
ccc300 CONTINUE
ccc CALL establishview(y-1, pv, el, slice, left, MAXXRES,
ccc * MAXYRES,MAXZRES,ngrid)
ccc dummy=left
ccc left=centre
ccc centre=right
ccc right=dummy
ccc CALL clearslice(slice,right, MAXXRES, MAXZRES)
ccc800 CONTINUE
ccc y=numY*resXY
ccc CALL establishview(y, pv, el, slice, left, MAXXRES,
ccc * MAXYRES, MAXZRES,ngrid)
ccc if(dilution) then
ccc CALL encodegraphics(21, 1, pv, el, MAXXRES, MAXYRES, MAXZRES)
ccc endif
ccc if(concentration) then
ccc CALL encodegraphics(22, 2, pv, el, MAXXRES, MAXYRES, MAXZRES)
ccc endif
ccc
ccc
ccc DO x=1, numX*resXY
ccc DO y=1, numY*resXY
ccc top_pv(x,y)=dilution_scale(toplayer(x,y))
ccc bot_pv(x,y)=dilution_scale(botlayer(x,y))
ccc END DO
ccc END DO
ccc CALL encodedoublegraphics(23, top_pv, bot_pv, MAXXRES, MAXYRES)
ccc
ccc if(fbindump) call fbindumproutine(toplayer, MAXXRES, MAXYRES)
ccc
print *,'Return from outputting.'
RETURN RETURN
END END
@ -777,8 +717,7 @@ c -------------------------------------------------------------------
c ---------------- Encode Graphics Subroutine ---------------------- c ---------------- Encode Graphics Subroutine ----------------------
c ------------------------------------------------------------------ c ------------------------------------------------------------------
SUBROUTINE encodegraphics(stream,nc,pv,el,MAXXRES,MAXYRES, SUBROUTINE encodegraphics(stream,nc,pv,el,MAXXRES,MAXYRES,MAXZRES)
1 MAXZRES)
INCLUDE '3dgeom.cb' INCLUDE '3dgeom.cb'
INCLUDE '3dpoints.cb' INCLUDE '3dpoints.cb'
BYTE pv(MAXXRES,MAXYRES,2) BYTE pv(MAXXRES,MAXYRES,2)

40
source/3drwalk.for
Unescape Escape View File

@ -28,18 +28,8 @@ c the same in the main block and the subroutines.
include '3dpoints.cb' include '3dpoints.cb'
include 'rise_fall.cb' include 'rise_fall.cb'
parameter (ntotal=250000000) parameter (ntotal=50000000)
character fname*100,buff*100,partl*100,fltnam*100 character fname*100,buff*100,partl*100,fltnam*100
cBMM2019
character header*1000
cBMM2019
character gridfname*100
integer gridnumber
common /bmhack/gridfname,gridnumber
character plmnam*100,cext*5 character plmnam*100,cext*5
real r real r
integer t,NUMPARTS,max_numps,v_time integer t,NUMPARTS,max_numps,v_time
@ -101,24 +91,15 @@ c Read in the user paramters
fname='920820' fname='920820'
call askstr(' Enter output graphics filename',fname) call askstr(' Enter output graphics filename',fname)
partl=fname partl=fname
gridfname=fname
gridnumber=1
dimopt='10' dimopt='10'
call askstr(' Enter RMA-2 or RMA-10 option',dimopt) call askstr(' Enter RMA-2 or RMA-10 option',dimopt)
if(dimopt(1:2).eq.'10') then if(dimopt(1:2).eq.'10') then
fname='bub003.3dg' fname='sydney.3dg'
call askstr(' Enter 3d geometry filename',fname) call askstr(' Enter 3d geometry filename',fname)
cBM2019 update for OUTBN3GE format open(9,file=fname,status='old',form='unformatted',err=10)
c open(9,file=fname,status='old',form='unformatted',err=10) fname='sydney.res'
OPEN(9,FILE=fname,STATUS='old',FORM='BINARY',err=10)
READ(9) header ! Skip the 1000 character header
cBM2019 fname='sydney.res'
fname='bub007.rma'
call askstr(' Enter RMA-10 velocity filename',fname) call askstr(' Enter RMA-10 velocity filename',fname)
cBM2019 open(2,file=fname,status='old',form='unformatted',err=10) open(2,file=fname,status='old',form='unformatted',err=10)
cBM2019 update for OUTBNRMA format
OPEN(2,FILE=fname,STATUS='old',FORM='BINARY',err=10)
READ(2) header ! Skip the 1000 character header
else else
fname='sydney.geo' fname='sydney.geo'
call askstr(' Enter 2D geometry filename',fname) call askstr(' Enter 2D geometry filename',fname)
@ -130,7 +111,6 @@ cBM2019 update for OUTBNRMA format
call askstr(' Enter RMA-2 velocity filename',fname) call askstr(' Enter RMA-2 velocity filename',fname)
open(2,file=fname,status='old',form='unformatted',err=10) open(2,file=fname,status='old',form='unformatted',err=10)
endif endif
fname='usrparam.rw' fname='usrparam.rw'
call askstr(' Enter user parameters filename',fname) call askstr(' Enter user parameters filename',fname)
open(3,file=fname,status='old',form='formatted',err=10) open(3,file=fname,status='old',form='formatted',err=10)
@ -228,12 +208,8 @@ c
write(*,'(/a/)') 'Loading the geometry of FE mesh ...' write(*,'(/a/)') 'Loading the geometry of FE mesh ...'
rewind 9 rewind 9
cBM2019
read(9) header
read(9) np,ne read(9) np,ne
cBM2019
rewind 9 rewind 9
read(9) header
call defini('nop ',nnop ,10,ne+4) call defini('nop ',nnop ,10,ne+4)
call defini('ilst',nilst,1,ne/2+4) call defini('ilst',nilst,1,ne/2+4)
call defini('imat',nimat,1,ne/2+4) call defini('imat',nimat,1,ne/2+4)
@ -703,9 +679,9 @@ c
resZ=reszz(m) resZ=reszz(m)
call outputgraphicscodes(m) call outputgraphicscodes(m)
c call addpointers(21) call addpointers(21)
c call addpointers(22) call addpointers(22)
c call addpointers(23) call addpointers(23)
close(21) close(21)
close(22) close(22)
close(23) close(23)

Write Preview
Loading…
Cancel
Save