subroutine matvec(u, v, ipar)

!  routine to apply the discrete differential operator A
!                                  v <-- Au
!  A is the compact representation of transformed Laplacian operator
!
!      Interior points:
!      V <- A[U] = 1/det_J*(Cxx U_x)_x + 1/y_eta*(Cyy U_y)_y + 1/det_J*(Czz U_z)_z
!
!     Boundary points:
!     V <- B[U]
!      B = boundary operator
!      for Dirichlet conditions B = I
!      for Neumann   conditions B = |d/dn|
!                    where n is outward facing unit normal
!**********************************************************

!grid _info contains Grid(0:nlevels-1),  current_grid_level
use grid_info
use intermediate_variables
use boundary_information
use mpi_parameters
use counters_flags_etc
implicit none
include 'mpif.h'

real  ::  U(Grid(current_grid_level)%n1_it,Grid(current_grid_level)%n2_it,Grid(current_grid_level)%n3_it)
real  ::  V(Grid(current_grid_level)%n1_it,Grid(current_grid_level)%n2_it,Grid(current_grid_level)%n3_it)

!Local variables
real, dimension(:,:,:), allocatable        :: dUdzeta,dUdxi,dUdeta
real, dimension(:,:,:,:), allocatable      :: scratch

integer                                    :: ipar(13),ii,flag,npts
integer                                    :: status(MPI_STATUS_SIZE)
real                                       :: h
real                                       :: end_values(2)
integer                                    :: nx,ny,nz,locnz
real                                       :: dx,dy,dz
real                                       :: xx, mean_value

V=0.0
ii = current_grid_level
nx = Grid(ii)%n1            
ny = Grid(ii)%n2            
nz = Grid(ii)%n3            
locnz = Grid(ii)%locn3      
dx = Grid(ii)%dxi          
dy = Grid(ii)%deta        
dz = Grid(ii)%dzeta      

flag=1               ! don't specify deriv values at ends
end_values=-999999.  ! obsolete but required


!***************************************************************************
!    Store approximation of A[U] in V:
!    V <-- 1/det_J*(Cxx U_x)_x + 1/y_eta*(Cyy U_y)_y + 1/det_J*(Czz U_z)_z
!***************************************************************************


!******************************************************************
!               Logic branches
!               depending on xy periodicity
!******************************************************************

 if( boundary_type2(1,1)=='periodic' .and. boundary_type4(1,1)=='periodic' ) then  

  allocate ( dUdzeta(1,1,2)  )

!***************************************************************************
! z part of the problem first, U & V are nz vectors
!***************************************************************************

 h=dz
 call compact_ddz(U,V,flag,end_values,1,1,nz)

! Save dU/dzeta 
 if( myid==0 .and. BC_pressure5(1,1) == 'Neumann' ) then   
  dUdzeta(1,1,1) = V(1,1,1)
 endif
 if( myid==nprocs-1 .and. BC_pressure6(1,1) == 'Neumann' ) then   
  dUdzeta(1,1,2) = V(1,1,nz)
 endif

!multiply U_z by Czz overwriting the result in held in V
 V(1,1,:) = Grid(ii)%Czz(1,:)*V(1,1,:)    ! Czz is indep of x

!differentiate (Czz U_z ) wrt z, overwrite results back into V
 call compact_ddz(V,V,flag,end_values,1,1,nz)

!multiply result by 1/det_J
 V(1,1,:)=V(1,1,:)/Grid(ii)%det_J(1,:)

! We've aready done the z part of the operator, now add on the algebraic, horizontal part.
! We're working here in xy transform space so U is the current xy transformed P(z) estimate
! at a fixed kx,ky location  (kx_val,ky_val; available via counters_flags_etc)

  V(1,1,:) = V(1,1,:)         &
           -( (kx_val)**2 * (Grid(ii)%Cxx(1,:)/Grid(ii)%det_J(1,:))   &
             +(ky_val)**2 * (Grid(ii)%Cyy(1)/Grid(ii)%y_eta(1))  ) * U(1,1,:)

! Neumann conditions at top&bottom
  V(1,1,1)=dUdzeta(1,1,1)
  V(1,1,nz)=dUdzeta(1,1,2)

  deallocate (dUdzeta )

 else    ! NONPERIODIC LOGIC

!***************************************************************************
! z part of the problem first
!***************************************************************************
allocate ( dUdzeta(nx,ny,2), dUdxi(ny,locnz,2), dUdeta(nx,locnz,2)  )

 h=dz
 call compact_ddz_mpi(U,V,flag,end_values,nx,ny,locnz)

! Save dU/dzeta 
 if( myid==0 .and. BC_pressure5(1,1) == 'Neumann' ) then   
  dUdzeta(:,:,1) = V(:,:,1)
 endif
 if( myid==nprocs-1 .and. BC_pressure6(1,1) == 'Neumann' ) then   
  dUdzeta(:,:,2) = V(:,:,locnz)
 endif

!multiply U_z by Czz overwriting the result in held in V
 do j=1,ny
   V(:,j,:) = Grid(ii)%Czz*V(:,j,:)
 enddo

!differentiate (Czz U_z ) wrt z, overwrite results back into V
 call compact_ddz_mpi(V,V,flag,end_values,nx,ny,locnz)

!multiply result by 1/det_J
 do j=1,ny
  V(:,j,:)=V(:,j,:)/Grid(ii)%det_J
 enddo

!******************************************************************
!                 do x part, compact differentiation
!******************************************************************

!Take x derivative of U(:,:,1:locnz) and store result in scratch(:,:,:,1),
 h=dx
 call compact_ddx(U,scratch(1,1,1,1),flag,end_values,nx,ny,locnz)

! Save dU/dxi on FACES 1&2
 dUdxi(:,:,1) = scratch(1,:,:,1)
 dUdxi(:,:,2) = scratch(nx,:,:,1)

!multiply U_x by Cxx, overwriting the result in scratch(:,:,:,1)
 do j=1,ny
  scratch(:,j,:,1) = Grid(ii)%Cxx*scratch(:,j,:,1)
 enddo

!differentiate (Cxx U_x) wrt x, writing results into scratch(:,:,:,2)
 call compact_ddx(scratch(1,1,1,1),scratch(1,1,1,2),flag,end_values,nx,ny,locnz)

!divide (Cxx U_x)_x by det_J  and add the result to V
!V then contains 1/det_J*{ (Cxx U_x)_x + (Czz U_z)_z  }
 do j=1,ny
  V(:,j,:) =  V(:,j,:) + scratch(:,j,:,2)/Grid(ii)%det_J
 enddo

!******************************************************************
!                 do y part, compact differentiation 
!******************************************************************

!Take y derivative of U(:,:,1:locnz) and store result in scratch(:,:,:,1).
 h=dy
 call compact_ddy(U,scratch(1,1,1,1),flag,end_values,nx,ny,locnz)

! Save dU/deta on FACES 3&4
 dUdeta(:,:,1) = scratch(:,1,:,1)
 dUdeta(:,:,2) = scratch(:,ny,:,1)

!multiply U_y by Cyy, overwriting the result in scratch
 do j=1,ny
  scratch(:,j,:,1) = Grid(ii)%Cyy(j)*scratch(:,j,:,1)
 enddo

!differentiate (Cyy U_y) wrt y, write results back into scratch(:,:,:,2)
  call compact_ddy(scratch(1,1,1,1),scratch(1,1,1,2),flag,end_values,nx,ny,locnz)

!add 1/y_eta (Cyy U_y)_y to V. Note Cyy=1/y_eta. V then contains:
!1/det_J*(Cxx U_x)_x + 1/y_eta*(Cyy U_y)_y + 1/det_J*(Czz U_z)_z
 do j=1,ny
  V(:,j,:) = V(:,j,:) + Grid(ii)%Cyy(j)*scratch(:,j,:,2)
 enddo

!******************************************************************
!        overwrite boundary locations w/ boundary operators
!******************************************************************

! USE IDENTITY OPERATOR FOR DIRICHLET BCs:*************************
 if( BC_pressure1(1,1) == 'Dirichlet' ) then   
  V(1,:,:)=U(1,:,:)
 endif
 if( BC_pressure2(1,1) == 'Dirichlet' ) then   
  V(nx,:,:)=U(nx,:,:)
 endif

 if( BC_pressure3(1,1) == 'Dirichlet' ) then   
  V(:,1,:)=U(:,1,:)
 endif
 if( BC_pressure4(1,1) == 'Dirichlet' ) then   
  V(:,ny,:)=U(:,ny,:)
 endif

 if( myid==0 .and. BC_pressure5(1,1) == 'Dirichlet' ) then   
  V(:,:,1)=U(:,:,1)
 endif
 if( myid==nprocs-1 .and. BC_pressure6(1,1) == 'Dirichlet' ) then   
  V(:,:,locnz)=U(:,:,locnz)
 endif

! USE NORMAL DERIV. OPERATOR FOR NEUMANN BCs:*************************
 if( BC_pressure1(1,1) == 'Neumann' ) then   
  V(1,:,:)=dUdxi(:,:,1)
 endif
 if( BC_pressure2(1,1) == 'Neumann' ) then   
  V(nx,:,:)=dUdxi(:,:,2)
 endif

 if( BC_pressure3(1,1) == 'Neumann' ) then   
  V(:,1,:)=dUdeta(:,:,1)
 endif
 if( BC_pressure4(1,1) == 'Neumann' ) then   
  V(:,ny,:)=dUdeta(:,:,2)
 endif

 if( myid==0 .and. BC_pressure5(1,1) == 'Neumann' ) then   
  V(:,:,1)=dUdzeta(:,:,1)
 endif
 if( myid==nprocs-1 .and. BC_pressure6(1,1) == 'Neumann' ) then   
  V(:,:,locnz)=dUdzeta(:,:,2)
 endif

 deallocate (dUdzeta, dUdxi,dUdeta)
endif    ! END XY PERIODIC IF BLOCK

end subroutine matvec



subroutine mg_solve(b,x,ipar)
 
! Routine to apply M ~ A_inverse to the rhs
! vector b, producing an approximate solution x.
! M is defined to be n_mg_cycles full multigrid
! cycles. The FMG cycles start and end at the
! finest (h=0) grid level.
!

use grid_info
use counters_flags_etc, only: n_mg_cycles,istep,istart,tol
use mpi_parameters, only: myid
implicit none

integer        :: i,N,ipar(13)
real           :: b(Grid(0)%n1,Grid(0)%n2,Grid(0)%locn3)
real           :: x(Grid(0)%n1,Grid(0)%n2,Grid(0)%locn3)
character(len=3) :: mg_cycle_type

  if( istep<istart ) then
   mg_cycle_type='FMG'
  else
   mg_cycle_type='V'
  endif
  N=n_mg_cycles

! (1)  store the input vector b in the "f" location of the fine grid:
   Grid(0)%f = b

! (2) MG cycles:
   if( mg_cycle_type=='FMG') then
    call FMG(0)    !  use Grid(0)%f as input and updates Grid(0)%v
   elseif(mg_cycle_type=='V') then
    Grid(0)%v=0.0
    do i=1,N
     call V_cycle(0)
     if(myid==0)write(0,*) i,Grid(0)%spar(2)
     if( Grid(0)%spar(2) < tol ) goto 999   ! convergence
    enddo
   else
    write(0,*) 'mg_cycle_type not properly defined'
    STOP
   endif
   
! (3) extract the solution from Grid(0), store in output vector x:
999 continue
   x = Grid(0)%v
   
   current_grid_level=0    ! return global variable to fine grid value
 end subroutine mg_solve


subroutine diagl(U,V,ipar)
!
! Routine to apply M ~ (1/diag) to the rhs
! vector b, producing an approximate solution x.
!

use grid_info
use boundary_information
use counters_flags_etc
use mpi_parameters
implicit none

integer        :: h,ipar(13)
real           :: diag_FD   ! approximate diagonal element in a FD discretization at level h
real           :: pi
real :: U(Grid(current_grid_level)%n1,Grid(current_grid_level)%n2,Grid(current_grid_level)%locn3)
real :: V(Grid(current_grid_level)%n1,Grid(current_grid_level)%n2,Grid(current_grid_level)%locn3)
real, allocatable                          :: kx(:),ky(:)
integer                                    :: N


h=current_grid_level

if( boundary_type2(1,1)=='periodic' .and. boundary_type2(1,1)=='periodic' ) then

!  create the x & y wavenumber arrays, align with the data arrays
  N=Grid(h)%n1-1         ! n1 includes both boundary points, N only 1
  pi=4.*atan(1.)
  allocate(kx(N+1))
  do i=1,N/2+1
   kx(i)=2.*pi*(i-1.)    ! FFTW real to half-complex storage
  enddo
  do i=2,N/2
   kx(N-i+2)=kx(i)       ! FFTW real to half-complex storage
  enddo
  kx(N+1)=0.             ! extra location, convenient for alignment of arrays
 

  N=Grid(h)%n2-1         ! n2 includes both boundary points, N only 1
  allocate(ky(N+1))
  do j=1,N/2+1
   ky(j)=2.*pi*(j-1.)    ! FFTW real to half-complex storage
  enddo
  do j=2,N/2
   ky(N-j+2)=ky(j)       ! FFTW real to half-complex storage
  enddo
  ky(N+1)=0.             ! extra location, convenient for alignment of arrays

