Friday, November 14, 2014

assnod.f --- First source for free

      SUBROUTINE assnod(inputf, nomodi)
      implicit none
      character*(*) inputf
      logical nomodi
      character*72 smthg
      double precision DVOL
      double precision SL
      double precision SB    
      double precision SD    
      double precision XG    
      double precision ZG    
      double precision ZB    
      double precision S0    
      double precision S1    
      double precision BMT   
      double precision BML   
      double precision KAPPAT
      double precision KAPPAL
      double precision P
      double precision FSPanW,FSPANWP

      double precision fFSL
      double precision aFSL
      double precision FSL, dFSL
      double precision FSW
      double precision eps
      double precision FSPB
      double precision sum
      integer fNPFSL
      integer aNPFSL
      integer NPFSW
      integer NPFSL
      integer NPFS
      integer NGRFSL
      integer NGRFSW
      integer NGRFS

      double precision FSP

      integer NET
      integer Nmax
      parameter (Nmax=2501)
      dimension P(Nmax,3)
      dimension FSP(Nmax,3)
      dimension FSPB(Nmax,3)
      dimension NET(Nmax,4)
      integer i, j
      integer NGRDL, NGRDH, NGRD, NPNLL, NPNLH, NPNL
      common /gen/DVOL,SL,SB,SD,XG,ZG,ZB,S0,S1,BMT,BML,KAPPAT,KAPPAL
      common /grid/NGRDL,NGRDH,NGRD,NPNLL,NPNLH,NPNL
     &            ,NGRFSL,NGRFSW,NGRFS,NPFSL,NPFSW,NPFS
      common /topo/P,NET
      common /fs/FSW, FSL, dFSL
      open (unit=10, FILE=inputf,STATUS = 'OLD')
      open (unit=11, FILE='.\nodes.dat',STATUS = 'UNKNOWN')

      do 100 i=1,7
        read(10, 101) smthg
100   continue
101   format (a72)

      read(10, 102) DVOL
      read(10, 102) SL
      read(10, 102) SB    
      read(10, 102) SD    
      read(10, 102) XG    
      read(10, 102) ZG    
      read(10, 102) ZB    
      read(10, 102) S0    
      read(10, 102) S1    
      read(10, 102) BMT   
      read(10, 102) BML   
      read(10, 102) KAPPAT
      read(10, 102) KAPPAL
102   format (21X,d12.5) 

      read(10, 103) smthg
      read(10, 103) smthg
103   format (a72)

      read(10, 104) NGRD, NGRDL, NGRDH
      read(10, 104) NPNL, NPNLL, NPNLH
104   format(7X,I6,1X,I6,1X,I6) 

      read(10, 105) smthg
      read(10, 105) smthg
      read(10, 105) smthg
105   format (a72)

      do 300 i=1, NGRD
        read (10, 301) P(i,1), P(i,2), P(i,3)
300   continue
301   format (6X,3(1X,D15.8))

      if (nomodi) then
        do 350 i=1, NGRD
          P(i,2)=(SB/2.d0)
     &          *(1.d0-(P(i,3)/SD)*(P(i,3)/SD))
     &          *(1.d0-(P(i,1)/(SL/2.d0))*(P(i,1)/(SL/2.d0)))
c     &          *(1.d0+0.2d0*(P(i,1)/(SL/2.d0))*(P(i,1)/(SL/2.d0)))
350     continue
      end if

      write (11, 11) 
      write (11, 12)

      write (11, 13) NGRDH, NGRDL
      do 3 i=1, NGRD
        write (11, 1) P(i, 1), P(i, 2), P(i, 3) 
3     continue

1     format (3(F20.15, 1X))
11    format ('TITLE = ', '"NODES"')
12    format ('VARIABLES = "x" , "y", "z"')
13    format ('ZONE ', 'I=', I3, 1H, 'J=', I3, 1H, 'F=POINT')
15    format (I4, 4(1X, I5))


      read(10, 106) smthg
      read(10, 106) smthg
      read(10, 106) smthg
106   format (a72)

      do 400 i=1, NPNL
        read (10, 401) NET(i,1), NET(i,2), NET(i,3), NET(i,4)
400   continue
401   format (6X,4(1X,I4))

      close(10)
c     ------------Free Surface
      fFSL=.2D0
      aFSL=.5D0
      FSW=.75D0

      eps=1.0D-6
      j=0      
      do 500 i=1, NGRD
        if (abs(P(i, 3)).lt.eps) then
          j=j+1
          FSPB(j, 1)=P(i, 1)
          FSPB(j, 2)=P(i, 2)
          FSPB(j, 3)=0.D0
        end if
500   continue

      fFSL=fFSL*SL
      aFSL=aFSL*SL
      FSW=FSW*SL

c     no of panels on FS
      fNPFSL=8
      aNPFSL=20
      NPFSW=20

      NPFSL=fNPFSL+NPNLL+aNPFSL
      NPFS=NPFSL*NPFSW

