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(*,42) resXY, resZ
WRITE(*,*)
WRITE(*,'(5A15)') 'newcondTS','outputTS','dt1','dt2',
1 'maxoutputs'
WRITE(*,'(5A15)') 'newcondTS','outputTS','dt1','dt2','maxoutputs'
WRITE(*,50) newcondTS, outputTS, dt,dtnew, maxoutputs
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(*,60) tot_dead_age, mass_pp, c_init
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'
i=1
DO WHILE (EDIFF(i).ne.-1)
@ -125,7 +129,8 @@ c *** Read in the new pollutant releases.
DO 60 w=1, num_pollwins
READ(4,*,ERR=200,end=205) d,(poll(i,d),i=1,4)
IF (d.gt.num_pollwins.or.d.lt.1) THEN
PRINT *,'ERROR - Pollutant window'
PRINT *,'ERROR - Pollutant window not as '//
1 'expected.'
STOP
ENDIF
masspB=mass_pp
@ -157,8 +162,6 @@ c-----------------------------------------------------getflowconditions
c---------------------------------------------------------------------c
c purpose: 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
INCLUDE '3duparms.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 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 cord(3,1),a(2,2)
real*8 ao(1),detj
@ -187,57 +190,23 @@ c *** Read in the velocity field
key=0
if(dimopt(1:2).eq.'10') then
cBM2019 updating reads to OUTBNRMA format
if(compact) then
c 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)
WRITE(*,*) 'No compact!'
STOP
read(2,err=210,end=215) TET,no3dn,NDG,no3de,
1 ((VEL(K,J),K=1,2),wsel(j),VEL(3,j),J=1,No3dn)
else
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)
cDRC 04/06/99 IYRR added for RMA10v6.5
read(2,err=210,end=215) TET,no3dn,NDG,no3de,IYRR,
1 ((VEL(K,J),K=1,2),depth,vel4,vel5,vel6,
2 VEL(3,j),wsel(j),J=1,no3dn),
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
1 ((VEL(K,J),K=1,2),wsel(j),vel4,vel5,vel6,
2 VEL(3,j),J=1,No3dn),(DFCT,J=1,No3dE)
print*,'Velocity file read at time',TET
endif
else
c 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)
WRITE(*,*) 'No RMA2!'
STOP
read(2,err=210,end=215) TET,no3dn,((vel(k,j),k=1,2),
1 wsel(j),j=1,no3dn),(vel4,j=1,no3dn)
endif
c
c-------For 2D case - START
c-------artificially create a 3D velocity field
c-------The 2D field could be from RMA2 or a depth averaged RMA10
c-------for 2D case
c
if(no3dn.ne.np) then
do i=1,np
ic(i)=0
@ -366,11 +335,7 @@ c
c
c----------calculate dhdt(i) for next time step
c
cBM2019
WRITE(*,*) 'Shouldnt be here'
STOP
c
c goto 40
if(dimopt(1:2).eq.'10') then
if(compact) then
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
endif
c-------For 2D case - END
40 key=1
c do i=1,No3dn
c do j=1,3
c vel(j,i)=0.0
c enddo
c enddo
RETURN
210 WRITE(*,*) ' *** Error in velocities! *** '
@ -419,55 +387,21 @@ c---------------------------------------------------------------------c
include '3dmesh.cb'
integer j,k
integer nfix,nfix1,nref
integer*2 imat(1),ncorn,iem,ndep,netyp
integer*4 nfixh !BM2019 changed from integer*2
integer*2 imat(1),ncorn,nfixh,iem,ndep,netyp
real width(1),paxis,alfa,th
cBM2019
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)
integer*2 nop(20,1),nsurf(1),ielist(1)
real cord(3,1),spec(3)
real*8 ao(1)
real diffs(2,1)
c
if(dimopt(1:2).eq.'10') then
read(9,err=10) NP,NE,NPM,NES,((CORD8(k,j),SPEC(k),K=1,3),
1 ALFA,NFIX,NFIX1,AO(J),NSURF4(J),J=1,NP),
read(9,err=10) NP,NE,NPM,NES,((CORD(k,j),SPEC(k),K=1,3),
1 ALFA,NFIX,NFIX1,AO(J),NSURF(J),J=1,NP),
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)
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
cBM2019
write(*,*) 'Stop: Geometry: No 2D'
stop
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)
do i=1,np
@ -512,9 +446,6 @@ c print *,i,(nop(k,i),k=1,nen)
c
c----------generate nodes below surface
c
cbm2019
write(*,*) 'Stop: no 2d!'
stop
np3=np
do n=1,np
cord(3,n)=0.0

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

@ -6,8 +6,6 @@ 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 --------------------------------------------------------------
@ -212,12 +210,6 @@ 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)
@ -247,19 +239,40 @@ cBMM2019
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
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
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)
pnl=partN(m)
pzl=partZ(m)
@ -269,198 +282,125 @@ cBM2019
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(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
IF (pvv.ne.0.0) THEN
IF (partV(m).eq.-9999.0) THEN
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
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)
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
@ -777,8 +717,7 @@ c -------------------------------------------------------------------
c ---------------- Encode Graphics Subroutine ----------------------
c ------------------------------------------------------------------
SUBROUTINE encodegraphics(stream,nc,pv,el,MAXXRES,MAXYRES,
1 MAXZRES)
SUBROUTINE encodegraphics(stream,nc,pv,el,MAXXRES,MAXYRES,MAXZRES)
INCLUDE '3dgeom.cb'
INCLUDE '3dpoints.cb'
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 'rise_fall.cb'
parameter (ntotal=250000000)
parameter (ntotal=50000000)
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
real r
integer t,NUMPARTS,max_numps,v_time
@ -101,24 +91,15 @@ c Read in the user paramters
fname='920820'
call askstr(' Enter output graphics filename',fname)
partl=fname
gridfname=fname
gridnumber=1
dimopt='10'
call askstr(' Enter RMA-2 or RMA-10 option',dimopt)
if(dimopt(1:2).eq.'10') then
fname='bub003.3dg'
fname='sydney.3dg'
call askstr(' Enter 3d geometry filename',fname)
cBM2019 update for OUTBN3GE format
c open(9,file=fname,status='old',form='unformatted',err=10)
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'
open(9,file=fname,status='old',form='unformatted',err=10)
fname='sydney.res'
call askstr(' Enter RMA-10 velocity filename',fname)
cBM2019 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
open(2,file=fname,status='old',form='unformatted',err=10)
else
fname='sydney.geo'
call askstr(' Enter 2D geometry filename',fname)
@ -130,7 +111,6 @@ cBM2019 update for OUTBNRMA format
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)
@ -228,12 +208,8 @@ c
write(*,'(/a/)') 'Loading the geometry of FE mesh ...'
rewind 9
cBM2019
read(9) header
read(9) np,ne
cBM2019
rewind 9
read(9) header
call defini('nop ',nnop ,10,ne+4)
call defini('ilst',nilst,1,ne/2+4)
call defini('imat',nimat,1,ne/2+4)
@ -703,9 +679,9 @@ c
resZ=reszz(m)
call outputgraphicscodes(m)
c call addpointers(21)
c call addpointers(22)
c call addpointers(23)
call addpointers(21)
call addpointers(22)
call addpointers(23)
close(21)
close(22)
close(23)

Write Preview
Loading…
Cancel
Save