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

update to a working version

Browse Source
master
Brett Miller 5 years ago
parent dc2c438b85
commit 4f6ce5a63d
1 changed files with 684 additions and 623 deletions
Show Stats Download Patch File Download Diff File

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

@ -6,6 +6,8 @@ c - FILE: output.f
c -
c - University of NSW - Water Research Laboratory.
c - Brett Miller.
c
c 2019 reworking to simply output a large array in text format
c -
c --------------------------------------------------------------
@ -210,6 +212,12 @@ c ------------------------------------------------------------------
INCLUDE '3dgeom.cb'
INCLUDE '3dpoints.cb'
cBMM2019
character gridfname*100
integer gridnumber
common /bmhack/gridfname,gridnumber
character fname*100
INTEGER*2 nop(20,1)
INTEGER itype(1)
INTEGER partI(1)
@ -239,40 +247,19 @@ c ------------------------------------------------------------------
BYTE bot_pv(MAXXRES,MAXYRES)
BYTE dilution_scale
n1=50
c print*,'Call 2'
call in4bytes(numviews,n1,21)
numviews=numviews+1
if(p_count.and.ngrid.eq.1) then
do i=1,no_poly
no_p(i)=0
ppm(i)=0.0
do k=1,tot_dead_age/3600+1
no_p_time(i,k)=0
enddo
cBM2019
REAL, DIMENSION(:,:,:),allocatable :: grid
allocate (grid(0:MAXXRES,0:MAXYRES,0:MAXZRES))
do x=1,numX
do y=1,numY
do z=1,numZ
grid(x,y,z)=0.0
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
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
DO 5 m=1, NUMPARTS
pel=partE(m)
pnl=partN(m)
pzl=partZ(m)
@ -282,125 +269,198 @@ c IF (partE(m).ne.-1.) THEN
pm=partM(m)
call calc_XYZ(nop,cord,itype,pel,pnl,pzl,pil,pe,pn,pz)
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)
IF (pvv.ne.0.0) THEN
IF (partV(m).eq.-9999.0) THEN
x=INT(a*resXY)+1
y=INT(b*resXY)+1
botlayer(x,y)=botlayer(x,y)+partM(m)
END IF
IF (pz.gt.-0.1.and.pvv.gt.0.0) THEN
x=INT(a*resXY)+1
y=INT(b*resXY)+1
toplayer(x,y)=toplayer(x,y)+partM(m)
nn0=nn0+1
END IF
ELSE
IF (c.ge.0.0.and.c.lt.float(numZ)) THEN
if(a.lt.0.0.or.a.ge.float(numX)) goto 5
if(b.lt.0.0.or.b.ge.float(numY)) goto 5
if(c.lt.0.0.or.c.ge.float(numY)) goto 5
x=INT(a*resXY)+1
y=INT(b*resXY)+1
z=INT(c*resZ)+1
table_x(table_next)=x
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)
z=INT(c*resXY)+1
grid(x,y,z)=grid(x,y,z)+pm
5 enddo
WRITE(fname,'(A,I5.5,A)') gridfname(1:lenstr(gridfname)),
1 gridnumber,'.grid'
OPEN(UNIT=93,file=fname)
do x=1,numX
do y=1,numY
do z=1,numZ
call calc_EN(a,b,float(x),float(y))
write(93,'(2F9.0,F6.2,F15.0)') a,b,z*dimZ,
1 grid(x,y,z)
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)
end do
end do
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
END
@ -717,7 +777,8 @@ c -------------------------------------------------------------------
c ---------------- Encode Graphics Subroutine ----------------------
c ------------------------------------------------------------------
SUBROUTINE encodegraphics(stream,nc,pv,el,MAXXRES,MAXYRES,MAXZRES)
SUBROUTINE encodegraphics(stream,nc,pv,el,MAXXRES,MAXYRES,
1 MAXZRES)
INCLUDE '3dgeom.cb'
INCLUDE '3dpoints.cb'
BYTE pv(MAXXRES,MAXYRES,2)

Write Preview
Loading…
Cancel
Save