### Blame view

highEkmanNLin/rhsZeta.f90 1.92 KB
 ```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71``` ``````! 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) ``````