update to a working version

source/3dcondits.for

@@ -1,5 +1,5 @@
 c----------------------------------------------------------getusrparams
-	SUBROUTINE getusrparams
+        SUBROUTINE getusrparams
 c---------------------------------------------------------------------c
 c purpose:                                                            c
 c         To read user parameters.                                    c
@@ -12,73 +12,69 @@ c---------------------------------------------------------------------c
         call getoutputflags
         call getconservparams
 c BMM 6/8/96
-	call getgridoutputs
+        call getgridoutputs
 
-	return
-	END
+        return
+        END
 
 c-----------------------------------------------------echoallparameters
-	SUBROUTINE echoallparameters
+        SUBROUTINE echoallparameters
 c---------------------------------------------------------------------c
 c purpose:                                                            c
 c         To display all user parameters.                             c
 c---------------------------------------------------------------------c
-	INCLUDE '3duparms.cb'
-	INCLUDE '3dgeom.cb'
+        INCLUDE '3duparms.cb'
+        INCLUDE '3dgeom.cb'
 
-	WRITE(*,*)
-	WRITE(*,*) '--------------------------------------------'
-	WRITE(*,*) 'User parameters are:-'
-	WRITE(*,*)
-	WRITE(*,'(2A15)') 'goE','goN'
-	WRITE(*,20) goE, goN
-	WRITE(*,*)
-	WRITE(*,'(2A15)') 'gaS','gaC'
-	WRITE(*,25) gaS, gaC
-	WRITE(*,*)
-	WRITE(*,'(3A15)') 'numX', 'numY', 'numZ'
-	WRITE(*,30) numX, numY, numZ
-	WRITE(*,*)
-	WRITE(*,'(2A15)') 'dimXY', 'dimZ'
-	WRITE(*,40) dimXY, dimZ
-	WRITE(*,*)
-	WRITE(*,'(2A15)') 'resXY', 'resZ'
-	WRITE(*,42) resXY, resZ
-	WRITE(*,*)
-	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)') 'Diffusivities: Element types not listed use the'
-     *			//' first'
-	WRITE(*,'(3A15)') 'Element type', 'HDIFF','VDIFF'
-	i=1
-	DO WHILE (EDIFF(i).ne.-1)
-		WRITE(*,'(I15,2F15.7)') EDIFF(i),HDIFF(i),VDIFF(i)
-		i=i+1
-	ENDDO
-	WRITE(*,*)
-	IF (dieoff) THEN
-		WRITE(*,'(A)') 'dieoff  = TRUE'
-	ELSE
-		WRITE(*,'(A)') 'dieoff  = FALSE'
-	ENDIF
-	IF (dieoffc) THEN
-		WRITE(*,'(A)') 'dieoffc = TRUE'
-	ELSE
-		WRITE(*,'(A)') 'dieoffc = FALSE'
-	ENDIF
-	WRITE(*,'(3A15)') 'T90_S','T90_D','TDPTH'
-	WRITE(*,'(3F15.3)') t90_s,t90_d,tdpth
-	WRITE(*,*)
-	WRITE(*,*) '--------------------------------------------'
-	WRITE(*,*)
+        WRITE(*,*)
+        WRITE(*,*) '--------------------------------------------'
+        WRITE(*,*) 'User parameters are:-'
+        WRITE(*,*)
+        WRITE(*,'(2A15)') 'goE','goN'
+        WRITE(*,20) goE, goN
+        WRITE(*,*)
+        WRITE(*,'(2A15)') 'gaS','gaC'
+        WRITE(*,25) gaS, gaC
+        WRITE(*,*)
+        WRITE(*,'(3A15)') 'numX', 'numY', 'numZ'
+        WRITE(*,30) numX, numY, numZ
+        WRITE(*,*)
+        WRITE(*,'(2A15)') 'dimXY', 'dimZ'
+        WRITE(*,40) dimXY, dimZ
+        WRITE(*,*)
+        WRITE(*,'(2A15)') 'resXY', 'resZ'
+        WRITE(*,42) resXY, resZ
+        WRITE(*,*)
+        WRITE(*,'(5A15)') 'newcondTS','outputTS','dt1','dt2',
+     1                    'maxoutputs'
+        WRITE(*,50) newcondTS, outputTS, dt,dtnew, maxoutputs
+        WRITE(*,*)
+        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(*,'(3A15)') 'Element type', 'HDIFF','VDIFF'
+        i=1
+        DO WHILE (EDIFF(i).ne.-1)
+               WRITE(*,'(I15,2F15.7)') EDIFF(i),HDIFF(i),VDIFF(i)
+               i=i+1
+        ENDDO
+        WRITE(*,*)
+        IF (dieoff) THEN
+                WRITE(*,'(A)') 'dieoff  = TRUE'
+        ELSE
+                WRITE(*,'(A)') 'dieoff  = FALSE'
+        ENDIF
+        IF (dieoffc) THEN
+                WRITE(*,'(A)') 'dieoffc = TRUE'
+        ELSE
+                WRITE(*,'(A)') 'dieoffc = FALSE'
+        ENDIF
+        WRITE(*,'(3A15)') 'T90_S','T90_D','TDPTH'
+        WRITE(*,'(3F15.3)') t90_s,t90_d,tdpth
+        WRITE(*,*)
+        WRITE(*,*) '--------------------------------------------'
+        WRITE(*,*)
 c
 c-------check whether the time steps meet the required relations
 c
