rhsZeta.f90 1.92 KB
! set specified boundaries for zeta

! Laplace 
A = 1/dx**2 
laplace(:,:)=0.0

do j = 2,ny,1
   jp = j+1
   jm = j-1
   do  i = 2,nx,1
      ip = i+1
      im = i-1                            
      laplace(i,j)= A*( psi(ip,j)+psi(im,j)+psi(i,jp)+psi(i,jm)-4.*psi(i,j) )
   enddo
end do
print*, '1'
call bndyNull(laplace,2)
call bndyNull(laplace,1)
call periodic(laplace)
! Jacobien
A = 1/(4.0*(dx**2)) 

do j = 2,ny-1,1
   jp = j+1
   jm = j-1
   do  i = 2,nx,1
     ip = i+1
     im = i-1     
     ! DO NOT FORGET TO CHECK SIGN
     jacoZ(i,j) = - A*(  ((laplace(ip,j)-laplace(im,j))*(psi(i,jp)-psi(i,jm))                               &
 &                    -(laplace(i,jp)-laplace(i,jm))*(psi(ip,j)-psi(im,j)))                              &
 &                    +(laplace(ip,j)*(psi(ip,jp)-psi(ip,jm))-laplace(im,j)*(psi(im,jp)-psi(im,jm))      &
 &                     -laplace(i,jp)*(psi(ip,jp)-psi(im,jp))+laplace(i,jm)*(psi(ip,jm)-psi(im,jm)))     &
 &                    +(laplace(ip,jp)*(psi(i,jp)-psi(ip,j))-laplace(im,jm)*(psi(im,j)-psi(i,jm))        &
 &                     -laplace(im,jp)*(psi(i,jp)-psi(im,j))+laplace(ip,jm)*(psi(ip,j)-psi(i,jm)))  )/3.0
   enddo
end do
call periodic(jacoZ)


! Making sure that PSI is periodic before making any calcul
call periodic(v)

! Set the right hans side
rhs(:,:) = 0.

! Calculating v_z
do j=2,ny-1
    rhs(:,j)   =  (v(1:nx,j+1)-v(1:nx,j-1))/(2*dy)  
end do

! Boundary ! For surface boundary Vz = tau/nu 
!rhs(:,1) = 0.0
!rhs(:,ny)  =  sin(2* Pi * x -  Pi/2 ) * 0.05
!0.75 * rhs(:,ny-1) + 0.25 * rhs(:,ny-3)

! Set the surface as Tau(x)
do i=1,nx
  x = float(i-1)*dx
  rhs(i,ny) = (v(i,ny)-v(i,ny-1))/dy!- cos(2* Pi * x / Lx ) * Vmean
  rhs(i,1)  =  (v(i,2)-v(i,1))/dy 
end do

! Multiply by constantes
rhs = (1/nu)  * (-f * rhs + epsi * jacoZ(1:nx,:))

! Boundary of Zeta 
!call bndyNull(zeta,2)
zeta(:,ny) = (realU(:,ny)-realU(:,ny-1))/dx
zeta(:,1) = (realU(:,2)-realU(:,1))/dx
call bndyNull(zeta,1)