text2.f90

fortran程序

       m=i-1
       do k=l,m; b(k)=b(k)-x*a(ki+k); end do
     end if
  end do
 b(1)=b(1)/a(1)
 return
end subroutine spabac               
subroutine checon(loads,oldlds,tol,converged)
! sets converged to .false. if relative change in loads and
! oldlds is greater than tol and updates oldlds
implicit none
real,intent(in)::loads(0:),tol;real,intent(in out)::oldlds(0:)
logical,intent(out)::converged
  converged=.true.
  converged=(maxval(abs(loads-oldlds))/maxval(abs(loads))<=tol)
  oldlds=loads
 return
end subroutine checon
subroutine formkb(kb,km,g)
! lower triangular global stiffness kb stored as kb(n,iw+1)
implicit none
real,intent(in)::km(:,:);real,intent(out)::kb(:,:)
integer,intent(in)::g(:);integer::iw,idof,i,j,icd
idof=size(km,1);  iw=size(kb,2)-1
   do i=1,idof
      if(g(i)>0) then 
         do j=1,idof
            if(g(j)>0) then
               icd=g(j)-g(i)+iw+1
               if(icd-iw-1<=0) kb(g(i),icd)= kb(g(i),icd) +km(i,j)
            end if
         end do
      end if
   end do
 return
end subroutine formkb
subroutine formtb(kb,km,g)
! assembles unsymmetrical band matrix kb from constituent km
  implicit none
  real,intent(in)::km(:,:); integer,intent(in)::g(:)
  real,intent(out)::kb(:,:); integer::i,j,idof,icd,iw
  idof=size(km,1); iw=(size(kb,2)-1)/2
  do i=1,idof
     if(g(i)/=0) then
        do j=1,idof
           if(g(j)/=0) then
              icd=g(j)-g(i)+iw+1
              kb(g(i),icd)=kb(g(i),icd)+km(i,j)
           end if
        end do
     end if
  end do
 return
end subroutine formtb
subroutine bantmul(kb,loads,ans)
! multiplies unsymmetrical band kb by vector loads
! could be much improved for vector processors
! look out for chance for zero-sized arrays
implicit none
  real,intent(in)::kb(:,:),loads(0:); real,intent(out)::ans(0:)
  integer::i,j,k,l,m,n,iw; real::x
  n=size(kb,1); l=size(kb,2); iw=(l-1)/2
    do i=1,n
       x=.0; k=iw+2
       do j=1,l
          k=k-1; m=i-k+1
          if(m<=n.and.m>=1)x=x+kb(i,j)*loads(m)
       end do
       ans(i)=x
    end do
 return
end subroutine bantmul 
subroutine linmul(bk,disps,loads)
! matrix-vector multiply for symmetric matrix bk
! stored in upper triangular form as a vector
implicit none
  real,intent(in)::bk(:),disps(0:); real,intent(out)::loads(0:)
   integer::i,j,n,iw; real::x; n=ubound(disps,1);iw=ubound(bk,1)/n-1 
    do i=1,n
       x=0.0
       do j=1,iw+1
          if(i+j<=n+1)    x=x+bk(n*(j-1)+i)*disps(i+j-1)
       end do
       do j=2,iw+1
          if(i-j+1>=1)    x=x+bk((n-1)*(j-1)+i)*disps(i-j+1)
       end do    
       loads(i)=x
    end do
  return
end subroutine linmul
subroutine linmul_sky(bp,disps,loads,kdiag)
! skyline product of symmetric matrix and a vector
implicit none
 real,intent(in)::bp(:),disps(0:);real,intent(out)::loads(0:)
 integer,intent(in)::kdiag(:); integer::n,i,j,low,lup,k; real::x
  n=ubound(disps,1)
  do i = 1 , n
     x = .0 ; lup=kdiag(i)
     if(i==1)low=lup; if(i/=1)low=kdiag(i-1)+1
     do j = low , lup
        x = x + bp(j) * disps(i + j - lup) 
     end do
     loads(i) = x
     if(i == 1) cycle   ; lup = lup - 1
     do j = low , lup
        k = i + j -lup - 1
        loads(k) = loads(k) + bp(j)*disps(i)        
     end do
  end do
 return