c     no of grid points on FS
      NGRFSL=fNPFSL+NPNLL+aNPFSL+1
      NGRFSW=NPFSW+1
      NGRFS=NGRFSL*NGRFSW

c     free surface length & length of the panel
      FSL=fFSL+SL+aFSL
      dFSL=FSL/dfloat(NPFSL)

c     -----------fore part

      do 600 j=1, fNPFSL
        sum=0.D0
        do 700 i=1, NGRFSW
          FSP((j-1)*NGRFSW+i, 1)=P(1, 1)+dfloat(fNPFSL-j+1)
     &                          *fFSL/dfloat(fNPFSL)
c          sum=sum+FSPanW(FSW, i-1, NPFSW)
c          FSP((j-1)*NGRFSW+i, 2)=FSW-sum
          FSP((j-1)*NGRFSW+i, 2)=FSPanWP(FSW, i, NGRFSW)
          FSP((j-1)*NGRFSW+i, 3)=0.D0
700     continue
600   continue

c     -----------ship side
      do 630 j=fNPFSL+1, fNPFSL+NGRDL
        sum=0.D0
        do 730 i=1, NGRFSW
          FSP((j-1)*NGRFSW+i, 1)=FSPB(j-fNPFSL, 1)
c          sum=sum+FSPanW(FSW-FSPB(j-fNPFSL,2), i-1, NPFSW)
c          FSP((j-1)*NGRFSW+i, 2)=FSW-sum
          FSP((j-1)*NGRFSW+i, 2)=FSPanWP(FSW-FSPB(j-fNPFSL,2),i,NGRFSW)
     &                          +FSPB(j-fNPFSL,2)
          FSP((j-1)*NGRFSW+i, 3)=0.D0
730     continue
630   continue

c     -----------aft part
      do 660 j=fNPFSL+NGRDL+1, fNPFSL+NGRDL+aNPFSL
        sum=0.D0
        do 760 i=1, NGRFSW
          FSP((j-1)*NGRFSW+i, 1)=P(NGRD, 1)-dfloat(j-fNPFSL-NGRDL)
     &                          *aFSL/dfloat(aNPFSL)
c          sum=sum+FSPanW(FSW, i-1, NPFSW)
c          FSP((j-1)*NGRFSW+i, 2)=FSW-sum
          FSP((j-1)*NGRFSW+i, 2)=FSPanWP(FSW, i, NGRFSW)
          FSP((j-1)*NGRFSW+i, 3)=0.D0
760     continue
660   continue

      do 800 i=1, NGRFS
        do 800 j=1, 3
          P(NGRD+i,j)=FSP(i,j)
800   continue

      do 900 j=1, NPFSL
        do 950 i=1, NPFSW
          NET((j-1)*NPFSW+i+NPNL,1)=NGRD+i+(j-1)*(NPFSW+1)
          NET((j-1)*NPFSW+i+NPNL,2)=NGRD+i+1+NPFSW+(j-1)*(NPFSW+1)
          NET((j-1)*NPFSW+i+NPNL,3)=NGRD+i+2+NPFSW+(j-1)*(NPFSW+1)
          NET((j-1)*NPFSW+i+NPNL,4)=NGRD+i+1+(j-1)*(NPFSW+1)
950     continue
900   continue

c     ------------ end Free Surface      

      write (11, 13) NGRFSW, NGRFSL
      do 4 i=1, NGRFS
        write (11, 1) P(i+NGRD, 1), P(i+NGRD, 2), P(i+NGRD, 3) 
4     continue

      close(11)

      open (unit=12, FILE='.\net-fs.dat', STATUS='UNKNOWN')
      do 5 i=NPNL+1, NPNL+NPFS
        write (12, 15) i, NET(i, 1), NET(i, 2), NET(i, 3), NET(i, 4)
5     continue
      close(12)

      return
      END


      FUNCTION FSPanWP(WW, i, NW)
      implicit none
      double precision FSPanWP, WW
      double precision aa, bb

      dimension aa(3,3), bb (3)
      double precision Cmin, Cmax, PanWm, PanWmin, PanWmax
      integer i, NW

      Cmin=3.D0/5.D0
      PanWm=WW/dfloat(NW)
      PanWmin=Cmin*PanWm

      if (i.eq.NW) then
        FSPanWP=0.D0
        return
      end if

      if (i.eq.1) then
        FSPanWP=WW
        return
      end if

      aa(1,1)=dfloat(NW*NW)
      aa(1,2)=dfloat(NW)
      aa(1,3)=1.D0
      aa(2,1)=1.D0
      aa(2,2)=1.D0
      aa(2,3)=1.D0
      aa(3,1)=dfloat(NW*NW)/4.D0
      aa(3,2)=dfloat(NW)/2.D0
      aa(3,3)=1.D0
      bb(1)=0.D0
      bb(2)=WW
      bb(3)=2.D0*WW/5.D0
      call gaussj(aa,3,3,bb,1,1)      
      FSPanWP=bb(1)*dfloat(i*i)+bb(2)*dfloat(i)+bb(3)
      return
      END

No comments: