subMuddx.f90 3.13 KB

subroutine computeRHS(stage, rhs)
! Set the ringht and side field 
  implicit none
  integer i,j 
  integer stage, nx,ny
  real dlx,dly,Pi , x, y
  real rhs(nx,ny)
  
   do i=1,nx
        do j=1,ny
            x = float(i-1)*dx
            y = float(j-1)*dy
            
            if (stage == 1)
              rhs(i,j)    =  -8* Pi * cos(2*Pi * x) * sin(2*Pi*y)
            else if (stage == 2)
              rhs(i,j)    =  -8* Pi * cos(2*Pi * x) * sin(2*Pi*y)
            else if (stage ==3)
              rhs(i,j)    =  -8* Pi * cos(2*Pi * x) * cos(Pi*y)
            end if
        end do
    end do
end

subroutine setBndy(stage,psi)
! set boundary condition for a specific stage
  implicit none
  integer i,j 
  integer stage, nx, ny
  real dlx,dly,Pi , x, y
  real psi(nx,ny)

! ---------------- zeta ----------------
  ! set boundary condition flags
  if (stage == 1) then 
    ! set specified boundaries in phi
    x = 0
    do j=1,ny
        y = float(j-1)*dy
        psi(1,j) = sin(2*Pi*y)  
    end do

    x = Lx
    do j=1,ny
        y = float(j-1)*dy
        psi(nx,j) = sin(2*Pi*y)
    end do

    y = 0
    do i=1,nx
        x = float(i-1)*dx
        psi(i,1) = 0.0
    end do

    y = Ly
    do i=1,nx
        x = float(i-1)*dx
        psi(i,ny) = 0.0
    end do

  ! ---------------- phi ----------------
  ! set boundary condition flags
  else if (stage == 2) then 
    ! set specified boundaries in phi
    x = xa
    do j=1,ny
        y = float(j-1)*dy
        psi(1,j) = sin(2*Pi*y)  
    end do

    x = xb
    do j=1,ny
        y = float(j-1)*dy
        psi(nx,j) = sin(2*Pi*y)
    end do

    y = yc
    do i=1,nx
        x = float(i-1)*dx
        psi(i,1) = 0.0
    end do

    y = yd
    do i=1,nx
        x = float(i-1)*dx
        psi(i,ny) = 0.0
    end do



! --------------------- v ------------------------
  ! set boundary condition flags 
  else if (stage ==3) then
    ! set specified boundaries in phi    
    x = xa
    do j=1,ny
        y = yc+float(j-1)*dly
        psi(1,j) = cos(Pi*y)  
    end do

    x = xb
    do j=1,ny
        y = yc+float(j-1)*dly
        psi(nx,j) = cos(Pi*y)
    end do

    y = yc
    do i=1,nx
      x = xa+float(i-1)*dlx
        psi(i,1) = cos(2*pi*x)
    end do

    y = yd
    do i=1,nx
        x = xa+float(i-1)*dlx
        psi(i,ny) =-cos(2*pi*x)
    end do
!------------ else error ---------------
  else
    write(*,*) 'Wrong stage for setBndy'
    call exit(2)
  end if 
end

subroutine bndc(kbdy,xory,alfa,gbdy)
! input mixed derivative b.c. to mud2
  implicit none
  integer kbdy
  real xory,alfa,gbdy,x,y,pe,px,py,pxx,pyy
  real xa,xb,yc,yd,tolmax,relmax
  common/ftmud2/xa,xb,yc,yd,tolmax,relmax
  real tau, nu
  
  tau  = 0.1
  nu = 1E-6

  if (kbdy.eq.3) then  ! y=yc tboundary
      x = xory
      y = yc
      alfa = 0
      gbdy = 0
      return
  end if
  if (kbdy.eq.4) then  ! y=yd boundary
      y = yd
      x = xory
      alfa = 0
      gbdy = tau/nu
      return
  end if
end


subroutine cof(x,y,cxx,cyy,cx,cy,ce)
! input pde coefficients at any grid point (x,y) in the solution region

  implicit none
  real x,y,cxx,cyy,cx,cy,ce
  cxx = 1
  cyy = 1
  cx = 0.
  cy = 0.
  ce = 0.
  return
end

!END MODULE