program f90_FEM_SEEP !****************************************************************************************** !2D Unsteady Seepage Flow Analysis with the function of saturation-unsaturation analysis !****************************************************************************************** implicit none integer::i,j,k,l,k1,k2,k3,n,ia,ja,ne,nnn character::strcom*256 !Comment integer::NODT !Number of nodes integer::NELT !Number of elements integer::MATEL !Number of material sets integer::KOH !Number of nodes with given total head integer::KOQ !Number of nodes with given discharge integer::KOU !Number of nodes with the condition of seepage point integer::KOP !Number of nodes with zero pressure head in seepage points integer::nt !Number of degrees of freedom of all nodes (NODT x nhen) integer::mm !Number of degrees of freedom of reduced FE equation integer::ib !Band width of matrix integer::nod !Number of nodes per element integer::nhen=1 !Number of degree of freedom per node integer::idan !Section (0: vertical, 1: horizontal) integer,allocatable::kakom(:,:) !Element connectivity real(8),allocatable::Ak0(:) !Saturated permeability coefficient of element real(8),allocatable::Akr(:) !Relative hydraulic conductivity (ratio of permeability: unsat/sat) real(8),allocatable::alpha(:) !Unsaturated property parameter 'alpha' of element real(8),allocatable::em(:) !Unsaturated property parameter 'm' of element integer,allocatable::matno(:) !Material set number real(8),allocatable::wm(:,:) !work for material properties real(8),allocatable::x(:) !x-coordinate of node real(8),allocatable::z(:) !z-coordinate of node integer,allocatable::nokh(:) !Node number with given total head integer,allocatable::nokq(:) !Node number with given discharge integer,allocatable::noku(:) !Node number with the condition of seepage point integer,allocatable::nokp(:) !Node number with zero pressure head in seepage points real(8),allocatable::Hinp(:) !Given total head real(8),allocatable::Qinp(:) !Given discharge real(8),allocatable::tak(:,:) !Permeability matrix real(8),allocatable::wak(:,:) !work for nodes with giventotal head real(8),allocatable::hvec(:) !Total head vector real(8),allocatable::h1vec(:) !Total head vector real(8),allocatable::h2vec(:) !Total head vector real(8),allocatable::qvec(:) !Discharge vector real(8),allocatable::q0vec(:) !Discharge vector real(8),allocatable::reac(:) !work vector real(8),allocatable::vx(:) !Element velocity vector in x-direction real(8),allocatable::vz(:) !Element velocity vector in z-direction real(8),allocatable::vm(:) !Element mean velocity vector real(8),allocatable::mindex(:) !Index for treatment of nodes with given total head real(8)::eak(1:4, 1:4) !Element permeability matrix real(8)::fact0=1.0E-30 !0 judgement real(8)::elarea !Half area of element integer::itmax=100 !Limit value of iterative calculation integer::icount !Number of converged degrees of freedom real(8)::QTi,QTo !Discharge of in-flow, discharge of out-flow real(8)::vmmax,xgm,zgm !Max.velocity, coordinates of centroid of element with max.velocity real(8)::vimax,xgi,zgi !Max.velocity in in-flow area, coordinates of centroid of element with max.velocity real(8)::vomax,xgo,zgo !Max.velocity in out-flow area, coordinates of centroid of element with max.velocity real(8)::xg,zg !Coordinates of centroid of element integer::kmmax,kimax,komax !Element number with max.velocity character::fnameR*50,fnameW*50 !Input file name, output file name character :: linebuf*1000 character :: fmt1*200,fmt2_3*200,fmt2_4*200,fmt3*200,fmt4*200 real(4) :: t1,t2 fmt1="(i5,2(',',e15.7),3(',',i2),2(',',e15.7))" !Node fmt2_3="(i5,3(',',i5),3(',',e15.7),',',i5)" !3 node element fmt2_4="(i5,4(',',i5),3(',',e15.7),',',i5)" !4 node element fmt3="(i5,2(',',e15.7),3(',',e15.7))" !head and discharge fmt4="(i5,2(',',e15.7),3(',',e15.7),2(',',e15.7),',',i5)" !Velocity and permeability call getarg(1,fnameR) !Input file name call getarg(2,fnameW) !Output file name !To start the measurement of calculation taime t1=DIFT() !*************************** !data input !*************************** open(11,file=fnameR,status='old') !!-------------------------------------------------------------------------------------- !!Comment !!-------------------------------------------------------------------------------------- read(11,'(a)') strcom !------------------------------------------------------------------------------------------ !Basic value for analysis !------------------------------------------------------------------------------------------ read(11,*) nod,NODT,NELT,MATEL,KOH,KOQ,KOU,idan !------------------------------------------------------------------------------------------ !To declare array size !------------------------------------------------------------------------------------------ nt=NODT*nhen allocate(kakom(1:NELT,1:nod)) allocate(matno(1:NELT)) allocate(wm(1:NELT,1:3)) allocate(Ak0(1:NELT)) allocate(Akr(1:NELT)) allocate(alpha(1:NELT)) allocate(em(1:NELT)) allocate(x(1:NODT)) allocate(z(1:NODT)) if(0{q0} !(In-flow is not allowed at seepage point. In-flow means positive discharge.) !(If discharge has positive value at seepage point, the discharge is set to zero) !----------------------------------------------------- do i=1,nt q0vec(i)=0.0D0 end do if(00)then mm=mm+1 mindex(mm)=i end if end do do i=1,mm ia=mindex(i) reac(i)=qvec(ia) do j=1,mm ja=mindex(j) tak(i,j)=tak(ia,ja) end do end do !----------------------------------------------------- !Calculation of band width & memory of band matrix !----------------------------------------------------- ib=IBAND(mm,tak,fact0,nt) call BAND(mm,ib,tak,nt) !----------------------------------------------------- !To verify diagonal elements & to stop if zero value !----------------------------------------------------- do i=1,mm if(abs(tak(i,1))m1)then mm=m1 else mm=n1-i+1 end if do j=1,mm z=array(i,j) if(j/=m1)then if(m1-j+11)then array(i,j)=z/array(i,1) else array(i,1)=sqrt(z) end if else array(i,j)=z/array(i,1) end if end do end do end subroutine CHBAND10 subroutine CHBAND11(n1,m1,array,vector,nt) integer,intent(in)::n1,m1,nt real(8),intent(in)::array(1:nt,1:nt) real(8),intent(out)::vector(nt) integer :: i,j,mm real(8) :: z vector(1)=vector(1)/array(1,1) do i=2,n1 z=vector(i) if(i>m1)then mm=m1 else mm=i end if do j=2,mm z=z-array(i-j+1,j)*vector(i-j+1) end do vector(i)=z/array(i,1) end do vector(n1)=vector(n1)/array(n1,1) do i=n1-1,1,-1 z=vector(i) do j=2,m1 if(i+j-1>n1)exit z=z-array(i,j)*vector(i+j-1) end do vector(i)=z/array(i,1) end do end subroutine CHBAND11 integer function IBAND(mm,array,fact0,nt) integer,intent(in)::mm,nt real(8),intent(in)::array(1:nt,1:nt),fact0 integer :: i,j !----------------------------------------------------- !Calculation of band width !----------------------------------------------------- IBAND=1 do i=1,mm-1 do j=mm,i+1,-1 if(fact0<=abs(array(i,j)))exit end do if(IBAND<=j-i+1)IBAND=j-i+1 end do if(mm