!
!  Main.f
!
!  Free-Format Fortran Source File
!  Generated by PGI Visual Fortran(R)
!  11/25/2013 3:44:30 PM
!
      program morison
	  implicit none
	  real Hs,Tz,d,Wmax,g,hw(1000),deltaw,deltat,deltaz
      real H_airy,T_airy,lambda_airy,depth,R2
      real Tmax,t,deltay,h,diam,F,deltaF,rho,Cd,Cm,pippo
      real pi,u,dudt,knum(1000),lambda(1000),y,fase(1000)
	  real Sw(1000),w(1000)
      real lambda0,lambda1,lambda2,dlambda,factor
	  real SEED,A,B,XX,RANDOM
      real x_cilindro(10),z
      real Fmorison,Ftot,omega,wave_num,Finer,Fdrag
	  integer nw,i,imode,nt,ny,j,k,JJ
      integer Ntime,Nstat

	  pi=acos(-1.e0)
	  rho=1025.e0
	  g=9.81e0
	  Cm=2.e0
	  Cd=1.6e0
	  diam=1.e0
	  h=10.e0
	  d=10.e0
	  write(*,*)'onda monocromatica (0) oppure spettro(1)'
	  read(*,*)imode

	  if (imode .eq. 0) then
! Parametri ambientali
        write(*,*)'Inserisci H altezza, T periodo e fondale'
        read(*,*)H_airy,T_airy,depth
! Parametri geometrici del pilone
        write(*,*)'Inserisci Diametro e Immersione'
        read(*,*)R2,h

        Ntime=200
        Nstat=100
        deltat=T_airy/Ntime
        deltaz=h/Nstat
		x_cilindro(1)=0.

		omega=2*pi/T_airy

        lambda0=g*T_airy**2/(2.*pi) !da mettere a posto
        dlambda=1.e0
        lambda1=lambda0
        do while(dlambda.ge.1.e-3)
	      lambda2=lambda0*tanh(2.e0*pi*depth/lambda1)   ! da mettere a posto
		  dlambda=abs(lambda2-lambda1)/lambda0
          lambda1=lambda2
        end do
        lambda_airy=lambda2
        wave_num=2*pi/lambda_airy
        open (unit=10,file='Fmorison_z.dat',status='unknown')
        open (unit=11,file='Ftot_t.dat',status='unknown')
        do i=1,Ntime+1
          t=deltat*(i-1)
          Ftot=0.
          Finer=0.
          Fdrag=0.
          write(10,*)'time=',t
		  do k=1,Nstat
            z=-deltaz*(2*k-1)/2
            if(depth/lambda_airy.gt.0.5e0) then   !DEEP WATER
	             u =    H_airy / 2.e0 * exp(wave_num*z) *
     >                  g * wave_num / omega *
     >                  cos(wave_num*x_cilindro(1)-omega*t)
	          dudt =    H_airy / 2.e0 * exp(wave_num*z) *
     >                  g * wave_num        *
     >                  sin(wave_num*x_cilindro(1)-omega*t)
              else    !INTERMEDIATE OR SHALLOW WATERS
	             u =    H_airy / 2.e0 * cosh(wave_num *(depth+z))/
     >                  cosh(wave_num*depth) * g * wave_num / omega *
     >                  cos(wave_num*x_cilindro(1)-omega*t)
	          dudt =    H_airy / 2.e0 * cosh(wave_num *(depth+z))/
     >                  cosh(wave_num*depth) * g * wave_num *
     >                  sin(wave_num*x_cilindro(1)-omega*t)
              end if
		    Fmorison=Cm*pi*R2**2/4.*rho*dudt+Cd/2.*rho*R2*u*abs(u)
            Finer = Finer + Cm*pi*R2**2/4.*rho*dudt*deltaz
            Fdrag = Fdrag + Cd/2.*rho*R2*u*abs(u)*deltaz
            Ftot=Ftot+Fmorison*deltaz
			write(10,*)-z,Fmorison
          end do
          write(10,*)
          write(11,*)t,Ftot,Finer,Fdrag
        end do
        close(10)
        close(11)
	  elseif (imode .eq. 1) then
        call Bretschneider(nw,Wmax,Hs,Tz,Sw,w,deltaw,d)
	    call Airy(nw,d,w,Sw,hw,deltaw)
        do i=1,nw
	      write(10,*)w(i),Sw(i),hw(i)
	    end do
	    close(10)
	    nt=2400