end subroutine linmul_sky  
function formxi(fsoil,fmax,rf,rm,ro) result(xi)
! soil spring stiffness in coupled pile-soil analysis
implicit none     ; real :: xi
 real,intent(in)::fsoil,fmax,rf,rm,ro;  real::phi
 phi = fsoil*ro*rf/fmax
 xi = log((rm-phi)/(ro-phi))+phi*(rm-ro)/((rm-phi)*(ro-phi))
return
end function formxi  
subroutine fkdiag(kdiag,g)
! finds the maximum bandwidth for each freedom
implicit none
 integer,intent(in)::g(:); integer,intent(out)::kdiag(:)
 integer::idof,i,iwp1,j,im,k
  idof=size(g)
  do i = 1,idof
     iwp1=1
     if(g(i)/=0) then
        do j=1,idof
           if(g(j)/=0) then
              im=g(i)-g(j)+1
              if(im>iwp1) iwp1=im
           end if
        end do
        k=g(i);   if(iwp1>kdiag(k))kdiag(k)=iwp1
     end if
  end do
 return
end subroutine fkdiag
subroutine invar(stress,sigm,dsbar,theta)
! forms the stress invariants in 2-d or 3-d
implicit none
  real,intent(in)::stress(:);real,intent(out)::sigm,dsbar,theta
  real::sx,sy,sz,txy,dx,dy,dz,xj3,sine,s1,s2,s3,s4,s5,s6,ds1,ds2,ds3,d2,d3,sq3
  integer :: nst ; nst = ubound(stress,1)
 select case (nst)
 case(4)
  sx=stress(1); sy=stress(2); txy=stress(3); sz=stress(4)
  sigm=(sx+sy+sz)/3.
  dsbar=sqrt((sx-sy)**2+(sy-sz)**2+(sz-sx)**2+6.*txy**2)/sqrt(2.)
  if(dsbar<1.e-10) then
     theta=.0
  else
     dx=(2.*sx-sy-sz)/3.; dy=(2.*sy-sz-sx)/3.; dz=(2.*sz-sx-sy)/3.
     xj3=dx*dy*dz-dz*txy**2
     sine=-13.5*xj3/dsbar**3
     if(sine>1.) sine=1.
     if(sine<-1.) sine=-1.
     theta=asin(sine)/3.
  end if
 case(6)
  sq3=sqrt(3.);  s1=stress(1)  ;  s2=stress(2)
  s3=stress(3) ;  s4=stress(4);  s5=stress(5);  s6=stress(6)
  sigm=(s1+s2+s3)/3.
  d2=((s1-s2)**2+(s2-s3)**2+(s3-s1)**2)/6.+s4*s4+s5*s5+s6*s6
  ds1=s1-sigm ;  ds2=s2-sigm  ;  ds3=s3-sigm
  d3=ds1*ds2*ds3-ds1*s5*s5-ds2*s6*s6-ds3*s4*s4+2.*s4*s5*s6
  dsbar=sq3*sqrt(d2)
  if(dsbar==0.)then
      theta=0.
    else
      sine=-3.*sq3*d3/(2.*sqrt(d2)**3)
      if(sine>1.)sine=1. ;  if(sine<-1.)sine=-1. ; theta=asin(sine)/3.
  end if
 case default
  print*,"wrong size for nst in invar"
 end select
 return
