rhsV.f90 1.19 KB
! set specified boundaries of V

! Jacobien
A = 1/(4.0*(dx**2))

do j = 2,ny,1
jp = j+1
jm = j-1
do  i = 2,nx,1
ip = i+1
im = i-1
! DO NOT FORGET TO CHECK SIGN
jacoV(i,j) = - A*(  ((v(ip,j)-v(im,j))*(psi(i,jp)-psi(i,jm))                               &
&                    -(v(i,jp)-v(i,jm))*(psi(ip,j)-psi(im,j)))                              &
&                    +(v(ip,j)*(psi(ip,jp)-psi(ip,jm))-v(im,j)*(psi(im,jp)-psi(im,jm))      &
&                     -v(i,jp)*(psi(ip,jp)-psi(im,jp))+v(i,jm)*(psi(ip,jm)-psi(im,jm)))     &
&                    +(v(ip,jp)*(psi(i,jp)-psi(ip,j))-v(im,jm)*(psi(im,j)-psi(i,jm))        &
&                     -v(im,jp)*(psi(i,jp)-psi(im,j))+v(ip,jm)*(psi(ip,j)-psi(i,jm)))  )/3.0

enddo
end do

call periodic(jacoV)

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

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

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

! boundary of RHS
rhs(:,1) = 0
rhs(:,ny)  = 0.75 * rhs(:,ny-1) + 0.25 * rhs(:,ny-3)
rhs = (1/nu) * ( f * rhs - epsi * jacoV(1:nx,:))

! Boundary of V
call bndyNull(v,2)
call bndyNull(v,1)