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)