From 4f6ce5a63d773d9a2fea1d1d59e087c7e75116e3 Mon Sep 17 00:00:00 2001 From: Brett Miller Date: Mon, 14 Oct 2019 16:03:42 +1100 Subject: [PATCH] update to a working version --- source/3doutput-block.for | 1307 +++++++++++++++++++------------------ 1 file changed, 684 insertions(+), 623 deletions(-) diff --git a/source/3doutput-block.for b/source/3doutput-block.for index 51cc5ef..feede24 100644 --- a/source/3doutput-block.for +++ b/source/3doutput-block.for @@ -6,17 +6,19 @@ 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 -------------------------------------------------------------- - SUBROUTINE outputgraphicscodes(ngrid) + SUBROUTINE outputgraphicscodes(ngrid) c--------------------------------------------------------------------c common mtot,npp,ia(1) common /dbsys/numa,next,idir,ipp(3) integer ngrid include '3dgeom.cb' - WRITE(*,*) 'Outputting...' + WRITE(*,*) 'Outputting...' nyy=numY*resXY+1 nxx=numX*resXY+1 nzz=numZ*resZ+1 @@ -46,22 +48,22 @@ c--------------------------------------------------------------------c 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 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), + 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('bpv ') + call delete('tpv ') + call delete('boly') + call delete('toly') call delete('el ') call delete('pv ') call delete('slic') @@ -71,57 +73,57 @@ c--------------------------------------------------------------------c call delete('tabx') call delete('ykup') - RETURN - END + RETURN + END c -------------------------------------------------------------- - SUBROUTINE addoutputheader(stream) - INCLUDE '3dgeom.cb' - INCLUDE '3dpolls.cb' - INCLUDE '3dpoints.cb' + 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) + 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 @@ -149,29 +151,29 @@ c add four bytes for negtive values of goE, goN, gE and gN WRITE(stream,REC=n) 0 endif - n=50 - CALL out4bytes(numviews, n, stream) + 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 + 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 + RETURN + END c ------------------------------------------------------------------ - SUBROUTINE addpointers(stream) - INCLUDE '3dpoints.cb' - INTEGER ptoffset, n - INTEGER stream + SUBROUTINE addpointers(stream) + INCLUDE '3dpoints.cb' + INTEGER ptoffset, n + INTEGER stream cdrc byte b byte b,c @@ -185,94 +187,79 @@ c print*,'Call 1, Stream =',stream c print*,'After in4bytes, n0=',n0 read(stream,rec=n0) c c print*,'After read' - INQUIRE(stream, NEXTREC=n0) - ptoffset=66 + 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) + 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 - + RETURN + END + c ------------------------------------------------------------------ c ------- Routines to prepare graphics views ----------------------- c ------------------------------------------------------------------ - SUBROUTINE outputplumes(NUMPARTS,MAXXRES,MAXYRES,MAXZRES,nop, + 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' - + 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) - 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 + 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 + 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 + +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 + end do + end do + DO 5 m=1, NUMPARTS pel=partE(m) pnl=partN(m) pzl=partZ(m) @@ -281,251 +268,324 @@ c IF (partE(m).ne.-1.) THEN 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 + CALL calc_ab(a, b, pe, pn) + CALL calc_c(c, pz) + 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*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 + 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 c ------------------------------------------------------------------ - SUBROUTINE addyslice(xp,yp,zp,mass,slice,lf,cn,rt,MAXXRES, + 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)) + 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)) + 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)) + 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)) + 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)) + 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)) + IF (slice(x,z,rt).gt.slice(x,0,rt)) 1 slice(x,0,rt)=slice(x,z,rt) - ENDIF + 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 + 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 + RETURN + END - SUBROUTINE establishview(y, pv, el, slice, left, MAXXRES, + 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 + 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 + DO x=1,numX*resXY mass=slice(x,0,left) - pv(x,y,1)=dilution_scale(mass) + pv(x,y,1)=dilution_scale(mass) enddo - DO z=1,numZ*resZ + DO z=1,numZ*resZ mass=slice(0,z,left) - el(y,z,1)=dilution_scale(mass) + el(y,z,1)=dilution_scale(mass) enddo endif if(concentration) then - DO x=1,numX*resXY + DO x=1,numX*resXY mass=slice(x,0,left) - pv(x,y,2)=concentration_scale(mass) + pv(x,y,2)=concentration_scale(mass) enddo - DO z=1,numZ*resZ + DO z=1,numZ*resZ mass=slice(0,z,left) - el(y,z,2)=concentration_scale(mass) + el(y,z,2)=concentration_scale(mass) enddo endif IF(pbindump(ngrid)) 1 CALL pbindumproutine(slice, left, MAXXRES, MAXZRES) - RETURN - END + RETURN + END - + c--------------------------------------------------------pbindumproutine - SUBROUTINE pbindumproutine(slice, left, MAXXRES, MAXZRES) + 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 + 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) + WRITE(37) ((slice(x,z,left),x=1,numX*resXY),z=1,numZ*resZ) c pos=1 c count=0 @@ -586,328 +646,329 @@ c8 print *,8,recpos-1,mass c stop - END + END c--------------------------------------------------------fbindumproutine - SUBROUTINE fbindumproutine(toplayer, MAXXRES, MAXYRES) + 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 + 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 + 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 + 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 + 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 + 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 - + 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 + 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 + 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 - + 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 + SUBROUTINE encodegraphics(stream,nc,pv,el,MAXXRES,MAXYRES, + 1 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)) + 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 + 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 + 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 + 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) + 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 + 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 + 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)) + 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 + 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 + 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 + 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) + 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 + 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 + 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-------- @@ -938,29 +999,29 @@ 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 + 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