diff --git a/source/3dcondits.for b/source/3dcondits.for new file mode 100644 index 0000000..68d950b --- /dev/null +++ b/source/3dcondits.for @@ -0,0 +1,1251 @@ +c----------------------------------------------------------getusrparams + SUBROUTINE getusrparams +c---------------------------------------------------------------------c +c purpose: c +c To read user parameters. c +c---------------------------------------------------------------------c + + call gettitle + call getgridparams + call getgeneralparams + call getanalysisoptn + call getoutputflags + call getconservparams +c BMM 6/8/96 + call getgridoutputs + + return + END + +c-----------------------------------------------------echoallparameters + SUBROUTINE echoallparameters +c---------------------------------------------------------------------c +c purpose: c +c To display all user parameters. c +c---------------------------------------------------------------------c + 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(*,*) +c +c-------check whether the time steps meet the required relations +c + if(mod(newcondTS,dt).ne.0.or.mod(outputTS,dt).ne.0) goto 70 + if(mod(newcondTS,outputTS).ne.0.and. + * 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) + + RETURN + +70 write(*,*) '** Time Steps of newcondTS, outputTS and dt'// + * ' not as required **' + stop + + END + +c-----------------------------------------------------getpollconditions + 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' + + 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 + 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 + 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 + ENDIF + masspB=mass_pp + numd=poll(1,d)*poll(2,d)*min(outputTS,newcondTS) + * /mass_pp + if(numd.gt.float(lnum_ps)) then + masspB=poll(1,d)*poll(2,d)* + * min(outputTS,newcondTS)/lnum_ps + endif + massp(d)=max(masspB,masspA(d)) + numd=poll(1,d)*poll(2,d)*min(outputTS,newcondTS) + * /massp(d)+25 + max_numps=max_numps+int(numd) + masspA(d)=masspB +60 CONTINUE + key=1 +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 + +c-----------------------------------------------------getflowconditions + 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---------------------------------------------------------------------c + INCLUDE '3duparms.cb' + INCLUDE '3dmesh.cb' + INCLUDE '3dgeom.cb' + INCLUDE '3dpolls.cb' + common /shape/shap(20),shpx(20),shpy(20),shpz(20) + + 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) + 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 + data prsm/0.,-1.,-1.,1.,0.,-1.,0.,1.,-1., + & 0.,-1.,1.,1.,0.,1.,0.,1.,1./ + data brck/-1.,-1.,-1.,1.,-1.,-1.,1.,1.,-1.,-1.,1.,-1., + & -1.,-1.,1.,1.,-1.,1.,1.,1.,1.,-1.,1.,1./ + +c *** Read in the velocity field + + key=0 + + if(dimopt(1:2).eq.'10') then + 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) + 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 + 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) + endif +c +c-------for 2D case +c + if(no3dn.ne.np) then + do i=1,np + ic(i)=0 + enddo + do ii=1,n1_2d + i=ielist(ii) + if(nop(20,i).ne.0) then + nen=20 + else + nen=15 + endif + do k=1,nen + vel(3,nop(k,i))=0.0 + enddo + enddo + do i=no3dn+1,np + wsel(i)=wsel(nsurf(i)) + enddo + if(restart) then + do i=1,np + dhdt(i)=(wsel(i)-wsel1(i))/float(newcondTS) + wsel0=wsel(i) + wsel(i)=wsel1(i) + wsel1(i)=wsel0 + enddo + endif + do i=no3dn+1,np + if(nsurf(i).eq.nsurf(i-1)) then + cord(3,i)=-wsel(nsurf(i))/2.0 + else + cord(3,i)=-wsel(nsurf(i)) + endif + do k=1,3 + vel(k,i)=vel(k,nsurf(i)) + enddo + enddo +c +c----------calculate vertical velocity at corner nodes +c + do ii=1,n1_2d + i=ielist(ii) + if(nop(6,i).eq.0) goto 10 + if(nop(20,i).ne.0) then + nen=20 + iad=12 + it=3 + do k=1,8 + xi(k)=brck(1,k) + eta(k)=brck(2,k) + zeta(k)=brck(3,k) + enddo + else + nen=15 + iad=9 + it=2 + do k=1,6 + xi(k)=prsm(1,k) + eta(k)=prsm(2,k) + zeta(k)=prsm(3,k) + enddo + endif + k=0 + do n=1,2*nen/5,2 + k=k+1 + call xn3x(it,nen,xi(k),eta(k),zeta(k)) +c +c----------------calculate the Jacobian matrix +c + do m=1,2 + a(1,m)=0.0 + a(2,m)=0.0 + dadx=0.0 + dady=0.0 + dhdx=0.0 + dhdy=0.0 + do j=1,nen + nod=nop(j,i) + a(1,m)=a(1,m)+shpx(j)*cord(m,nod) + a(2,m)=a(2,m)+shpy(j)*cord(m,nod) + dadx=dadx+shpx(j)*ao(nod) + dady=dady+shpy(j)*ao(nod) + dhdx=dhdx+shpx(j)*(wsel(nod)+ao(nod)) + dhdy=dhdy+shpy(j)*(wsel(nod)+ao(nod)) + enddo + enddo + detj=a(1,1)*a(2,2)-a(1,2)*a(2,1) + vel(3,nop(n,i))=vel(3,nop(n,i))+ + & vel(1,nop(n,i))*(dadx*a(2,2)-dady*a(1,2))/detj+ + & vel(2,nop(n,i))*(dady*a(1,1)-dadx*a(2,1))/detj + vel(3,nop(n+iad,i))=vel(3,nop(n+iad,i))+ + & vel(1,nop(n+iad,i))*(dhdx*a(2,2)-dhdy*a(1,2))/detj+ + & vel(2,nop(n+iad,i))*(dhdy*a(1,1)-dhdx*a(2,1))/detj+ + & dhdt(nop(n+iad,i)) + ic(nop(n,i))=ic(nop(n,i))+1 + ic(nop(n+iad,i))=ic(nop(n+iad,i))+1 + enddo + 10 enddo + do i=1,np + if(ic(i).gt.1) then + vel(3,i)=vel(3,i)/float(ic(i)) + endif + enddo +c +c----------interpolate the mid-side nodal values +c + do i=1,ne + if(nop(6,i).eq.0) goto 20 + if(nop(20,i).ne.0) then + nen=8 + iad=12 + else + nen=6 + iad=9 + endif + do n=2,nen,2 + vel(3,nop(n,i))=(vel(3,nop(n-1,i))+ + & vel(3,nop(mod(n+1,nen),i)))/2.0 + vel(3,nop(n+iad,i))=(vel(3,nop(n+iad-1,i))+ + & vel(3,nop(mod(n+1,nen)+iad,i)))/2.0 + enddo + do n=1,nen/2 + vel(3,nop(nen+n,i))=(vel(3,nop(2*n-1,i))+ + & vel(3,nop(2*n-1+iad,i)))/2.0 + enddo + 20 enddo +c +c----------calculate dhdt(i) for next time step +c +c goto 40 + if(dimopt(1:2).eq.'10') then + if(compact) then + 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, + 1 (vel1,vel2,wsel1(j),vel4,vel5,vel6, + 2 VEL1,J=1,No3dn),(DFCT,J=1,No3dE) + endif + else + read(2,err=30,end=30) TET,no3dn,(vel1,vel2,wsel1(j), + 1 j=1,no3dn),(vel1,j=1,no3dn) + endif + backspace 2 + do i=no3dn+1,np + wsel1(i)=wsel1(nsurf(i)) + enddo + do i=1,np + dhdt(i)=(wsel1(i)-wsel(i))/float(newcondTS) + enddo + + 30 continue + endif + + 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! *** ' + STOP +215 WRITE(*,*) ' *** Run out of velocity data *** ' + return + END + +c-----------------------------------------------------------getgeometry + subroutine getgeometry(cord,nop,nsurf,imat,width,ao, + * ielist,diffs) +c---------------------------------------------------------------------c +c purpose: c +c To read geometric data of FE mesh. c +c---------------------------------------------------------------------c + include '3duparms.cb' + include '3dmesh.cb' + integer j,k + integer nfix,nfix1,nref + integer*2 imat(1),ncorn,nfixh,iem,ndep,netyp + real width(1),paxis,alfa,th + 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,((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), +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) + else + 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 + ao(i)=cord(3,i) + cord(1,i)=cord(1,i)*xscale + cord(2,i)=cord(2,i)*yscale + enddo + endif + print *,'np,npm=',np,npm + print *,'ne,nes=',ne,nes + do i=1,NE + if(nop(20,i).ne.0) then + nen=20 + goto 11 + else if(nop(15,i).ne.0.and.nop(16,i).eq.0) then + nen=15 + goto 11 + else if(nop(13,i).ne.0.and.nop(14,i).eq.0) then + nen=13 + goto 11 + else if(nop(10,i).ne.0.and.nop(11,i).eq.0) then + nen=10 + goto 11 + else if(nop(8,i).ne.0.and.nop(9,i).eq.0) then + nen=8 + goto 11 + else if(nop(6,i).ne.0.and.nop(7,i).eq.0) then + nen=6 + goto 11 + else if(nop(3,i).ne.0.and.nop(4,i).eq.0) then + nen=3 + goto 11 + endif +c print *,i,imat(i),nen +c print *,i,(nop(k,i),k=1,20) +c 11 print *,i,imat(i),nen +c print *,i,(nop(k,i),k=1,nen) + 11 enddo + + n1_2d=0 + if(np.eq.npm.or.dimopt(1:1).eq.'2') then +c +c----------generate nodes below surface +c + np3=np + do n=1,np + cord(3,n)=0.0 + do i=1,ne + if(nop(6,i).eq.0) goto 15 !1D or 1D to 2D junction elements + if(nop(7,i).eq.0) then + nen=6 + else + nen=8 + endif + do k=1,nen-1,2 + if(n.eq.nop(k,i)) goto 20 !corner nodes + enddo + do k=2,nen,2 + if(n.eq.nop(k,i)) goto 30 !mid-side nodes + enddo + 15 enddo + 20 nn=2 + nsurf(n)=n + goto 40 + 30 nn=1 + nsurf(n)=-n + 40 do k=1,nn + np3=np3+1 + cord(1,np3)=cord(1,n) + cord(2,np3)=cord(2,n) + cord(3,np3)=ao(n)/float(k) + ao(np3)=ao(n) + nsurf(np3)=n + enddo + enddo +c +c----------form 3D element +c + do i=1,ne + if(nop(6,i).eq.0) goto 60 + if(nop(7,i).eq.0) then + nen=6 + iad=9 + else + nen=8 + iad=12 + endif + do k=1,nen + nop(k+iad,i)=nop(k,i) + enddo + do k=1,nen + do n=np+1,np3 + if(nsurf(n).eq.nop(k+iad,i)) goto 50 + enddo + write(*,*) ' ** Error in code **' + stop + 50 nop(k,i)=n + if(mod(k,2).eq.1) then + ik=nen+(k+1)/2 + nop(ik,i)=n+1 + endif + enddo + n1_2d=n1_2d+1 + ielist(n1_2d)=i +c print *,'nop=',i,(nop(k,i),k=1,nen+iad) + 60 enddo + npm=np + np=np3 + else +c +c----------check whether there are any 1D & 2D elements in the 3D mesh +c + np3=np + do i=1,ne + if(imat(i).lt.900) then + if(nop(10,i).ne.0) goto 65 + nen=0 + if(nop(8,i).ne.0.and.nop(9,i).eq.0) then + nen=8 + iad=12 + else if(nop(6,i).ne.0.and.nop(7,i).eq.0) then + nen=6 + iad=9 + else if(nop(3,i).ne.0.and.nop(4,i).eq.0) then + nen=3 +c Modify 22/1/96 ver 2.2 +c Works for meshes that don't have all element numbers. + else if(nop(1,i).eq.0.and.nop(2,i).eq.0) then + nen=0 + else + nen=5 + endif + if(nen.eq.0) goto 65 + n1_2d=n1_2d+1 + ielist(n1_2d)=i + if(nen.gt.5) then !2D elements + do k=1,nen + nop(k+iad,i)=nop(k,i) + nsurf(nop(k+iad,i))=nop(k,i) + cord(3,nop(k+iad,i))=cord(3,1) + ao(nop(k+iad,i))=ao(nop(k,i)) + enddo + else if(nen.le.5) then !1D or 1D to 2D junction elements + print *,'TT',i,nen,nop(1,i),nop(2,i) + iad=12 + aa=sqrt((cord(2,nop(3,i))-cord(2,nop(1,i)))* + & (cord(2,nop(3,i))-cord(2,nop(1,i)))+ + & (cord(1,nop(3,i))-cord(1,nop(1,i)))* + & (cord(1,nop(3,i))-cord(1,nop(1,i)))) + asin=-(cord(2,nop(3,i))-cord(2,nop(1,i)))/aa + acos=(cord(1,nop(3,i))-cord(1,nop(1,i)))/aa + wd=width(nop(1,i)) + x=cord(1,nop(1,i))+asin*wd/2.0 + y=cord(2,nop(1,i))+acos*wd/2.0 + z=cord(3,1) + ao0=ao(nop(1,i)) + call existnod(cord,ao,nsurf,nop,np,np3,x,y,z,ao0, + & 1,0,i,l) + nop(1+iad,i)=l + x=cord(1,nop(1,i))-asin*wd/2.0 + y=cord(2,nop(1,i))-acos*wd/2.0 + z=cord(3,1) + ao0=ao(nop(1,i)) + call existnod(cord,ao,nsurf,nop,np,np3,x,y,z,ao0, + & 1,0,i,l) + nop(3+iad,i)=l + nop(2+iad,i)=nop(1,i) + nsurf(nop(2+iad,i))=nop(1,i) + if(nen.eq.5) then !1D to 2D junction elements + nop(6+iad,i)=nop(3,i) + nsurf(nop(6+iad,i))=nop(3,i) + x13=cord(1,nop(4,i))-cord(1,nop(1,i)) + y13=cord(2,nop(4,i))-cord(2,nop(1,i)) + x23=cord(1,nop(5,i))-cord(1,nop(1,i)) + y23=cord(2,nop(5,i))-cord(2,nop(1,i)) + area=x13*y23-x23*y13 + if(area.gt.0) then + nop(5+iad,i)=nop(4,i) + nsurf(nop(5+iad,i))=nop(4,i) + nop(7+iad,i)=nop(5,i) + nsurf(nop(7+iad,i))=nop(5,i) + else + nop(5+iad,i)=nop(5,i) + nsurf(nop(5+iad,i))=nop(5,i) + nop(7+iad,i)=nop(4,i) + nsurf(nop(7+iad,i))=nop(4,i) + endif + else !1D elements + wd=width(nop(3,i)) + x=cord(1,nop(3,i))-asin*wd/2.0 + y=cord(2,nop(3,i))-acos*wd/2.0 + z=cord(3,1) + ao0=ao(nop(3,i)) + call existnod(cord,ao,nsurf,nop,np,np3,x,y,z, + & ao0,3,0,i,l) + nop(5+iad,i)=l + x=cord(1,nop(3,i))+asin*wd/2.0 + y=cord(2,nop(3,i))+acos*wd/2.0 + z=cord(3,1) + ao0=ao(nop(3,i)) + call existnod(cord,ao,nsurf,nop,np,np3,x,y,z, + & ao0,3,0,i,l) + nop(7+iad,i)=l + nop(6+iad,i)=nop(3,i) + nsurf(nop(6+iad,i))=nop(3,i) + endif + np3=np3+1 + nop(4+iad,i)=np3 + nsurf(nop(4+iad,i))=nop(2,i) + cord(1,np3)=(cord(1,nop(3+iad,i)) + & +cord(1,nop(5+iad,i)))/2.0 + cord(2,np3)=(cord(2,nop(3+iad,i)) + & +cord(2,nop(5+iad,i)))/2.0 + cord(3,np3)=(cord(3,nop(3+iad,i)) + & +cord(3,nop(5+iad,i)))/2.0 + ao(np3)=(ao(nop(3+iad,i))+ao(nop(5+iad,i)))/2.0 + np3=np3+1 + nop(8+iad,i)=np3 + nsurf(nop(8+iad,i))=nop(2,i) + cord(1,np3)=(cord(1,nop(1+iad,i)) + & +cord(1,nop(7+iad,i)))/2.0 + cord(2,np3)=(cord(2,nop(1+iad,i)) + & +cord(2,nop(7+iad,i)))/2.0 + cord(3,np3)=(cord(3,nop(1+iad,i)) + & +cord(3,nop(7+iad,i)))/2.0 + ao(np3)=(ao(nop(1+iad,i))+ao(nop(7+iad,i)))/2.0 + nen=8 + endif +c +c----------------generate nodes below surface +c + if(nen.gt.5) then + do k=1,nen + if(mod(k,2).eq.0) then + nn=1 + else + nn=2 + endif + do m=1,nn + x=cord(1,nop(k+iad,i)) + y=cord(2,nop(k+iad,i)) + if(m.eq.1) then + z=ao(nop(k+iad,i)) + else + z=(cord(3,nop(k+iad,i)) + & +ao(nop(k+iad,i)))/2.0 + endif + ao0=ao(nop(k+iad,i)) +c +c-------------------------check whether it exists +c + call existnod(cord,ao,nsurf,nop,np,np3,x,y,z, + & ao0,k,iad,i,l) + if(m.eq.1) then + nop(k,i)=l + else + nop(k/2+nen+1,i)=l + endif + enddo + enddo + endif + endif + 65 enddo + np=np3 + print *,'npm,np=',npm,np + endif +c print *,'After generation ....' +c do i=1,ne +c if(nop(20,i).ne.0) then +c nen=20 +c else if(nop(15,i).ne.0.and.nop(16,i).eq.0) then +c nen=15 +c else if(nop(13,i).ne.0.and.nop(14,i).eq.0) then +c nen=13 +c else if(nop(10,i).ne.0.and.nop(11,i).eq.0) then +c nen=10 +c else if(nop(8,i).ne.0.and.nop(9,i).eq.0) then +c nen=8 +c else if(nop(6,i).ne.0.and.nop(7,i).eq.0) then +c nen=6 +c else if(nop(3,i).ne.0.and.nop(4,i).eq.0) then +c nen=3 +c else +c nen=5 +c endif +c print *,i,imat(i),nen +c print *,i,(nop(k,i),k=1,nen) +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 + + return + + 10 write(*,'(/a/)') ' *** Error in geometry file ***' + stop + end + +c---------------------------------------------------------------existnod + subroutine existnod(cord,ao,nsurf,nop,np,np3,x,y,z,ao0,k,iad,i, + & l) +c----------------------------------------------------------------------c +c purpose: c +c To check whether the created node exists in the extension c +c list of nodes. c +c----------------------------------------------------------------------c + real cord(3,1) + real*8 ao(1) + integer*2 nsurf(1),nop(20,1) +c + do l=np+1,np3 + if(x.eq.cord(1,l).and.y.eq.cord(2,l). + & and.z.eq.cord(3,l)) return + enddo + np3=np3+1 + l=np3 + cord(1,np3)=x + cord(2,np3)=y + cord(3,np3)=z + ao(np3)=ao0 + nsurf(np3)=nop(k+iad,i) + return + end + +c-----------------------------------------------------------getpollwins + subroutine getpollwins(pollwin,masspA) +c---------------------------------------------------------------------c +c purpose: c +c To read pollutant windows. c +c---------------------------------------------------------------------c + INCLUDE '3duparms.cb' + include '3dpolls.cb' + real pollwin(4,1),masspA(1) +c + do i=1,num_pollwins + read(4,*,err=10) (pollwin(j,i),j=1,4) + masspA(i)=mass_pp + enddo + return + + 10 write(*,'(/a/)') ' *** Error in pollutant file ***' + stop + end + +c-------------------------------------------------------------poll_skip + 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' + + 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 + + key=0 + +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 not as '// + 1 'expected.' + STOP + ENDIF +60 CONTINUE + key=1 +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 + +c----------------------------------------------------------getopenbound + subroutine getopenbound(no_open,opbd) +c---------------------------------------------------------------------c +c purpose: c +c To read open boundaries. c +c---------------------------------------------------------------------c +c no_open -- number of open boundaries specified. c +c opbd(1,i) -- x-component of the first point c +c opbd(2,i) -- y-component of the first point c +c opbd(3,i) -- x-component of the second point c +c opbd(4,i) -- y-component of the second point c +c---------------------------------------------------------------------c + integer no_open + real opbd(4,1) +c + do i=1,no_open + read(7,*) (opbd(k,i),k=1,4) + enddo + return + end +c---------------------------------------------------------getareaparams + subroutine getareaparams(nply,nump,cord,np,no_poly,ii) +c---------------------------------------------------------------------c +c purpose: c +c To read parameters of areas in which particles are counted. c +c---------------------------------------------------------------------c + include '3duparms.cb' + dimension nump(1),nply(2,1),cord(3,1),xy(2,500) + character*1 line + common /cline/line(160) + common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 + * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 +c + npadd=np + ntmp=nin + nin=3 + rewind nin + ii=1 + no_poly=0 + call find('AREAS',5,key) + if(key.eq.0) then + p_count=.false. + return + endif + p_count=.true. + 10 call free + if(line(1).eq.' ') goto 20 + call freei('N',nns,1) + call freei('R',nump(ii),nns) + if(nump(ii).eq.0) then + call freer('Y',xy(1,1),2*nns) + do k=1,nns + npadd=npadd+1 + if(npadd-np.gt.500) goto 30 + nump(ii+k-1)=npadd + cord(1,npadd)=xy(1,k) + cord(2,npadd)=xy(2,k) + enddo + endif + no_poly=no_poly+1 + nply(1,no_poly)=nns + nply(2,no_poly)=ii + ii=ii+nns + goto 10 + 20 ii=ii-1 + nin=ntmp + print *,'nply,nump=',nply(1,1),nply(2,1),(nump(i),i=1,ii) + return + 30 write(*,'(a)') ' *** Total no of polygon nodes > 500 ***' + stop + end +c--------------------------------------------------------getoutputflags + subroutine getoutputflags +c---------------------------------------------------------------------c +c purpose: c +c To read output flags. c +c---------------------------------------------------------------------c + include '3duparms.cb' + character*1 line + character*4 flag + common /cline/line(160) + common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 + * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 + data dlevel/100.0,200.0,500.0,1000.0,2000.0/ + data clevel/100.0,500.0,1000.0,5000.0,20000.0/ +c + ntmp=nin + nin=3 + rewind nin + call find('OUTPUT',6,key) + if(key.eq.0) then + dilution=.true. + write(*,*) 'Dilution' + concentration=.true. + write(*,*) 'Concentration' + return + endif + 10 call free + if(line(1).eq.' ') goto 20 + if(line(1).eq.'P') then + call freeh('N',flag,4,1) + if(flag(1:4).eq.'CONC') then + dilution=.false. + write(*,*) 'No Dilution' + concentration=.true. + write(*,*) 'Concentration' + elseif(flag(1:4).eq.'BOTH') then + dilution=.true. + write(*,*) 'Dilution' + concentration=.true. + write(*,*) 'Concentration' + else + dilution=.true. + write(*,*) 'Dilution' + concentration=.false. + write(*,*) 'No Concentration' + endif + elseif(line(1).eq.'D') then + call freer('S',dlevel,5) + write(*,'(a,5f12.1)') 'Input Dilution Levels=',dlevel + elseif(line(1).eq.'C') then + call freer('S',clevel,5) + write(*,'(a,5f12.1)') 'Input Concentration Levels=',clevel + endif + goto 10 + 20 nin=ntmp + return + end +c------------------------------------------------------getconservparams + subroutine getconservparams +c---------------------------------------------------------------------c +c purpose: c +c To read conservative option & parameters. c +c---------------------------------------------------------------------c + include '3duparms.cb' + character*1 line + character*4 flag + common /cline/line(160) + common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 + * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 +c + ntmp=nin + nin=3 + rewind nin + call find('CONSERVATION',12,key) + if(key.eq.0) then + dieoff=.false. + dieoffc=.false. + return + endif + call free + call freeh('F',flag,1,1) + if(flag(1:1).eq.'Y') then + dieoff=.true. + dieoffc=.false. + else if(flag(1:1).eq.'C') then + dieoff=.true. + dieoffc=.true. + endif + call freer('W',t90_s,1) + call freer('P',t90_d,1) + call freer('H',tdpth,1) + if(tdpth.eq.0.0) tdpth=4.0 + nin=ntmp + return + end + +c---------------------------------------------------------getgridouputs + subroutine getgridoutputs +c---------------------------------------------------------------------c +c purpose: c +c To read which grids should have bindump files c +c added BMM 6/8/96 c +c---------------------------------------------------------------------c + include '3duparms.cb' + character*1 line + character*4 flag + common /cline/line(160) + common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 + * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 +c + ntmp=nin + nin=3 + rewind nin + call find('GRID OUTPUT',11,key) + if(key.eq.0) then + return + endif + + do i=1,50 + pbindump(i)=.FALSE. + enddo +10 call free + if(line(1).eq.' ') goto 90 + call freei('P',n,1) + pbindump(n)=.TRUE. + goto 10 + +90 nin=ntmp + return + end + + +c--------------------------------------------------------------gettitle + subroutine gettitle +c---------------------------------------------------------------------c +c purpose: c +c To read the title of the problem. c +c---------------------------------------------------------------------c + include '3duparms.cb' + CHARACTER*80 dummstr + character*1 line + common /cline/line(160) + common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 + * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 +c + ntmp=nin + nin=3 + rewind nin + call find('TITLE',5,key) + if(key.eq.0) return + call free + call freeh(' ',dummstr,1,80) + print *,dummstr(1:lenstr(dummstr)) + nin=ntmp + return + end +c---------------------------------------------------------getgridparams + subroutine getgridparams +c---------------------------------------------------------------------c +c purpose: c +c To read grid parameters. c +c---------------------------------------------------------------------c + include '3duparms.cb' + INCLUDE '3dgeom.cb' + character*1 line + integer numxx(50),numyy(50),numzz(50),resxyy(50),reszz(50) + real goee(50),gonn(50),gass(50),gacc(50),dimxyy(50),dimzz(50) + common /mgird/goee,gonn,gass,gacc,numxx,numyy,numzz,dimxyy, + * dimzz,resxyy,reszz + common /cline/line(160) + common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 + * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 +c + ntmp=nin + nin=3 + rewind nin + no_grids=1 + call find('GRID PARAM',10,key) + if(key.eq.0) then + write(*,*) ' ** Grid parameters are not defined **' + stop + endif + 10 call free + if(line(1).eq.' ') goto 90 + call freei('R',n,1) + if(n.ne.0) then + call free + else + n=1 + endif + no_grids=max(no_grids,n) + call freer('E',goE,1) + call freer('N',goN,1) + call freer('S',gaS,1) + call freer('C',gaC,1) + call free + call freei('X',numX,1) + call freei('Y',numY,1) + call freei('Z',numZ,1) + call free + call freer('Y',dimXY,1) + call freer('Z',dimZ,1) + call free + call freei('Y',resXY,1) + call freei('Z',resZ,1) + print *,'goE,goN=',goE,goN + print *,'gaS,gaC=',gaS,gaC + print *,'numX,numY,numZ=',numX,numY,numZ + print *,'dimXY,dimZ=',dimXY,dimZ + print *,'resXY,resZ=',resXY,resZ + goee(n)=goE + gonn(n)=goN + gass(n)=gaS + gacc(n)=gaC + numxx(n)=numX + numyy(n)=numY + numzz(n)=numZ + dimxyy(n)=dimXY + dimzz(n)=dimZ + resxyy(n)=resXY + reszz(n)=resZ + goto 10 + 90 nin=ntmp + return + end +c------------------------------------------------------getgeneralparams + subroutine getgeneralparams +c---------------------------------------------------------------------c +c purpose: c +c To read general parameters. c +c---------------------------------------------------------------------c + include '3duparms.cb' + character*1 line + common /cline/line(160) + common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 + * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 +c + ntmp=nin + nin=3 + rewind nin + call find('GENERAL PARAM',13,key) + if(key.eq.0) then + write(*,*) ' ** General parameters are not defined **' + stop + endif + call free + call freei('S',maxoutputs,1) + call free + call freei('S',newcondTS,1) + call free + call freei('S',outputTS,1) + call free + call freei('T',dt,1) + if(dt.eq.0) then + call freei('1',dt,1) + call freei('2',dtnew,1) + else + dtnew=dt + endif + call free + call freei('E',tot_dead_age,1) + +cBMM additional lines can be read in for varying diffusivity ver 2.2 + do i=1,10 + EDIFF(i)=-1 + end do + call free + call freer('Y',HDIFF(1),1) + if (HDIFF(1).ge.0.0) then + 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 + call freer('Y',HDIFF(i),1) + call free + call freer('Y',VDIFF(i),1) + enddo + endif + + call free + call freer('E',mass_pp,1) + call freei('R',lnum_ps,1) + if(lnum_ps.eq.0) lnum_ps=2000000000 + call free + call freer('N',c_init,1) + nin=ntmp + return + end +c-------------------------------------------------------getanalysisoptn + subroutine getanalysisoptn +c---------------------------------------------------------------------c +c purpose: c +c To read analysis option. c +c---------------------------------------------------------------------c + include '3duparms.cb' + include 'rise_fall.cb' + character*1 line + character*40 flag + common /cline/line(160) + common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 + * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 +c + hsteadystate=.true. + psteadystate=.true. + plumeonly=.false. + floatableonly=.false. + euler=.true. + compact=.false. + ntmp=nin + nin=3 + rewind nin + call find('ANALYSIS OPT',12,key) + if(key.eq.0) then + return + endif + call free + if(line(1).eq.' ') goto 90 + call freeh('W',flag,1,8) + if(flag(1:8).eq.'UNSTEADY') then + hsteadystate=.false. + else if(flag(1:7).eq.'COMPACT') then + compact=.true. + hsteadystate=.false. + endif + call freeh('E',flag,1,8) + if(flag(1:8).eq.'UNSTEADY') then + psteadystate=.false. + endif + call freeh('S',flag,1,5) + if(flag(1:5).eq.'PLUME'.and..not.settles) then + plumeonly=.true. + rr=-1.0 + rf=-1.0 + elseif(flag(1:5).eq.'FLOAT'.and..not.settles) then + floatableonly=.true. + rr=100.0 + rf=-1.0 + y0r=y0r*100.0/5.0 + y1r=y1r*100.0/5.0 + endif + call free + if(line(1).eq.' ') goto 90 + call freeh('N',flag,1,5) + if(flag(1:5).eq.'SECAN') euler=.false. + if(.not.plumeonly) then + call free + if(line(1).eq.' ') goto 90 + if(.not.floatableonly) then + call freer('F',rf,1) + if(rf.eq.0.0) rf=-1.0 + call freer('A',y0f,1) + call freer('B',y1f,1) + endif + if(.not.settles) then + call freer('R',rr,1) + if(rr.eq.0.0) rr=-1.0 + call freer('C',y0r,1) + call freer('D',y1r,1) + endif + if(y1r.eq.0.0.or.y1f.eq.0.0) then + write(*,*) ' *** b and d cannot be zero. ***' + stop + endif + if(rr.gt.100.0.or.rf.gt.100.0) then + write(*,*) ' *** rf and rr having a range of 0 - 100 ***' + stop + endif + endif + 90 nin=ntmp + return + end diff --git a/source/3dgeneral.for b/source/3dgeneral.for new file mode 100644 index 0000000..1349892 --- /dev/null +++ b/source/3dgeneral.for @@ -0,0 +1,70 @@ +c---------------------------------------------------------------calc_ab + SUBROUTINE calc_ab(a, b, E, N) +c---------------------------------------------------------------------c +c purpose: c +c To rotate and translate coordinates from E-N to grid. c +c---------------------------------------------------------------------c + INCLUDE '3dgeom.cb' + REAL E, N, a, b + REAL nt, et + + nt= N-goN + et= E-goE + a = (nt*gaS+et*gaC)/dimXY + b = (nt*gaC-et*gaS)/dimXY + RETURN + END + +c---------------------------------------------------------------calc_EN + subroutine calc_EN(E,N,x,y) +c---------------------------------------------------------------------c +c Purpose: c +c To rotate and translate coordinates from grid to E-N. c +c---------------------------------------------------------------------c + INCLUDE '3dgeom.cb' + real E,N,x,y +c + N=goN+x*dimXY*gaS+y*dimXY*gaC + E=goE+x*dimXY*gaC-y*dimXY*gaS + return + end + +c----------------------------------------------------------------calc_c + SUBROUTINE calc_c(c, Z) +c---------------------------------------------------------------------c +c purpose: c +c To translate z-coordinates from E-N to grid. c +c---------------------------------------------------------------------c + INCLUDE '3dgeom.cb' + REAL c, Z + + c=-Z/dimZ + RETURN + END + +c--------------------------------------------------------------calc_XYZ + subroutine calc_XYZ(nop,cord,itype,pel,pnl,pzl,pil,pe,pn,pz) +c---------------------------------------------------------------------c +c purpose: c +c To interplate the global coordinates. c +c---------------------------------------------------------------------c + common /etype/inode(4) + common /shape/shap(20),shpx(20),shpy(20),shpz(20) + real pe,pn,pz,pel,pnl,pzl,cord(3,1),shap + integer*2 nop(20,1) + integer pil,it,nen,i,itype(1) +c + it=itype(pil) + nen=inode(it) + pe=0.0 + pn=0.0 + pz=0.0 + call xn3(it,nen,pel,pnl,pzl) + do i=1,nen + pe=pe+shap(i)*cord(1,nop(i,pil)) + pn=pn+shap(i)*cord(2,nop(i,pil)) + pz=pz+shap(i)*cord(3,nop(i,pil)) + enddo + return + end + diff --git a/source/3dgeom.cb b/source/3dgeom.cb new file mode 100644 index 0000000..b4557f7 --- /dev/null +++ b/source/3dgeom.cb @@ -0,0 +1,13 @@ +c GEOMETRY COMMON BLOCK. +c Used as part of the random walk model. +c BMM 17/5/93. + + REAL goN, goE, gaS, gaC + INTEGER numX, numY, numZ + REAL dimXY, dimZ + INTEGER resXY, resZ + + COMMON /GEOM/ bathym, goN, goE, gaS, gaC, numX, numY, numZ, + * dimXY, dimZ, resXY, resZ + + diff --git a/source/3dmesh.cb b/source/3dmesh.cb new file mode 100644 index 0000000..a14f536 --- /dev/null +++ b/source/3dmesh.cb @@ -0,0 +1,2 @@ + common /mesh/np,ne,npm,nes,n1_2d + integer np,ne,npm,nes,n1_2d diff --git a/source/3doutput-block.for b/source/3doutput-block.for new file mode 100644 index 0000000..51cc5ef --- /dev/null +++ b/source/3doutput-block.for @@ -0,0 +1,1014 @@ +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 diff --git a/source/3dpoints.cb b/source/3dpoints.cb new file mode 100644 index 0000000..64bac31 --- /dev/null +++ b/source/3dpoints.cb @@ -0,0 +1,10 @@ + +c POINTS COMMON BLOCK. +c Used as part of the random walk model. +c BMM 17/5/93. + + INTEGER numviews + INTEGER viewpoint(3) + + COMMON /POINT/ numviews, viewpoint + diff --git a/source/3dpolls.cb b/source/3dpolls.cb new file mode 100644 index 0000000..d601a01 --- /dev/null +++ b/source/3dpolls.cb @@ -0,0 +1,14 @@ +c POLLUTANTS COMMON BLOCK +c Used as part of the random walk model. + + INTEGER poll_type + INTEGER num_pollwins + + COMMON /POLLS/ poll_type, num_pollwins + +c Descriptions +c ============ +c poll_type : +c 1 - Plume particles (no die off) +c 2 - Floatable particles +c diff --git a/source/3drwalk.for b/source/3drwalk.for new file mode 100644 index 0000000..d5953c8 --- /dev/null +++ b/source/3drwalk.for @@ -0,0 +1,790 @@ +c----------------------------------------------------------------RWALK_3D + program RWALK_3D +c-----------------------------------------------------------------------c +c Purpose: c +c RWALK_3D - Random Walk Model c +c-----------------------------------------------------------------------c + include '3duparms.cb' + include '3dpolls.cb' + include '3dgeom.cb' + include '3dmesh.cb' + include '3dpoints.cb' + include 'rise_fall.cb' + parameter (ntotal=50000000) + character fname*100,buff*100,partl*100,fltnam*100 + character plmnam*100,cext*5 + real r + integer t,NUMPARTS,max_numps,v_time + logical flipflag,restout,existing + integer numxx(50),numyy(50),numzz(50),resxyy(50),reszz(50) + real goee(50),gonn(50),gass(50),gacc(50),dimxyy(50),dimzz(50) + common /mgird/goee,gonn,gass,gacc,numxx,numyy,numzz,dimxyy, + * dimzz,resxyy,reszz + common /dbsys/numa,next,idir,ipp(3) + common /etype/inode(4) + common /randm/isd, isd1 + common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10, + * nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 + common mtot,npp,ia(ntotal) + integer isd, isd1 + integer ADAY,AYEAR,AMONTH,ATIME,EGDS,EGDF + data inode/10,15,20,13/ + data isd/0/, isd1/0/ + data y0f/10.6/, y1f/-11.34/, y0r/1.5951/, y1r/-0.97525/ + data rf/30.0/, rr/5.0/ +c + write(*,'(a)') '-----------------------------------------------' + write(*,'(a)') '--- ---' + write(*,'(a)') '--- 3D_RWALK Version 2.2 ---' + write(*,'(a)') '--- 06/08/96 ---' + write(*,'(a)') '--- ---' + write(*,'(a)') '--- Y.C. Wang & B. Miller ---' + write(*,'(a)') '--- Water Research Lab. U.N.S.W. ---' + write(*,'(a)') '--- Compiled with ifort 12/7/07 ---' + write(*,'(a)') '--- ---' + write(*,'(a)') '-----------------------------------------------' + +c------------------------------------------------------------------ +c Use Defined variables to echo where and when this program was compiled +c added to code by AF 12/8/04 +c +c WRITE(*,'(2A)') 'This program was compiled in: ' +c 1 ,PATH +c WRITE(*,'(2A)') ' with filename : ' +c 1 ,__FILE__ +c WRITE(*,'(2A)') 'It was compiled on : ' +c 1 ,__DATE__ +c WRITE(*,'(2A)') ' at time : ' +c 1 ,__TIME__ +c +c IF (INDEX +c 1 (PATH,'/wrlnmr/') +c 2 .eq.0) THEN +c PRINT *,'This is not being run from the QA area' +c ELSE +c PRINT*, 'This is being run from the QA area' +c ENDIF +c +c print *,'Running the program' +c---------------------------------------------------------------- + mtot=ntotal + idir=mtot + call setcst +c +c-------open required files +c +c 21 : output stream for RWD graphics file +c 22 : output stream for RWC graphics file +c 23 : output stream for RWZ graphics file +c 9 : geometry file +c 2 : RMA10 velocity file +c 3 : input user parameters file +c 4 : input pollutants file +c 7 : input open-boundary file +c 8 : input hard-boundary file +c 10 : input restart file +c 11 : output restart file +c 19 : Floatable Bin Dump file +c 20 : Plume Bin Dump file +c 30 : settling particles file +c 35 : particle tracking path file +c 36 : particle counting file +c + WRITE(*,*) ' Enter start-stamp (yyyymmddhhmm): ' + READ(*,'(A12)') start_stamp + WRITE(*,'(A18,A12,A1)') ' START STAMP is "',start_stamp,'"' + WRITE(*,*) + + fname='920820' + call askstr(' Enter output graphics filename',fname) + partl=fname + dimopt='10' + call askstr(' Enter RMA-2 or RMA-10 option',dimopt) + if(dimopt(1:2).eq.'10') then + fname='sydney.3dg' + call askstr(' Enter 3d geometry filename',fname) + open(9,file=fname,status='old',form='unformatted',err=10) + fname='sydney.res' + call askstr(' Enter RMA-10 velocity filename',fname) + open(2,file=fname,status='old',form='unformatted',err=10) + else + fname='sydney.geo' + call askstr(' Enter 2D geometry filename',fname) + open(9,file=fname,status='old',form='unformatted',err=10) + buff='1.0,1.0' + call askstr(' Enter scales in x and y directions',buff) + read(buff,*) xscale,yscale + fname='sydney.res' + 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) + fname='polls' + call askstr(' Enter pollutants filename',fname) + open(4,file=fname,status='old',form='formatted',err=10) + fname='no' + call askstr(' Enter open boundary filename',fname) + if(fname(1:1).ne.'n'.and.fname(1:1).ne.'N') then + open(7,file=fname,status='old',form='formatted',err=10) + rewind 7 + openbd=.true. + else + openbd=.false. + endif + fname='no' + call askstr(' Enter input restart filename',fname) + if(fname(1:1).ne.'n'.and.fname(1:1).ne.'N') then + open(10,file=fname,status='unknown',form='unformatted',err=10) + restart=.true. + else + restart=.false. + endif + fname='no' + call askstr(' Enter output restart filename',fname) + if(fname(1:1).ne.'n'.and.fname(1:1).ne.'N') then + open(11,file=fname,status='unknown',form='unformatted',err=10) + restout=.true. + else + restout=.false. + endif + fname='no' + CALL askstr(' Enter Floatables Bin-Dump filename',fname) + fltnam=fname + if(fname(1:1).ne.'n'.and.fname(1:1).ne.'N') then + fbindump=.TRUE. + else + fbindump=.false. + endif + fname='no' + CALL askstr(' Enter Plumes Bin-Dump filename',fname) + plmnam=fname + if(fname(1:1).ne.'n'.and.fname(1:1).ne.'N') then + do itemp=1,50 + pbindump(itemp)=.TRUE. + enddo + else + do itemp=1,50 + pbindump(itemp)=.false. + enddo + endif +c print*,plmnam,'pbindump=',pbindump(1) + fname='no' + CALL askstr(' Enter settleable particles filename',fname) + if(fname(1:1).ne.'n'.and.fname(1:1).ne.'N') then + inquire(file=fname(1:lenstr(fname)),EXIST=existing) + if(existing) then + open(30,file=fname(1:lenstr(fname)),status='old') + close(30,status='delete') + endif + OPEN(30,FILE=fname, ACCESS='direct', + 1 STATUS='new',FORM='unformatted',RECL=8, ERR=10) + settles=.TRUE. + rr=-1.0 + rf=100.0 + y0f=y0f*100.0/30.0 + y1f=y1f*100.0/30.0 + else + settles=.false. + endif + fname='no' + CALL askstr(' Enter particle tracking path filename',fname) + if(fname(1:1).ne.'n'.and.fname(1:1).ne.'N') then + inquire(file=fname(1:lenstr(fname)),EXIST=existing) + if(existing) then + open(35,file=fname(1:lenstr(fname)),status='old') + close(35,status='delete') + endif +ctmp OPEN(35,FILE=fname, ACCESS='direct', +ctmp 1 STATUS='new',FORM='unformatted',RECL=8, ERR=10) + OPEN(35,FILE=fname, + 1 STATUS='new',FORM='formatted',ERR=10) + tracking=.TRUE. + else + tracking=.false. + endif + goto 20 + 10 len0=lnblnk(fname) + write(6,'(a)') '*** File: '''//fname(1:len0)//''' not found ***' + stop +c +c-------read data +c +20 write(*,'(/a/)') 'Loading the user parameters ...' + call getusrparams + print *,'plumeonly,floatableonly=',plumeonly,floatableonly + call echoallparameters + + write(*,'(/a/)') 'Loading the geometry of FE mesh ...' + rewind 9 + read(9) np,ne + rewind 9 + call defini('nop ',nnop ,10,ne+4) + call defini('ilst',nilst,1,ne/2+4) + call defini('imat',nimat,1,ne/2+4) + call definr('cord',ncord,3,7*np) !leave more space for 2D case + call definr('ao ',nao ,2,7*np+12) + call defini('nsur',nnsur,1,np*7/2+12) + call definr('widt',nwidt,1,np*7+12) +c BMM extra line for diffusivities + call definr('difs',ndifs,2,ne+4) + print *,'np=',np,npm,ne + call getgeometry(ia(ncord),ia(nnop),ia(nnsur),ia(nimat), + & ia(nwidt),ia(nao),ia(nilst),ia(ndifs)) + print *,'np=',np,npm,ne + call refinr('cord',ncord,3,np+800) !leave some space for polygon + call refinr('ao ',nao ,2,np+4) + call refini('nsur',nnsur,1,np/2+4) + call refinr('widt',nwidt,1,np+4) + call refinr('difs',ndifs,2,ne+4) +c +c-------calculate the day of year DAYOFY +c + if(dieoffc) then + READ(start_stamp,'(I4,2I2,I4)') ISY, ISM, ISD, IST + RMIN=-outputTS/60 + CALL TINC2(IST,ISD,ISM,ISY,ATIME,ADAY,AMONTH,AYEAR,RMIN,IER,0) + call gdate(1,AMONTH,ADAY,AYEAR,IDWK,IDYR,IDMON,EGDF,IER,0) + call gdate(1,1,1,AYEAR,IDWK,IDYR,IDMON,EGDS,IER,0) + DAYOFY=EGDF-EGDS+1 + TOFDAY=int(ATIME/100)+mod(ATIME,100)/60.0 + call DIEOFFS(0.0,DAYOFY,T90_D) +c print *,'T90_D=',T90_D + T90_D=3600.0*T90_D + endif + print *,'Euler=',euler + + call defini('nply',nnply,2,1000) + call defini('nump',nnump,1,10000) + call getareaparams(ia(nnply),ia(nnump),ia(ncord),np,no_poly, + & nlengthp) + print *,'no_poly,nlengthp=',no_poly,nlengthp + call refini('nump',nnump,1,nlengthp) + call refini('nply',nnply,2,no_poly) + call defini('no_p',nno_p,1,no_poly) + call definr('ppm ',nppm ,1,no_poly) +c +c-------create an array storing number of particles in specified area with +c time tag +c + mlp=tot_dead_age/3600+2 + call defini('nopt',nnopt,no_poly,mlp) + + if(p_count) then + fname=partl(1:lenstr(partl))//'.par' + CALL askstr(' Enter particle counting filename',fname) + inquire(file=fname(1:lenstr(fname)),EXIST=existing) + if(existing) then + open(36,file=fname(1:lenstr(fname)),status='old') + close(36,status='delete') + endif + OPEN(36,FILE=fname, + 1 STATUS='new',FORM='formatted', ERR=10) + + write(36,'(a)') ' Start_stamp='//start_stamp + write(36,'(a,i15)') ' outputTS=',outputTS + write(36,'(a,i15)') 'tot_dead_age=',tot_dead_age + write(36,'(a,i15)') 'No. of Areas=',no_poly + + endif + + write(*,'(/a/)') 'Loading the pollutant windows ...' + read(4,*) poll_type,num_pollwins + call definr('pwin',npwin,4,num_pollwins) + call definr('masA',nmasA,1,num_pollwins) + call getpollwins(ia(npwin),ia(nmasA)) +c +c-------divide FE elements into blocks +c + nx=2 + ny=5 + nz=1 + nxyz=nx*ny*nz + call defini('ityp',nityp,1,ne) + call definr('ebox',nebox,6,ne) + call definr('egrp',negrp,6,nxyz) + call defini('nbeg',nnbeg,1,nxyz) + call defini('nend',nnend,1,nxyz) + call defini('nelg',nnelg,3,ne) + call definr('egps',negps,6,2000) + call defini('nbb ',nnbb ,1,2000) + call defini('nee ',nnee ,1,2000) + call defini('nblk',nnblk,2,2000) + call defini('nelb',nnelb,1,1) + call blocks(nx,ny,nz,np,ne,ia(nnop),ia(ncord),nxyz,ia(negrp), + * ia(nnbeg),ia(nnend),ia(nnelg),ia(nnbb),ia(nnee),ia(nnelb), + * ia(nnblk),ia(negps),ia(nityp),ia(nebox),m,num) + call refini('nelb',nnelb,1,num) + call delete('nelg') + call delete('nbeg') + call delete('nend') + call refinr('egrp',negrp,6,nxyz) + call refinr('egps',negps,6,m) + call refini('nbb ',nnbb ,1,m) + call refini('nee ',nnee ,1,m) + call refini('nblk',nnblk,2,m) + call locate('nelb',nnelb,mlt,num) +c +c-------read open boundaries +c + if(openbd) then + read(7,*) no_open + call definr('opbd',nopbd,4,no_open) + call getopenbound(no_open,ia(nopbd)) + endif + call bcname +c +c-------setup boundary conditions +c + +c +c-------main program loop +c + write(*,*) 'Start the walking ...' + + no_settles=1 + NUMPARTS=0 + max_numpsA=0 + max_numpsB=0 + write(*,'(/a/)') 'Loading velocity field ...' + call definr('Avel',nAvel,3,np) + call definr('Bvel',nBvel,3,np) + call definr('vel ',nvel ,3,np) + call definr('dudt',ndudt,3,np) + call definr('dhdt',ndhdt,1,np) + call definr('wsel',nwsel,1,np) + call definr('wsl1',nwsl1,1,np) + call defini('ic ',nic ,1,np/2+4) + call definr('Apol',nApol,4,num_pollwins) + call definr('Bpol',nBpol,4,num_pollwins) + call definr('poll',npoll,4,num_pollwins) + call definr('pdot',npdot,1,num_pollwins) + call definr('masp',nmasp,1,num_pollwins) + call clearr(ia(npdot),1,num_pollwins) + call clearr(ia(ndhdt),1,np) +c +c-------initialise dudt(3,np) +c + if(.not.Euler) then + call clearr(ia(ndudt),3,np) + endif +c +c-------skip the first pollutant release for hot start +c + if(.not.psteadystate.and.restart) then + call getpollconditions(ia(nApol),ia(nmasp),ia(nmasA), + * NUMPARTS,max_numpsA,keyp) + if(keyp.eq.0) goto 50 + endif + + call getflowconditions(ia(nnsur),ia(nAvel),ia(ncord), + * ia(nao),ia(nwsel),ia(nwsl1),ia(nic), + * ia(ndhdt),ia(nnop),ia(nilst),keyf) + call getpollconditions(ia(nApol),ia(nmasp),ia(nmasA), + * NUMPARTS,max_numpsA,keyp) + if(keyf.eq.0.or.keyp.eq.0) goto 50 + if(.not.restart) then + if(.not.hsteadystate) then + call getflowconditions(ia(nnsur), + * ia(nBvel),ia(ncord),ia(nao),ia(nwsel),ia(nwsl1), + * ia(nic),ia(ndhdt),ia(nnop),ia(nilst),keyf) + if(keyf.eq.0) goto 50 + else + nBvel=nAvel + endif + if(.not.psteadystate) then + call getpollconditions(ia(nBpol),ia(nmasp),ia(nmasA), + * NUMPARTS,max_numpsB,keyp) + if(keyp.eq.0) goto 50 + else + nBpol=nApol + max_numpsB=max_numpsA + endif + endif + flipflag=.FALSE. + max_numps=max(max_numpsA,max_numpsB) + + print *,'no_grids=',no_grids + do m=1,no_grids + if(no_grids.eq.1) then + cext=' ' + else + write(cext,'(i2.2)') m + endif + inquire(file=partl(1:lenstr(partl))//'.rwd'//cext(1:2), + * EXIST=existing) + if(existing) then + open(21,file=partl(1:lenstr(partl))//'.rwd'//cext(1:2), + * status='old') + close(21,status='delete') + endif + open(21,file=partl(1:lenstr(partl))//'.rwd'//cext(1:2), + * access='direct',status='new',form='unformatted',recl=1) + inquire(file=partl(1:lenstr(partl))//'.rwc'//cext(1:2), + * EXIST=existing) + if(existing) then + open(22,file=partl(1:lenstr(partl))//'.rwc'//cext(1:2), + * status='old') + close(22,status='delete') + endif + open(22,file=partl(1:lenstr(partl))//'.rwc'//cext(1:2), + * access='direct',status='new',form='unformatted',recl=1) + inquire(file=partl(1:lenstr(partl))//'.rwz'//cext(1:2), + * EXIST=existing) + if(existing) then + open(23,file=partl(1:lenstr(partl))//'.rwz'//cext(1:2), + * status='old') + close(23,status='delete') + endif + open(23,file=partl(1:lenstr(partl))//'.rwz'//cext(1:2), + * access='direct',status='new',form='unformatted',recl=1) + goE=goee(m) + goN=gonn(m) + gaS=gass(m) + gaC=gacc(m) + numX=numxx(m) + numY=numyy(m) + numZ=numzz(m) + dimXY=dimxyy(m) + dimZ=dimzz(m) + resXY=resxyy(m) + resZ=reszz(m) + call addoutputheader(21) + call addoutputheader(22) + call addoutputheader(23) + close(21) + close(22) + close(23) + + IF (pbindump(m)) THEN + inquire(file=plmnam(1:lenstr(plmnam))//cext(1:2), + * EXIST=existing) + if(existing) then + open(20,file=plmnam(1:lenstr(plmnam))//cext(1:2), + * status='old') + close(20,status='delete') + endif + +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') + + OPEN(20,FILE=plmnam(1:lenstr(plmnam))//cext(1:2), +cdrc 1 ACCESS='direct',STATUS='new',FORM='unformatted',RECL=1) + 1 ACCESS='direct',STATUS='new',FORM='unformatted',RECL=4) + r=FLOAT(numX*resXY) + WRITE(20,REC=1) r + r=FLOAT(numY*resXY) + WRITE(20,REC=2) r + r=FLOAT(numZ*resZ) + WRITE(20,REC=3) r + r=float(4+maxoutputs) + WRITE(20,REC=3) int(r) + WRITE(20,REC=4) 0.0 + WRITE(20,REC=4+maxoutputs) 0.0 + close(20) + + ENDIF + + IF (fbindump) THEN + inquire(file=fltnam(1:lenstr(fltnam))//cext(1:2), + * EXIST=existing) + if(existing) then + open(19,file=fltnam(1:lenstr(fltnam))//cext(1:2), + * status='old') + close(19,status='delete') + endif + OPEN(19,FILE=fltnam(1:lenstr(fltnam))//cext(1:2), +cdrc 1 ACCESS='direct',STATUS='new',FORM='unformatted',RECL=1) + 1 ACCESS='direct',STATUS='new',FORM='unformatted',RECL=4) + r=FLOAT(numX*resXY) + WRITE(19,REC=1) r + r=FLOAT(numY*resXY) + WRITE(19,REC=2) r + r=float(4+maxoutputs) + WRITE(19,REC=3) r + WRITE(19,REC=4) 0.0 + WRITE(19,REC=4+maxoutputs) 0.0 + close(19) + ENDIF + enddo + + if(settles) write(30,rec=1) no_settles,mass_pp + + if(restart) then + rewind 10 + read(10) NUMPARTS + endif + + call defini('prtI',nprtI,1,NUMPARTS+1) + call definr('prtE',nprtE,1,NUMPARTS+1) + call definr('prtN',nprtN,1,NUMPARTS+1) + call definr('prtZ',nprtZ,1,NUMPARTS+1) + call definr('prtM',nprtM,1,NUMPARTS+1) + call definr('prtA',nprtA,1,NUMPARTS+1) + call definr('prtV',nprtV,1,NUMPARTS+1) + +c +c-------read data for restart +c + if(restart) then + call inputRS(NUMPARTS,max_numps,ia(nprtI),ia(nprtE), + * ia(nprtN),ia(nprtZ),ia(nprtM),ia(nprtA),ia(nprtV),ia(nBvel), + * ia(nwsl1),ia(nBpol),ia(nmasA)) + flipflag=.true. + endif + + v_time=0 !sum of dtnew until newcondTS + + t_p=0 !sum of time until termination + + do i=1,maxoutputs + + t=0 + +c +c----------estimate the maximum number of particles released +c + max_num=max_numps+NUMPARTS + call refini('prtI',nprtI,1,max_num) + call refinr('prtE',nprtE,1,max_num) + call refinr('prtN',nprtN,1,max_num) + call refinr('prtZ',nprtZ,1,max_num) + call refinr('prtM',nprtM,1,max_num) + call refinr('prtA',nprtA,1,max_num) + call refinr('prtV',nprtV,1,max_num) + +c +c----------unsteady case +c + if(.not.hsteadystate.or.(.not.psteadystate)) then + + 30 if(flipflag) then + call walkem(ia(ncord),ia(nnop),ia(nityp),ia(nebox), + * ia(nBvel),ia(nAvel),ia(nvel),ia(ndudt),ia(nao),nxyz, + * ia(negrp),ia(negps),ia(nnblk),ia(nnelb),ia(nnbb), + * ia(nnee),ia(nBpol),ia(nApol),ia(npoll),ia(npdot), + * ia(nmasp),ia(npwin),ia(nprtE),ia(nprtN),ia(nprtZ), + * ia(nprtI),ia(nprtM),ia(nprtA),ia(nprtV),ia(nopbd), + * no_open,NUMPARTS,t,t_p,v_time,ia(ndifs)) + if(v_time.eq.newcondTS) then + v_time=0 + if(.not.hsteadystate) then + call getflowconditions(ia(nnsur), + * ia(nBvel),ia(ncord),ia(nao),ia(nwsel),ia(nwsl1), + * ia(nic),ia(ndhdt),ia(nnop),ia(nilst),keyf) + endif + if(.not.psteadystate) then + call getpollconditions(ia(nBpol),ia(nmasp), + * ia(nmasA),NUMPARTS,max_numpsB,keyp) + endif + flipflag=.FALSE. + if(max_numps.lt.max_numpsB) then + max_num=max_num + max_numpsB - max_numps + max_numps=max_numpsB + call refini('prtI',nprtI,1,max_num) + call refinr('prtE',nprtE,1,max_num) + call refinr('prtN',nprtN,1,max_num) + call refinr('prtZ',nprtZ,1,max_num) + call refinr('prtM',nprtM,1,max_num) + call refinr('prtA',nprtA,1,max_num) + call refinr('prtV',nprtV,1,max_num) + endif + endif + else + call walkem(ia(ncord),ia(nnop),ia(nityp),ia(nebox), + * ia(nAvel),ia(nBvel),ia(nvel),ia(ndudt),ia(nao),nxyz, + * ia(negrp),ia(negps),ia(nnblk),ia(nnelb),ia(nnbb), + * ia(nnee),ia(nApol),ia(nBpol),ia(npoll),ia(npdot), + * ia(nmasp),ia(npwin),ia(nprtE),ia(nprtN),ia(nprtZ), + * ia(nprtI),ia(nprtM),ia(nprtA),ia(nprtV),ia(nopbd), + * no_open,NUMPARTS,t,t_p,v_time,ia(ndifs)) + if(v_time.eq.newcondTS) then + v_time=0 + if(.not.hsteadystate) then + call getflowconditions(ia(nnsur), + * ia(nAvel),ia(ncord),ia(nao),ia(nwsel),ia(nwsl1), + * ia(nic),ia(ndhdt),ia(nnop),ia(nilst),keyf) + endif + if(.not.psteadystate) then + call getpollconditions(ia(nApol),ia(nmasp), + * ia(nmasA),NUMPARTS,max_numpsA,keyp) + endif + flipflag=.TRUE. + if(max_numps.lt.max_numpsA) then + max_num=max_num + max_numpsA - max_numps + max_numps=max_numpsA + call refini('prtI',nprtI,1,max_num) + call refinr('prtE',nprtE,1,max_num) + call refinr('prtN',nprtN,1,max_num) + call refinr('prtZ',nprtZ,1,max_num) + call refinr('prtM',nprtM,1,max_num) + call refinr('prtA',nprtA,1,max_num) + call refinr('prtV',nprtV,1,max_num) + endif + endif + endif + + if(t.lt.outputTS) then + if(.not.hsteadystate.and.keyf.eq.0) goto 35 + if(.not.psteadystate.and.keyp.eq.0) goto 35 + goto 30 + endif + + +c +c----------steady-state case +c + else + call walkem(ia(ncord),ia(nnop),ia(nityp),ia(nebox), + * ia(nAvel),ia(nAvel),ia(nAvel),ia(ndudt),ia(nao),nxyz, + * ia(negrp),ia(negps),ia(nnblk),ia(nnelb),ia(nnbb), + * ia(nnee),ia(nApol),ia(nApol),ia(nApol),ia(npdot), + * ia(nmasp),ia(npwin),ia(nprtE),ia(nprtN),ia(nprtZ), + * ia(nprtI),ia(nprtM),ia(nprtA),ia(nprtV),ia(nopbd), + * no_open,NUMPARTS,t,t_p,v_time,ia(ndifs)) + endif +c +c----------output graphics +c + 35 call refini('prtI',nprtI,1,NUMPARTS) + call refinr('prtE',nprtE,1,NUMPARTS) + call refinr('prtN',nprtN,1,NUMPARTS) + call refinr('prtZ',nprtZ,1,NUMPARTS) + call refinr('prtM',nprtM,1,NUMPARTS) + call refinr('prtA',nprtA,1,NUMPARTS) + call refinr('prtV',nprtV,1,NUMPARTS) + if(.not.settles.and.t.eq.outputTS) then + do m=1,no_grids + if(no_grids.eq.1) then + cext=' ' + else + write(cext,'(i2.2)') m + endif + open(21,file=partl(1:lenstr(partl))//'.rwd'//cext(1:2), + * access='direct',status='old',form='unformatted',recl=1) + open(22,file=partl(1:lenstr(partl))//'.rwc'//cext(1:2), + * access='direct',status='old',form='unformatted',recl=1) + open(23,file=partl(1:lenstr(partl))//'.rwz'//cext(1:2), + * access='direct',status='old',form='unformatted',recl=1) + if(fbindump) then + OPEN(19,FILE=fltnam(1:lenstr(fltnam))//cext(1:2), +cdrc 1 ACCESS='direct',STATUS='old',FORM='unformatted',RECL=1) + 1 ACCESS='direct',STATUS='old',FORM='unformatted',RECL=4) + endif +c print*,plmnam,'pbindump=',pbindump(m) + if(pbindump(m)) then + OPEN(20,FILE=plmnam(1:lenstr(plmnam))//cext(1:2), +crdc 1 ACCESS='direct',STATUS='old',FORM='unformatted',RECL=1) + 1 ACCESS='direct',STATUS='old',FORM='unformatted',RECL=4) + endif + goE=goee(m) + goN=gonn(m) + gaS=gass(m) + gaC=gacc(m) + numX=numxx(m) + numY=numyy(m) + numZ=numzz(m) + dimXY=dimxyy(m) + dimZ=dimzz(m) + resXY=resxyy(m) + resZ=reszz(m) + call outputgraphicscodes(m) +c print*,'here' + call addpointers(21) +c print*,'here2' + call addpointers(22) + call addpointers(23) + close(21) + close(22) + close(23) + if(fbindump) close(19) + if(pbindump(m)) close(20) + enddo + endif + if(.not.hsteadystate.and.keyf.eq.0) goto 40 + if(.not.psteadystate.and.keyp.eq.0) goto 40 + + enddo + + 40 if(settles) write(30,rec=1) no_settles,mass_pp + call bcname +c +c-------output restart file +c + if(restout) then + if(hsteadystate.and.psteadystate) then + call outputRS(NUMPARTS,max_numps,ia(nprtI),ia(nprtE), + * ia(nprtN),ia(nprtZ),ia(nprtM),ia(nprtA),ia(nprtV),ia(nAvel), + * ia(nwsl1),ia(nApol),ia(nmasA)) + elseif(flipflag) then + call outputRS(NUMPARTS,max_numps,ia(nprtI),ia(nprtE), + * ia(nprtN),ia(nprtZ),ia(nprtM),ia(nprtA),ia(nprtV),ia(nAvel), + * ia(nwsl1),ia(nApol),ia(nmasA)) + else + call outputRS(NUMPARTS,max_numps,ia(nprtI),ia(nprtE), + * ia(nprtN),ia(nprtZ),ia(nprtM),ia(nprtA),ia(nprtV),ia(nBvel), + * ia(nwsl1),ia(nBpol),ia(nmasA)) + endif + endif + + 50 close(2) + close(3) + close(4) + close(9) + if(restart) close(10) + if(restout) close(11) + if(settles) close(30) + stop + end + +c----------------------------------------------------------------outputRS + subroutine outputRS(NUMPARTS,max_numps,partI,partE,partN,partZ, + * partM,partA,partV,vel,wsel1,poll,masspA) +c-----------------------------------------------------------------------c +c Purpose: c +c To output neccessary data for restart. c +c-----------------------------------------------------------------------c + include '3dpolls.cb' + include '3dmesh.cb' + integer NUMPARTS,max_numps,i,j + integer partI(1) + real partE(1),partN(1),partZ(1),partM(1),partA(1),partV(1) + real poll(4,1),vel(3,1),wsel1(1),masspA(1) +c + rewind 11 + write(11) NUMPARTS,max_numps,(partI(i),i=1,NUMPARTS),(partE(i), + 1 i=1,NUMPARTS),(partN(i),i=1,NUMPARTS),(partZ(i),i=1, + 2 NUMPARTS),(partM(i),i=1,NUMPARTS),(partA(i),i=1, + 3 NUMPARTS),(partV(i),i=1,NUMPARTS),((vel(j,i),j=1,3), + 4 wsel1(i),i=1,np),((poll(j,i),j=1,4),i=1,num_pollwins), + 5 (masspA(i),i=1,num_pollwins) + + return + end +c-----------------------------------------------------------------inputRS + subroutine inputRS(NUMPARTS,max_numps,partI,partE,partN,partZ, + * partM,partA,partV,vel,wsel1,poll,masspA) +c-----------------------------------------------------------------------c +c Purpose: c +c To input neccessary data for restart. c +c-----------------------------------------------------------------------c + include '3dpolls.cb' + include '3dmesh.cb' + integer NUMPARTS,max_numps,i,j + integer partI(1) + real partE(1),partN(1),partZ(1),partM(1),partA(1),partV(1) + real poll(4,1),vel(3,1),wsel1(1),masspA(1) +c + rewind 10 + read(10) NUMPARTS,max_numps,(partI(i),i=1,NUMPARTS),(partE(i), + 1 i=1,NUMPARTS),(partN(i),i=1,NUMPARTS),(partZ(i),i=1, + 2 NUMPARTS),(partM(i),i=1,NUMPARTS),(partA(i),i=1, + 3 NUMPARTS),(partV(i),i=1,NUMPARTS),((vel(j,i),j=1,3), + 4 wsel1(i),i=1,np),((poll(j,i),j=1,4),i=1,num_pollwins), + 5 (masspA(i),i=1,num_pollwins) + + return + end diff --git a/source/3duparms.cb b/source/3duparms.cb new file mode 100644 index 0000000..85456b7 --- /dev/null +++ b/source/3duparms.cb @@ -0,0 +1,35 @@ +c USER PARAMETERS COMMON BLOCK +c Used as part of the random walk model + + INTEGER newcondTS + INTEGER outputTS + INTEGER dt,dtnew + INTEGER tot_dead_age + INTEGER maxoutputs + INTEGER lnum_ps + INTEGER no_grids + REAL hdiff(1:10), vdiff(1:10) + INTEGER ediff(1:10) + REAL mass_pp,c_init + real xscale,yscale,clevel,dlevel + real t90_s,t90_d,tdpth,tofday,dayofy + character*12 start_stamp + character*5 dimopt + logical hsteadystate,psteadystate,restart,fbindump + logical pbindump(50),settles,openbd + logical p_count,dilution,concentration,dieoff,dieoffc + logical plumeonly,floatableonly + logical tracking,euler + logical compact + + COMMON /UPARMS/newcondTS,outputTS,dt,dtnew,tot_dead_age, + 1 maxoutputs,hdiff,vdiff,ediff,mass_pp,c_init,lnum_ps,no_grids, + 2 start_stamp,hsteadystate,psteadystate,restart,fbindump, + 3 pbindump,settles,openbd,p_count,dilution,concentration, + 4 plumeonly,floatableonly,tracking,euler,compact,dimopt + common /settle/no_settles + common /dimen2/xscale,yscale + common /conserv/t90_s,t90_d,tdpth,tofday,dayofy,dieoff,dieoffc + common /cdlvl/clevel(5),dlevel(5) + + diff --git a/source/3dwalking.for b/source/3dwalking.for new file mode 100644 index 0000000..2d7d1a1 --- /dev/null +++ b/source/3dwalking.for @@ -0,0 +1,772 @@ +c---------------------------------------------------------------walkem + subroutine walkem(cord,nop,itype,ebox,vel1,vel2,vel,dudt,ao, + * nxyz,egrp,egps,nblk,nelb,nbb,nee,poll1,poll2, + * poll,PDOT,massp,pollwin,partE,partN,partZ, + * partI,partM,partA,partV,opbd,no_open,NUMPARTS, + * t,t_p,v_time,diffs) +c---------------------------------------------------------------------c +c purpose: c +c To perform random walking. c +c---------------------------------------------------------------------c + include '3duparms.cb' + include '3dpolls.cb' + include '3dgeom.cb' + +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 + + common /shape/shap(20),shpx(20),shpy(20),shpz(20) + common /etype/inode(4) + common /randm/isd, isd1 + LOGICAL dead, walkjustone + INTEGER isd, isd1, p, pr, m, w, i, j, key, it, nen, num + INTEGER oldparts, newparts, deadparts, NUMPARTS, nxyz, pil + REAL r, ran3, tZ, tR, q, c, pel, pnl, pzl, pa, pv + REAL pe,pn,pz,xi,eta,zeta,shap + INTEGER t,t_p,v_time,no_open + INTEGER*2 nop(20,1) + INTEGER partI(1),nblk(2,1),nbb(1),nee(1),nelb(1),itype(1) + REAL partE(1),partN(1),partZ(1),partM(1),partA(1),vel(3,1) + REAL cord(3,1),egrp(6,1),egps(6,1),poll(4,1),pollwin(4,1) + REAL ebox(6,1),vel1(3,1),vel2(3,1),poll1(4,1),poll2(4,1) + REAL partV(1),opbd(4,1),dudt(3,1),PDOT(1),massp(1) + REAL TOFDAY1,DAYOFY1 + REAL*8 ao(1) + REAL diffs(2,1) + + oldparts=NUMPARTS + newparts=0 + deadparts=0 + num_steps=outputTS/dt + if(.not.hsteadystate.or.(.not.psteadystate)) then + num_steps=min(num_steps,newcondTS/dt) + endif + ibd=0 + + do i=1,num_steps +c +c----------calculate dieoff +c + if(dieoffc) then + TOFDAY=TOFDAY+dt/3600.0 + if(TOFDAY.ge.24.0) then + TOFDAY=TOFDAY-24.0 + DAYOFY=DAYOFY+1.0 + if(DAYOFY.gt.366.0) DAYOFY=DAYOFY-366.0 + endif + DELT=float(dt) + TOFDAY1=TOFDAY + DAYOFY1=DAYOFY + call DIEOFFS(TOFDAY1,DAYOFY1,T90_S) +cdrc + call DIEOFFS(0.0,DAYOFY1,T90_D) +c print *,'T90_S=',T90_S +crdc +c print *,'T90_D=',T90_D + T90_S=3600.0*T90_S +crdc + T90_D=3600.0*T90_D + endif + +c +c----------Temporal interpolation of velocities +c + if(.not.hsteadystate.or.(.not.psteadystate)) then + call interp_vels(vel1,vel2,vel,dudt,poll1, + * poll2,poll,v_time) + endif +c +c----------Old particles +c + m=1 + do while (m.le.NUMPARTS) +c +c---------------The particle is "active" and in the system. +c + pel=partE(m) + pnl=partN(m) + pzl=partZ(m) + pil=partI(m) + pa=partA(m) + pv=partV(m) + pm=partM(m) + it=itype(pil) + nen=inode(it) + call xn3(it,nen,pel,pnl,pzl) + pe=0.0 + pn=0.0 + pz=0.0 + do j=1,nen + pe=pe+shap(j)*cord(1,nop(j,pil)) + pn=pn+shap(j)*cord(2,nop(j,pil)) + pz=pz+shap(j)*cord(3,nop(j,pil)) + enddo + dead=walkjustone(pel,pnl,pzl,pil,pa,pv,pm,cord,itype, + * ebox,nop,vel,dudt,ao,egrp,egps,nblk,nelb,nbb, + * nee,nxyz,pe,pn,pz,no_open,opbd,dt,diffs,ibd) + partE(m)=pel + partN(m)=pnl + partZ(m)=pzl + partI(m)=pil + partA(m)=pa + partV(m)=pv + partM(m)=pm + + IF (dead.and.(.not.tracking)) THEN +c +c-----------------------particle is dead or out of system +c + NUMPARTS=NUMPARTS-1 + do j=m,NUMPARTS + partI(j)=partI(j+1) + partE(j)=partE(j+1) + partN(j)=partN(j+1) + partZ(j)=partZ(j+1) + partM(j)=partM(j+1) + partA(j)=partA(j+1) + partV(j)=partV(j+1) + enddo + deadparts=deadparts+1 + m=m-1 + END IF + m=m+1 + enddo +c +c----------New Particles. +c + NUMPARTS0=NUMPARTS + do ns=1,int(dt/dtnew) +c +c------------Temporal interpolation of velocities +c + if(ns.ne.1) then + if(.not.hsteadystate.or.(.not.psteadystate)) then + call interp_vels(vel1,vel2,vel,dudt,poll1, + * poll2,poll,v_time) + endif + endif + +ccc num_s=num_steps-ns+1 + +c +c------------Old particles related to time increment dtnew +c + m=1+NUMPARTS0 + do while (m.le.NUMPARTS) +c +c---------------The particle is "active" and in the system. +c + pel=partE(m) + pnl=partN(m) + pzl=partZ(m) + pil=partI(m) + pa=partA(m) + pv=partV(m) + pm=partM(m) + it=itype(pil) + nen=inode(it) + call xn3(it,nen,pel,pnl,pzl) + pe=0.0 + pn=0.0 + pz=0.0 + do j=1,nen + pe=pe+shap(j)*cord(1,nop(j,pil)) + pn=pn+shap(j)*cord(2,nop(j,pil)) + pz=pz+shap(j)*cord(3,nop(j,pil)) + enddo + dead=walkjustone(pel,pnl,pzl,pil,pa,pv,pm,cord,itype, + * ebox,nop,vel,dudt,ao,egrp,egps,nblk,nelb,nbb, + * nee,nxyz,pe,pn,pz,no_open,opbd,dtnew,diffs,ibd) + partE(m)=pel + partN(m)=pnl + partZ(m)=pzl + partI(m)=pil + partA(m)=pa + partV(m)=pv + partM(m)=pm + + IF (dead.and.(.not.tracking)) THEN +c +c-----------------------particle is dead or out of system +c + NUMPARTS=NUMPARTS-1 + do j=m,NUMPARTS + partI(j)=partI(j+1) + partE(j)=partE(j+1) + partN(j)=partN(j+1) + partZ(j)=partZ(j+1) + partM(j)=partM(j+1) + partA(j)=partA(j+1) + partV(j)=partV(j+1) + enddo + deadparts=deadparts+1 + m=m-1 + END IF + m=m+1 + enddo + + do w=1,num_pollwins +c +c-------------Establish the number of particles to be released. +c + q=poll(1,w) + c=poll(2,w) + PDOT(w)=PDOT(w)+(c*q*dtnew)/massp(w) + P=INT(PDOT(w)) + PDOT(w)=PDOT(w)-P +c p=100000 +c if(i.gt.1.or.ns.gt.1) p=0 + do pr=1,P + NUMPARTS=NUMPARTS+1 + m=NUMPARTS + r=ran3(isd) + pe=pollwin(1,w)+r*(pollwin(3,w)-pollwin(1,w)) + pn=pollwin(2,w)+r*(pollwin(4,w)-pollwin(2,w)) + partA(m)=0. + pm=massp(w) + tZ=-poll(3,w) + tR=poll(4,w) + r=2.0*(ran3(isd)-0.5) + pz=tZ+r*tR + pz_keep=pz + +c +c---------------Check if the particle is in the water column. +c + IF (pz.gt.0.0) pz=0.0 +c i.e. if it is above the water surface then put it on the surface + num=-1 + call trinvs3(nop,cord,itype,ebox,egrp,nblk,egps,nelb, + * nbb,nee,nxyz,nen,it,pe,pn,pz,xi,eta,zeta,num,key) + num_keep1=num + if(key.eq.1) then +c i.e. the pe,pn,pz coordinate was inside of element num + partE(m)=xi + partN(m)=eta + partZ(m)=zeta + partI(m)=num + else +c try and see if the plan form (pe,pn) is inside an element + num=-1 + call trinvs3(nop,cord,itype,ebox,egrp,nblk,egps, + * nelb,nbb,nee,nxyz,nen,it,pe,pn,0.0,xi,eta, + * zeta,num,key) + num_keep2=num + if(key.eq.1) then +c the (pe, pn) point was in element num, but we knew that (pe,pn,pz) wasn't so the +c old pz value must have been below the bed. + pz=0.0 + call xn3(it,nen,xi,eta,zeta) + do j=1,nen + pz=pz+shap(j)*ao(nop(j,num)) + enddo + pz=0.99*pz +c set the pz value to be at 99% of the depth and all is well. + else +c the (pe,pn) pair was not found to be in an element. +c This may happen if the (pe, pn) pair falls directly onto a element boundary. + write(*,'(a)') ' *** Window for discharging '// + * 'plumes out of the modelling domain (1)***' + write(*,'(i5,3f15.5)') 0,pe,pn,pz + stop + endif + num=-1 + call trinvs3(nop,cord,itype,ebox,egrp,nblk,egps, + * nelb,nbb,nee,nxyz,nen,it,pe,pn,pz,xi,eta, + * zeta,num,key) + if(key.eq.1) then + partE(m)=xi + partN(m)=eta + partZ(m)=zeta + partI(m)=num + else + write(*,'(a)') ' *** Window for discharging '// + * 'plumes out of the modelling domain (2)***' + write(*,'(i5,3f15.5)') 1,pe,pn,pz + write(*,*) 'orig pz, num_keep1, num_keep2',pz_keep + * ,num_keep1, num_keep2 + stop + endif + endif + + pel=partE(m) + pnl=partN(m) + pzl=partZ(m) + pil=partI(m) + pa=partA(m) + it=itype(pil) + nen=inode(it) + call xn3(it,nen,pel,pnl,pzl) + call rise_fall_vel(pv) +c following line for plume only + if(plumeonly) pv=0.0 + +c following few lines for floatables only + if(floatableonly) then + if(pv.le.0.0) then + NUMPARTS = NUMPARTS - 1 + goto 70 + endif + endif +c------------------------------------------ + + if(settles.and.pv.ge.0.0) then + NUMPARTS = NUMPARTS - 1 + goto 70 + endif + + dead=walkjustone(pel,pnl,pzl,pil,pa,pv,pm,cord,itype, + * ebox,nop,vel,dudt,ao,egrp,egps,nblk,nelb,nbb, + * nee,nxyz,pe,pn,pz,no_open,opbd,dtnew,diffs,ibd) + + newparts=newparts+1 + partE(m)=pel + partN(m)=pnl + partZ(m)=pzl + partI(m)=pil + partA(m)=pa + partV(m)=pv + partm(m)=pm + + IF (dead.and.(.not.tracking)) THEN + NUMPARTS=NUMPARTS-1 + deadparts=deadparts+1 + END IF + + 70 enddo + enddo + v_time=v_time+dtnew + + enddo +c--------------------- + +c +c----------output tracking path +c + if(tracking.and.t+dt.eq.outputTS) then +ctmp inquire(35,NEXTREC=ntrack) +ctmp write(35,rec=ntrack) NUMPARTS,t_p + write(35,*) NUMPARTS,t_p + ntrack=ntrack+1 + do m=1,NUMPARTS + pel=partE(m) + pnl=partN(m) + pzl=partZ(m) + pil=partI(m) + pa=partA(m) + pv=partV(m) + pm=partM(m) + it=itype(pil) + nen=inode(it) + call xn3(it,nen,pel,pnl,pzl) + pe=0.0 + pn=0.0 + pz=0.0 + do j=1,nen + pe=pe+shap(j)*cord(1,nop(j,pil)) + pn=pn+shap(j)*cord(2,nop(j,pil)) + pz=pz+shap(j)*cord(3,nop(j,pil)) + enddo + call calc_ab(aE,bN,pe,pn) +ctmp write(35,rec=ntrack) aE*resXY,bN*resXY + write(35,*) pe,pn,pz,pm + ntrack=ntrack+1 + enddo + endif + + t=t+dt + t_p=t_p+dt + + enddo + print *,'ibd=',ibd + + if(.not.hsteadystate.or.(.not.psteadystate)) then + WRITE(*,200) min(outputTS,newcondTS), oldparts, + * newparts,deadparts + else + WRITE(*,200) outputTS, oldparts, + * newparts,deadparts + endif +200 FORMAT('During the last',I5,' seconds:- [ Old:',I8, + 1 ' ; New:',I7,' ; Dead:',I7,' ]') + RETURN + END + +c-----------------------------------------------------------walkjustone + FUNCTION walkjustone(pel,pnl,pzl,pil,pa,pv,pm,cord,itype,ebox, + * nop,vel,dudt,ao,egrp,egps,nblk,nelb,nbb,nee,nxyz,pe, + * pn,pz,no_open,opbd,dtt,diffs,ibd) +c---------------------------------------------------------------------c +c purpose: c +c To walk one step. c +c---------------------------------------------------------------------c + include '3duparms.cb' + common /shape/shap(20),shpx(20),shpy(20),shpz(20) + common /etype/inode(4) + common /randm/isd, isd1 + LOGICAL walkjustone + + REAL cord(3,1),vel(3,1),egrp(6,1),egps(6,1),ebox(6,1) + REAL opbd(4,1),dudt(3,1) + REAL pe,pn,pz,pa,pel,pnl,pzl,xi,eta,zeta + REAL peold,pnold,pzold,peold1,pnold1,pzold1 + REAL*8 ao(1) + INTEGER*2 nop(20,1) + INTEGER nblk(2,1),nelb(1),nbb(1),nee(1),itype(1) + INTEGER num,pil,key,nen,nxyz,it,no_open + REAL ran3, normdis + LOGICAL boundaries + INTEGER isd, isd1 + REAL len,ang + REAL vE,vN,vZ + REAL pv + real*8 a(3,3),b(3),ai(3,3),du(3,3),a1,a2,a3,dudt0,dvdt0,dwdt0 + real*8 detj,ass + integer dtt + real diffs(2,1) + + walkjustone=.FALSE. +ccc do i=1,num_steps + pa=pa+dtt + if(pa.gt.tot_dead_age.or.pm/mass_pp.le.0.001) then +c if(pa.gt.tot_dead_age) then + walkjustone=.TRUE. + return + + endif +c +c----------gradually die-off +c + if(dieoff.and.pv.eq.0.0) then +c print*,tdpth + if(abs(pz).le.tdpth) then +c print*,'surface dieoff=',10.0**(-dtt/t90_s) +c print*,pm + pm=pm*10.0**(-dtt/t90_s) +c print*,pm + else +c print*,'deep dieoff=',10.0**(-dtt/t90_d) +c print*,pm + pm=pm*10.0**(-dtt/t90_d) +c print*,pm + endif + endif +c +c----------settleable particles on bed -- return +c + if(pv.eq.-9999.0) return +c +c----------interpolate velocity from FE mesh +c + vE=0.0 + vN=0.0 + vZ=0.0 + it=itype(pil) + nen=inode(it) +ccc call xn3(it,nen,pel,pnl,pzl) + do j=1,nen + vE=vE+shap(j)*vel(1,nop(j,pil)) + vN=vN+shap(j)*vel(2,nop(j,pil)) + vZ=vZ+shap(j)*vel(3,nop(j,pil)) +c print *,'pil,vel=',pil,j,vel(3,nop(j,pil)) + enddo + +c print *,'vE,vN,vZ1=',vE,vN,vZ +c +c----------predict the trajectory of the particle +c + if(.not.euler) then + call xn3x(it,nen,pel,pnl,pzl) +c +c----------calculate the Jacobian matrix +c + do j=1,3 + a(j,1)=0.0 + a(j,2)=0.0 + a(j,3)=0.0 + du(1,j)=0.0 + du(2,j)=0.0 + du(3,j)=0.0 + do k=1,nen + nele=nop(k,pil) + a(j,1)=a(j,1)+shpx(k)*cord(j,nele) + a(j,2)=a(j,2)+shpy(k)*cord(j,nele) + a(j,3)=a(j,3)+shpz(k)*cord(j,nele) + du(1,j)=du(1,j)+shpx(k)*vel(j,nele) + du(2,j)=du(2,j)+shpy(k)*vel(j,nele) + du(3,j)=du(3,j)+shpz(k)*vel(j,nele) + enddo + enddo + a1=a(2,2)*a(3,3)-a(2,3)*a(3,2) + a2=a(3,2)*a(1,3)-a(3,3)*a(1,2) + a3=a(1,2)*a(2,3)-a(1,3)*a(2,2) + detj=a(1,1)*a1+a(2,1)*a2+a(3,1)*a3 +c if(detj.lt.1.0e-10) then +c print *,'detj,pil,nen,pel,pnl,pzl=',detj,pil,nen,pel,pnl,pzl +c stop +c endif +c +c----------calculate the inverse Jacobi matrix +c + ai(1,1)=(a(2,2)*a(3,3)-a(2,3)*a(3,2))/detj + ai(1,2)=(a(3,1)*a(2,3)-a(2,1)*a(3,3))/detj + ai(1,3)=(a(2,1)*a(3,2)-a(3,1)*a(2,2))/detj + ai(2,1)=(a(3,2)*a(1,3)-a(1,2)*a(3,3))/detj + ai(2,2)=(a(1,1)*a(3,3)-a(1,3)*a(3,1))/detj + ai(2,3)=(a(3,1)*a(1,2)-a(1,1)*a(3,2))/detj + ai(3,1)=(a(1,2)*a(2,3)-a(2,2)*a(1,3))/detj + ai(3,2)=(a(2,1)*a(1,3)-a(1,1)*a(2,3))/detj + ai(3,3)=(a(1,1)*a(2,2)-a(1,2)*a(2,1))/detj + + do m=1,3 + do j=1,3 + a(m,j)=0.0 + do k=1,3 + a(m,j)=a(m,j)+dtt*ai(j,k)*du(k,m) + enddo +c a(m,j)=dtt*a(m,j) + enddo +c print *,'a=',(a(m,j),j=1,3) + enddo + do k=1,3 + a(k,k)=a(k,k)-2.0 + enddo + dudt0=0.0 + dvdt0=0.0 + dwdt0=0.0 + do k=1,nen + dudt0=dudt0+shap(k)*dudt(1,nop(k,pil)) + dvdt0=dvdt0+shap(k)*dudt(2,nop(k,pil)) + dwdt0=dwdt0+shap(k)*dudt(3,nop(k,pil)) + enddo + b(1)=-2.0*vE-dudt0*dtt + b(2)=-2.0*vN-dvdt0*dtt + b(3)=-2.0*vZ-dwdt0*dtt + +c do i=1,3 +c ass=max(dabs(a(i,1)),dabs(a(i,2)),dabs(a(i,3))) +c do j=1,3 +c a(i,j)=a(i,j)/ass +c enddo +c b(i)=b(i)/ass +c 25 enddo + a1=a(2,2)*a(3,3)-a(2,3)*a(3,2) + a2=a(3,2)*a(1,3)-a(3,3)*a(1,2) + a3=a(1,2)*a(2,3)-a(1,3)*a(2,2) + detj=a(1,1)*a1+a(2,1)*a2+a(3,1)*a3 +c detj=detj*dtt + vE=(a1*b(1)+a2*b(2)+a3*b(3))/detj + a1=a(3,1)*a(2,3)-a(3,3)*a(2,1) + a2=a(1,1)*a(3,3)-a(1,3)*a(3,1) + a3=a(2,1)*a(1,3)-a(2,3)*a(1,1) + vN=(a1*b(1)+a2*b(2)+a3*b(3))/detj + a1=a(2,1)*a(3,2)-a(2,2)*a(3,1) + a2=a(3,1)*a(1,2)-a(3,2)*a(1,1) + a3=a(1,1)*a(2,2)-a(1,2)*a(2,1) + vZ=(a1*b(1)+a2*b(2)+a3*b(3))/detj +c print *,'vE,vN,vZ2=',vE,vN,vZ,dtt +c stop + + endif +c------------------------------- + + peold=pe + pnold=pn + pzold=pz + + len=sqrt(6*diffs(1,pil)*dtt)*ran3(isd) + ang=6.2831853*ran3(isd) + angz=6.2831853*ran3(isd) +c if(vdiff.ne.0.0) then +c angz=6.2831853*ran3(isd) +c else +c angz=0.0 +c endif + peold1 = pe + vE*dtt + pnold1 = pn + vN*dtt + pe = peold1 + len*cos(ang)*cos(angz) + pn = pnold1 + len*sin(ang)*cos(angz) + if(pv.eq.0.0) then + pzold1 = pz + vZ*dtt + pz = pzold1 + ran3(isd)*sqrt(6*diffs(2,pil)*dtt)*sin(angz) + else + pzold1 = pz + (vZ+pv)*dtt + pz = pzold1 + endif + IF(boundaries(pe,pn,pz,peold,pnold,pzold,no_open,opbd)) THEN + walkjustone=.TRUE. + RETURN + + ENDIF +c --> + IF (pz.gt.0.0) pz=0.0 + num=pil + call trinvs3(nop,cord,itype,ebox,egrp,nblk,egps,nelb,nbb, + * nee,nxyz,nen,it,pe,pn,pz,xi,eta,zeta,num,key) + if(key.eq.1) then + pel=xi + pnl=eta + pzl=zeta + pil=num + else +c +c-------------neglect the random displacements +c + num=pil + if(pnold1.gt.0.0) pnold1=0.0 + call trinvs3(nop,cord,itype,ebox,egrp,nblk,egps,nelb,nbb, + * nee,nxyz,nen,it,peold1,pnold1,pzold1,xi,eta,zeta,num,key) + if(key.eq.1) then + pe=peold1 + pn=pnold1 + pz=pzold1 + pel=xi + pnl=eta + pzl=zeta + pil=num +c +c-------------neutral bouyant and floatable particales +c + else if(pv.ge.0.0) then + pe=peold + pn=pnold + pz=pzold + ibd=ibd+1 +c +c-------------settleable particles +c + else + num=pil + call trinvs3(nop,cord,itype,ebox,egrp,nblk,egps,nelb, + * nbb,nee,nxyz,nen,it,pe,pn,0.0,xi,eta,zeta,num,key) + if(key.eq.1) then + if(settles) then + no_settles=no_settles+1 + write(30,rec=no_settles) pe,pn + walkjustone=.true. + else + pz=0.0 + call xn3(it,nen,xi,eta,zeta) + do j=1,nen + pz=pz+shap(j)*ao(nop(j,num)) + enddo + pv=-9999.0 + endif + else + pe=peold + pn=pnold + pz=pzold + endif + endif + endif +ccc enddo + RETURN + END + +c------------------------------------------------------------boundaries + FUNCTION boundaries(pe,pn,pz,peold,pnold,pzold,no_open,opbd) +c---------------------------------------------------------------------c +c purpose: c +c To check whether the particle within the caculation region. c +c---------------------------------------------------------------------c + include '3duparms.cb' + include '3dgeom.cb' + + REAL pe, pn, pz, peold, pnold, pzold + REAL a, b + INTEGER no_open + LOGICAL boundaries + real opbd(4,1),s,r + real*8 deta + +c This is very rudimentary. It must be modified later. +c A particle can leave the domain in two ways:- +c i) Leave the (x,y) range. +c ii) Enter a part of the (x,y) range that is "land". +c Any particle leaving by mode (i) are considered to be "gone" +c Particles leaving by (ii) are returned to their old position. +c +c To particle is considered to entered a land area if:- +c a) It is beneath the bathymetry. +c b) It is in a grid cell that doesn't have all four +c corners within the ocean area. +c +c The function returns true only if the particle is "gone". + + boundaries=.FALSE. +c CALL calc_ab(a, b, pe, pn) + +c Check if particle has left the array region. Case i + +c IF ((a.le.0.0).or.(a.ge.float(numX)).or. +c 1 (b.le.0.0).or.(b.ge.float(numY))) THEN +c boundaries=.TRUE. +c RETURN +c ENDIF + +c-------check if particle has left from specified open boundaries + + if(openbd) then + do i=1,no_open + deta=(pnold-pn)*(opbd(3,i)-opbd(1,i)) + * -(peold-pe)*(opbd(4,i)-opbd(2,i)) + if(dabs(deta).lt.1.0e-8) goto 10 + r=((pnold-opbd(2,i))*(opbd(3,i)-opbd(1,i)) + * -(peold-opbd(1,i))*(opbd(4,i)-opbd(2,i)))/deta + if(r.lt.-1.0e-5.or.r.gt.1.0+1.0e-5) goto 10 + s=((pnold-pn)*(peold-opbd(1,i)) + * -(peold-pe)*(pnold-opbd(2,i)))/deta + if(s.le.1.0+1.0e-5.and.s.ge.-1.0e-5) then + boundaries=.TRUE. + return + endif + 10 enddo + endif + +c Check for particles below the bathymetry. (Case (iia)) + + RETURN + END + +c-----------------------------------------------------------interp_vels + subroutine interp_vels(vel1,vel2,vel,dudt,poll1,poll2,poll, + * v_time) +c---------------------------------------------------------------------c +c purpose: c +c To perform temporal interpolation of velocities. c +c---------------------------------------------------------------------c + INCLUDE '3duparms.cb' + INCLUDE '3dpolls.cb' + INCLUDE '3dmesh.cb' + dimension vel1(3,1),vel2(3,1),vel(3,1),poll1(4,1),poll2(4,1), + * poll(4,1),dudt(3,1) + REAL r + INTEGER i,j,v_time +c + r=float(v_time)/float(newcondTS) + do i=1,np + do j=1,3 + vel(j,i)=vel1(j,i)+r*(vel2(j,i)-vel1(j,i)) + enddo + enddo + do i=1,num_pollwins + do j=1,4 + poll(j,i)=poll1(j,i)+r*(poll2(j,i)-poll1(j,i)) + enddo + enddo +c +c-------calculate dudt +c + if(.not.euler) then + do i=1,np + do j=1,3 + dudt(j,i)=(vel2(j,i)-vel1(j,i))/float(newcondTS) + enddo + enddo + endif + + return + end + diff --git a/source/askstr.for b/source/askstr.for new file mode 100644 index 0000000..b18c0f6 --- /dev/null +++ b/source/askstr.for @@ -0,0 +1,25 @@ + SUBROUTINE ASKSTR(PROMPT, STRING) +C +C Prompt standard input for a string, and return that string. +C If just is pressed, the passed string is not changed. + CHARACTER*(*) PROMPT,STRING +C + INTEGER STRLEN + CHARACTER*80 TMPSTR +C + STRLEN = LEN(STRING) + DO WHILE(STRING(STRLEN:STRLEN) .EQ. ' ') + IF (STRLEN .EQ. 1) GOTO 150 + STRLEN = STRLEN - 1 + END DO + 150 CONTINUE +cdrc TYPE 10,PROMPT,STRING(1:STRLEN) + PRINT 10,PROMPT,STRING(1:STRLEN) + 10 FORMAT(1X,A,' [',A,']: ',$) + READ(*,'(A)') TMPSTR + IF (TMPSTR .NE. ' ') THEN + STRING = TMPSTR + END IF +C + RETURN + END diff --git a/source/blocks.for b/source/blocks.for new file mode 100644 index 0000000..239a9e6 --- /dev/null +++ b/source/blocks.for @@ -0,0 +1,245 @@ +c----------------------------------------------------------------blocks + subroutine blocks(nx,ny,nz,np,ne,nop,cord,nxyz,egrp,nbeg,nend, + * nelg,nbb,nee,nelb,nblk,egps,itype,ebox,m,num) +c---------------------------------------------------------------------c +c Purpose: c +c To divide FE mesh into blocks. c +c---------------------------------------------------------------------c + common /etype/inode(4) + integer*2 nop + dimension nop(20,1),cord(3,1),egrp(6,1),nbeg(1),nend(1), + * nelg(1),nee(1),nbb(1),nelb(1),nblk(2,1),egps(6,1), + * itype(1),ebox(6,1) +c + xmax=cord(1,1) + xmin=cord(1,1) + ymax=cord(2,1) + ymin=cord(2,1) + zmax=cord(3,1) + zmin=cord(3,1) + disd=0.0 + do i=1,ne + if(nop(9,i).eq.0) goto 10 +c +c----------store element type +c + if(nop(20,i).ne.0) then + itype(i)=3 + else if(nop(15,i).ne.0) then + itype(i)=2 + else if(nop(13,i).ne.0) then + itype(i)=4 + else + itype(i)=1 + endif +c +c----------find B.L. and T.R. conners of each element +c + ebox(1,i)=cord(1,nop(1,i)) + ebox(4,i)=cord(1,nop(1,i)) + ebox(2,i)=cord(2,nop(1,i)) + ebox(5,i)=cord(2,nop(1,i)) + ebox(3,i)=cord(3,nop(1,i)) + ebox(6,i)=cord(3,nop(1,i)) + do k=2,inode(itype(i)) + ebox(1,i)=min(ebox(1,i),cord(1,nop(k,i))) + ebox(4,i)=max(ebox(4,i),cord(1,nop(k,i))) + ebox(2,i)=min(ebox(2,i),cord(2,nop(k,i))) + ebox(5,i)=max(ebox(5,i),cord(2,nop(k,i))) + ebox(3,i)=min(ebox(3,i),cord(3,nop(k,i))) + ebox(6,i)=max(ebox(6,i),cord(3,nop(k,i))) + enddo + ebox(1,i)=ebox(1,i)-1.0e-5*abs(ebox(1,i)) + ebox(4,i)=ebox(4,i)+1.0e-5*abs(ebox(4,i)) + ebox(2,i)=ebox(2,i)-1.0e-5*abs(ebox(2,i)) + ebox(5,i)=ebox(5,i)+1.0e-5*abs(ebox(5,i)) + ebox(3,i)=ebox(3,i)-1.0e-5*abs(ebox(3,i)) + ebox(6,i)=ebox(6,i)+1.0e-5*abs(ebox(6,i)) +c +c----------calculate the maximum diagonal distance +c + disd0=sqrt((ebox(4,i)-ebox(1,i))*(ebox(4,i)-ebox(1,i))+ + * (ebox(5,i)-ebox(2,i))*(ebox(5,i)-ebox(2,i))) + disd=max(disd,disd0) +c +c-------find B.L. and T.R. conners of the FE mesh +c + xmax=max(xmax,ebox(4,i)) + xmin=min(xmin,ebox(1,i)) + ymax=max(ymax,ebox(5,i)) + ymin=min(ymin,ebox(2,i)) + zmax=max(zmax,ebox(6,i)) + zmin=min(zmin,ebox(3,i)) + 10 enddo + + m=1 + dx=(xmax-xmin)/float(nx) + if(dx.lt.disd) then + nx=(xmax-xmin)/disd + if(nx.lt.1) nx=1 + dx=(xmax-xmin)/float(nx) + endif + dy=(ymax-ymin)/float(ny) + if(dy.lt.disd) then + ny=(ymax-ymin)/disd + if(ny.lt.1) ny=1 + dy=(ymax-ymin)/float(ny) + endif + dz=(zmax-zmin)/float(nz) + x=xmin +c +c-------divide into blocks +c + do i=1,nx + y=ymin + do j=1,ny + z=zmin + do k=1,nz + egrp(1,m)=x-dx*0.01 + egrp(2,m)=y-dy*0.01 + egrp(3,m)=z-dz*0.01 + egrp(4,m)=x+dx*1.01 + egrp(5,m)=y+dy*1.01 + egrp(6,m)=z+dz*1.01 + z=z+dz + m=m+1 + enddo + y=y+dy + enddo + x=x+dx + enddo + nxyz=nx*ny*nz + nxyz0=nxyz + num=1 + do i=1,nxyz + nbeg(i)=num + do j=1,ne + if(nop(9,j).ne.0) then + do k=1,inode(itype(j)) + if(cord(1,nop(k,j)).gt.egrp(4,i).or. + * cord(1,nop(k,j)).lt.egrp(1,i)) goto 30 + if(cord(2,nop(k,j)).gt.egrp(5,i).or. + * cord(2,nop(k,j)).lt.egrp(2,i)) goto 30 + if(cord(3,nop(k,j)).gt.egrp(6,i).or. + * cord(3,nop(k,j)).lt.egrp(3,i)) goto 30 + nelg(num)=j + num=num+1 + goto 40 + 30 enddo + endif + 40 enddo + nend(i)=num-1 + enddo + i=1 + 45 if(nend(i)-nbeg(i).lt.0) then + nxyz=nxyz-1 + do j=i,nxyz + nbeg(j)=nbeg(j+1) + nend(j)=nend(j+1) + do k=1,6 + egrp(k,j)=egrp(k,j+1) + enddo + enddo + else + i=i+1 + endif + if(i.le.nxyz) goto 45 +c +c-------divide into sub-blocks +c + m=1 + num=1 + do n=1,nxyz + num_elem=nend(n)-nbeg(n)+1 + if(num_elem.gt.nxyz+nxyz) then + a_multiplier=sqrt(float(num_elem))/float(nx*ny) + nxs=int(nx*a_multiplier) + if(nxs.lt.1) nxs=1 + nys=int(ny*a_multiplier) + if(nys.lt.1) nys=1 + dx=(egrp(4,n)-egrp(1,n))/float(nxs) + dy=(egrp(5,n)-egrp(2,n))/float(nys) + if(dx.lt.disd) then + nxs=(egrp(4,n)-egrp(1,n))/disd + if(nxs.lt.1) nxs=1 + dx=(egrp(4,n)-egrp(1,n))/float(nxs) + endif + if(dy.lt.disd) then + dy=disd + nys=(egrp(5,n)-egrp(2,n))/disd + if(nys.lt.1) nys=1 + dy=(egrp(5,n)-egrp(2,n))/float(nys) + endif + nxyzs=nxs*nys + if(nxyzs.le.1) goto 70 + dz=(egrp(6,n)-egrp(3,n)) + + nblk(1,n)=nxyzs + nblk(2,n)=m + x=egrp(1,n) + do i=1,nxs + y=egrp(2,n) + do j=1,nys + z=egrp(3,n) + egps(1,m)=x-0.01*dx + egps(2,m)=y-0.01*dy + egps(3,m)=z-0.01*dz + egps(4,m)=x+dx*1.01 + egps(5,m)=y+dy*1.01 + egps(6,m)=z+dz*1.01 + y=y+dy + m=m+1 + enddo + x=x+dx + enddo + do i=nblk(2,n),nblk(2,n)+nblk(1,n)-1 + nbb(i)=num + do j=nbeg(n),nend(n) + if(nop(9,nelg(j)).ne.0) then + do k=1,inode(itype(nelg(j))) + if(cord(1,nop(k,nelg(j))).gt.egps(4,i).or. + * cord(1,nop(k,nelg(j))).lt.egps(1,i)) goto 50 + if(cord(2,nop(k,nelg(j))).gt.egps(5,i).or. + * cord(2,nop(k,nelg(j))).lt.egps(2,i)) goto 50 + if(cord(3,nop(k,nelg(j))).gt.egps(6,i).or. + * cord(3,nop(k,nelg(j))).lt.egps(3,i)) goto 50 + nelb(num)=nelg(j) + num=num+1 + goto 60 + 50 enddo + endif + 60 enddo + nee(i)=num-1 + enddo + i=nblk(2,n) + 65 if(nee(i)-nbb(i).lt.0) then + nblk(1,n)=nblk(1,n)-1 + m=m-1 + do j=i,nblk(2,n)+nblk(1,n)-1 + nbb(j)=nbb(j+1) + nee(j)=nee(j+1) + do k=1,6 + egps(k,j)=egps(k,j+1) + enddo + enddo + else + i=i+1 + endif + if(i.le.nblk(2,n)+nblk(1,n)-1) goto 65 + else + 70 nblk(1,n)=1 + nblk(2,n)=m + nbb(m)=num + do i=nbeg(n),nend(n) + nelb(num)=nelg(i) + num=num+1 + enddo + do j=1,6 + egps(j,m)=egrp(j,n) + enddo + nee(m)=num-1 + m=m+1 + endif + enddo + return + end diff --git a/source/dieoffs.for b/source/dieoffs.for new file mode 100644 index 0000000..df3f700 --- /dev/null +++ b/source/dieoffs.for @@ -0,0 +1,43 @@ + SUBROUTINE DIEOFFS(TOFDAY,DAYOFY,DELT,T90) + +C Calculates faecal coliform dieoff T90 based on solar radiation during +C the time interval DELT(sec). +C Input Time of day (hrs) TOFDAY, Day of year (1-366) DAYOFY. +C Output T90 (hrs) + + REAL TOFDAY,DAYOFY,LAT,LONG,ELEVAT,DELT,DELTH,STANM + REAL TSOLHR,T90 + + COMMON/INPUTS/LAT,LONG,ELEVAT,STANM + COMMON/MET/WETBLB,ATMPR,DRYBLB,CLOUD,DAT + + LAT=-33.8 + LONG=151.3 + ELEVAT=0.0 ! Elevation (m) + DELTH=DELT/3600. + STANM=150.0 ! Longitude corresponding to standard time + WETBLB=20.0 ! Dry bulb temp + DRYBLB=17.5 ! Wet bulb temp (depends on humidity) + CLOUD=0.0 ! Cloud cover (0-1) + ATMPR=1015. ! Atm pressure + DAT=0.07 ! Dust attenuation + +C Dieoff-Solar Radiation Parameters +C Ref. Caldwell-Connell Engineers (1980) Fig. 7.6 + A=3.3 + B=-0.7 + + CALL HEATEX(TOFDAY,DAYOFY,DELT,TSOLHR) + +C Convert radiation in DELT to hourly solar radiation in MJ/m2 + TSOLHR=TSOLHR/1000./DELTH +C PRINT*,TSOLHR +C Calculate the T90 + IF (TSOLHR.LT.0.08) THEN + T90=20. + ELSE + T90=A*(TSOLHR**(B)) + ENDIF + + RETURN + END diff --git a/source/flib.for b/source/flib.for new file mode 100644 index 0000000..44e27b5 --- /dev/null +++ b/source/flib.for @@ -0,0 +1,875 @@ +c----------------------------------------------------------------setcst + subroutine setcst +c---------------------------------------------------------------------c +c purpose: c +c To set consats. c +c---------------------------------------------------------------------c +ccc implicit real*8 (a-h,o-z) + common mtot,np,ia(1) + common /dbsys/numa,next,idir,ipp(3) + common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 + * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 +c + ipp(1)=4 + ipp(2)=8 + ipp(3)=1 +c mtot=100000 + numa=0 + next=1 +c idir=mtot + ntm=6 + ntr=5 + nin=1 + not=2 + nsp=3 + nfl=4 + nt7=7 + nt8=8 + nt9=9 + nt10=10 + nt11=11 + nt12=12 + nt13=13 + nt14=14 + nt15=15 + nt16=16 + nt17=17 + nt18=18 + nt19=19 +c + return + end +c-----------------------------------------------------------------upper + subroutine upper +c---------------------------------------------------------------------c +c purpose: c +c This routine converts lower case to upper case c +c---------------------------------------------------------------------c +ccc implicit real*8 (a-h,o-z) + character*1 line,sp,sc +c + common /cline/line(160) + common /iline/ii +c + data sp/' '/,sc/':'/ +c +c------check end of line and covert to upper case +c + isp=ichar(sp) + do 80 i=1,ii + if(line(i).eq.sc) goto 90 + nn=ichar(line(i)) + if(nn.gt.96) line(i)=char(nn-isp) + 80 continue +c +c------remove blanks at end of line +c + jj=0 + do 50 j=1,ii + if(line(j).ne.sp) goto 70 + jj=jj+1 + 50 continue + 70 if(jj.eq.0) goto 85 + if(jj.eq.ii) goto 100 + kk=jj+1 + do 60 j=1,ii-jj + line(j)=line(kk) + kk=kk+1 + 60 continue + ii=ii-jj-1 + 85 if(line(ii).ne.sp) goto 95 + ii=ii-1 + goto 85 + 90 ii=i +c + 95 return + 100 ii=1 + return + end +c----------------------------------------------------------------bcmane + subroutine bcname +c---------------------------------------------------------------------c +c purpose: c +c To print arry names in 'ia'. c +c---------------------------------------------------------------------c +ccc implicit real*8 (a-h,o-z) + character*1 name(4) + common mtot,np,ia(1) + common /dbsys/numa,next,idir,ipp(3) + common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 + * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 +c + k=idir + do 10 i=1,numa + do 20 j=1,4 + name(j)=char(ia(k+j-1)) + 20 continue +ccc write(not,1000) i,name,ia(k+4),ia(k+5),ia(k+7),ia(k+8) + write(ntm,1000) i,name,ia(k+4),ia(k+5),ia(k+7),ia(k+8) + 1000 format(1x,i7,5x,4a1,',nr=',i7,',nc=',i7,',na=',i12, + * ',size=',i12) + k=k+10 + 10 continue + return + end +c----------------------------------------------------------------definr + subroutine definr(name,na,nr,nc) +c---------------------------------------------------------------------c +c purpose: c +c To reserve storage for REAL data. c +c---------------------------------------------------------------------c +ccc implicit real*8 (a-h,o-z) + character*1 name(4) + common mtot,np,ia(1) +c + np=1 + call defin(name,na,nr,nc) + return + end +c----------------------------------------------------------------defini + subroutine defini(name,na,nr,nc) +c---------------------------------------------------------------------c +c purpose: c +c To reserve storage for INTEGER data. c +c---------------------------------------------------------------------c +ccc implicit real*8 (a-h,o-z) + character*1 name(4) + common mtot,np,ia(1) +c + np=1 + call defin(name,na,nr,nc) + return + end +c-----------------------------------------------------------------defin + subroutine defin(name,na,nr,nc) +c---------------------------------------------------------------------c +c purpuse: c +c To define and reserve storage for array. c +c---------------------------------------------------------------------c +ccc implicit real*8 (a-h,o-z) + character*1 name(4) + common mtot,np,ia(1) + common /dbsys/numa,next,idir,ipp(3) + common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 + * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 +c +c---------------------------------------------------------------------c +c where name = name of array--4 logicals maxmum c +c na = location of array if in blank common c +c nr = number of rows c +c nc = number of columns c +c mtot = end of directory c +c numa = number of arrays in data base c +c next = next avialable storage location c +c idir = start of directory in data base c +c ipp = number of logicals contained in data type c +c lenr = number of logicals in physical record c +c np = type of data c +c = 1 integer data c +c = 2 real data c +c = 3 logical data c +c------directory definition for core or sequential files--------------c +c idir(i,n) = name of array --iname(4,char) +c idir(5,n) = number of rows -- nr c +c idir(6,n) = number of colomns -- nc c +c idir(7,n) = type of data -- np c +c idir(8,n) = incore address -- na c +c = -1 if sequential file on disk c +c = -2 if direct access on disk c +c idir(9,n) = size of array c +c idir(10,n) = 0 if in core storage c +c------directory defintion for direct access files--------------------c +c idir(5,n) = number of integers c +c idir(6,n) = number of real words c +c idir(7,n) = number of logicals c +c idir(8,n) = number of logical records c +c idir(9,n) = logical record number c +c idir(10,n) = lun if on logical unit lun c +c---------------------------------------------------------------------c +c +c------evaluate storage requirements +c + if((np.ne.1).and.(np.ne.3)) np=2 + nsize=(nr*nc*ipp(np)-1)/(ipp(1)*2) + nsize=nsize*2+2 + na=next + next=next+nsize +c +c------set up new directory +c + numa=numa+1 + idir=idir-10 + i=idir +c +c------check storage limits +c + if(i.ge.next) goto 100 + i=next-i+mtot-1 + write(ntm,2000) i,mtot,(name(k),k=1,4) + stop + 100 call icon(name,ia(i)) + ia(i+4)=nr + ia(i+5)=nc + ia(i+6)=np + ia(i+7)=na + ia(i+8)=nsize + ia(i+9)=0 + 900 return +c---------------------------------------------------------------------c + 2000 format(/,'storage fault in defin',/,'storage required =',i9, + * /,'storage available =',i9, + * /,'name = ',4a) + end +c +c----------------------------------------------------------------locate + subroutine locate(name,na,nr,nc) +c---------------------------------------------------------------------c +c purpose: c +c To locate array NAME in IA. c +c---------------------------------------------------------------------c +ccc implicit real*8 (a-h,o-z) + character*1 name + dimension name(4),iname(4) + common mtot,np,ia(1) + common /dbsys/numa,next,idir,ipp(3) + common /iolist/ntm,ntr,nun,not,nsp,nfl,nt7,nt8,nt9,nt10 + * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 +c +c------locate and return properties on array +c + na=0 + call icon(name,iname) + i=ifind(iname,0) + if(i.eq.0) goto 900 +c +c------return array properties +c + na=ia(i+7) + nr=ia(i+4) + nc=ia(i+5) + np=ia(i+6) + return +c---------------------------------------------------------------------c + 900 write(ntm,1000) name + 1000 format(/,' NO SUCH ARRAY ',4a1,' IN IA ',/) + return + end +c +c----------------------------------------------------------------delete + subroutine delete(name) +c---------------------------------------------------------------------c +c purpose: c +c To delete array NAME. c +c---------------------------------------------------------------------c +ccc implicit real*8 (a-h,o-z) + character*1 name + dimension name(4),iname(4) + common mtot,np,ia(1) + common /dbsys/numa,next,idir,ipp(3) + common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 + * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 +c +c------delete array from storage +c + 100 call icon(name,iname) + i=ifind(iname,0) + if(i.eq.0) goto 900 +c +c------check on storage location +c + 200 nsize=ia(i+8) +c +c------set size of array +c + next=next-nsize + numa=numa-1 + na=ia(i+7) +c +c------check if out of core or direct access +c + if(na.gt.0) goto 500 + write(ntm,1000) name + goto 800 + 500 if(na.eq.next) goto 800 +c +c------compact storage +c + ii=na+nsize + nnxt=next-1 + do 700 j=na,nnxt + ia(j)=ia(ii) + ii=ii+1 + 700 continue +c +c------compact and update directory +c + 800 na=i-idir + idir=idir+10 + if(na.eq.0) goto 900 + na=na/10 + do 860 k=1,na + ii=i+9 + do 850 j=1,10 + ia(ii)=ia(ii-10) + ii=ii-1 + 850 continue + if(ia(i+7).le.0) goto 860 + if(ia(i+9).eq.0) ia(i+7)=ia(i+7)-nsize + i=i-10 + 860 continue +c + 900 return +c---------------------------------------------------------------------c + 1000 format(/,' NAME ',4a1,' IS BEING USED FOR AN OUT OF CORE FILE'/) + end +c +c------------------------------------------------------------------icon + subroutine icon(name,iname) +c---------------------------------------------------------------------c +c purpose: c +c To convert logicals to integer data. c +c---------------------------------------------------------------------c +ccc implicit real*8 (a-h,o-z) + character*1 name(4) + dimension iname(4) +c + do 100 i=1,4 + iname(i)=ichar(name(i)) + 100 continue +c + return + end +c-----------------------------------------------------------------ifind + function ifind(iname,lun) +c---------------------------------------------------------------------c +c purpose: c +c This function is used to find location of array. c +c---------------------------------------------------------------------c +ccc implicit real*8 (a-h,o-z) + dimension iname(4) + common mtot,np,ia(1) + common /dbsys/numa,next,idir,ipp(3) + common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 + * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 +c + i=idir + do 100 n=1,numa + if(lun.ne.ia(i+9)) goto 300 + if(iname(1).ne.ia(i )) goto 300 + if(iname(2).ne.ia(i+1)) goto 300 + if(iname(3).ne.ia(i+2)) goto 300 + if(iname(4).eq.ia(i+3)) goto 200 + 300 i=i+10 + 100 continue + i=0 + 200 ifind=i +c + return + end +c----------------------------------------------------------------refinr + subroutine refinr(name,na,nr,nc) +c---------------------------------------------------------------------c +c purpose: c +c To reserve storage for Real array. c +c---------------------------------------------------------------------c +ccc implicit real*8 (a-h,o-z) + character*1 name(4) + common mtot,np,ia(1) +c + np=1 + call refin(name,na,nr,nc) + return + end +c----------------------------------------------------------------refini + subroutine refini(name,na,nr,nc) +c---------------------------------------------------------------------c +c purpose: c +c To reserve storage for Integer array. c +c---------------------------------------------------------------------c +ccc implicit real*8 (a-h,o-z) + character*1 name(4) + common mtot,np,ia(1) +c + np=1 + call refin(name,na,nr,nc) + return + end +c-----------------------------------------------------------------refin + subroutine refin(name,na,nr,nc) +c---------------------------------------------------------------------c +c purpose: c +c To refinde and reserve new storage for array. c +c---------------------------------------------------------------------c +ccc implicit real*8 (a-h,o-z) + character*1 name(4) + dimension iname(4) + common mtot,np,ia(1) + common /dbsys/numa,next,idir,ipp(3) + common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 + * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 +c +c------locate origin array +c + call icon(name,iname) + i=ifind(iname,0) + nr1=ia(i+4) + nc1=ia(i+5) + na=ia(i+7) + nsize1=ia(i+8) + nsize=(nr*nc*ipp(np)-1)/(ipp(1)*2) + nsize=nsize*2+2 + nszicr=nsize-nsize1 + nrincr=nr-nr1 + ncincr=nc-nc1 + if(nszicr.eq.0) goto 450 +c +c------set size of array +c + nexto=next-1 + next=next+nszicr + if(i.gt.next) goto 450 + i=next-i+mtot-1 + print *,'name,nr,nc=',name,nr,nc + write(ntm,1000) i,mtot,(name(k),k=1,4) + stop + 450 ia(i+4)=nr + ia(i+5)=nc + ia(i+8)=nsize + if(nszicr.eq.0) goto 900 +c +c------compact storage +c + nextn=next-1 + nnew=na+nsize + nold=na+nsize1 + if(np.eq.3) goto 500 + if(nrincr.eq.0) goto 500 + nrincr=nrincr*np + nexto=nexto-nrincr + nnew=na + nold=na + do 10 k=1,nc1 + nnew=nnew+nr*np + nold=nnew-nrincr + nexto=nexto+nrincr + if(i.gt.nexto) goto 480 + i=nexto-i+mtot + print *,'name,nr,nc=',name,nr,nc + write(ntm,1000) i,mtot,(name(m),m=1,4) + stop + 480 if(nrincr.gt.0) goto 510 + ij=nnew + do 20 j=nold,nexto + ia(ij)=ia(j) + ij=ij+1 + 20 continue + goto 10 + 510 ij1=nexto+nrincr + do 50 j=nold,nexto + ii1=nexto-j+nold + ia(ij1)=ia(ii1) + ij1=ij1-1 + 50 continue + 10 continue + nnew=na+nsize + nexto=nexto+nrincr + nold=nnew-(nextn-nexto) + 500 if(i.eq.idir) goto 900 + if(nextn.eq.nexto) goto 600 + if(nnew.gt.nold) goto 520 + ij=nold + do 30 j=nnew,nextn + ia(j)=ia(ij) + ij=ij+1 + 30 continue + goto 600 + 520 ij1=nexto + do 60 j=nnew,nextn + ii1=nextn-j+nnew + ia(ii1)=ia(ij1) + ij1=ij1-1 + 60 continue +c +c------compact and update directory +c + 600 ij=(i-idir)/10 + if(ij.eq.0) goto 900 + i=i-10 + do 40 j=1,ij + if(ia(i+7).le.0) goto 40 + if(ia(i+9).eq.0) ia(i+7)=ia(i+7)+nszicr + i=i-10 + 40 continue +c + 900 return +c---------------------------------------------------------------------c + 1000 format(/,' STORAGE FAULT IN REFIN',/,' STORAGE REQUIRED =',I9, + * /,' STORAGE AVIABLE =',I9, + * /,' NAME = ',4a) + end +c------------------------------------------------------------------find + subroutine find(sep,nostr,nn) +c---------------------------------------------------------------------c +c input parameters: c +c sep--separator c +c nn--flag whether there is separator or not c +c---------------------------------------------------------------------c +c implicit real*8 (a-h,o-z) + character line*160,sep*80 + common /cline/line + common /iline/ii + common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 + * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 +c +c------find separator in input file +c + nn=0 + 50 read(nin,1000,err=900,end=900) line(1:80) + ii=80 +c +c------convert to upper case +c + call upper + if(sep(1:nostr).ne.line(1:nostr)) goto 50 + nn=1 + 900 return +c---------------------------------------------------------------------c + 1000 format(a80) + end +c------------------------------------------------------------------free + subroutine free +c---------------------------------------------------------------------c +c purpose: c +c To read a line of free field data. c +c---------------------------------------------------------------------c +c implicit real*8 (a-h,o-z) + character*1 line,sp,bs,sc,c + common /cline/line(160) + common /iline/ ii + common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 + * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 + data sp/' '/,c/'%'/,sc/':'/,bs/'\\'/ +c +c------read line of free field data +c + do 40 i=1,160 + line(i)=sp + 40 continue + 50 i=1 + ii=80 + 60 read(nin,1000,err=100) (line(k),k=i,ii) +c +c------check for additional line +c + 70 do 80 k=i,ii + if(line(k).ne.bs) goto 80 + i=k + ii=k+79 + if(ii.gt.160) ii=160 + goto 60 + 80 continue + call upper +c +c------check for comment +c + if(line(1).ne.c) goto 900 + write(ntm,2000) (line(i),i=1,ii) + goto 50 +c +c------error in data +c + 100 write(ntm,2100) + call freept +c + 900 return +c---------------------------------------------------------------------c + 1000 format(80a1) + 2000 format(1x,80a1) + 2100 format(/'* Error in reading input line *'/) + end +c----------------------------------------------------------------preept + subroutine freept +c---------------------------------------------------------------------c +c purpose: c +c To write record to file NOT and NTM. c +c---------------------------------------------------------------------c +c implicit real*8 (a-h,o-z) + character*1 line +c + common /iline/ii + common /cline/line(160) + common /iolist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 + * ,nt11,nt12,nt13,nt14,nt15,nt16,nt17,nt18,nt19 +c +c------echo of free field information +c + write(ntm,1000) (line(i),i=1,ii) + return + 1000 format(1x,80a1) + end +c-----------------------------------------------------------------freeh + subroutine freeh(ic,idata,nc,num) +c---------------------------------------------------------------------c +c purpose: c +c To extract hollerith data. c +c---------------------------------------------------------------------c +c num -- negative read strings separated by blank space or comma c +c num -- positive read nc number of characters c +c---------------------------------------------------------------------c +c implicit real*8 (a-h,o-z) + character*1 line,ic,blk,comma,io,eqs,colin,idata + dimension idata(nc,num) + common /cline/line(160) + common /iline/ii + data blk/' '/,comma/','/,io/'0'/,eqs/'='/,colin/':'/ +c +c------find hollerith string +c + num0=abs(num) + 90 i=0 + do 210 j=1,num0 + do 210 n=1,nc + idata(n,j)=blk + 210 continue + if(ic.eq.blk) goto 200 + do 100 i=2,ii + if((line(i-1).eq.ic).and.(line(i).eq.eqs)) goto 200 + 100 continue + return +c +c------extract hollerith string +c + 200 do 300 j=1,num0 + 260 i=i+1 + if(i.gt.ii) goto 400 + do 290 n=1,nc + if(line(i).eq.eqs) goto 400 + if(line(i).eq.colin) goto 300 + if(num.gt.0) goto 280 + if(line(i).eq.blk) goto 300 + if(line(i).eq.comma) goto 300 + 280 idata(n,j)=line(i) + if(n.eq.nc) goto 290 + i=i+1 + 290 continue + 300 continue +c + 400 return + end +c-----------------------------------------------------------------freei + subroutine freei(ic,idata,num) +c---------------------------------------------------------------------c +c purpose: c +c To extract integer data. c +c---------------------------------------------------------------------c +c input parameter: c +c ic--flag of data ic=' ' unidentifier c +c ic=else identifier c +c num--number of data c +c---------------------------------------------------------------------c +c output parameter: c +c idata--integer array c +c---------------------------------------------------------------------c +c implicit real*8 (a-h,o-z) + character*1 line,ic,blk,comma,io,eqs,neg,colin,lne,pls + dimension idata(num) + common /cline/line(160) + common /iline/ii + common /ilist/ntm,ntr,nin,not,nsp,nfl,nt7,nt8,nt9,nt10 + data blk/' '/,comma/','/,io/'0'/,eqs/'='/,neg/'-'/,colin/':'/ + data pls/'+'/ +c +c------find integer string +c + 90 i=0 + do 210 j=1,num + idata(j)=0 + 210 continue + if(ic.eq.blk) goto 200 + do 100 i=1,ii + if(line(i).eq.ic.and.line(i+1).eq.eqs) goto 200 + 100 continue + return +c +c------zero integer string +c + 200 if(line(i+1).eq.eqs) i=i+1 + do 250 j=1,num + isign=1 +c +c------skip blanks between integers +c + 215 if(line(i+1).ne.blk) goto 220 + i=i+1 + if(i.gt.ii) goto 900 + goto 215 + 220 i=i+1 + if(i.gt.ii) goto 230 +c +c------check for sign +c + if(line(i).ne.pls) goto 235 + i=i+1 + if(i.gt.ii) goto 230 + 235 lne=line(i) + if(lne.ne.neg) goto 225 + isign=-1 + goto 220 +c +c------extract integer +c + 225 if(lne.eq.blk) goto 230 + if(lne.eq.comma) goto 230 + if(lne.eq.colin) goto 230 + nn=ichar(lne)-ichar(io) + if((nn.lt.0).or.(nn.gt.9)) goto 900 + idata(j)=10*idata(j)+nn + goto 220 +c +c------set sign +c + 230 idata(j)=idata(j)*isign + 250 continue + 900 return + end +c-----------------------------------------------------------------freer + subroutine freer(ic,data,num) +c---------------------------------------------------------------------c +c purpose: c +c To extract real data. c +c---------------------------------------------------------------------c +c implicit real*8 (a-h,o-z) + character*1 line,blk,ic,mul,div,add,sub,eqs,e + dimension data(num) + common /cline/line(160) + common /iline/ii +c + data blk/' '/,mul/'*'/,div/'\/'/,add/'+'/,sub/'-'/,e/'E'/, + * eqs/'='/ +c +c------find real string +c + 90 i=0 + do 260 j=1,num + data(j)=0.0 + 260 continue + if(ic.eq.blk) goto 250 + do 100 i=1,ii + if((line(i).eq.ic).and.(line(i+1).eq.eqs)) goto 250 + 100 continue + return +c +c------extract real data +c + 250 do 300 j=1,num + jj=0 + 270 if(i.gt.ii) goto 300 + call rdata(i,xx,nn) + if(jj.ne.0) goto 275 + data(j)=xx + goto 290 +c +c------arithmetric statement +c + 275 if(jj.eq.1) data(j)=data(j)*xx + if(jj.eq.2) data(j)=data(j)/xx + if(jj.eq.3) data(j)=data(j)+xx + if(jj.eq.4) data(j)=data(j)-xx + if(jj.ne.5) goto 290 +c +c------exponential data +c + jj=abs(xx) + if(jj.eq.0) goto 290 + do 280 k=1,jj + if(xx.lt.0.0) data(j)=data(j)/10.0 + if(xx.gt.0.0) data(j)=data(j)*10.0 + 280 continue +c +c------set type of statement +c + 290 jj=0 + if(line(i).eq.mul) jj=1 + if(line(i).eq.div) jj=2 + if(line(i).eq.add) jj=3 + if(line(i).eq.sub) jj=4 + if(line(i).eq.e) jj=5 + if(line(i+1).eq.eqs) jj=0 + if(jj.ne.0) goto 270 + if(nn.gt.9) return + 300 continue +c---------------------------------------------------------------------c + return + end +c-----------------------------------------------------------------rdata + subroutine rdata(i,xx,nn) +c---------------------------------------------------------------------c +c purpose: c +c To convert string to real floating point number. c +c---------------------------------------------------------------------c +ccc implicit real*8 (a-h,o-z) + character*1 line,blk,comma,io,dot,neg,eqs,pls + common /cline/line(160) + common /iline/ii + data blk/' '/,comma/','/,io/'0'/,dot/'.'/,neg/'-'/,eqs/'='/ + data pls/'+'/ +c +c------converts string to real floating point number +c + if(line(i+1).eq.eqs) i=i+1 + y=0 + is=1 + xx=0.0 + 250 if(line(i+1).ne.blk) goto 255 + i=i+1 + if(i.gt.ii) goto 300 + goto 250 + 255 if(line(i+1).ne.neg) goto 260 + is=-1 + i=i+1 + if(i.gt.ii) goto 300 + 260 if(line(i+1).ne.pls) goto 265 + i=i+1 + if(i.gt.ii) goto 300 + 265 if(line(i+1).ne.blk) goto 270 + i=i+1 + if(i.gt.ii) goto 300 + goto 265 + 270 i=i+1 + if(i.gt.ii) goto 300 + if(line(i).eq.blk.and.line(i+1).eq.blk) goto 270 + nn=ichar(line(i))-ichar(io) + xn=isign(nn,is) + if(line(i).ne.dot) goto 275 + y=1.0 + goto 270 + 275 if(line(i).eq.blk) goto 300 + if(line(i).eq.comma) goto 300 + if(nn.lt.0.or.nn.gt.9) goto 300 + if(y.eq.0) goto 280 + y=y/10.0 + xn=xn*y + xx=xx+xn + goto 270 + 280 xx=10.0*xx+xn + goto 270 + 300 return + end +c +c---------------------------------------------------------------clearr + subroutine clearr(adata,nr,nc) +c--------------------------------------------------------------------c +c purpose: c +c This routine is used to initialize real array ADATA. c +c--------------------------------------------------------------------c +ccc implicit real*8 (a-h,o-z) + dimension adata(nr,nc) +c + do 10 i=1,nr + do 20 j=1,nc + adata(i,j)=0.0 + 20 continue + 10 continue + return + end + diff --git a/source/gdate.for b/source/gdate.for new file mode 100644 index 0000000..aff02db --- /dev/null +++ b/source/gdate.for @@ -0,0 +1,209 @@ +C +C ANALYST : D.GORHAM +C INDEX NO. : UT 95.00 +C TITLE : EXTENDED GREGORIAN DATE SUBROUTINE FROM YEAR 1/1/1 +C DATE : 08/03/82 +C REV.01 : 19-SEP-83 BY DJG: MODIFIED FOR VAX +C +C THIS SUBROUTINE IS ADAPTED BY DJG FROM THAT IN HP COMMUNICATOR +C 1981 VOL. 5 ISSUE 3 WITH INTEGER*4 TYPE BEING CHANGED TO REAL +C PLUS ASSOCIATED CHANGES. +C +C PARAMETER LIST : +C OPT = 1 TO CONVERT CALENDAR DATE TO EXTENDED GREGORIAN DATE (EGD) +C = 2 TO CONVERT EXTENDED GREGORIAN DATE (EGD) TO CALENDAR DATE +C MON = MONTH (1-12) +C DAY = DAY OF MONTH +C YEAR= CALENDAR YEAR EG. 1982 +C DWK = DAY OF WEEK (1=SUNDAY, 7=SATURDAY) +C DYR = DAY OF YEAR +C DMON= DAYS IN MONTH +C EGD = EXTENDED GREGORIAN DATE SINCE 1/1/1 +C IER = 0 FOR NO ERRORS +C 1 IF ERRORS DETECTED +C LU = 0 FOR NO ERROR MESSAGE PRINTOUT +C = N FOR DISPLAY OF ERROR MESSAGE ON LOGICAL UNIT 'N' +C +C TITLE TRUE GEORGIAN DATE ROUTINE ISS. 2 801126 R.A.G. (MWR) +C + SUBROUTINE GDATE (OPT,MON,DAY,YEAR,DWK,DYR,DMON,EGD,IER,LU) +C +C GREGORIAN DATE ROUTINE 801126 +C + IMPLICIT INTEGER*4 (A-Z) +C +C +C -------------------------------------------------------------- +C +C REVISION LIST +C +C --DATE-- ---BY--- --DESCRIPTION-- +C +C 19/03/80 R.A.G. -ORIGINAL ISSUE +C 26/11/80 R.A.G. -CONVERTED REAL NUMBERS TO DOUBLE INTEGER AND +C ADDED DAY OF THE WEEK, DAY OF THE YEAR (JULIAN +C DATE), AND NUMBER OF DAYS IN THE MONTH. +C +C --------------------------------------------------------------- +C +C + INTEGER*4 M(12) +C + LOGICAL LYEAR +C + DATA M/31,28,31,30,31,30,31,31,30,31,30,31/ +C +C INITIALIZE SOME VARIABLES +C + M(2)=28 + DYR=0 + IER=0 +C +C CHECK FOR OPTION +C 1=CALLER SUPPLIES DAY MONTH YEAR +C >1=CALLER SUPPLIES EXTENDED GREGORIAN DATE +C + IF (OPT.GT.1) GO TO 120 +C +C THIS SECTION CONVERTS TO A GREGORIAN DATE +C +C TEST ARGUMENTS FOR VALIDITY +C + IF (MON.LT.1.OR.MON.GT.12.OR.DAY.LT.1.OR.DAY.GT.31. + 1 OR.YEAR.LT.1) THEN +cdrc TYPE 5000 + PRINT 5000 +5000 FORMAT(' **GDATE ERROR') +cdrc TYPE '(A,I6)',' MONTH =',MON +cdrc TYPE '(A,I6)',' DAY =',DAY +cdrc TYPE '(A,I6)',' YEAR =',YEAR + PRINT '(A,I6)',' MONTH =',MON + PRINT '(A,I6)',' DAY =',DAY + PRINT '(A,I6)',' YEAR =',YEAR + IER=1 + GO TO 180 + END IF +C + IF (LYEAR(YEAR)) M(2)=29 +C + IF (DAY.GT.M(MON)) THEN +cdrc TYPE 5000 +cdrc TYPE '(A,I6)',' DAY =',DAY + PRINT 5000 + PRINT '(A,I6)',' DAY =',DAY + IER=2 + GO TO 180 + END IF +C +C CALCULATE GREG. DATE TO 1ST OF REQUESTED YEAR +C + Y=YEAR-1 + EGD=GDATS(Y) +C +C CALCULATE TO CURRENT GREGORIAN DATE +C + J=MON-1 + IF (J.EQ.0) GO TO 110 + DO 100 I=1,J + DYR=DYR+M(I) + 100 CONTINUE + EGD=EGD+DYR + 110 EGD=EGD+DAY + DYR=DYR+DAY + GO TO 170 +C +C THIS SECTION CONVERTS FROM A GREGORIAN DATE +C + 120 IF (EGD.LT.1) THEN + IER=3 +cdrc TYPE 5000 +cdrc TYPE '(A,I6)',' EGD =',EGD + PRINT 5000 + PRINT '(A,I6)',' EGD =',EGD + GO TO 180 + END IF +C +C CALCULATE CURRENT DATE (DD/MM/YYYY) +C + YEAR=(EGD/366)-1 + 130 YEAR=YEAR+1 + EEGD=GDATS(YEAR) + IF (EGD-EEGD-368) 140,140,130 + 140 YEAR=YEAR+1 + DYR=EGD-EEGD +C + IF (LYEAR(YEAR)) M(2)=29 +C + DO 150 MON=1,12 + EEGD=EEGD+M(MON) + IF (EGD.LE.EEGD) GO TO 160 + 150 CONTINUE + M(2)=28 + GO TO 140 +C +C CALCULATE THE REMAINING ARGUMENTS +C + 160 DAY=EGD+M(MON)-EEGD + 170 DMON=M(MON) + DWK=MOD(EGD,7)+1 + RETURN +C +180 RETURN + END +C +C +C TITLE FUNCTION LYEAR(LEAP YEAR) ISS. 1 801126 R.A.G. (MWR) + FUNCTION LYEAR(YEAR) +C + IMPLICIT INTEGER*4 (A-Z) +C +C THIS FUNCTION WILL TEST A GIVEN YEAR AND RETURN A TRUE/FALSE +C INDICATION +C +C -------------------------------------------------------------- +C +C REVISION LIST +C +C --DATE-- ---BY--- --DESCRIPTION-- +C +C 26/11/80 R.A.G. -ORIGINAL ISSUE +C +C -------------------------------------------------------------- +C + LYEAR=0 + IF (MOD(YEAR,4).EQ.0.AND.MOD(YEAR,100).NE.0.OR. + & MOD(YEAR,400).EQ.0) LYEAR= -1 + RETURN + END +C TITLE FUNCTION GDATS(GDATE) ISS. 1 801126 R.A.G. (MWR) + FUNCTION GDATS(YEAR) +C +C + IMPLICIT INTEGER*4 (A-Z) + +C +C +C THIS FUNCTION IS PART OF THE GDATE SUBROUTINE, AND RETURNS A +C GREGORIAN DATE TO THE 1ST OF THE YEAR BASED ON THAT YEAR +C +C ----------------------------------------------------------------- +C +C REVISION LIST +C +C --DATE-- ---BY--- --DESCRIPTION +C +C 26/11/80 R.A.G. -ORIGINAL ISSUE +C +C -------------------------------------------------------------------- +C +C **NOTE** +C +C GDATS MUST BE DECLARED AS AN INTEGER*4 FUNCTION!!!!!!! +C + CON=365 + GDATS=CON*YEAR + GDATS=GDATS+(24*(YEAR/100)) + GDATS=GDATS+(YEAR/400) + GDATS=GDATS+(MOD(YEAR,100)/4) + RETURN + END diff --git a/source/heatex.for b/source/heatex.for new file mode 100644 index 0000000..8a1f9cf --- /dev/null +++ b/source/heatex.for @@ -0,0 +1,364 @@ + SUBROUTINE HEATEX(TOFDAY,DAYOFY,DELT,TSOLHR) +C +C HEATEX COMPUTES THE NET AMOUNT OF HEAT +C RADIATION FLUX BEING TRANSFERRED ACROSS +C THE AIR-WATER INTERFACE BASED ON AN +C ENERGY BUDGET WHICH CONSIDERS SOLAR +C RADIATION, ATMOSPHERIC RADIATION, BACK +C RADIATION, CONDUCTION, AND EVAPORATION. +C SAVE + + INTEGER ITIME,INWDAY,METRIC + REAL CON1,CON2,CON3,CON4,CON5,CON6,DELTSL + REAL SOLCON,TSOLHR,HA,CS,SOLAR,ALPHT + REAL ELEVAT,LAT,LONG + REAL DELTH,DELT,TOFDAY,DAYOFY,TLEFT + REAL DECLIN,REARTH,EQTIME,DECLON + REAL PI,STR,STB,STE,STS,ST1,ST2 + + COMMON/INPUTS/LAT,LONG,ELEVAT,STANM + COMMON/MET/WETBLB,ATMPR,DRYBLB,CLOUD,DAT + +C For metric units, radiation is in kJ/m2*hr + METRIC=1 + LOUT=6 +CC +CCC + +CCC NCASI Commentary, HEATEX Section A. (QUAL2 Step 1-0) +CCC A. Compute and/or define required constants. +CCC +C +CCC NCASI Commentary, HEATEX Section F. +CCC F. Test for beginning of a new day. +CCC +CJFD cheap fix to prevent incrementing time for each material type ... +CJFD jump passed time calcs, this only works because local values are saved +C +C DATA ITIME/0/ +C +CIPK MAY94 CHANGED ORDER +C + ITIME=0 + IF(ITIME .EQ. 0) THEN + PI=4.0*atan(1.) +c PI=3.141628 + CON1=2.0*PI/365.0 + CON2=PI/180.0*LAT + CON3=180.0/PI + CON4=23.45*PI/180.0 + CON5=PI/12.0 + CON6=12.0/PI + DELTSL=(LONG-STANM)/15.0 +CIPK MAR94 + IF(METRIC .EQ. 0) THEN + SOLCON=438.0 + ELEXP=EXP(-ELEVAT/2532.0) + ELSE + SOLCON=4974.4 +c SOLCON=4870.8 !Original value incorrect DRC +c ELEXP=EXP(-ELEV*3.2808/2532.0) + ELEXP=EXP(-ELEVAT/771.76) + ENDIF +CIPK MAR94 END CHANGES + ENDIF + TLEFT=0. + TSOLHR=0. + 55 CONTINUE + ITIME=ITIME+1 + IF(ITIME .EQ. 1) THEN + INWDAY=1 + ELSE + INWDAY=0 + ENDIF +C INWDAY =1 signals a new day or start of simulation +C +C Get Delta Time in hours +C + DELTH = DELT/3600.0 + TOFDAY = TOFDAY+DELTH +C + IF(DELTH .EQ. 0.) THEN + DELTH=1.0 + TOFDAY=0.00 + TLEFT=24. + ENDIF +C +C write(*,*) TOFDAY, dayofy +c WRITE(IOT,*) 'TOFDAY DAYOFY',TOFDAY,DAYOFY + 56 CONTINUE + IF(DELT .EQ. 0.) THEN + TLEFT=TLEFT-DELTH*1.0001 + TOFDAY=TOFDAY+DELTH + GO TO 85 + ENDIF + IF(TLEFT .GT. 0.0) THEN + TOFDAY=TLEFT + DAYOFY=DAYOFY + 1.0 + INWDAY=1 + IF (TOFDAY .GT. 24.001) THEN + TLEFT=TOFDAY-24.00 + TOFDAY=24.00 + DELTH=24.00 + ELSE + DELTH=TLEFT + TLEFT=0.0 + ENDIF + GO TO 85 + ENDIF + IF (TOFDAY .GT. 24.001) THEN + TLEFT=TOFDAY-24.00 + TOFDAY=24.00 + DELTH=DELTH-TLEFT + INWDAY=0 + ENDIF + 85 CONTINUE +CDRC write(LOUT,*)'at 85 tleft,tofday,delth,inwday',tleft,tofday,delth +CDRC + ,inwday +cDRC CALL GETMET(TOFDAY) + IF(INWDAY .EQ. 1) THEN + +CCC +CCC NCASI Commentary, HEATEX Section B. (QUAL2 Step 2-0) +CCC B. Begin computations for calculating the +CCC net solar radiation term. +CCC +CCC B.1 Test for beginning of a new day. +CCC +CCC +CCC B.1a Calculate seasonal and daily position +CCC of the sun relative to the location +CCC of the basin on the earth's +CCC surface. (QUAL2 Step 2-1) +CCC + REARTH=1.0+0.017*COS(CON1*(186.0-DAYOFY)) + DECLIN=CON4*COS(CON1*(172.0-DAYOFY)) + RR=REARTH**2 + EQTIME=0.000121-0.12319*SIN(CON1*(DAYOFY-1.0)-0.07014) + * -0.16549*SIN(2.0*CON1*(DAYOFY-1.0)+0.3088) + DECLON=ABS(DECLIN) +CC +CC Replace TAN function with SIN/COS. +CC + TANA = SIN(CON2)/COS(CON2) + TANB = SIN(DECLON)/COS(DECLON) + ACS = TANA*TANB +CC + IF (ACS .NE. 0.0) THEN + XX=SQRT(1.0-ACS*ACS) + XX=XX/ACS +cipk oct94 ACS=ATAN(XX) + if(xx .gt. 0.) then + ACS=abs(ATAN(XX)) + IF (DECLIN.GT.0.0) ACS=PI-ACS + else + acs=abs(atan(xx)) + if (declin.lt.0.0) acs=pi-acs + endif + ELSE + ACS=PI/2.0 + ENDIF +CCC +CCC B.1a Calculate the standard time of +CCC sunrise (STR) and sunset (STS). +CCC (QUAL2 Step 2-2) +CCC +C + STR=12.0-CON6*ACS+DELTSL + STS=24.0-STR+2.0*DELTSL + DAYLEN=STS-STR + STB=TOFDAY-DELTH + STE=STB+DELTH +cDRC WRITE(LOUT,*) 'AFTER 2-2 STR,STS,STB,STE',STR,STS,STB,STE +CCC +CCC B.2 Increment the variables that define the +CCC time of the beginning(STB) and the +CCC end (STE) of the time interval. +CCC + ELSE + STB=STB+DELTH + STE=STB+DELTH + ENDIF +CCC +CCC B.3 Test if time to read in local +CCC climatological data. (QUAL2 Step 2-3) +CCC +CJFD IF (TRLCD.NE.1.0) GO TO 82 +CCC +CCC B.7 Compute vapor pressures (VPWB and +CCC VPAIR), dew point (DEWPT), AND +CCC dampening effect of clouds (CNS +CCC and CNL). (QUAL2 Step 2-4) +CCC +CIPK MAR94 + IF(METRIC .EQ. 0) THEN + VPWB=0.1001*EXP(0.03*WETBLB)-0.0837 + VPAIR=VPWB-0.000367*ATMPR*(DRYBLB-WETBLB) + * *(1.0+(WETBLB-32.0)/1571.0) + DEWPT=ALOG((VPAIR+0.0837)/0.1001)/0.03 + ELSE +c VPWB=(0.1001*EXP(0.03*(WETBLB(NN)*1.8+32.))-0.0837)/29.53*1000. +c VPAIR=VPWB-0.000367*ATMPR(NN)*(DRYBLB(NN)-WETBLB(NN))*1.8 +c * *(1.0+(WETBLB(NN)*1.8)/1571.0) +c DEWPT=((ALOG((VPAIR/1000.*29.53+0.0837)/0.1001)/0.03)-32.0)/1.8 +c IF(NN .EQ. 1) WRITE(75,*) 'VPWB',VPWB,VPAIR,DEWPT + vpwb=8.8534*exp(0.054*wetblb)-2.8345 + VPAIR=VPWB-0.0006606*ATMPR*(DRYBLB-WETBLB) + * *(1.0+WETBLB/872.78) +c vpwb in millibars +c vpair in millibars + DEWPT=ALOG((VPAIR+2.8345)/8.8534)/0.054 +c dewpt in deg C +c IF(NN .EQ. 1) WRITE(75,*) 'VPWB',VPWB,VPAIR,DEWPT + ENDIF +CIPK MAR94 END CHANGES + CS=1.0-0.65*CLOUD**2 + IF (CLOUD.GT.0.9) CS=0.50 + CNL=CLOUD*10.0+1.0 + NL=CNL +82 CONTINUE + IF (STS.LE.STB.OR.STR.GE.STE) GO TO 35 +C IF(STR.GT.STB.AND.STR.LT.STE) GO TO 41 +C IF (STS.LT.STE.AND.STS.GT.STB) GO TO 42 + ST1=STB + ST2=STE + IF(STB .LT. STR) ST1=STR + IF(STE .GT. STS) ST2=STS +CCC +CCC NCASI Commentary, HEATEX Section C. (QUAL2 Step 2-5) +CCC C. Continue with calculations for solar +CCC radiation. +CCC +CCC C.1 Calculate hour angles (TB and TE). +CCC +C TB=STB-12.0-DELTSL+EQTIME +C TE=STE-12.0-DELTSL+EQTIME +C WRITE(LOUT,*) '40 ,DELTSL,EQTIME,STB,STE,TB,TE', +C +DELTSL,EQTIME,STB,STE,TB,TE +C GO TO 43 +C 41 TB=STR-12.0-DELTSL+EQTIME +C TE=STE-12.0-DELTSL+EQTIME +C WRITE(LOUT,*) '41 ,DELTSL,EQTIME,STR,STE,TB,TE', +C +DELTSL,EQTIME,STR,STE,TB,TE +C GO TO 43 +C 42 TB=STB-12.0-DELTSL+EQTIME +C TE=STS-12.0-DELTSL+EQTIME +C WRITE(LOUT,*) '42 ,DELTSL,EQTIME,STB,STS,TB,TE', +C +DELTSL,EQTIME,STB,STR,TB,TE +C 43 CONTINUE + TB=ST1-12.0-DELTSL+EQTIME + TE=ST2-12.0-DELTSL+EQTIME +CDRC WRITE(LOUT,*) '43 ,DELTSL,EQTIME,ST1,ST1,TB,TE', +CDRC +DELTSL,EQTIME,ST1,ST2,TB,TE + + TALT=(TB+TE)/2.0 +CDRC WRITE(LOUT,*) 'TB,TE',TB,TE +CCC +CCC C.2 Compute amount of clear sky, solar +CCC radiation(SOLAR), and altitude of +CCC the sun (ALPHT). (QUAL2 Step 2-6) +CCC + + SOLAR=SOLCON/RR*(SIN(CON2)*SIN(DECLIN)*(TE-TB)+CON6*COS(CON2)* + * COS(DECLIN)*(SIN(CON5*TE)-SIN(CON5*TB))) + ALPHT=SIN(CON2)*SIN(DECLIN)+COS(CON2)*COS(DECLIN)*COS(CON5*TALT) + IF (ABS(ALPHT).EQ.1.0) GO TO 4 + Y=SQRT(1.0-ALPHT*ALPHT) + Y=ALPHT/Y + ALPHT=ATAN(Y) + GO TO 5 + 4 IF (ALPHT.EQ.-1.0) GO TO 6 + ALPHT=PI/2.0 + GO TO 5 + 6 ALPHT=-PI/2.0 + 5 CONTINUE +CDRC write(LOUT,*) 'alpht',alpht + IF (ALPHT.LT.0.01) GO TO 35 +CCC +CCC C.3 Compute absorption and scattering due +CCC to atmospheric conditions. (QUAL2 +CCC Step 2-7) +CCC +CIPK MAR94 + IF(METRIC .EQ. 0) THEN + PWC=0.00614*EXP(0.0489*DEWPT) + ELSE +c PWC=0.00614*EXP(0.0489*(DEWPT*1.8+32.0)) +c IF(NN.EQ. 1) WRITE(75,*) 'PWC',PWC + PWC=0.02936*EXP(0.08802*DEWPT) +c IF(NN.EQ. 1) WRITE(75,*) 'PWC',PWC + ENDIF +CIPK MAR94 END CHANGES + OAM=ELEXP/(SIN(ALPHT)+0.15*(ALPHT*CON3+3.885)**(-1.253)) + A1=EXP(-(0.465+0.0408*PWC)*(0.129+0.171*EXP(-0.880*OAM))*OAM) + A2=EXP(-(0.465+0.0408*PWC)*(0.179+0.421*EXP(-0.721*OAM))*OAM) +CCC +CCC C.4 Compute reflectivity coefficient (RS). +CCC (QUAL2 Step 2-8) +CCC + GO TO (30,31,31,31,31,31,32,32,32,32,33), NL + 30 AR=1.18 + BR=-0.77 + GO TO 34 + 31 AR=2.20 + BR=-0.97 + GO TO 34 + 32 AR=0.95 + BR=-0.75 + GO TO 34 + 33 AR=0.35 + BR=-0.45 + 34 CONTINUE + RS=AR*(CON3*ALPHT)**BR +CDRC write(LOUT,*) 'rs',rs +CC +CC Add test for RS greater than 1.0. +CC + IF(RS.GE.1.0) GO TO 35 +CC +CCC +CCC C.5 Compute atmospheric transmission term (ATC). +CCC + ATC=(A2+0.5*(1.0-A1-DAT))/(1.0-0.5*RS*(1.0-A1+DAT)) +CCC +CCC C.6 Compute net solar radiaiont for the time +CCC interval delta t. (QUAL2 Step 2-9) +CCC +CDRC write(LOUT,*) 'solar,atc,cs',solar,atc,cs + TSOLHR = TSOLHR+SOLAR*ATC*CS*(1.0-RS) + GO TO 36 + 35 TSOLHR = TSOLHR+0.0 + 36 CONTINUE + CLC=1.0+0.17*CLOUD**2 +CCC +CCC NCASI Commentary, HEATEX Section D. (QUAL2 Step 3-0) +CCC D. Compute heat fluxes from other terms. +CCC +CCC D.1 Long wave atmospheric radiation (HA). +CCC +CIPK MAR94 + IF(METRIC .EQ. 0) THEN + HA = HA+ + + 0.97*1.73E-09*2.89E-06*(DRYBLB+460.0)**6*CLC*DELTH + ELSE + +c +c HA(NN) = HA(NN) + +c + 0.97*1.73E-09*2.89E-06*(DRYBLB(NN)*1.8+32.+460.0)**6 +c + *CLC*DELTH*4870.8/438. +c IF(NN.EQ. 1) WRITE(75,*) 'HA',HA(NN) + HA = HA + + + 0.97*9.37e-06*2.0412E-07*(DRYBLB+273.0)**6 + + *CLC*DELTH +c IF(NN.EQ. 1) WRITE(75,*) 'HA',HA + ENDIF +CIPK MAR94 END CHANGES +C WRITE(LOUT,*) 'tofday,tsolhr,ha', +C + TOFDAY, TSOLHR, HA, STB, STE, STR, STS + IF(TLEFT .GT. 0.) GO TO 56 + +CJFD cheap fix to prevent incrementing time for each material type ... +CC + DELTH = DELT/3600. + RETURN + END diff --git a/source/lenstr.for b/source/lenstr.for new file mode 100644 index 0000000..6059a58 --- /dev/null +++ b/source/lenstr.for @@ -0,0 +1,18 @@ + + integer function lenstr(string) + character*(*) string + integer val,min,max + + min=ichar('!') + max=ichar('z') + lenstr = len(string) + val=ichar(string(lenstr:lenstr)) + do while(((val.lt.min).or.(val.gt.max)).and.(lenstr.gt.0)) + lenstr = lenstr - 1 + val=ichar(string(lenstr:lenstr)) + end do + return + end + + + diff --git a/source/random.for b/source/random.for new file mode 100644 index 0000000..1251fa0 --- /dev/null +++ b/source/random.for @@ -0,0 +1,133 @@ +c FILE : random.f +c ---------------------------------------------------------------------- +c ------ Uniform Distribution Random Number Generator ------------------ +c ---------------------------------------------------------------------- +c From Numerical Recipies. + + FUNCTION ran3(idum) + INTEGER*2 i,ii,k,inext,inextp,iff + INTEGER*4 idum,mbig,mseed,mz,mj,ma(55),mk + REAL*4 fac, ran3 + parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9) + save inext,inextp,ma,iff + data iff /0/ +C + if(idum.lt.0.or.iff.eq.0)then + iff=1 + mj=mseed-iabs(idum) + mj=mod(mj,mbig) + ma(55)=mj + mk=1 + do i=1,54 + ii=mod(21*i,55) + ma(ii)=mk + mk=mj-mk + if(mk.lt.mz) mk=mk+mbig + mj=ma(ii) + end do + do k=1,4 + do i=1,55 + ma(i)=ma(i)-ma(1+mod(i+30,55)) + if(ma(i).lt.mz)ma(i)=ma(i)+mbig + end do + end do + inext=0 + inextp=31 + idum=1 + endif + inext=inext+1 + if(inext.eq.56)inext=1 + inextp=inextp+1 + if(inextp.eq.56)inextp=1 + mj=ma(inext)-ma(inextp) + if(mj.lt.mz)mj=mj+mbig + ma(inext)=mj + ran3=real(mj)*fac + return + end + + + FUNCTION ran4(idum) + INTEGER*2 i,ii,k,inext,inextp,iff + INTEGER*4 idum,mbig,mseed,mz,mj,ma(55),mk + REAL*4 fac, ran4 + parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9) + save inext,inextp,ma,iff + data iff /0/ +C + if(idum.lt.0.or.iff.eq.0)then + iff=1 + mj=mseed-iabs(idum) + mj=mod(mj,mbig) + ma(55)=mj + mk=1 + do i=1,54 + ii=mod(21*i,55) + ma(ii)=mk + mk=mj-mk + if(mk.lt.mz) mk=mk+mbig + mj=ma(ii) + end do + do k=1,4 + do i=1,55 + ma(i)=ma(i)-ma(1+mod(i+30,55)) + if(ma(i).lt.mz)ma(i)=ma(i)+mbig + end do + end do + inext=0 + inextp=31 + idum=1 + endif + inext=inext+1 + if(inext.eq.56)inext=1 + inextp=inextp+1 + if(inextp.eq.56)inextp=1 + mj=ma(inext)-ma(inextp) + if(mj.lt.mz)mj=mj+mbig + ma(inext)=mj + ran4=real(mj)*fac + return + end + + +c --------------------------------------------------------------------- +c -------- Normally Distributed Random number generator --------------- +c --------------------------------------------------------------------- +c From numerical recipies +c Mean=0, Stdev=1 + + FUNCTION normdis(idum) + INTEGER iset, idum + REAL v1, v2, r, fac, normdis, ran4, gset + SAVE gset, iset + DATA iset/0/ +c idum=0 + IF (iset.eq.0) THEN +1 v1=2.*RAN4(idum)-1. + v2=2.*RAN4(idum)-1. + r=v1**2+v2**2 + IF (r.ge.1.) GOTO 1 + fac=SQRT(-2.*LOG(r)/r) + gset=v1*fac + normdis=v2*fac + iset=1 + ELSE + normdis=gset + iset=0 + END IF + RETURN + END + +c --------------------------------------------------------------- + + FUNCTION rbit() + INTEGER rbit, isd + REAL r, ran3 + r=ran3(isd) + IF (r.le.0.5) THEN + rbit=1 + ELSE + rbit=-1 + ENDIF + RETURN + END diff --git a/source/rise_fall.cb b/source/rise_fall.cb new file mode 100644 index 0000000..1127c08 --- /dev/null +++ b/source/rise_fall.cb @@ -0,0 +1,2 @@ + real y0r,y1r,y0f,y1f,rr,rf + common /risefall/y0r,y1r,y0f,y1f,rr,rf diff --git a/source/rise_fall_vel.for b/source/rise_fall_vel.for new file mode 100644 index 0000000..dc1604d --- /dev/null +++ b/source/rise_fall_vel.for @@ -0,0 +1,44 @@ +c--------------------------------------------------------rise_fall_vel + subroutine rise_fall_vel(vv) +c--------------------------------------------------------------------c +c vv > 0 : rise and vv < 0 : fall vv = 0.0 : nuetral c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c This program is to assign a velocity to one of the c +c sewage particles which have various sizes and densities and c +c paricle size and density are uniformly distributed among c +c the particles c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c +c yo and y1 are coefficients of speed regression equation +c which has form 10.0**((m-y0)/y1) +c particles which fall have possitive velocity +c particles which rise have negantive velocity +c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c + include 'rise_fall.cb' + REAL rnumb, vv, ran3, fspe, rspe + common /randm/isd,isd1 + INTEGER isd,isd1 +c + rnumb=100.0*ran3(isd) + vv=0.0 +c + if(rnumb.le.rf) vv=fspe(rnumb) + if(rnumb.ge.(100.0-rr)) vv=rspe(rnumb) + return + end + + function fspe(rnumb) + include 'rise_fall.cb' + real fspe, rnumb + fspe=-10.0**((rnumb-y0f)/y1f)/3600.0 + return + end + + function rspe(rnumb) + include 'rise_fall.cb' + real rspe, rnumb + rspe=10.0**(((100.0-rnumb)-y0r)/y1r)/3600.0 + return + end diff --git a/source/tinc2.for b/source/tinc2.for new file mode 100644 index 0000000..e34403e --- /dev/null +++ b/source/tinc2.for @@ -0,0 +1,49 @@ +C + SUBROUTINE TINC2(IHRMNS,IDAYS,MONS,IYEARS,IHRMNF,IDAYF,MONF, + & IYEARF,RMIN,IER,LU) +C +C AUTHOR : D.GORHAM +C DATE : 4-MAR-83 +C TITLE : GIVEN START TIME/DATE AND NO. OF MINUTE INTERVALS +C FINISH TIME/DATE +C REV.01 : 22-SEP-83 BY DJG: MODIFIED FOR VAX +C REV.02 : 21-FEB-84 BY DJG: MODIFIED SUCH THAT IF 'RMIN' IS +C NEGATIVE CORRECT TIME IS ALWAYS RETURNED +C + INTEGER*4 EGDS,EGDF + CALL GDATE(1,MONS,IDAYS,IYEARS,IDWK,IDYR,IDMON,EGDS,IER,LU) + IF(IER.EQ.0) GO TO 10 + IF(LU.EQ.0) GO TO 10 + WRITE(LU,1000) +1000 FORMAT(' **ERROR CALLING START DATE/TIME IN ROUTINE TINC2') +C +10 IHRS=IHRMNS/100 + MINS=IHRMNS-IHRS*100 + TMINS=RMIN+MINS+IHRS*60 + ITMINS=TMINS + NDAYS=ITMINS/1440 + IREM=ITMINS-NDAYS*1440 + NHRS=IREM/60 + NMINS=IREM-NHRS*60 +C + IF(NMINS.LT.0) THEN + NHRS=NHRS-1 + NMINS=60+NMINS + END IF +C + IF(NHRS.LT.0) THEN + NDAYS=NDAYS-1 + NHRS=24+NHRS + END IF +C + EGDF=EGDS+NDAYS + IHRMNF=NHRS*100+NMINS + print *,rmin,nhrs,nmins +C + CALL GDATE(2,MONF,IDAYF,IYEARF,IDWK,IDYR,IDMON,EGDF,IER,LU)!M.W.S. + IF(IER.EQ.0) GO TO 40 +C + WRITE(LU,1100) +1100 FORMAT(' **ERROR CALLING FINISH DATE/TIME IN ROUTINE TINC2') +40 RETURN + END diff --git a/source/trinvs.for b/source/trinvs.for new file mode 100644 index 0000000..f3f1a68 --- /dev/null +++ b/source/trinvs.for @@ -0,0 +1,713 @@ +c----------------------------------------------------------------trinvs3 + subroutine trinvs3(nel,cno,itype,ebox,egrp,nblk,egps,nelb,nbb, + * nee,nxyz,nen,it,x,y,z,xi,eta,zeta,num,nn) +c----------------------------------------------------------------------c +c purpose: c +c To locate the element the coordinate belongs to and to c +c calculate the shape functions. c +c----------------------------------------------------------------------c +ccc implicit real*8 (a-h,o-z) + common /etype/inode(4) + common /shape/shap(20),shpx(20),shpy(20),shpz(20) + integer*2 nel(20,1) + dimension pos(3),cno(3,1),nbb(1),nee(1),itype(1),ebox(6,1), + * nelb(1),egps(6,1),nblk(2,1),egrp(6,1) + real shap,xi,eta,zeta +c + pos(1)=x + pos(2)=y + pos(3)=z + nn=0 + if(num.gt.0) then +c +c----------box test +c + if(pos(1).lt.ebox(1,num).or.pos(1).gt.ebox(4,num)) goto 1 + if(pos(2).lt.ebox(2,num).or.pos(2).gt.ebox(5,num)) goto 1 + if(pos(3).lt.ebox(3,num).or.pos(3).gt.ebox(6,num)) goto 1 + it=itype(num) + nen=inode(it) + call invs_el3(nel,cno,pos,xi,eta,zeta,nen,it,num,key,ii) + if(key.eq.1) goto 30 + endif +c call trinvs1(nel,cno,egrp,nblk,egps,nelb,nbb,nee,nxyz,nen, +c * it,x,y,z,xi,eta,zeta,num,nn) +c if(nn.eq.1) goto 30 +c goto 35 + 1 do n=1,nxyz + if(pos(1).gt.egrp(4,n).or.pos(1).lt.egrp(1,n)) goto 10 + if(pos(2).gt.egrp(5,n).or.pos(2).lt.egrp(2,n)) goto 10 + if(pos(3).gt.egrp(6,n).or.pos(3).lt.egrp(3,n)) goto 10 + if(nblk(1,n).gt.1) then + do j=nblk(2,n),nblk(2,n)+nblk(1,n)-1 + if(pos(1).gt.egps(4,j).or.pos(1).lt.egps(1,j)) goto 5 + if(pos(2).gt.egps(5,j).or.pos(2).lt.egps(2,j)) goto 5 + if(pos(3).gt.egps(6,j).or.pos(3).lt.egps(3,j)) goto 5 + ng=j + goto 20 + 5 enddo + else + ng=nblk(2,n) + goto 20 + endif + 10 enddo +c +c------out of the finite element mesh +c + goto 900 + + 20 do 35 i=nbb(ng),nee(ng) + num=nelb(i) +c +c----------box test +c + if(pos(1).lt.ebox(1,num).or.pos(1).gt.ebox(4,num)) goto 35 + if(pos(2).lt.ebox(2,num).or.pos(2).gt.ebox(5,num)) goto 35 + if(pos(3).lt.ebox(3,num).or.pos(3).gt.ebox(6,num)) goto 35 + it=itype(num) + nen=inode(it) + call invs_el3(nel,cno,pos,xi,eta,zeta,nen,it,num,key,ii) + if(key.eq.1) goto 30 + 35 continue +c write(*,'(/a/)') 'The input coordinate is out of the FE mesh.' + goto 900 + 30 nn=1 + 900 return + 1000 format(2i10) + 1100 format(i10,3f10.0) + 1200 format(21i5) + 1300 format(i6,6e12.4,i6) + end +c--------------------------------------------------------------invs_el3 + subroutine invs_el3(nel,cno,pos,xi,eta,zeta,nen,it,num,key,ii) +c---------------------------------------------------------------------c +c Purpose: c +c To perform inverse mapping for element i. c +c---------------------------------------------------------------------c + common /shape/shap(20),shpx(20),shpy(20),shpz(20) + integer*2 nel(20,1) + dimension nele(20),cnod(3,20),pos(3),cno(3,1) + real shap,xi,eta,zeta + real a(3,3),b(3),a1,a2,a3 + real*8 detj,ass +c + do 40 j=1,nen + nele(j)=nel(j,num) + cnod(1,j)=cno(1,nele(j)) + cnod(2,j)=cno(2,nele(j)) + cnod(3,j)=cno(3,nele(j)) + 40 continue +c +c------Newton-Raphson iteration method +c + ii=0 + xi0=0.0 + eta0=0.0 + zeta0=0.0 +c +c------calculate the values of shape funcations and their derivatives +c + 1 ii=ii+1 + call xn3x(it,nen,xi0,eta0,zeta0) +c +c------calculate the Jacobian matrix +c + do 20 j=1,3 + a(j,1)=0.0 + a(j,2)=0.0 + a(j,3)=0.0 + b(j)=-pos(j) + do 20 i=1,nen + a(j,1)=a(j,1)+shpx(i)*cnod(j,i) + a(j,2)=a(j,2)+shpy(i)*cnod(j,i) + a(j,3)=a(j,3)+shpz(i)*cnod(j,i) + b(j)=b(j)+shap(i)*cnod(j,i) + 20 continue + + if(abs(b(1)).lt.1.0e-5.and.abs(b(2)).lt.1.0e-5.and. + & abs(b(3)).lt.1.0e-5) goto 30 + + do i=1,3 + ass=max(abs(a(i,1)),abs(a(i,2)),abs(a(i,3))) + if(dabs(ass).lt.1.0e-10) goto 25 + do j=1,3 + a(i,j)=a(i,j)/ass + enddo + b(i)=b(i)/ass + 25 enddo + a1=a(2,2)*a(3,3)-a(2,3)*a(3,2) + a2=a(3,2)*a(1,3)-a(3,3)*a(1,2) + a3=a(1,2)*a(2,3)-a(1,3)*a(2,2) + detj=a(1,1)*a1+a(2,1)*a2+a(3,1)*a3 + if(dabs(detj).lt.1.0e-10) goto 900 + xi=xi0-(a1*b(1)+a2*b(2)+a3*b(3))/detj + a1=a(3,1)*a(2,3)-a(3,3)*a(2,1) + a2=a(1,1)*a(3,3)-a(1,3)*a(3,1) + a3=a(2,1)*a(1,3)-a(2,3)*a(1,1) + eta=eta0-(a1*b(1)+a2*b(2)+a3*b(3))/detj + a1=a(2,1)*a(3,2)-a(2,2)*a(3,1) + a2=a(3,1)*a(1,2)-a(3,2)*a(1,1) + a3=a(1,1)*a(2,2)-a(1,2)*a(2,1) + zeta=zeta0-(a1*b(1)+a2*b(2)+a3*b(3))/detj + if(abs(eta-eta0).le.0.5e-4.and.abs(xi-xi0).le.0.5e-4.and. + * abs(zeta-zeta0).le.0.5e-4) goto 30 + if(ii.gt.11) goto 30 + xi0=xi + eta0=eta + zeta0=zeta + goto 1 +c 900 write(6,*) 'Jacobian Determinant = 0',num,nen,it,xi0,eta0,zeta0, +c * detj +c stop + 900 goto 35 +c goto 35 +c +c------determine whether it is inside the element +c + 30 goto (50,60,70,80), it + 50 if(xi.lt.-5.0e-3.or.eta.lt.-5.0e-3.or.zeta.lt.-5.0e-3) + * goto 35 + if((xi+eta+zeta).gt.1.005) goto 35 + goto 90 + 60 if(xi.lt.-5.0e-3.or.xi.gt.1.005) goto 35 + if(eta.lt.-1.005.or.eta.gt.1.005) goto 35 + if(zeta.lt.-1.005.or.zeta.gt.1.005) goto 35 + if((xi+abs(eta)).gt.1.005) goto 35 + goto 90 + 70 if(xi.lt.-1.005.or.xi.gt.1.005) goto 35 + if(eta.lt.-1.005.or.eta.gt.1.005) goto 35 + if(zeta.lt.-1.005.or.zeta.gt.1.005) goto 35 + goto 90 + 80 if(xi.lt.-5.0e-3.or.xi.gt.1.005) goto 35 + if(eta.lt.-1.005.or.eta.gt.1.005) goto 35 + if(zeta.lt.-1.005.or.zeta.gt.1.005) goto 35 + if((xi+abs(eta)).gt.1.005) goto 35 + if((xi+abs(zeta)).gt.1.005) goto 35 + 90 key=1 + return + 35 key=0 + return + end + +c-------------------------------------------------------------------xn3 + subroutine xn3(it,nen,x,y,z) +c---------------------------------------------------------------------c +c +c +c subroutine to evaluate shape function +c +c IT is element type 1 10 point tetrahedron +c type 2 15 point prism +c type 3 20 point parallelipiped +c type 4 13 point rectangular base pyramid +c nen is number of shape functions +c X, Y, Z are cordinates of point to be evaluated in local coord +c +c---------------------------------------------------------------------c +ccc implicit real*8 (a-h,o-z) + save + common /shape/shap(20),shpx(20),shpy(20),shpz(20) + dimension irf(10,2) + dimension a(5),b(5),c(5),xt(3),yt(3),jx(3),kx(3),xm(5) + 1 ,sn(2),shpp(15,4),ilokup(13) +c + data ilokup/5,6,1,2,3,4,9,7,14,15,10,11,13/ + data irf/1,1,2,2,3,3,1,2,3,4 + 1 ,1,2,2,3,3,1,4,4,4,4/ + data xt/0.0,1.0,0.0/,yt/-1.0,0.0,1.0/,jx/2,3,1/,kx/3,1,2/ + data shpp/0.0,0.0,1.0,12*0.0, + + 0.25,-1.0,3.0,-1.0,0.25,4*0.0,0.25,-1.0,0.0,-1.0,0.25,0.0, + + 0.25,-1.0,0.0,1.0,-0.25,4*0.0,0.25,-1.0,0.0,1.0,-0.25,0.0, + + 0.25,-1.0,0.0,-1.0,0.25,4*0.0,-0.25,1.0,0.0,1.0,-0.25,0.0/ + data ncall/0/ +c + if(it .eq. 4) go to 80 + if( it - 2 ) 500,80,300 +c- +c-----shape functions for right prism and pyramid..... +c- + 80 do j0=1,nen + if(it.eq.2) then + i=j0 + else + i=ilokup(j0) + endif + ncall = ncall + 1 + if( ncall .gt. 1 ) go to 125 +c- +c-----calculate invarient triangular functions..... +c- + n = 0 + do 100 j = 1,5,2 + n = n + 1 + jj = jx(n) + kk = kx(n) + a(j) = xt(jj)*yt(kk) - xt(kk)*yt(jj) + b(j) = yt(jj) - yt(kk) + c(j) = xt(kk) - xt(jj) + 100 continue +c- +c-.....shape function calculations..... +c- + 125 do 130 k = 1, 5, 2 + xm(k) = 0.5*(a(k)+b(k)*x+c(k)*y) + 130 continue +c- +c-.....set indexes and normalize coordinates..... +c- + if( i .le. 6 ) then + zooo = -z + else + zooo = z + endif + if(it .eq. 2) go to 132 + IF(X .EQ. 1.0) GO TO 250 +c if(x .eq. 1.0) x=0.999 + if(abs(1.0-x) .le. 1.0e-6) then +c print *,'x=',x + goto 250 + endif + zooo=zooo/(1.-x) +c if(x .eq. 1.) zooo=sign(x,zooo) + if(i .eq. 3) go to 200 + 132 if( i .gt. 6 .and. i .lt. 10 ) go to 170 + l = i + if( i .gt. 9 ) l = i - 9 + if( mod(l,2) .eq. 0 ) go to 150 +c- +c-.....corner nodes..... +c- + shap(j0)=xm(l)*((2.0*xm(l)-1.0)*(1.0+zooo)-(1.0-zooo**2))/2.0 + if(it .eq. 4) shap(j0)=shap(j0)+xm(l)/2.0*x*(1.-zooo**2) + goto 240 +c- +c-.....mid side node..... +c- + 150 n1 = l - 1 + n2 = mod(l+1,6) + shap(j0) = 2.0*xm(n1)*xm(n2)*(1.0+zooo) + goto 240 +c- +c-.....mid side rectangle..... +c- + 170 n1 = 1 + if( i .eq. 8 ) n1 = 3 + if( i .eq. 9 ) n1 = 5 + shap(j0) = xm(n1)*(1.0-zooo**2) + if(it .eq. 4) shap(j0)=shap(j0)*(1.-x) + 180 goto 240 +c- +c......special case for no 3 shape function of pyramid +c- + 200 l=i + shap(j0)=xm(l)*(2.*xm(l)-1.0) + goto 240 +c +c-----Special case when x exactly equals 1.0 +c + 250 continue + shap(j0)=shpp(i,1) + 240 enddo + return +c- +c-.....shape functions for cube..... +c- + 300 zm1=0.125*(1.0-z) + zp1=0.125*(1.0+z) +c- +c......corner nodes +c- + shap(1)=(1.-x)*(1.-y)*zm1*(-x-y-z-2.) + shap(3)=(1.+x)*(1.-y)*zm1*(x-y-z-2.) + shap(5)=(1.+x)*(1.+y)*zm1*(x+y-z-2.) + shap(7)=(1.-x)*(1.+y)*zm1*(-x+y-z-2.) + shap(13)=(1.-x)*(1.-y)*zp1*(-x-y+z-2.) + shap(15)=(1.+x)*(1.-y)*zp1*(x-y+z-2.) + shap(17)=(1.+x)*(1.+y)*zp1*(x+y+z-2.) + shap(19)=(1.-x)*(1.+y)*zp1*(-x+y+z-2.) +c- +c......mid-side nodes +c- + x21=0.25*(1.0-x*x) + y21=0.25*(1.0-y*y) + z21=0.25*(1.0-z*z) + shap(2)=x21*(1.-y)*(1.-z) + shap(4)=(1.+x)*y21*(1.-z) + shap(6)=x21*(1.+y)*(1.-z) + shap(8)=(1.-x)*y21*(1.-z) + shap(9)=(1.-x)*(1.-y)*z21 + shap(10)=(1.+x)*(1.-y)*z21 + shap(11)=(1.+x)*(1.+y)*z21 + shap(12)=(1.-x)*(1.+y)*z21 + shap(14)=x21*(1.-y)*(1.+z) + shap(16)=(1.+x)*y21*(1.+z) + shap(18)=x21*(1.+y)*(1.+z) + shap(20)=(1.-x)*y21*(1.+z) + return + 500 continue +c- +c......section for tetrahedron +c- + do i=1,nen + i1=irf(i,1) + i2=irf(i,2) + ia=i1 + do 550 k=1,2 + go to (505,515,525,535),ia + 505 sn(k)=1.-x-y-z + go to 550 + 515 sn(k)=z + go to 550 + 525 sn(k)=x + go to 550 + 535 sn(k)=y + 550 ia=i2 +c- +c......test for corner nodes +c- + if(i1 .ne. i2) then +c- +c......for midsides evaluate function etc +c- + shap(i)=4.*sn(1)*sn(2) + else +c- +c......corner node evaluation +c- + shap(i)=sn(1)*(2.*sn(1)-1.) + endif + enddo +c- +c......final step +c- + return + end + +c------------------------------------------------------------------xn3x + subroutine xn3x(it,nen,x,y,z) +c---------------------------------------------------------------------c +c +c +c subroutine to evaluate shape function and its derivatives +c for three dimensions +c +c IT is element type 1 10 point tetrahedron +c type 2 15 point prism +c type 3 20 point parallelipiped +c type 4 13 point rectangular base pyramid +c I is shape function number +c X, Y, Z are cordinates of point to be evaluated in local coord +c +c---------------------------------------------------------------------c +ccc implicit real*8 (a-h,o-z) + save + common /shape/shap(20),shpx(20),shpy(20),shpz(20) + dimension irf(10,2) + dimension a(5),b(5),c(5),xt(3),yt(3),jx(3),kx(3),xm(5) + 1 ,sn(2),snx(2),sny(2),snz(2),shpp(15,4),ilokup(13) +c + data ilokup/5,6,1,2,3,4,9,7,14,15,10,11,13/ + data irf/1,1,2,2,3,3,1,2,3,4 + 1 ,1,2,2,3,3,1,4,4,4,4/ + data xt/0.0,1.0,0.0/,yt/-1.0,0.0,1.0/,jx/2,3,1/,kx/3,1,2/ + data shpp/0.0,0.0,1.0,12*0.0, + + 0.25,-1.0,3.0,-1.0,0.25,4*0.0,0.25,-1.0,0.0,-1.0,0.25,0.0, + + 0.25,-1.0,0.0,1.0,-0.25,4*0.0,0.25,-1.0,0.0,1.0,-0.25,0.0, + + 0.25,-1.0,0.0,-1.0,0.25,4*0.0,-0.25,1.0,0.0,1.0,-0.25,0.0/ + data ncall/0/ +c + if(it .eq. 4) go to 80 + if( it - 2 ) 500,80,300 +c- +c-----shape functions for right prism and pyramid..... +c- + 80 do j0=1,nen + if(it.eq.2) then + i=j0 + else + i=ilokup(j0) + endif + ncall = ncall + 1 + if( ncall .gt. 1 ) go to 125 +c- +c-----calculate invarient triangular functions..... +c- + n = 0 + do 100 j = 1,5,2 + n = n + 1 + jj = jx(n) + kk = kx(n) + a(j) = xt(jj)*yt(kk) - xt(kk)*yt(jj) + b(j) = yt(jj) - yt(kk) + c(j) = xt(kk) - xt(jj) + 100 continue +c- +c-.....shape function calculations..... +c- + 125 do 130 k = 1, 5, 2 + xm(k) = 0.5*(a(k)+b(k)*x+c(k)*y) + 130 continue +c- +c-.....set indexes and normalize coordinates..... +c- + if( i .le. 6 ) then + zooo = -z + else + zooo = z + endif + if(it .eq. 2) go to 132 + IF(X .EQ. 1.0) GO TO 250 +c if(x .eq. 1.0) x=0.999 + if(abs(1.0-x) .le. 1.0e-6) then +c print *,'x=',x + goto 250 + endif + zooo=zooo/(1.-x) +c if(x .eq. 1.) zooo=sign(x,zooo) + if(i .eq. 3) go to 200 + 132 if( i .gt. 6 .and. i .lt. 10 ) go to 170 + l = i + if( i .gt. 9 ) l = i - 9 + if( mod(l,2) .eq. 0 ) go to 150 +c- +c-.....corner nodes..... +c- + shap(j0)=xm(l)*((2.0*xm(l)-1.0)*(1.0+zooo)-(1.0-zooo**2))/2.0 + if(it .eq. 4) shap(j0)=shap(j0)+xm(l)/2.0*x*(1.-zooo**2) + shpx(j0)=0.5*b(l)*((1.+zooo)*(2.*xm(l)-0.5)-0.5*(1.0-zooo**2)) + if(it .eq. 4) then + shpx(j0)=shpx(j0)+xm(l)/2.*((2.*xm(l)-1.)*zooo + * /(1.-x)+zooo**2+1.)+b(l)/4.*(1.-zooo**2)*x + endif + shpy(j0)=0.5*c(l)*((1.+zooo)*(2.*xm(l)-0.5)-0.5*(1.0-zooo**2)) + if(it .eq. 4) shpy(j0)=shpy(j0)+c(l)/4.*(1-zooo**2)*x + shpz(j0) = xm(l)*(xm(l)+zooo-0.5) + if(it .eq. 4) then + shpz(j0)=(shpz(j0)-zooo*xm(l))/(1.-x)+xm(l)*zooo + endif + if(i .le. 6) shpz(j0) = -shpz(j0) + goto 240 +c- +c-.....mid side node..... +c- + 150 n1 = l - 1 + n2 = mod(l+1,6) + shap(j0) = 2.0*xm(n1)*xm(n2)*(1.0+zooo) + shpx(j0) = (1.0+zooo)*(xm(n1)*b(n2) + xm(n2)*b(n1)) + if(it .eq. 4) then + shpx(j0)=shpx(j0)+2.*xm(n1)*xm(n2)*zooo/(1.-x) + endif + shpy(j0) = (1.0+zooo)*( xm(n1)*c(n2)+xm(n2)*c(n1)) + shpz(j0) = 2.0*xm(n1)*xm(n2) + if(i .le. 6) shpz(j0) = -shpz(j0) + if(it .eq. 2) goto 240 + shpz(j0)=shpz(j0)/(1.-x) + goto 240 +c- +c-.....mid side rectangle..... +c- + 170 n1 = 1 + if( i .eq. 8 ) n1 = 3 + if( i .eq. 9 ) n1 = 5 + shap(j0) = xm(n1)*(1.0-zooo**2) + if(it .eq. 4) shap(j0)=shap(j0)*(1.-x) + shpx(j0) = 0.5*(1.0-zooo**2)*b(n1) + if(it .eq. 4) shpx(j0)=shpx(j0)*(1.-x)-xm(n1)*(zooo**2+1.) + shpy(j0) = 0.5*(1.0-zooo**2)*c(n1) + if(it .eq. 4) shpy(j0)=shpy(j0)*(1.-x) + shpz(j0) = -2.0*xm(n1)*zooo + 180 goto 240 +c- +c......special case for no 3 shape function of pyramid +c- + 200 l=i + shap(j0)=xm(l)*(2.*xm(l)-1.0) + shpx(j0)=b(l)/2.*(4.*xm(l)-1.0) + shpy(j0)=c(l)/2.*(4.*xm(l)-1.0) + shpz(j0)=0.0 + goto 240 +c +c--------Special case when x exactly equals 1.0 +c + 250 continue + shap(j0)=shpp(i,1) + shpx(j0)=shpp(i,2) + shpy(j0)=shpp(i,3) + shpy(j0)=shpp(i,4) + 240 enddo + return +c- +c-.....shape functions for cube..... +c- + 300 zm1=0.125*(1.0-z) + zp1=0.125*(1.0+z) + ym1=0.125*(1.0-y) + yp1=0.125*(1.0+y) +c- +c......corner nodes +c- + shap(1)=(1.-x)*(1.-y)*zm1*(-x-y-z-2.) + shpx(1)=(1.-y)*zm1*(2.*x+y+z+1.) + shpy(1)=(1.-x)*zm1*(x+2.*y+z+1.) + shpz(1)=(1.-x)*ym1*(x+y+2.*z+1.) + + shap(3)=(1.+x)*(1.-y)*zm1*(x-y-z-2.) + shpx(3)=(1.-y)*zm1*(2.*x-y-z-1.) + shpy(3)=(1.+x)*zm1*(-x+2.*y+z+1.) + shpz(3)=(1.+x)*ym1*(-x+y+2.*z+1.) + + shap(5)=(1.+x)*(1.+y)*zm1*(x+y-z-2.) + shpx(5)=(1.+y)*zm1*(2.*x+y-z-1.) + shpy(5)=(1.+x)*zm1*(x+2.*y-z-1.) + shpz(5)=(1.+x)*yp1*(-x-y+2.*z+1.) + + shap(7)=(1.-x)*(1.+y)*zm1*(-x+y-z-2.) + shpx(7)=(1.+y)*zm1*(2.*x-y+z+1.) + shpy(7)=(1.-x)*zm1*(-x+2.*y-z-1.) + shpz(7)=(1.-x)*yp1*(x-y+2.*z+1.) + + shap(13)=(1.-x)*(1.-y)*zp1*(-x-y+z-2.) + shpx(13)=(1.-y)*zp1*(2.*x+y-z+1.) + shpy(13)=(1.-x)*zp1*(x+2.*y-z+1.) + shpz(13)=(1.-x)*ym1*(-x-y+2.*z-1.) + + shap(15)=(1.+x)*(1.-y)*zp1*(x-y+z-2.) + shpx(15)=(1.-y)*zp1*(2.*x-y+z-1.) + shpy(15)=(1.+x)*zp1*(-x+2.*y-z+1.) + shpz(15)=(1.+x)*ym1*(x-y+2.*z-1.) + + shap(17)=(1.+x)*(1.+y)*zp1*(x+y+z-2.) + shpx(17)=(1.+y)*zp1*(2.*x+y+z-1.) + shpy(17)=(1.+x)*zp1*(x+2.*y+z-1.) + shpz(17)=(1.+x)*yp1*(x+y+2.*z-1.) + + shap(19)=(1.-x)*(1.+y)*zp1*(-x+y+z-2.) + shpx(19)=(1.+y)*zp1*(2.*x-y-z+1.) + shpy(19)=(1.-x)*zp1*(-x+2.*y+z-1.) + shpz(19)=(1.-x)*yp1*(-x+y+2.*z-1.) +c- +c......mid-side nodes +c- + x21=0.25*(1.0-x*x) + y21=0.25*(1.0-y*y) + z21=0.25*(1.0-z*z) + x2=0.5*x + y2=0.5*y + z2=0.5*z + shap(2)=x21*(1.-y)*(1.-z) + shpx(2)=-x2*(1.-y)*(1.-z) + shpy(2)=-x21*(1.-z) + shpz(2)=-x21*(1.-y) + + shap(4)=(1.+x)*y21*(1.-z) + shpx(4)=y21*(1.-z) + shpy(4)=-y2*(1.+x)*(1.-z) + shpz(4)=-y21*(1.+x) + + shap(6)=x21*(1.+y)*(1.-z) + shpx(6)=-x2*(1.+y)*(1.-z) + shpy(6)=x21*(1.-z) + shpz(6)=-x21*(1.+y) + + shap(8)=(1.-x)*y21*(1.-z) + shpx(8)=-y21*(1.-z) + shpy(8)=-y2*(1.-x)*(1.-z) + shpz(8)=-y21*(1.-x) + + shap(9)=(1.-x)*(1.-y)*z21 + shpx(9)=-(1.-y)*z21 + shpy(9)=-(1.-x)*z21 + shpz(9)=-z2*(1.-x)*(1.-y) + + shap(10)=(1.+x)*(1.-y)*z21 + shpx(10)=(1.-y)*z21 + shpy(10)=-(1.+x)*z21 + shpz(10)=-z2*(1.+x)*(1.-y) + + shap(11)=(1.+x)*(1.+y)*z21 + shpx(11)=(1.+y)*z21 + shpy(11)=(1.+x)*z21 + shpz(11)=-z2*(1.+x)*(1.+y) + + shap(12)=(1.-x)*(1.+y)*z21 + shpx(12)=-(1.+y)*z21 + shpy(12)=(1.-x)*z21 + shpz(12)=-z2*(1.-x)*(1.+y) + + shap(14)=x21*(1.-y)*(1.+z) + shpx(14)=-x2*(1.-y)*(1.+z) + shpy(14)=-x21*(1.+z) + shpz(14)=x21*(1.-y) + + shap(16)=(1.+x)*y21*(1.+z) + shpx(16)=y21*(1.+z) + shpy(16)=-y2*(1.+x)*(1.+z) + shpz(16)=y21*(1.+x) + + shap(18)=x21*(1.+y)*(1.+z) + shpx(18)=-x2*(1.+y)*(1.+z) + shpy(18)=x21*(1.+z) + shpz(18)=x21*(1.+y) + + shap(20)=(1.-x)*y21*(1.+z) + shpx(20)=-y21*(1.+z) + shpy(20)=-y2*(1.-x)*(1.+z) + shpz(20)=y21*(1.-x) + + return + 500 continue +c- +c......section for tetrahedron +c- + do i=1,nen + i1=irf(i,1) + i2=irf(i,2) + ia=i1 + do 550 k=1,2 + go to (505,515,525,535),ia + 505 sn(k)=1.-x-y-z + snx(k)=-1. + sny(k)=-1. + snz(k)=-1. + go to 550 + 515 sn(k)=z + snx(k)=0. + sny(k)=0. + snz(k)=1. + go to 550 + 525 sn(k)=x + snx(k)=1. + sny(k)=0. + snz(k)=0. + go to 550 + 535 sn(k)=y + snx(k)=0. + sny(k)=1. + snz(k)=0. + 550 ia=i2 +c- +c......test for corner nodes +c- + if(i1 .ne. i2) then +c- +c......for midsides evaluate function etc +c- + shap(i)=4.*sn(1)*sn(2) + shpx(i)=4.*(snx(1)*sn(2)+sn(1)*snx(2)) + shpy(i)=4.*(sny(1)*sn(2)+sn(1)*sny(2)) + shpz(i)=4.*(snz(1)*sn(2)+sn(1)*snz(2)) + else +c- +c......corner node evaluation +c- + shap(i)=sn(1)*(2.*sn(1)-1.) + shpx(i)=snx(1)*(4.*sn(1)-1.) + shpy(i)=sny(1)*(4.*sn(1)-1.) + shpz(i)=snz(1)*(4.*sn(1)-1.) + endif + enddo +c- +c......final step +c- + return + end