end subroutine invar
subroutine formm(stress,m1,m2,m3)
! forms the derivatives of the invariants with respect to stress 2- or 3-d
 implicit none
  real,intent(in)::stress(:); real,intent(out)::m1(:,:),m2(:,:),m3(:,:)
  real::sx,sy,txy,tyz,tzx,sz,dx,dy,dz,sigm ; integer::nst , i , j
  nst=ubound(stress,1)
  select case (nst)
  case(4)
  sx=stress(1); sy=stress(2); txy=stress(3); sz=stress(4)
  dx=(2.*sx-sy-sz)/3.; dy=(2.*sy-sz-sx)/3.; dz=(2.*sz-sx-sy)/3.
  sigm=(sx+sy+sz)/3.
  m1=.0; m2=.0; m3=.0
  m1(1,1:2)=1.; m1(2,1:2)=1.; m1(4,1:2)=1.
  m1(1,4)=1.; m1(4,4)=1.; m1(2,4)=1.
  m1=m1/9./sigm
  m2(1,1)=.666666666666666; m2(2,2)=.666666666666666; m2(4,4)=.666666666666666
  m2(2,4)=-.333333333333333;m2(4,2)=-.333333333333333;m2(1,2)=-.333333333333333
  m2(2,1)=-.333333333333333;m2(1,4)=-.333333333333333;m2(4,1)=-.333333333333333
  m2(3,3)=2.; m3(3,3)=-dz
  m3(1:2,3)=txy/3.; m3(3,1:2)=txy/3.; m3(3,4)=-2.*txy/3.; m3(4,3)=-2.*txy/3.
  m3(1,1)=dx/3.; m3(2,4)=dx/3.; m3(4,2)=dx/3.
  m3(2,2)=dy/3.; m3(1,4)=dy/3.; m3(4,1)=dy/3.
  m3(4,4)=dz/3.; m3(1,2)=dz/3.; m3(2,1)=dz/3.
 case(6)
  sx=stress(1); sy=stress(2)    ;   sz=stress(3)
  txy=stress(4)  ;   tyz=stress(5) ;   tzx=stress(6)
  sigm=(sx+sy+sz)/3.
  dx=sx-sigm  ;   dy=sy-sigm ;  dz=sz-sigm
  m1 = .0; m2 = .0; m1(1:3,1:3) = 1./(3.*sigm)
  do  i=1,3 ; m2(i,i)=2. ;  m2(i+3,i+3)=6. ; end do
  m2(1,2)=-1.; m2(1,3)=-1. ; m2(2,3)=-1.; m3(1,1)=dx
  m3(1,2)=dz ; m3(1,3)=dy ; m3(1,4)=txy  ;  m3(1,5)=-2.*tyz
  m3(1,6)=tzx ; m3(2,2)=dy ; m3(2,3)=dx ; m3(2,4)=txy
  m3(2,5)=tyz ; m3(2,6)=-2.*tzx ;  m3(3,3)=dz
  m3(3,4)=-2.*txy; m3(3,5)=tyz ;  m3(3,6)=tzx
  m3(4,4)=-3.*dz ;  m3(4,5)=3.*tzx;  m3(4,6)=3.*tyz
  m3(5,5)=-3.*dx;  m3(5,6)=3.*txy ;  m3(6,6)=-3.*dy
  do  i=1,6 ;  do  j=i,6
      m1(i,j)=m1(i,j)/3.;  m1(j,i)=m1(i,j) ;  m2(i,j)=m2(i,j)/3.
      m2(j,i)=m2(i,j)   ;  m3(i,j)=m3(i,j)/3. ; m3(j,i)=m3(i,j)
  end do; end do
 case default
  print*,"wrong size for nst in formm"
 end select
 return   
 end subroutine formm              
  subroutine fmkdke(km,kp,c,ke,kd,theta)
  ! builds up 'coupled' stiffnesses ke and kd from 'elastic'
  ! stiffness km , fluid stiffness kp and coupling matrix c
  implicit none
  real,intent(in)::km(:,:),kp(:,:),c(:,:),theta
  real,intent(out)::ke(:,:),kd(:,:)
  integer::idof; idof=size(km,1)
  ke(1:idof,1:idof)=theta*km; ke(1:idof,idof+1:)=theta*c
  ke(idof+1:,1:idof)=theta*transpose(c); ke(idof+1:,idof+1:)=-theta**2*kp
  kd(1:idof,1:idof)=(theta-1.)*km; kd(1:idof,idof+1:)=(theta-1.)*c
  kd(idof+1:,1:idof)=ke(idof+1:,1:idof)
  kd(idof+1:,idof+1:)=theta*(1.-theta)*kp
  end subroutine fmkdke
subroutine banmul(kb,loads,ans)
implicit none
real,intent(in)::kb(:,:),loads(0:); real,intent(out):: ans(0:)
integer :: neq,nband,i,j ; real :: x
neq = ubound(kb,1) ; nband = ubound(kb,2) -1
  do i = 1 , neq
     x = .0
     do j = nband + 1 , 1 , -1
        if(i+j>nband+1) x = x +kb(i,j)*loads(i+j-nband-1)
     end do
     do j = nband , 1 , -1
        if(i-j<neq - nband) x = x + kb(i-j+nband+1,j)*loads(i-j+nband+1)
     end do
     ans(i) = x
  end do
 return