! precondition using the diagonal elements of an approximated operator
  do i=1,Grid(h)%n1
   do j=1,Grid(h)%n2
    do k=1,Grid(h)%locn3
     diag_FD = -2.*( Grid(h)%Czz(i,k)/(Grid(h)%dzeta)**2 )  &  ! diagonal term for FD approx of z operator
                   - (kx(i)**2) * (Grid(h)%Cxx(i,k)/(Grid(h)%det_J(i,k)))           &   ! true diag term in F space
                   - (ky(j)**2) * (Grid(h)%Cyy(j)/(Grid(h)%y_eta(j)))                   ! true diag term in F space

     if( j==1 .and. k==1 ) then
      write(0,*)'****** diag:  ',i, -2.*( Grid(h)%Czz(i,k)/(Grid(h)%dzeta)**2 ), - (kx(i)**2) * (Grid(h)%Cxx(i,k)/(Grid(h)%det_J(i,k)))
     endif
     V(i,j,k) = U(i,j,k)/diag_FD   ! interior points
    enddo
   enddo
  enddo
  STOP

else

! precondition using the diagonal elements of an approximated operator
  do i=1,Grid(h)%n1
   do j=1,Grid(h)%n2
    do k=1,Grid(h)%locn3
     diag_FD = -2.*( Grid(h)%Cxx(i,k)/(Grid(h)%dxi)**2       &
                   + Grid(h)%Cyy(j  )/(Grid(h)%deta)**2      &
                   + Grid(h)%Czz(i,k)/(Grid(h)%dzeta)**2 )
     V(i,j,k) = U(i,j,k)/diag_FD   ! interior points
    enddo
   enddo
  enddo

endif

! Redo Neumann boundaries: i.e. overwrite locations where the boundary operator differs from the interior operator

if( BC_pressure1(1,1) == 'Neumann' ) then
 do k=1,Grid(h)%locn3
  do j=1,Grid(h)%n2
   diag_FD = 1./Grid(h)%dxi
   V(1,j,k)= U(1,j,k)/diag_FD
  enddo
 enddo
endif

 if( BC_pressure2(1,1) == 'Neumann' ) then
 do k=1,Grid(h)%locn3
  do j=1,Grid(h)%n2
   diag_FD = 1./Grid(h)%dxi
   V(Grid(h)%n1,j,k)= U(Grid(h)%n1,j,k)/diag_FD
  enddo
 enddo
 endif

 if( BC_pressure3(1,1) == 'Neumann' ) then 
 do k=1,Grid(h)%locn3
  do i=1,Grid(h)%n1
   diag_FD = 1./Grid(h)%deta
   V(i,1,k)= U(i,1,k)/diag_FD
  enddo
 enddo
 endif

 if( BC_pressure4(1,1) == 'Neumann' ) then  
 do k=1,Grid(h)%locn3
  do i=1,Grid(h)%n1
   diag_FD = 1./Grid(h)%deta
   V(i,Grid(h)%n2,k)= U(i,Grid(h)%n2,k)/diag_FD
  enddo
 enddo
 endif

 if( BC_pressure5(1,1) == 'Neumann' .and. myid==0 ) then
 do i=1,Grid(h)%n1
  do j=1,Grid(h)%n2
   diag_FD = 1./Grid(h)%dzeta
   V(i,j,1)= U(i,j,1)/diag_FD
  enddo
 enddo
 endif

 if( BC_pressure6(1,1) == 'Neumann' .and. myid==nprocs-1) then
 do i=1,Grid(h)%n1
  do j=1,Grid(h)%n2
   diag_FD = 1./Grid(h)%dzeta
   V(i,j,Grid(h)%locn3)= U(i,j,Grid(h)%locn3)/diag_FD
  enddo
 enddo
 endif

 if( boundary_type2(1,1)=='periodic' ) then
 do k=1,Grid(h)%locn3
  do j=1,Grid(h)%n2
   diag_FD = 1.                                    ! Identity operator
   V(Grid(h)%n1,j,k)= U(Grid(h)%n1,j,k)/diag_FD    ! only "far" endpts have I[P]=0
  enddo
 enddo
 endif

 if( boundary_type4(1,1)=='periodic' ) then
 do k=1,Grid(h)%locn3
  do i=1,Grid(h)%n1
   diag_FD = 1
   V(i,Grid(h)%n2,k)= U(i,Grid(h)%n2,k)/diag_FD
  enddo
 enddo
 endif

 return
end