From dc2c438b85a4bbdb42a67fc83e86edc884a6e3c2 Mon Sep 17 00:00:00 2001 From: Brett Miller Date: Mon, 14 Oct 2019 16:03:27 +1100 Subject: [PATCH] update to a working version --- source/3dcondits.for | 395 +++++++++++++++++++++++++------------------ 1 file changed, 232 insertions(+), 163 deletions(-) diff --git a/source/3dcondits.for b/source/3dcondits.for index 68d950b..9b954fa 100644 --- a/source/3dcondits.for +++ b/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