end subroutine banmul
subroutine chobk1(kb,loads)
implicit none
real,intent(in)::kb(:,:); real,intent(in out)::loads(0:)
integer :: iw,n,i,j,k   ;   real :: x
n = ubound(kb,1); iw = ubound(kb,2) - 1
  loads(1) = loads(1) / kb(1,iw+1)
    do i = 2 , n
       x = .0; k = 1
       if(i <= iw + 1 ) k = iw - i + 2   ! zero sized arrays ?
       do j = k , iw ; x=x+kb(i,j)*loads(i+j-iw-1) ; end do
       loads(i) = ( loads(i) -x)/kb(i,iw+1)
    end do
  return
end subroutine chobk1
subroutine chobk2(kb,loads)
implicit none
real,intent(in)::kb(:,:); real,intent(in out)::loads(0:)
integer :: iw,n,i,j,l,m  ;  real :: x
 n = ubound(kb,1); iw = ubound(kb,2) - 1
 loads(n) = loads(n)/kb(n,iw+1)
 do i = n-1,1,-1
    x = .0; l = i + iw
    if(i > n-iw) l = n ; m = i + 1
    do j = m , l; x=x+kb(j,iw+i-j+1) * loads(j) ; end do
    loads(i) = (loads(i) - x)/ kb(i,iw+1)
 end do
return
end subroutine chobk2     
 subroutine fmbeam(der2,fun,points,i,ell)
!
! this subroutine forms the beam shape functions
! and their 2nd derivatives in local coordinates
!
 implicit none
 real,intent(in)::points(:,:),ell
 integer,intent(in)::i
 real,intent(out)::der2(:),fun(:)
 real::xi,xi2,xi3
 xi=points(i,1); xi2=xi*xi; xi3=xi2*xi
 fun(1)=.25*(xi3-3.*xi+2.); fun(2)=.125*ell*(xi3-xi2-xi+1.)
 fun(3)=.25*(-xi3+3.*xi+2.); fun(4)=.125*ell*(xi3+xi2-xi-1.)
 der2(1)=1.5*xi; der2(2)=.25*ell*(3.*xi-1.)
 der2(3)=-1.5*xi; der2(4)=.25*ell*(3.*xi+1.)
 end subroutine fmbeam                                                       
 subroutine beam_km(km,ei,ell)
!
!      this subroutine forms the stiffness matrix of a
!      line beam element(bending only)
!
 implicit none
   real,intent(in)::ei,ell; real,intent(out):: km(:,:)
   km(1,1)=12.*ei/(ell*ell*ell) ;km(3,3)=km(1,1)
   km(1,2)=6.*ei/(ell*ell) ;  km(2,1)=km(1,2) ; km(1,4)=km(1,2)
   km(4,1)=km(1,4) ;  km(1,3)=-km(1,1) ; km(3,1)=km(1,3) ; km(3,4)=-km(1,2)
   km(4,3)=km(3,4) ; km(2,3)=km(3,4) ;  km(3,2)=km(2,3); km(2,2)=4.*ei/ell
   km(4,4)=km(2,2) ;km(2,4)=2.*ei/ell ; km(4,2)=km(2,4)
  return
 end subroutine beam_km                                                      
 subroutine beam_mm(mm,fs,ell)
 implicit none
 real,intent(in)::fs,ell
 real,intent(out)::mm(:,:)
 real::fac
!
!     this subroutine forms the mass matrix of a beam element
!
 fac=(fs*ell)/420. 
 mm(1,1)=156.*fac; mm(3,3)=mm(1,1); mm(1,2)=22.*ell*fac; mm(2,1)=mm(1,2)
 mm(3,4)=-mm(1,2); mm(4,3)=mm(3,4); mm(1,3)=54.*fac; mm(3,1)=mm(1,3)
 mm(1,4)=-13.*ell*fac; mm(4,1)=mm(1,4); mm(2,3)=-mm(1,4)
 mm(3,2)=mm(2,3); mm(2,2)=4.*(ell**2)*fac; mm(4,4)=mm(2,2)
 mm(2,4)=-3.*(ell**2)*fac; mm(4,2)=mm(2,4)
 return
 end subroutine beam_mm  
subroutine beam_kp(kp,coord,pax)
!
!      this subroutine forms the terms of the beam stiffness
!      m