!	    Tmax=1200.e0
        Tmax=2*pi/w(1)
	    deltat=Tmax/nt
	    ny=100
	    deltay=h/ny

	    do j=1,nw
          factor=0.
          fase(j)=2*pi*factor
        end do

	    open(20,file='FORZA.dat',status='unknown',access='sequential')
	    do i=1,nt
	      t=float(i-1)*deltat
	      F=0.e0
	      do k=1,ny
	        y=-(deltay/2+(k-1)*deltay)
            dudt=0.e0
	        u=0.e0
	        do j=1,nw
	          lambda0=g*2*pi/w(j)**2   ! da mettere a posto
			  dlambda=1.e0
              lambda1=lambda0
              do while(dlambda.ge.1.e-3)
	            lambda2=lambda0*tanh(2.e0*pi*d/lambda1)   ! da mettere a posto
			    dlambda=abs(lambda2-lambda1)/lambda0
                lambda1=lambda2
              end do
              lambda(j)=lambda2
	          knum(j)=2*pi/lambda(j)

              if(d/lambda(j).gt.0.5e0) then
	             u =    u + hw(j) / 2.e0 * exp(knum(j)*y) *
     >                  g * knum(j) / w(j) * cos(-w(j)*t+fase(j))
	          dudt = dudt + hw(j) / 2.e0 * exp(knum(j)*y) *
     >                  g * knum(j)        * sin(-w(j)*t+fase(j))
              else
	             u =    u + hw(j) / 2.e0 * cosh(knum(j)*(d+y))/
     >              cosh(knum(j)*d) * g * knum(j) / w(j) *
     >              cos(-w(j)*t+fase(j))
	          dudt = dudt + hw(j) / 2.e0 * cosh(knum(j)*(d+y))/
     >              cosh(knum(j)*d) * g * knum(j)        *
     >              sin(-w(j)*t+fase(j))

              end if
	        end do
	        deltaF=(Cm*rho*pi*diam**2/4.e0*dudt+
     >            0.5*Cd*rho*diam*u*abs(u))*deltay
	        F=F+deltaF
	      end do
          write(*,*)i,t,F
	      write(20,*)t,F
	    end do
        close(20)
	  end if

	  stop
	  end

	  subroutine Bretschneider(nw,Wmax,Hs,Tz,Sw,w,deltaw,d)
	  implicit none
	  real w(1000),Sw(1000),Wmax,Hs,Tz,A,B,g,deltaw,d
	  integer i,nw
	  write (*,*)'scrivere i valori di Hs e Tz'
	  read(*,*)Hs,Tz
	  open(10,file='spettro.dat',status='unknown',access='sequential')
	  g=9.81E0
	  Wmax=SQRT(2*ACOS(-1.)*g*tanh(2*ACOS(-1.)*d))
	  nw=200
	  deltaw=Wmax/nw
	  do i=1,nw
	  w(i)=Wmax*i/nw
	  A=173.E0*Hs**2.E0/Tz**4.E0
	  B=691.E0/Tz**4.E0
	  Sw(i)=A/w(i)**5.E0*exp(-B/w(i)**4.E0)
	  end do
	  return
	  end

	  subroutine Airy(nw,d,w,Sw,hw,deltaw)
	  implicit none
	  real w(1000),Sw(1000),hw(1000),deltaw,d
	  integer i,nw

	  do i=1,nw
	  hw(i)=2.*sqrt(2.*deltaw*Sw(i))

	  end do

	  return
	  end