@@ -87,51 +83,50 @@ c
      *     mod(outputTS,newcondTS).ne.0) goto 70
         if(mod(dt,dtnew).ne.0) goto 70
 
-20	FORMAT(2F15.3)
-25	FORMAT(2F15.9) 
-30	FORMAT(3I15)
-40	FORMAT(2F15.3)
-42	FORMAT(2I15)
-50	FORMAT(5I15)
-60	FORMAT(I15,2E15.3)
+20      FORMAT(2F15.3)
+25      FORMAT(2F15.9) 
+30      FORMAT(3I15)
+40      FORMAT(2F15.3)
+42      FORMAT(2I15)
+50      FORMAT(5I15)
+60      FORMAT(I15,2E15.3)
 
-	RETURN
+        RETURN
 
 70      write(*,*) '** Time Steps of newcondTS, outputTS and dt'//
      *             ' not as required **'
         stop
 
-	END
+        END
 
 c-----------------------------------------------------getpollconditions
-	SUBROUTINE getpollconditions(poll,massp,masspA,NUMPARTS,
+        SUBROUTINE getpollconditions(poll,massp,masspA,NUMPARTS,
      *                                            max_numps,key)
 c---------------------------------------------------------------------c
 c purpose:                                                            c
 c         To input a new set of pollution conditions.                 c
 c---------------------------------------------------------------------c
         INCLUDE '3duparms.cb'
-	INCLUDE '3dpolls.cb'
+        INCLUDE '3dpolls.cb'
 
-	REAL	poll(4,1)
+        REAL    poll(4,1)
 c		poll(1,i) = Vol of discharge at window i
 c		poll(2,i) = Conc of discharge at window i
 c		poll(3,i) = Depth of discharge at window i
 c		poll(4,i) = Radius of discharge at window i
 
-	INTEGER	i,d,w,max_numps,NUMPARTS,key
+        INTEGER i,d,w,max_numps,NUMPARTS,key
         REAL    massp(1),masspA(1),masspB,numd
 
 c  ***	Read in the new pollutant releases.
 
         key=0
         max_numps=NUMPARTS
-	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)
               IF (d.gt.num_pollwins.or.d.lt.1) THEN
-		PRINT *,'ERROR - Pollutant window not as '//
-     1					'expected.'
-		STOP
+                 PRINT *,'ERROR - Pollutant window'
+                 STOP
               ENDIF
               masspB=mass_pp
               numd=poll(1,d)*poll(2,d)*min(outputTS,newcondTS)
@@ -145,38 +140,40 @@ c  ***	Read in the new pollutant releases.
      *                  /massp(d)+25
               max_numps=max_numps+int(numd)
               masspA(d)=masspB
-60	CONTINUE
+60      CONTINUE
         key=1
-65	FORMAT(I2, F10.5, E10.1, 2F10.3)
-	RETURN
+65      FORMAT(I2, F10.5, E10.1, 2F10.3)
+        RETURN
 
-200	WRITE(*,*) 'Error in pollutants!'
-	STOP
-205	WRITE(*,*) 'Run out of pollutant data'
-	return
-	END
+200     WRITE(*,*) 'Error in pollutants!'
+        STOP
+205     WRITE(*,*) 'Run out of pollutant data'
+        return
+        END
 
 c-----------------------------------------------------getflowconditions
-	SUBROUTINE getflowconditions(nsurf,vel,cord,ao,wsel,wsel1,ic,
+        SUBROUTINE getflowconditions(nsurf,vel,cord,ao,wsel,wsel1,ic,
      &                               dhdt,nop,ielist,key)
 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'
-	INCLUDE '3dgeom.cb'
-	INCLUDE '3dpolls.cb'
+        INCLUDE '3dgeom.cb'
+        INCLUDE '3dpolls.cb'
         common /shape/shap(20),shpx(20),shpy(20),shpz(20)
 
-	REAL    vel(3,1)
+        REAL    vel(3,1)
 c		vel(1,i) = Velocity to E at node i
 c		vel(2,i) = Velocity to N at node i
 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,vel4,vel5,vel6,tet,dfct,wsel(1),wsel1(1)
+        INTEGER i,k,j,no3dn,ndg,no3de,key
+        REAL    vel1,vel2,depth,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
@@ -190,23 +187,57 @@ c  ***	Read in the velocity field
         key=0
 
         if(dimopt(1:2).eq.'10') then
+cBM2019 updating reads to OUTBNRMA format
           if(compact) then
-	   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)
+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
           else
-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),wsel(j),vel4,vel5,vel6,
-     2        VEL(3,j),J=1,No3dn),(DFCT,J=1,No3dE)
-	  print*,'Velocity file read at time',TET
+
+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,
+     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
+
+          print*,'Velocity file read at time',TET
           endif
         else
-           read(2,err=210,end=215) TET,no3dn,((vel(k,j),k=1,2),
-     1        wsel(j),j=1,no3dn),(vel4,j=1,no3dn)
+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
         endif
 c
-c-------for 2D case
-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
+
         if(no3dn.ne.np) then
            do i=1,np
               ic(i)=0
@@ -335,15 +366,19 @@ c
 c
 c----------calculate dhdt(i) for next time step
 c
-c          goto 40
+
+cBM2019
+           WRITE(*,*) 'Shouldnt be here'
+           STOP
+c
            if(dimopt(1:2).eq.'10') then
              if(compact) then
-	      read(2,err=30,end=30) TET,no3dn,NDG,no3de,
+              read(2,err=30,end=30) TET,no3dn,NDG,no3de,
      1           (vel1,vel2,wsel1(j),
      2           VEL1,J=1,No3dn)
              else
 cDRC 04/06/99 IYRR added for RMA10v6.5
-	      read(2,err=30,end=30) TET,no3dn,NDG,no3de,IYRR,
+              read(2,err=30,end=30) TET,no3dn,NDG,no3de,IYRR,
      1           (vel1,vel2,wsel1(j),vel4,vel5,vel6,
      2           VEL1,J=1,No3dn),(DFCT,J=1,No3dE)
              endif
@@ -362,19 +397,16 @@ 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
+        RETURN
 
-210	WRITE(*,*) ' *** Error in velocities! *** '
-	STOP
-215	WRITE(*,*) ' *** Run out of velocity data *** '
-	return
-	END
+210     WRITE(*,*) ' *** Error in velocities! *** '
+        STOP
+215     WRITE(*,*) ' *** Run out of velocity data *** '
+        return
+        END
 
 c-----------------------------------------------------------getgeometry
         subroutine getgeometry(cord,nop,nsurf,imat,width,ao,
@@ -387,21 +419,55 @@ c---------------------------------------------------------------------c
         include '3dmesh.cb'
         integer j,k
         integer nfix,nfix1,nref
-        integer*2 imat(1),ncorn,nfixh,iem,ndep,netyp
+        integer*2 imat(1),ncorn,iem,ndep,netyp
+        integer*4 nfixh  !BM2019 changed from integer*2
         real width(1),paxis,alfa,th
-        integer*2 nop(20,1),nsurf(1),ielist(1)
-        real cord(3,1),spec(3)
+
+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)
         real*8 ao(1)
-	real diffs(2,1)
+        real diffs(2,1)
 c
         if(dimopt(1:2).eq.'10') then
-           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),
+
+           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),
      2             (NDEP,NREF,J=1,NPM),
-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,
+     3             ((NOP4(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
@@ -446,6 +512,9 @@ 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
@@ -693,24 +762,24 @@ c       enddo
 c       stop
 c
 c Creating array for diffusivities.
-	k=1
-	DO WHILE(EDIFF(k).ne.-1)
-		k=k+1
-	ENDDO
-	DO i=1,ne
-		j=1
-		DO WHILE ((imat(i).ne.ediff(j)).and.(j.lt.k))
-			j=j+1
-		ENDDO
-		if (imat(i).eq.ediff(j)) then
-			diffs(1,i)=hdiff(j)
-			diffs(2,i)=vdiff(j)
-		else
-			diffs(1,i)=hdiff(1)
-			diffs(2,i)=vdiff(1)
-		endif
-	ENDDO
-		
+        k=1
+        DO WHILE(EDIFF(k).ne.-1)
+        	k=k+1
+        ENDDO
+        DO i=1,ne
+        	j=1
+        	DO WHILE ((imat(i).ne.ediff(j)).and.(j.lt.k))
+        		j=j+1
+        	ENDDO
+        	if (imat(i).eq.ediff(j)) then
+        		diffs(1,i)=hdiff(j)
+        		diffs(2,i)=vdiff(j)
+        	else
+        		diffs(1,i)=hdiff(1)
+        		diffs(2,i)=vdiff(1)
+        	endif
+        ENDDO
+        	
         return
 
    10   write(*,'(/a/)') ' *** Error in geometry file ***'
@@ -764,43 +833,43 @@ c
         end
 
 c-------------------------------------------------------------poll_skip
-	SUBROUTINE poll_skip(poll,key)
+        SUBROUTINE poll_skip(poll,key)
 c---------------------------------------------------------------------c
 c purpose:                                                            c
 c         To skip the first pollutant for hot start in unsteady case. c
 c---------------------------------------------------------------------c
         INCLUDE '3duparms.cb'
-	INCLUDE '3dpolls.cb'
+        INCLUDE '3dpolls.cb'
 
-	REAL	poll(4,1)
+        REAL	poll(4,1)
 c		poll(1,i) = Vol of discharge at window i
 c		poll(2,i) = Conc of discharge at window i
 c		poll(3,i) = Depth of discharge at window i
 c		poll(4,i) = Radius of discharge at window i
 
-	INTEGER	d, w, key
+        INTEGER	d, w, key
 
         key=0
 
 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)
               IF (d.gt.num_pollwins.or.d.lt.1) THEN
-		PRINT *,'ERROR - Pollutant window not as '//
+        	PRINT *,'ERROR - Pollutant window not as '//
      1					'expected.'
-		STOP
+        	STOP
               ENDIF
 60	CONTINUE
         key=1
 65	FORMAT(I2, F10.5, E10.1, 2F10.3)
-	RETURN
+        RETURN
 
 200	WRITE(*,*) 'Error in pollutants!'
-	STOP
+        STOP
 205	WRITE(*,*) 'Run out of pollutant data'
-	return
-	END
+        return
+        END
 
 c----------------------------------------------------------getopenbound
         subroutine getopenbound(no_open,opbd)
@@ -988,7 +1057,7 @@ c
        rewind nin
        call find('GRID OUTPUT',11,key)
        if(key.eq.0) then 
-		return
+        	return
        endif
 
        do i=1,50
@@ -1142,19 +1211,19 @@ cBMM additional lines can be read in for varying diffusivity ver 2.2
        call free
        call freer('Y',HDIFF(1),1)
        if (HDIFF(1).ge.0.0) then
-		EDIFF(1)=1
+        	EDIFF(1)=1
        		call free
        		call freer('Y',VDIFF,1)
        else
-		j=-1*int(HDIFF(1))
-		do i=1,j
-			call free 
-			call freei('E',EDIFF(i),1)
-			call free
+        	j=-1*int(HDIFF(1))
+        	do i=1,j
+        		call free 
+        		call freei('E',EDIFF(i),1)
+        		call free
        			call freer('Y',HDIFF(i),1)
-			call free
+        		call free
        			call freer('Y',VDIFF(i),1)
-		enddo
+        	enddo
        endif
 
        call free

source/3drwalk.for

@@ -28,8 +28,18 @@ c the same in the main block and the subroutines.
         include '3dpoints.cb'
         include 'rise_fall.cb'
 
-        parameter (ntotal=50000000)
+        parameter (ntotal=250000000)
         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
@@ -91,15 +101,24 @@ 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='sydney.3dg'
+           fname='bub003.3dg'
            call askstr(' Enter 3d geometry filename',fname)
-           open(9,file=fname,status='old',form='unformatted',err=10)
-           fname='sydney.res'
+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'
            call askstr(' Enter RMA-10 velocity filename',fname)
-           open(2,file=fname,status='old',form='unformatted',err=10)
+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
         else
            fname='sydney.geo'
            call askstr(' Enter 2D geometry filename',fname)
@@ -111,6 +130,7 @@ c Read in the user paramters
            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)
@@ -208,8 +228,12 @@ 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)
@@ -458,7 +482,7 @@ c
 c BMM 070810 Addition of another file type for storing the bin file data
 
              OPEN(37,FILE=plmnam(1:lenstr(plmnam))//'-block'//cext(1:2),
-     1		  STATUS='unknown',FORM='unformatted')
+     1               STATUS='unknown',FORM='unformatted')
 
              OPEN(20,FILE=plmnam(1:lenstr(plmnam))//cext(1:2),
 cdrc     1        ACCESS='direct',STATUS='new',FORM='unformatted',RECL=1)
@@ -475,7 +499,7 @@ cdrc     1        ACCESS='direct',STATUS='new',FORM='unformatted',RECL=1)
              WRITE(20,REC=4+maxoutputs) 0.0
              close(20)
 
-   	   ENDIF
+           ENDIF
 
            IF (fbindump) THEN
              inquire(file=fltnam(1:lenstr(fltnam))//cext(1:2),
@@ -679,9 +703,9 @@ c
                resZ=reszz(m)
 
                call outputgraphicscodes(m)
-               call addpointers(21)
-               call addpointers(22)
-               call addpointers(23)
+c              call addpointers(21)
+c              call addpointers(22)
+c              call addpointers(23)
                close(21)
                close(22)
                close(23)