You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1015 lines
26 KiB
Fortran
1015 lines
26 KiB
Fortran
c --------------------------------------------------------------
|
|
c -
|
|
c - RWALK - Random Walk Model.
|
|
c - (OUTPUT ENCODING ROUTINES)
|
|
c - FILE: output.f
|
|
c -
|
|
c - University of NSW - Water Research Laboratory.
|
|
c - Brett Miller.
|
|
c -
|
|
c --------------------------------------------------------------
|
|
|
|
SUBROUTINE outputgraphicscodes(ngrid)
|
|
c--------------------------------------------------------------------c
|
|
common mtot,npp,ia(1)
|
|
common /dbsys/numa,next,idir,ipp(3)
|
|
integer ngrid
|
|
include '3dgeom.cb'
|
|
|
|
WRITE(*,*) 'Outputting...'
|
|
nyy=numY*resXY+1
|
|
nxx=numX*resXY+1
|
|
nzz=numZ*resZ+1
|
|
call locate('nop ',nnop ,mlt,mlt)
|
|
call locate('ityp',nityp,mlt,mlt)
|
|
call locate('cord',ncord,mlt,mlt)
|
|
call locate('prtE',nprtE,mlt,NUMPARTS)
|
|
call locate('prtN',nprtN,mlt,NUMPARTS)
|
|
call locate('prtZ',nprtZ,mlt,NUMPARTS)
|
|
call locate('prtI',nprtI,mlt,NUMPARTS)
|
|
call locate('prtV',nprtV,mlt,NUMPARTS)
|
|
call locate('prtM',nprtM,mlt,NUMPARTS)
|
|
call locate('prtA',nprtA,mlt,NUMPARTS)
|
|
call locate('nply',nnply,mlt,no_poly)
|
|
call locate('nump',nnump,mlt,mlt)
|
|
call locate('no_p',nno_p,mlt,no_poly)
|
|
call locate('ppm ',nppm ,mlt,no_poly)
|
|
call locate('nopt',nnopt,no_poly,mlt)
|
|
call defini('ykup',nykup,1,nyy)
|
|
call defini('tabx',ntabx,1,NUMPARTS)
|
|
call defini('tabz',ntabz,1,NUMPARTS)
|
|
call defini('tabm',ntabm,1,NUMPARTS)
|
|
call defini('tabp',ntabp,1,NUMPARTS)
|
|
call definr('slic',nslic,nxx+2,3*nzz+6)
|
|
nxx4=nxx
|
|
nyy4=nyy
|
|
nzz4=nzz
|
|
call defini('pv ',npv ,nxx4,nyy4*2)
|
|
call defini('el ',nel ,nyy4,nzz4*2)
|
|
call definr('toly',ntly ,nxx4,nyy4)
|
|
call definr('boly',nbly ,nxx4,nyy4)
|
|
call defini('tpv ',ntpv ,nxx4,nyy4)
|
|
call defini('bpv ',nbpv ,nxx4,nyy4)
|
|
|
|
CALL outputplumes(NUMPARTS,nxx,nyy,nzz,ia(nnop),ia(nityp),
|
|
* ia(ncord),ia(nprtE),ia(nprtN),ia(nprtZ),ia(nprtI),
|
|
* ia(nprtM),ia(nprtA),ia(nprtV),ia(nykup),ia(ntabx),
|
|
* ia(ntabz),ia(ntabm),ia(ntabp),ia(nslic),ia(npv),
|
|
* ia(nel),ia(ntly),ia(nbly),ia(ntpv),ia(nbpv),ia(nnply),
|
|
* ia(nnump),ia(nno_p),ia(nnopt),ia(nppm),no_poly,ngrid)
|
|
|
|
call delete('bpv ')
|
|
call delete('tpv ')
|
|
call delete('boly')
|
|
call delete('toly')
|
|
call delete('el ')
|
|
call delete('pv ')
|
|
call delete('slic')
|
|
call delete('tabp')
|
|
call delete('tabm')
|
|
call delete('tabz')
|
|
call delete('tabx')
|
|
call delete('ykup')
|
|
|
|
RETURN
|
|
END
|
|
|
|
c --------------------------------------------------------------
|
|
|
|
SUBROUTINE addoutputheader(stream)
|
|
INCLUDE '3dgeom.cb'
|
|
INCLUDE '3dpolls.cb'
|
|
INCLUDE '3dpoints.cb'
|
|
INCLUDE '3duparms.cb'
|
|
INTEGER i, n
|
|
BYTE bt
|
|
CHARACTER*1 bc
|
|
REAL gE, gN, a, b
|
|
INTEGER stream
|
|
|
|
n=1
|
|
DO 10 i=1,12
|
|
bc=start_stamp(i:i)
|
|
WRITE(stream,REC=n) bc
|
|
n=n+1
|
|
10 CONTINUE
|
|
i=outputTS
|
|
CALL out4bytes(i, n, stream)
|
|
bt=poll_type
|
|
WRITE(stream,REC=n) bt
|
|
n=18
|
|
i=NINT(abs(goE))
|
|
CALL out4bytes(i, n, stream)
|
|
i=NINT(abs(goN))
|
|
CALL out4bytes(i, n, stream)
|
|
a=numX
|
|
b=numY
|
|
CALL calc_EN(gE, gN, a, b)
|
|
i=NINT(abs(gE))
|
|
CALL out4bytes(i, n, stream)
|
|
i=NINT(abs(gN))
|
|
CALL out4bytes(i, n, stream)
|
|
i=NINT(numZ*dimZ)
|
|
CALL out4bytes(i, n, stream)
|
|
i=numX*resXY
|
|
CALL out4bytes(i, n, stream)
|
|
i=numY*resXY
|
|
CALL out4bytes(i, n, stream)
|
|
|
|
IF (stream.eq.21.or.stream.eq.22) THEN
|
|
i=numZ*resZ
|
|
ELSE
|
|
i=numX*resXY
|
|
END IF
|
|
CALL out4bytes(i, n, stream)
|
|
|
|
c add four bytes for negtive values of goE, goN, gE and gN
|
|
n=58
|
|
if(goE.lt.0.0) then
|
|
WRITE(stream,REC=n) 1
|
|
else
|
|
WRITE(stream,REC=n) 0
|
|
endif
|
|
n=n+1
|
|
if(goN.lt.0.0) then
|
|
WRITE(stream,REC=n) 1
|
|
else
|
|
WRITE(stream,REC=n) 0
|
|
endif
|
|
n=n+1
|
|
if(gE.lt.0.0) then
|
|
WRITE(stream,REC=n) 1
|
|
else
|
|
WRITE(stream,REC=n) 0
|
|
endif
|
|
n=n+1
|
|
if(gN.lt.0.0) then
|
|
WRITE(stream,REC=n) 1
|
|
else
|
|
WRITE(stream,REC=n) 0
|
|
endif
|
|
|
|
n=50
|
|
CALL out4bytes(numviews, n, stream)
|
|
c
|
|
c-------Leave some room for numviews and point table offset.
|
|
c
|
|
n=66+4*maxoutputs+1
|
|
WRITE(stream,REC=n) bt
|
|
numviews=0
|
|
c
|
|
c-------Extra line for recording the last record of the file
|
|
c
|
|
n1=62
|
|
call out4bytes(n,n1,stream)
|
|
|
|
RETURN
|
|
END
|
|
|
|
c ------------------------------------------------------------------
|
|
|
|
SUBROUTINE addpointers(stream)
|
|
INCLUDE '3dpoints.cb'
|
|
INTEGER ptoffset, n
|
|
INTEGER stream
|
|
cdrc byte b
|
|
byte b,c
|
|
|
|
c stream=21 is for the RWD file
|
|
c stream=22 is for the RWC file
|
|
c stream=23 is for the RWZ file
|
|
|
|
n=62
|
|
c print*,'Call 1, Stream =',stream
|
|
call in4bytes(n0,n,stream)
|
|
c print*,'After in4bytes, n0=',n0
|
|
read(stream,rec=n0) c
|
|
c print*,'After read'
|
|
INQUIRE(stream, NEXTREC=n0)
|
|
ptoffset=66
|
|
n=66
|
|
CALL out4bytes(viewpoint(stream-20), n+4*(numviews-1), stream)
|
|
n=50
|
|
CALL out4bytes(numviews, n, stream)
|
|
n=54
|
|
CALL out4bytes(ptoffset, n, stream)
|
|
read(stream, rec=n0-1) b
|
|
RETURN
|
|
END
|
|
|
|
c ------------------------------------------------------------------
|
|
c ------- Routines to prepare graphics views -----------------------
|
|
c ------------------------------------------------------------------
|
|
|
|
SUBROUTINE outputplumes(NUMPARTS,MAXXRES,MAXYRES,MAXZRES,nop,
|
|
* itype,cord,partE,partN,partZ,partI,partM,partA,partV,
|
|
* ylookup,table_x,table_z,table_M,table_point,slice,pv,el,
|
|
* toplayer,botlayer,top_pv,bot_pv,nply,nump,no_p,no_p_time,
|
|
* ppm,no_poly,ngrid)
|
|
INCLUDE '3duparms.cb'
|
|
INCLUDE '3dgeom.cb'
|
|
INCLUDE '3dpoints.cb'
|
|
|
|
INTEGER*2 nop(20,1)
|
|
INTEGER itype(1)
|
|
INTEGER partI(1)
|
|
INTEGER ylookup(MAXYRES)
|
|
INTEGER table_x(NUMPARTS)
|
|
INTEGER table_z(NUMPARTS)
|
|
REAL table_M(NUMPARTS)
|
|
INTEGER table_point(NUMPARTS)
|
|
INTEGER table_next
|
|
INTEGER point,ngrid
|
|
|
|
REAL partE(1),partN(1),partZ(1),partM(1),partV(1),partA(1)
|
|
REAL cord(3,1)
|
|
INTEGER nply(2,1),nump(1),no_p(1),ppm(1),no_p_time(no_poly,1)
|
|
REAL slice(0:MAXXRES, 0:MAXZRES, 3)
|
|
INTEGER left, centre, right, dummy
|
|
BYTE pv(MAXXRES,MAXYRES,2)
|
|
BYTE el(MAXYRES,MAXZRES,2)
|
|
|
|
INTEGER x, y, z, pil, no_poly
|
|
REAL mass
|
|
REAL a, b, c, pel,pnl,pzl,pe,pn,pz,pa,pm
|
|
INTEGER m, n1
|
|
REAL toplayer(MAXXRES, MAXYRES)
|
|
REAL botlayer(MAXXRES, MAXYRES)
|
|
BYTE top_pv(MAXXRES,MAXYRES)
|
|
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
|
|
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
|
|
pel=partE(m)
|
|
pnl=partN(m)
|
|
pzl=partZ(m)
|
|
pil=partI(m)
|
|
pa=partA(m)
|
|
pvv=partV(m)
|
|
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
|
|
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)
|
|
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)
|
|
|
|
RETURN
|
|
END
|
|
|
|
c ------------------------------------------------------------------
|
|
|
|
SUBROUTINE addyslice(xp,yp,zp,mass,slice,lf,cn,rt,MAXXRES,
|
|
* MAXZRES)
|
|
INCLUDE '3dgeom.cb'
|
|
INTEGER xp,yp,zp
|
|
REAL mass
|
|
REAL slice(0:MAXXRES,0:MAXZRES,3)
|
|
INTEGER lx,hx,lz,hz
|
|
INTEGER ymult, divider
|
|
INTEGER x,z
|
|
INTEGER lf, cn, rt
|
|
|
|
IF ((yp.eq.numY*resXY).or.(yp.eq.1)) THEN
|
|
ymult=2
|
|
ELSE
|
|
ymult=3
|
|
ENDIF
|
|
lx=-1
|
|
hx=1
|
|
lz=-1
|
|
hz=1
|
|
IF (xp.eq.1) lx=0
|
|
IF (xp.eq.numX*resXY) hx=0
|
|
IF (zp.eq.1) lz=0
|
|
IF (zp.eq.numZ*resZ) hz=0
|
|
divider=(hx-lx+1)*(hz-lz+1)*ymult
|
|
DO 400 x=xp+lx, xp+hx
|
|
DO 300 z=zp+lz,zp+hz
|
|
slice(x,z,cn)=slice(x,z,cn)+mass/divider
|
|
IF (slice(x,z,cn).gt.slice(0,z,cn))
|
|
1 slice(0,z,cn)=slice(x,z,cn)
|
|
IF (slice(x,z,cn).gt.slice(x,0,cn))
|
|
1 slice(x,0,cn)=slice(x,z,cn)
|
|
IF (yp.ne.1) THEN
|
|
slice(x,z,lf)=slice(x,z,lf)+mass/divider
|
|
IF (slice(x,z,lf).gt.slice(0,z,lf))
|
|
1 slice(0,z,lf)=slice(x,z,lf)
|
|
IF (slice(x,z,lf).gt.slice(x,0,lf))
|
|
1 slice(x,0,lf)=slice(x,z,lf)
|
|
ENDIF
|
|
IF (yp.ne.numY*resXY) THEN
|
|
slice(x,z,rt)=slice(x,z,rt)+mass/divider
|
|
IF (slice(x,z,rt).gt.slice(0,z,rt))
|
|
1 slice(0,z,rt)=slice(x,z,rt)
|
|
IF (slice(x,z,rt).gt.slice(x,0,rt))
|
|
1 slice(x,0,rt)=slice(x,z,rt)
|
|
ENDIF
|
|
300 CONTINUE
|
|
400 CONTINUE
|
|
RETURN
|
|
END
|
|
|
|
SUBROUTINE clearslice(slice, y, MAXXRES, MAXZRES)
|
|
INCLUDE '3dgeom.cb'
|
|
REAL slice(0:MAXXRES, 0:MAXZRES, 3)
|
|
INTEGER x,y,z
|
|
DO 20 x=0,numX*resXY
|
|
DO 10 z=0,numZ*resZ
|
|
slice(x,z,y)=0.0
|
|
10 CONTINUE
|
|
20 CONTINUE
|
|
RETURN
|
|
END
|
|
|
|
SUBROUTINE establishview(y, pv, el, slice, left, MAXXRES,
|
|
* MAXYRES,MAXZRES,ngrid)
|
|
INCLUDE '3dgeom.cb'
|
|
INCLUDE '3duparms.cb'
|
|
INTEGER y, left
|
|
BYTE dilution_scale
|
|
BYTE concentration_scale
|
|
REAL slice(0:MAXXRES, 0:MAXZRES, 3)
|
|
BYTE pv(MAXXRES,MAXYRES,2)
|
|
BYTE el(MAXYRES,MAXZRES,2)
|
|
INTEGER x,z,ngrid
|
|
REAL mass
|
|
if(dilution) then
|
|
DO x=1,numX*resXY
|
|
mass=slice(x,0,left)
|
|
pv(x,y,1)=dilution_scale(mass)
|
|
enddo
|
|
DO z=1,numZ*resZ
|
|
mass=slice(0,z,left)
|
|
el(y,z,1)=dilution_scale(mass)
|
|
enddo
|
|
endif
|
|
if(concentration) then
|
|
DO x=1,numX*resXY
|
|
mass=slice(x,0,left)
|
|
pv(x,y,2)=concentration_scale(mass)
|
|
enddo
|
|
DO z=1,numZ*resZ
|
|
mass=slice(0,z,left)
|
|
el(y,z,2)=concentration_scale(mass)
|
|
enddo
|
|
endif
|
|
IF(pbindump(ngrid))
|
|
1 CALL pbindumproutine(slice, left, MAXXRES, MAXZRES)
|
|
RETURN
|
|
END
|
|
|
|
|
|
c--------------------------------------------------------pbindumproutine
|
|
SUBROUTINE pbindumproutine(slice, left, MAXXRES, MAXZRES)
|
|
c----------------------------------------------------------------------c
|
|
c purpose: c
|
|
c To dump a 3D plume file. c
|
|
c----------------------------------------------------------------------c
|
|
INCLUDE '3duparms.cb'
|
|
INCLUDE '3dgeom.cb'
|
|
INCLUDE '3dpoints.cb'
|
|
INTEGER left
|
|
REAL slice(0:MAXXRES, 0:MAXZRES, 3)
|
|
INTEGER x, z, pos, count, recpos, n0, numv
|
|
REAL mass, out, sview, rn
|
|
|
|
c BMM 070810
|
|
c Change of the format of the bin file into a standard binary file.
|
|
c The file has been opened in 3drwalk.F on stream 37
|
|
|
|
WRITE(37) ((slice(x,z,left),x=1,numX*resXY),z=1,numZ*resZ)
|
|
|
|
c pos=1
|
|
c count=0
|
|
c read(20,REC=3) n0
|
|
cc read(20,REC=3) rn
|
|
cc n0=int(rn)
|
|
c read(20,REC=n0,err=1) mass
|
|
c INQUIRE(UNIT=20, NEXTREC=n0,err=2)
|
|
c recpos=n0
|
|
c DO 300 WHILE (pos.le.(numX*resXY*numZ*resZ))
|
|
c x=mod(pos-1, (numX*resXY))+1
|
|
c z=((pos-1)/(numX*resXY))+1
|
|
c mass=slice(x,z,left)
|
|
c IF (mass.eq.0.0) THEN
|
|
c count=count+1
|
|
c ELSE
|
|
c IF (count.ne.0) THEN
|
|
c out=(-1.0)*count
|
|
c WRITE(20, REC=recpos,err=3) out
|
|
c recpos=recpos+1
|
|
c count=0
|
|
c ENDIF
|
|
c WRITE(20, REC=recpos,err=4) mass
|
|
c recpos=recpos+1
|
|
c END IF
|
|
c pos=pos+1
|
|
c300 END DO
|
|
c IF (count.ne.0) THEN
|
|
c out=(-1.0)*count
|
|
c WRITE(20, REC=recpos,err=5) out
|
|
c recpos=recpos+1
|
|
c END IF
|
|
c read(20, REC=4,err=6) sview
|
|
c numv=int(sview)
|
|
c if(numv.ne.numviews) then
|
|
c write(20, REC=4) float(numviews)
|
|
c write(20, REC=4+numviews) float(n0)
|
|
c endif
|
|
c write(20, REC=3,err=7) recpos-1
|
|
cc write(20, REC=3,err=7) float(recpos-1)
|
|
c read(20, REC=recpos-1,err=8) mass
|
|
c RETURN
|
|
c1 print *,1,n0
|
|
c stop
|
|
c2 print *,2,n0
|
|
c stop
|
|
c3 print *,3,recpos,out
|
|
c stop
|
|
c4 print *,4,recpos,mass
|
|
c stop
|
|
c5 print *,5,recpos,out
|
|
c stop
|
|
c6 print *,6,4,sview
|
|
c stop
|
|
c7 print *,7,recpos-1
|
|
c stop
|
|
c8 print *,8,recpos-1,mass
|
|
c stop
|
|
|
|
|
|
END
|
|
|
|
c--------------------------------------------------------fbindumproutine
|
|
SUBROUTINE fbindumproutine(toplayer, MAXXRES, MAXYRES)
|
|
c----------------------------------------------------------------------c
|
|
c purpose: c
|
|
c To dump a floatable file. c
|
|
c----------------------------------------------------------------------c
|
|
INCLUDE '3duparms.cb'
|
|
INCLUDE '3dgeom.cb'
|
|
INCLUDE '3dpoints.cb'
|
|
REAL toplayer(MAXXRES, MAXYRES)
|
|
INTEGER x, y, pos, count, recpos, n0
|
|
REAL mass, out, rn
|
|
|
|
pos=1
|
|
count=0
|
|
read(19,REC=3) rn
|
|
n0=int(rn)
|
|
read(19,REC=n0) mass
|
|
INQUIRE(UNIT=19, NEXTREC=n0)
|
|
recpos=n0
|
|
DO 300 WHILE (pos.le.(numX*resXY*numY*resXY))
|
|
x=mod(pos-1, (numX*resXY))+1
|
|
y=((pos-1)/(numX*resXY))+1
|
|
mass=toplayer(x,y)
|
|
IF (mass.eq.0.0) THEN
|
|
count=count+1
|
|
ELSE
|
|
IF (count.ne.0) THEN
|
|
out=(-1.0)*count
|
|
WRITE(19, REC=recpos) out
|
|
recpos=recpos+1
|
|
count=0
|
|
ENDIF
|
|
WRITE(19, REC=recpos) mass
|
|
recpos=recpos+1
|
|
END IF
|
|
pos=pos+1
|
|
300 END DO
|
|
IF (count.ne.0) THEN
|
|
out=(-1.0)*count
|
|
WRITE(19, REC=recpos) out
|
|
recpos=recpos+1
|
|
END IF
|
|
write(19, REC=4) float(numviews)
|
|
write(19, REC=4+numviews) float(n0)
|
|
write(19, REC=3) float(recpos-1)
|
|
read(19, REC=recpos-1) mass
|
|
RETURN
|
|
END
|
|
|
|
|
|
c -----------------------------------------------------------------------------
|
|
|
|
FUNCTION dilution_scale(mass)
|
|
INCLUDE '3duparms.cb'
|
|
INCLUDE '3dgeom.cb'
|
|
REAL mass, dil
|
|
BYTE dilution_scale
|
|
|
|
IF (abs(mass).le.1.0e-6) THEN
|
|
dilution_scale=0
|
|
RETURN
|
|
END IF
|
|
dil=c_init*((dimXY/resXY)**2*(dimZ/resZ))/mass
|
|
IF (dil.lt.dlevel(1)) THEN
|
|
dilution_scale=1
|
|
ELSE
|
|
IF (dil.lt.dlevel(2)) THEN
|
|
dilution_scale=2
|
|
ELSE
|
|
IF (dil.lt.dlevel(3)) THEN
|
|
dilution_scale=3
|
|
ELSE
|
|
IF (dil.lt.dlevel(4)) THEN
|
|
dilution_scale=4
|
|
ELSE
|
|
IF (dil.lt.dlevel(5)) THEN
|
|
dilution_scale=5
|
|
ELSE
|
|
dilution_scale=6
|
|
END IF
|
|
END IF
|
|
END IF
|
|
END IF
|
|
END IF
|
|
RETURN
|
|
END
|
|
|
|
c -----------------------------------------------------------------------------
|
|
|
|
FUNCTION concentration_scale(mass)
|
|
INCLUDE '3duparms.cb'
|
|
INCLUDE '3dgeom.cb'
|
|
REAL mass, con
|
|
BYTE concentration_scale
|
|
|
|
con=1.0e-4*mass/((dimXY/resXY)**2*(dimZ/resZ))
|
|
IF (con.le.1.0e-5) THEN
|
|
concentration_scale=0
|
|
RETURN
|
|
END IF
|
|
IF (con.gt.clevel(5)) THEN
|
|
concentration_scale=1
|
|
ELSE
|
|
IF (con.gt.clevel(4)) THEN
|
|
concentration_scale=2
|
|
ELSE
|
|
IF (con.gt.clevel(3)) THEN
|
|
concentration_scale=3
|
|
ELSE
|
|
IF (con.gt.clevel(2)) THEN
|
|
concentration_scale=4
|
|
ELSE
|
|
IF (con.gt.clevel(1)) THEN
|
|
concentration_scale=5
|
|
ELSE
|
|
concentration_scale=6
|
|
END IF
|
|
END IF
|
|
END IF
|
|
END IF
|
|
END IF
|
|
RETURN
|
|
END
|
|
|
|
c -------------------------------------------------------------------
|
|
c ---------------- Encode Graphics Subroutine ----------------------
|
|
c ------------------------------------------------------------------
|
|
|
|
SUBROUTINE encodegraphics(stream,nc,pv,el,MAXXRES,MAXYRES,MAXZRES)
|
|
INCLUDE '3dgeom.cb'
|
|
INCLUDE '3dpoints.cb'
|
|
BYTE pv(MAXXRES,MAXYRES,2)
|
|
BYTE el(MAXYRES,MAXZRES,2)
|
|
INTEGER x, y, z, count, stream, nc, n0, n
|
|
BYTE c, c1
|
|
|
|
n=62
|
|
c print*,'Call 3, Stream =',stream
|
|
call in4bytes(n0,n,stream)
|
|
read(stream,rec=n0) c
|
|
INQUIRE(stream, NEXTREC=viewpoint(stream-20))
|
|
|
|
c Planview
|
|
x=1
|
|
y=1
|
|
c=pv(x,y,nc)
|
|
count=0 !should be 1 BUG
|
|
x=x+1
|
|
DO 200 WHILE (y.le.numY*resXY)
|
|
c1=pv(x,y,nc)
|
|
IF (c.eq.c1) THEN
|
|
count=count+1
|
|
ELSE
|
|
CALL grout(stream,c,count)
|
|
count=0 !should be 1 BUG
|
|
c=c1
|
|
END IF
|
|
x=x+1
|
|
IF (x.gt.numX*resXY) THEN
|
|
x=1
|
|
y=y+1
|
|
END IF
|
|
200 END DO
|
|
IF (count.ne.0) THEN
|
|
CALL grout(stream,c,count)
|
|
END IF
|
|
|
|
c Elevation.
|
|
y=1
|
|
z=1
|
|
c=el(y,z,nc)
|
|
count=0
|
|
z=z+1
|
|
DO 300 WHILE (y.le.numY*resXY)
|
|
c1=el(y,z,nc)
|
|
IF (c.eq.c1) THEN
|
|
count=count+1
|
|
ELSE
|
|
CALL grout(stream,c,count)
|
|
count=0
|
|
c=c1
|
|
END IF
|
|
z=z+1
|
|
IF (z.gt.numZ*resZ) THEN
|
|
z=1
|
|
y=y+1
|
|
END IF
|
|
300 END DO
|
|
IF (count.ne.0) THEN
|
|
CALL grout(stream,c,count)
|
|
END IF
|
|
INQUIRE(stream, NEXTREC=n0)
|
|
n=62
|
|
call out4bytes(n0-1,n,stream)
|
|
RETURN
|
|
END
|
|
|
|
c -------------------------------------------------------------------
|
|
c ----------- Encode Double Graphics Subroutine ---------------------
|
|
c -------------------------------------------------------------------
|
|
|
|
SUBROUTINE encodedoublegraphics(stream,tpv,bpv,MAXXRES,MAXYRES)
|
|
INCLUDE '3dgeom.cb'
|
|
INCLUDE '3dpoints.cb'
|
|
BYTE tpv(MAXXRES,MAXYRES)
|
|
BYTE bpv(MAXXRES,MAXYRES)
|
|
INTEGER x, y, count, stream, n0, n
|
|
BYTE c, c1
|
|
|
|
n=62
|
|
c print*,'Call 4, Stream =',stream
|
|
call in4bytes(n0,n,stream)
|
|
read(stream,rec=n0) c
|
|
INQUIRE(stream, NEXTREC=viewpoint(stream-20))
|
|
|
|
c Top Planview
|
|
x=1
|
|
y=1
|
|
c=tpv(x,y)
|
|
count=0
|
|
x=x+1
|
|
DO 200 WHILE (y.le.numY*resXY)
|
|
c1=tpv(x,y)
|
|
IF (c.eq.c1) THEN
|
|
count=count+1
|
|
ELSE
|
|
CALL grout(stream,c,count)
|
|
count=0
|
|
c=c1
|
|
END IF
|
|
x=x+1
|
|
IF (x.gt.numX*resXY) THEN
|
|
x=1
|
|
y=y+1
|
|
END IF
|
|
200 END DO
|
|
IF (count.ne.0) THEN
|
|
CALL grout(stream,c,count)
|
|
END IF
|
|
|
|
c Bottom Planview
|
|
x=1
|
|
y=1
|
|
c=bpv(x,y)
|
|
count=0
|
|
x=x+1
|
|
DO 300 WHILE (y.le.numY*resXY)
|
|
c1=bpv(x,y)
|
|
IF (c.eq.c1) THEN
|
|
count=count+1
|
|
ELSE
|
|
CALL grout(stream,c,count)
|
|
count=0
|
|
c=c1
|
|
END IF
|
|
x=x+1
|
|
IF (x.gt.numX*resXY) THEN
|
|
x=1
|
|
y=y+1
|
|
END IF
|
|
300 END DO
|
|
IF (count.ne.0) THEN
|
|
CALL grout(stream,c,count)
|
|
END IF
|
|
INQUIRE(stream, NEXTREC=n0)
|
|
n=62
|
|
call out4bytes(n0-1,n,stream)
|
|
RETURN
|
|
END
|
|
|
|
c -------------------------------------------------------------------
|
|
|
|
SUBROUTINE grout(stream,c,count)
|
|
BYTE c, b
|
|
INTEGER count, stream
|
|
INTEGER i1, i2, i3, i4, s, n
|
|
|
|
IF (count.le.15) THEN
|
|
s=128+8*count+c
|
|
b=s-128
|
|
INQUIRE(stream,NEXTREC=n)
|
|
WRITE(stream,REC=n) b
|
|
ELSE
|
|
i4=count/(256**3)
|
|
i3=(count-i4*256**3)/256**2
|
|
i2=(count-i4*256**3-i3*256**2)/256
|
|
i1=count-i4*256**3-i3*256**2-i2*256
|
|
s=c
|
|
if (i4.gt.0) s=s+64
|
|
if (i3.gt.0) s=s+32
|
|
if (i2.gt.0) s=s+16
|
|
if (i1.gt.0) s=s+8
|
|
INQUIRE(stream,NEXTREC=n)
|
|
b=s-128
|
|
WRITE(stream,REC=n) b
|
|
n=n+1
|
|
if (i4.gt.0) then
|
|
b=i4-128
|
|
WRITE(stream,REC=n) b
|
|
n=n+1
|
|
end if
|
|
if (i3.gt.0) then
|
|
b=i3-128
|
|
WRITE(stream,REC=n) b
|
|
n=n+1
|
|
end if
|
|
if (i2.gt.0) then
|
|
b=i2-128
|
|
WRITE(stream,REC=n) b
|
|
n=n+1
|
|
end if
|
|
if (i1.gt.0) then
|
|
b=i1-128
|
|
WRITE(stream,REC=n) b
|
|
n=n+1
|
|
end if
|
|
END IF
|
|
RETURN
|
|
END
|
|
|
|
c --------------------------------------------------------------------
|
|
c ---Subroutine for reading a 4 byte (pos) integer to stream 1--------
|
|
c --------------------------------------------------------------------
|
|
|
|
SUBROUTINE in4bytes(i,n,stream)
|
|
INTEGER i,n,stream
|
|
BYTE b
|
|
INTEGER i1, i2, i3, i4
|
|
|
|
READ(stream,REC=n) b
|
|
n=n+1
|
|
i1=b+128
|
|
READ(stream,REC=n) b
|
|
n=n+1
|
|
i2=b+128
|
|
READ(stream,REC=n) b
|
|
n=n+1
|
|
i3=b+128
|
|
READ(stream,REC=n) b
|
|
n=n+1
|
|
i4=b+128
|
|
i=(i4*256**3)+(i3*256**2)+(i2*256)+i1
|
|
RETURN
|
|
END
|
|
|
|
c --------------------------------------------------------------------
|
|
c ---Subroutine for writing a 4 byte (pos) integer to stream 1--------
|
|
c --------------------------------------------------------------------
|
|
|
|
SUBROUTINE out4bytes(i,n,stream)
|
|
INTEGER stream, i, n
|
|
BYTE b
|
|
INTEGER i1, i2, i3, i4
|
|
|
|
i4=i/(256**3)
|
|
i3=(i-i4*256**3)/256**2
|
|
i2=(i-i4*256**3-i3*256**2)/256
|
|
i1=i-i4*256**3-i3*256**2-i2*256
|
|
b=i1-128
|
|
WRITE(stream,REC=n) b
|
|
n=n+1
|
|
b=i2-128
|
|
WRITE(stream,REC=n) b
|
|
n=n+1
|
|
b=i3-128
|
|
WRITE(stream,REC=n) b
|
|
n=n+1
|
|
b=i4-128
|
|
WRITE(stream,REC=n) b
|
|
n=n+1
|
|
RETURN
|
|
END
|
|
|
|
|
|
|
|
c-----------------------------------------------------------particles_c
|
|
subroutine particles_c(cord,nply,nump,no_p,no_p_time,ppm,
|
|
* no_poly,pe,pn,pa,pm)
|
|
c---------------------------------------------------------------------c
|
|
c purpose: c
|
|
c To count particles in a convex polygon. c
|
|
c---------------------------------------------------------------------c
|
|
c input parameters: c
|
|
c no_poly -- number of polygons c
|
|
c nply(1,i) -- number of nodes in polygon i c
|
|
c nply(2,i) -- starting point of connectivity data c
|
|
c in array nump(k) c
|
|
c nump -- polygon connectivity data c
|
|
c cord(3,1) -- finite element nodal coordinates c
|
|
c pe,pn -- current location of the particle c
|
|
c pa -- life of the particle c
|
|
c pm -- mass of the particle c
|
|
c---------------------------------------------------------------------c
|
|
c output parameters: c
|
|
c no_p(i) -- number of particles in polygon i c
|
|
c ppm(i) -- total mass in polygon i c
|
|
c no_p_time(i,k) -- number of particles in polygon i and hour k c
|
|
c---------------------------------------------------------------------c
|
|
include '3duparms.cb'
|
|
real cord(3,1),ppm(1)
|
|
real pe,pn,pa,pm
|
|
integer nply(2,1),nump(1),no_p(1),no_p_time(no_poly,1)
|
|
real*8 area,ass
|
|
c
|
|
do i=1,no_poly
|
|
do j=1,nply(1,i)-1
|
|
i1=nply(2,i)+j-1
|
|
i2=nply(2,i)+j
|
|
x13=cord(1,nump(i1))-pe
|
|
y13=cord(2,nump(i1))-pn
|
|
x23=cord(1,nump(i2))-pe
|
|
y23=cord(2,nump(i2))-pn
|
|
ass=x13*y23
|
|
area=ass-x23*y13
|
|
if(area.lt.-0.1e-2*dabs(ass)) goto 20
|
|
enddo
|
|
no_p(i)=no_p(i)+1
|
|
ppm(i)=ppm(i)+pm
|
|
k=int((pa-dt)/3600)+1
|
|
no_p_time(i,k)=no_p_time(i,k)+1
|
|
20 enddo
|
|
return
|
|
end
|