text:
module ml_nd_module

  use bl_constants_module
  use bl_prof_module
  use mg_module
  use fabio_module
  use ml_layout_module
  use bndry_reg_module

  implicit none

!  private :: grid_res, grid_laplace_1d, grid_laplace_2d, grid_laplace_3d

contains

  subroutine ml_nd(mla,mgt,rh,full_soln,fine_mask,one_sided_ss,ref_ratio, &
                   do_diagnostics,rel_eps,abs_eps_in)

    use ml_norm_module        , only : ml_norm_inf
    use ml_restriction_module , only : ml_restriction, periodic_add_copy
    use ml_prolongation_module, only : ml_nodal_prolongation

    type(ml_layout), intent(in   ) :: mla
    type(mg_tower ), intent(inout) :: mgt(:)
    type( multifab), intent(inout) :: rh(:)
    type( multifab), intent(inout) :: full_soln(:)
    type(lmultifab), intent(in   ) :: fine_mask(:)
    type( multifab), intent(in   ) :: one_sided_ss(2:)
    integer        , intent(in   ) :: ref_ratio(:,:)
    integer        , intent(in   ) :: do_diagnostics
    real(dp_t)     , intent(in   ) :: rel_eps

    real(dp_t)     , intent(in   ), optional :: abs_eps_in

    integer    :: nlevs

    type(multifab), allocatable  ::      soln(:)
    type(multifab), allocatable  ::        uu(:)
    type(multifab), allocatable  ::   uu_hold(:)
    type(multifab), allocatable  ::       res(:)
    type(multifab), allocatable  ::  temp_res(:)

    ! This is just a holder for a zero array at the same size as rh --
    !   we only need it to pass in to crse_fine_residual
    type(multifab), allocatable  ::   zero_rh(:)

    type(bndry_reg), allocatable :: brs_flx(:)

    type(box   ) :: pd,pdc
    type(layout) :: la, lac
    integer :: n, dm
    integer :: mglev, mglev_crse, iter
    logical :: fine_converged

    real(dp_t) :: Anorm, bnorm, fac, tres, ttres, tres0, abs_eps, t1(3), t2(3)
    real(dp_t) :: stime, bottom_solve_time

    character(len=3)          :: number
    character(len=20)         :: filename

    logical nodal(get_dim(rh(1)))

    type(bl_prof_timer), save :: bpt

    call build(bpt, "ml_nd")

    dm                = get_dim(rh(1))
    stime             = parallel_wtime()
    nodal             = .True.
    nlevs             = mla%nlevel
    bottom_solve_time = zero

    if (nghost(rh(nlevs)) .ne. 1) &
        call bl_error("ml_nd requires one ghost cell for the RHS");

    if ( present(abs_eps_in) ) then
       abs_eps = abs_eps_in
    else
       abs_eps = mgt(nlevs)%abs_eps
    end if

    allocate(soln(nlevs), uu(nlevs), uu_hold(2:nlevs-1), res(nlevs))
    allocate(temp_res(nlevs))
    allocate(brs_flx(2:nlevs))
    allocate(zero_rh(2:nlevs))

    do n = 2,nlevs-1
       la = mla%la(n)
       call multifab_build( uu_hold(n), la, 1, 1, nodal)
       call setval( uu_hold(n), ZERO,all=.true.)
    end do

    do n = nlevs, 1, -1

       la = mla%la(n)
       call multifab_build(    soln(n), la, 1, 1, nodal)
       call multifab_build(      uu(n), la, 1, 1, nodal)
       call multifab_build(     res(n), la, 1, 1, nodal)
       call multifab_build(temp_res(n), la, 1, 1, nodal)
       call setval(    soln(n), ZERO,all=.true.)
       call setval(      uu(n), ZERO,all=.true.)
       call setval(     res(n), ZERO,all=.true.)
       call setval(temp_res(n), ZERO,all=.true.)

       if ( n == 1 ) exit

       call multifab_build(zero_rh(n), la, 1, nghost(rh(n)), nodal)
       call setval(zero_rh(n), ZERO,all=.true.)

       ! Build the (coarse resolution) flux registers to be used in computing
       !  the residual at a non-finest AMR level.

       pdc = layout_get_pd(mla%la(n-1))
       lac = mla%la(n-1)
       call bndry_reg_rr_build(brs_flx(n), la, lac, ref_ratio(n-1,:), pdc, nodal = nodal, other = .false.)

    end do
    !
    ! Let's elide some reductions by doing these reductions together.
    !
    bnorm = ml_norm_inf(rh,fine_mask,local=.true.)

    Anorm = stencil_norm(mgt(nlevs)%ss(mgt(nlevs)%nlevels),local=.true.)
    do n = 1, nlevs-1
       Anorm = max(stencil_norm(mgt(n)%ss(mgt(n)%nlevels),fine_mask(n),local=.true.), Anorm)
    end do

    do n = nlevs,1,-1
       mglev = mgt(n)%nlevels
       call compute_defect(mgt(n)%ss(mglev),res(n),rh(n),full_soln(n),mgt(n)%mm(mglev), &
                           mgt(n)%stencil_type, mgt(n)%lcross, mgt(n)%uniform_dh)
    end do

    do n = 1,nlevs
       call multifab_copy(rh(n),res(n),ng=nghost(rh(n)))
    end do

    t1(1) = bnorm
    t1(2) = Anorm
    t1(3) = ml_norm_inf(rh,fine_mask,local=.true.)

    call parallel_reduce(t2, t1, MPI_MAX)

    bnorm = t2(1)
    Anorm = t2(2)
    tres0 = t2(3)

    if ( parallel_IOProcessor() .and. mgt(nlevs)%verbose > 0 ) then
       write(unit=*, &
             fmt='("F90mg: Initial rhs                  = ",g15.8)') bnorm
       write(unit=*, &
             fmt='("F90mg: Initial residual (resid0)    = ",g15.8)') tres0
    end if

    bnorm = max(bnorm, tres0)
    ! wqz.  Otherwise, we may sometimes have trouble if tres0 is much larger than bnorm

    ! ****************************************************************************

    fine_converged = .false.

    if ( ml_converged(res, fine_mask, bnorm, &
                      rel_eps, abs_eps, mgt(nlevs)%verbose) ) then
       if ( parallel_IOProcessor() .and. mgt(nlevs)%verbose > 0 ) &
            write(unit=*, fmt='("F90mg: No iterations needed ")')

    else   ! Not already converged

     do iter = 1, mgt(nlevs)%max_iter

       if ( (iter .eq. 1) .or. fine_converged ) then
          if ( ml_converged(res, fine_mask, bnorm, rel_eps, abs_eps, mgt(nlevs)%verbose) ) exit
       end if

       ! Set: uu = 0
       do n = 1,nlevs
          call setval(uu(n), ZERO, all=.true.)
       end do

       ! Set: uu_hold = 0
       do n = 2,nlevs-1
          call setval(uu_hold(n), ZERO, all=.true.)
       end do

       !   Down the V-cycle
       do n = nlevs,1,-1

          mglev = mgt(n)%nlevels

          if ( do_diagnostics == 1 ) then
             tres = norm_inf(res(n))
             if ( parallel_ioprocessor() ) then
                print *,'DWN: RES BEFORE GSRB AT LEVEL ',n, tres
             end if
          end if

          ! Relax ...
          if (n > 1) then
             call mini_cycle(mgt(n), mglev, mgt(n)%ss(mglev), &
                  uu(n), res(n), mgt(n)%mm(mglev), mgt(n)%nu1, mgt(n)%nu2)
          else
             call mg_tower_cycle(mgt(n), mgt(n)%cycle_type, mglev, mgt(n)%ss(mglev), &
                  uu(n), res(n), mgt(n)%mm(mglev), mgt(n)%nu1, mgt(n)%nu2, &
                  bottom_solve_time = bottom_solve_time)
          end if

          ! Add: soln += uu
          call plus_plus(soln(n),uu(n))

          if (n > 1) then
             mglev_crse = mgt(n-1)%nlevels

             ! Compute COARSE Res = Rh - Lap(Soln)
             call compute_defect(mgt(n-1)%ss(mglev_crse),res(n-1), &
                  rh(n-1),soln(n-1),mgt(n-1)%mm(mglev_crse), &
                  mgt(n-1)%stencil_type, mgt(n-1)%lcross, mgt(n-1)%uniform_dh)

             if ( dm .eq. 3 ) then
                fac = (8.0_dp_t)**(ref_ratio(n-1,1)/2)
                call multifab_mult_mult_s(res(n-1),fac,nghost(res(n-1)))
             else if ( dm .eq. 2 ) then
                fac = (4.0_dp_t)**(ref_ratio(n-1,1)/2)
                call multifab_mult_mult_s(res(n-1),fac,nghost(res(n-1)))
             end if

             ! Compute FINE Res = Res - Lap(uu)
             mglev = mgt(n)%nlevels
             call compute_defect(mgt(n)%ss(mglev), temp_res(n), &
                  res(n),uu(n),mgt(n)%mm(mglev), &
                  mgt(n)%stencil_type, mgt(n)%lcross, mgt(n)%uniform_dh)
             call multifab_copy(res(n),temp_res(n),ng=nghost(res(n)))

             if ( do_diagnostics == 1 ) then
                tres = norm_inf(res(n))
                if ( parallel_ioprocessor() ) &
                   print *,'DWN: RES AFTER  GSRB AT LEVEL ',n, tres
             end if

             ! Restrict FINE Res to COARSE Res
             call ml_restriction(res(n-1), res(n), mgt(n)%mm(mglev),&
                                 mgt(n-1)%mm(mglev_crse), ref_ratio(n-1,:))

             ! Compute CRSE-FINE Res = Rh - Lap(Soln)
             pdc = layout_get_pd(mla%la(n-1))
             call crse_fine_residual_nodal(n,mgt,brs_flx(n),res(n-1),zero_rh(n),temp_res(n),temp_res(n-1), &
                  soln(n-1),soln(n),one_sided_ss(n),ref_ratio(n-1,:),pdc)

             ! Copy u_hold = uu
             if (n < nlevs) call multifab_copy(uu_hold(n),uu(n),ng=nghost(uu(n)))

             ! Set: uu = 0
             call setval(uu(n),ZERO,all=.true.)

             if ( dm .eq. 3 ) then
                fac = 1.0_dp_t / (8.0_dp_t)**(ref_ratio(n-1,1)/2)
                call multifab_mult_mult_s(res(n-1),fac,nghost(res(n-1)))
             else if ( dm .eq. 2 ) then
                fac = 1.0_dp_t / (4.0_dp_t)**(ref_ratio(n-1,1)/2)
                call multifab_mult_mult_s(res(n-1),fac,nghost(res(n-1)))
             end if

          else

             if (do_diagnostics == 1 ) then
                call compute_defect(mgt(n)%ss(mglev),temp_res(n), res(n),uu(n),mgt(n)%mm(mglev), &
                                    mgt(n)%stencil_type, mgt(n)%lcross, mgt(n)%uniform_dh)
                tres = norm_inf(temp_res(n))
                if ( parallel_ioprocessor() ) then
                   print *,'DWN: RES AFTER  GSRB AT LEVEL ',n, tres
                end if
             end if

          end if

        end do

        !   Back up the V-cycle
        do n = 2, nlevs

          pd = layout_get_pd(mla%la(n))
          mglev = mgt(n)%nlevels

          ! Interpolate uu from coarser level
          if (iter == 1) call plus_plus(uu(n-1),  full_soln(n-1))
          call ml_nodal_prolongation(uu(n), uu(n-1), ref_ratio(n-1,:))
          if (iter == 1) call sub_sub(uu(n-1), full_soln(n-1))

          ! Subtract: uu -= full_soln
          !     Must do this in order to remove interpolated full_soln...
          if (iter == 1) call sub_sub(uu(n),full_soln(n))

          ! Add: soln += uu
          call plus_plus(soln(n), uu(n), nghost(uu(n)))

          ! Add: uu_hold += uu
          if (n < nlevs) call plus_plus(uu_hold(n), uu(n), nghost(uu(n)))

          ! Compute Res = Res - Lap(uu)
          call compute_defect(mgt(n)%ss(mglev),temp_res(n),res(n),uu(n),mgt(n)%mm(mglev), &
                              mgt(n)%stencil_type, mgt(n)%lcross, mgt(n)%uniform_dh)
          call multifab_copy(res(n),temp_res(n),ng=nghost(res(n)))

          if ( do_diagnostics == 1 ) then
             tres = norm_inf(res(n))
             if ( parallel_ioprocessor() ) then
                print *,'UP : RES BEFORE GSRB AT LEVEL ',n, tres
             end if
          end if

          ! Set: uu = 0
          call setval(uu(n),ZERO,all=.true.)

          ! Relax ...
          call mini_cycle(mgt(n), mglev, mgt(n)%ss(mglev), &
               uu(n), res(n), mgt(n)%mm(mglev), mgt(n)%nu1, mgt(n)%nu2)

          ! Compute Res = Res - Lap(uu)
          call compute_defect(mgt(n)%ss(mglev),temp_res(n),res(n),uu(n),mgt(n)%mm(mglev), &
                              mgt(n)%stencil_type, mgt(n)%lcross, mgt(n)%uniform_dh)
          call multifab_copy(res(n),temp_res(n),ng=nghost(res(n)))

          if ( do_diagnostics == 1 ) then
             tres = norm_inf(res(n))
             if ( parallel_ioprocessor() ) then
                print *,'UP : RES AFTER  GSRB AT LEVEL ',n, tres
                if (n == nlevs) print *,' '
             end if
          end if

          ! Add: soln += uu
          call plus_plus(soln(n), uu(n), nghost(uu(n)))

          ! Add: uu += uu_hold so that it will be interpolated too.
          if (n < nlevs) call plus_plus(uu(n), uu_hold(n), nghost(uu(n)))

       end do

       !    Inject the solution to the coarser grids.
       do n = nlevs,2,-1
          mglev      = mgt(n)%nlevels
          mglev_crse = mgt(n-1)%nlevels
          call ml_restriction(soln(n-1), soln(n), mgt(n)%mm(mglev), &
               mgt(n-1)%mm(mglev_crse), ref_ratio(n-1,:), inject = .true.)
       end do

        do n = 1,nlevs
           call multifab_fill_boundary(soln(n), cross = mgt(n)%lcross)
        end do

       !    Optimization so don't have to do multilevel convergence test each time

       !    Compute the residual on just the finest level
        n = nlevs
        mglev = mgt(n)%nlevels
        call compute_defect(mgt(n)%ss(mglev),res(n),rh(n),soln(n),mgt(n)%mm(mglev), &
                            mgt(n)%stencil_type, mgt(n)%lcross, mgt(n)%uniform_dh)

        if ( ml_fine_converged(res, bnorm, rel_eps, abs_eps) ) then

          fine_converged = .true.

          !      Compute the residual on every level
          do n = 1,nlevs-1
             mglev = mgt(n)%nlevels
             call compute_defect(mgt(n)%ss(mglev),res(n),rh(n),soln(n),mgt(n)%mm(mglev), &
                                 mgt(n)%stencil_type, mgt(n)%lcross, mgt(n)%uniform_dh)
          end do

          do n = nlevs,2,-1
             !  Restrict the finer residual onto the coarser grid
             mglev      = mgt(n  )%nlevels
             mglev_crse = mgt(n-1)%nlevels
             if ( dm .eq. 3 ) then
                fac = (8.0_dp_t)**(ref_ratio(n-1,1)/2)
                call multifab_mult_mult_s(res(n-1),fac,nghost(res(n-1)))
             else if ( dm .eq. 2 ) then
                fac = (4.0_dp_t)**(ref_ratio(n-1,1)/2)
                call multifab_mult_mult_s(res(n-1),fac,nghost(res(n-1)))
             end if
             call ml_restriction(res(n-1), res(n), mgt(n)%mm(mglev),&
                  mgt(n-1)%mm(mglev_crse), ref_ratio(n-1,:))

             !  Compute the coarse-fine residual at coarse-fine nodes
             pdc = layout_get_pd(mla%la(n-1))
             call crse_fine_residual_nodal(n,mgt,brs_flx(n),res(n-1), &
                  zero_rh(n),temp_res(n),temp_res(n-1), &
                  soln(n-1),soln(n),one_sided_ss(n),ref_ratio(n-1,:),pdc)
             if ( dm .eq. 3 ) then
                fac = 1.0_dp_t / (8.0_dp_t)**(ref_ratio(n-1,1)/2)
                call multifab_mult_mult_s(res(n-1),fac,nghost(res(n-1)))
             else if ( dm .eq. 2 ) then
                fac = 1.0_dp_t / (4.0_dp_t)**(ref_ratio(n-1,1)/2)
                call multifab_mult_mult_s(res(n-1),fac,nghost(res(n-1)))
             end if
          end do

          if ( mgt(nlevs)%verbose > 1 ) then
             do n = 1,nlevs
                ttres = norm_inf(res(n),local=.true.)
                call parallel_reduce(tres, ttres, MPI_MAX, proc = parallel_IOProcessorNode())
                if ( parallel_IOProcessor() ) then
                   write(unit=*, fmt='("F90mg: Iteration   ",i3," Lev ",i1," resid/resid0 = ",g15.8)') &
                        iter,n,tres/tres0
                end if
             end do
          end if

       else

          fine_converged = .false.
          if ( mgt(nlevs)%verbose > 1 ) then

             if ( .false. ) then
                !
                ! Some debugging code I want to keep around for a while.
                ! I don't want to have to recreate this all the time :-)
                !
                write(number,fmt='(i3.3)') iter
                filename = 'res_fine_iter=' // number
                call fabio_write(res(nlevs), 'debug', trim(filename))
             end if

             ttres = norm_inf(res(nlevs),local=.true.)
             call parallel_reduce(tres, ttres, MPI_MAX, proc = parallel_IOProcessorNode())
             if ( parallel_IOProcessor() ) then
                write(unit=*, fmt='("F90mg: Iteration   ",i3," Fine  resid/resid0 = ",g15.8)') iter,tres/tres0
             end if
          end if

       end if

     end do

     iter = iter-1
     if (iter < mgt(nlevs)%max_iter) then
        if ( mgt(nlevs)%verbose > 0 ) then
           !
           ! Consolidate these reductions.
           !
           t1(1) = ml_norm_inf(res,fine_mask,local=.true.)
           t1(2) = (parallel_wtime() - stime)
           t1(3) = bottom_solve_time

           call parallel_reduce(t2, t1, MPI_MAX, proc = parallel_IOProcessorNode())

           if ( parallel_IOProcessor() ) then
              if ( tres0 .gt. 0.0_dp_t) then
                 write(unit=*, fmt='("F90mg: Final Iter. ",i3," resid/resid0 = ",g15.8)') iter,t2(1)/tres0
                 write(unit=*, fmt='("F90mg: Solve time: ",g13.6, " Bottom Solve time: ", g13.6)') t2(2), t2(3)
                 write(unit=*, fmt='("")')
              else
                 write(unit=*, fmt='("F90mg: Final Iter. ",i3," resid/resid0 = ",g15.8)') iter,0.0_dp_t
                 write(unit=*, fmt='("F90mg: Solve time: ",g13.6, " Bottom Solve time: ", g13.6)') t2(2), t2(3)
                 write(unit=*, fmt='("")')
              end if
           end if
        end if
     else
        if (mgt(nlevs)%abort_on_max_iter) then
           call bl_error("Multigrid Solve: failed to converge in max_iter iterations")
        end if
     end if

    ! Add: full_soln += soln
    do n = 1,nlevs
       call plus_plus(full_soln(n),soln(n))
    end do

    end if

    ! ****************************************************************************

    do n = 2,nlevs-1
       call multifab_destroy(uu_hold(n))
    end do

    do n = nlevs, 1, -1
       call multifab_destroy(    soln(n))
       call multifab_destroy(      uu(n))
       call multifab_destroy(     res(n))
       call multifab_destroy(temp_res(n))
       if ( n == 1 ) exit
       call multifab_destroy( zero_rh(n))
       call bndry_reg_destroy(brs_flx(n))
    end do

    call destroy(bpt)

  contains

    subroutine crse_fine_residual_nodal(n,mgt,brs_flx,crse_res,fine_rhs,temp_res,temp_crse_res, &
         crse_soln,fine_soln,one_sided_ss,ref_ratio,pdc)

      use nodal_interface_stencil_module , only : ml_crse_contrib, ml_fine_contrib

      integer        , intent(in   ) :: n
      type(mg_tower) , intent(inout) :: mgt(:)
      type(bndry_reg), intent(inout) :: brs_flx
      type(multifab) , intent(inout) :: crse_res
      type(multifab) , intent(in   ) :: fine_rhs
      type(multifab) , intent(inout) :: temp_res
      type(multifab) , intent(inout) :: temp_crse_res
      type(multifab) , intent(inout) :: crse_soln
      type(multifab) , intent(inout) :: fine_soln
      type(multifab) , intent(in   ) :: one_sided_ss
      integer        , intent(in   ) :: ref_ratio(:)
      type(box)      , intent(in   ) :: pdc

      integer :: i,dm,mglev_crse,mglev_fine

      mglev_crse = mgt(n-1)%nlevels
      mglev_fine = mgt(n  )%nlevels
      dm        = get_dim(temp_res)

      !    Compute the fine contributions at faces, edges and corners.

      !    First compute a residual which only takes contributions from the
      !       grid on which it is calculated.

      if (mgt(n)%lcross) then
        call grid_res(one_sided_ss,temp_res, &
             fine_rhs,fine_soln,mgt(n)%mm(mglev_fine),mgt(n)%face_type, &
             mgt(n)%lcross, mgt(n)%uniform_dh)
      else
        call grid_res(mgt(n)%ss(mglev_fine),temp_res, &
             fine_rhs,fine_soln,mgt(n)%mm(mglev_fine),mgt(n)%face_type, &
             mgt(n)%lcross, mgt(n)%uniform_dh)
      end if

      !    Zero out the flux registers which will hold the fine contributions
      call bndry_reg_setval(brs_flx, ZERO, all = .true.)

      do i = 1,dm
         call ml_fine_contrib(brs_flx%bmf(i,0), &
              temp_res,mgt(n)%mm(mglev_fine),ref_ratio,pdc,-i)
         call ml_fine_contrib(brs_flx%bmf(i,1), &
              temp_res,mgt(n)%mm(mglev_fine),ref_ratio,pdc,+i)
      end do

!     Compute the crse contributions at edges and corners and add to fine contributions
!        in temp_crse_res (need to do this in a temporary for periodic issues)
      call setval(temp_crse_res,ZERO,all=.true.)
      do i = 1,dm
         call ml_crse_contrib(temp_crse_res, brs_flx%bmf(i,0), crse_soln, &
              mgt(n-1)%ss(mgt(n-1)%nlevels), &
              mgt(n-1)%mm(mglev_crse), &
              mgt(n  )%mm(mglev_fine), &
              pdc,ref_ratio, -i)
         call ml_crse_contrib(temp_crse_res, brs_flx%bmf(i,1), crse_soln, &
              mgt(n-1)%ss(mgt(n-1)%nlevels), &
              mgt(n-1)%mm(mglev_crse), &
              mgt(n  )%mm(mglev_fine), &
              pdc,ref_ratio, +i)
      end do

!     Add to res(n-1).
      call plus_plus(crse_res,temp_crse_res)

      call periodic_add_copy(crse_res,temp_crse_res,synced=.true.)

!     Clear temp_crse_res (which is temp_res(n-1) from calling routine) just in case...
      call setval(temp_crse_res,ZERO,all=.true.)

    end subroutine crse_fine_residual_nodal

    function ml_fine_converged(res, bnorm, rel_eps, abs_eps) result(r)
      logical :: r
      type(multifab), intent(in) :: res(:)
      real(dp_t), intent(in) :: rel_eps, abs_eps, bnorm
      real(dp_t) :: ni_res
      integer    :: nlevs
      nlevs = size(res)
      ni_res = norm_inf(res(nlevs))
      r = ( ni_res <= rel_eps*(bnorm) .or. ni_res <= abs_eps )
    end function ml_fine_converged

    function ml_converged(res, mask, bnorm, rel_eps, abs_eps, verbose) result(r)
      use ml_norm_module, only : ml_norm_inf
      logical :: r
      integer :: verbose
      type(multifab), intent(in) :: res(:)
      type(lmultifab), intent(in) :: mask(:)
      real(dp_t), intent(in) :: rel_eps, abs_eps, bnorm
      real(dp_t) :: ni_res
      ni_res = ml_norm_inf(res, mask)
      r = ( ni_res <= rel_eps*(bnorm) .or. ni_res <= abs_eps )
      if ( r .and. parallel_IOProcessor() .and. verbose > 1) then
         if (ni_res <= rel_eps*bnorm) then
            print *,'Converged res < rel_eps*bnorm '
         else if (ni_res <= abs_eps) then
            print *,'Converged res < abs_eps '
         end if
      end if
    end function ml_converged

  end subroutine ml_nd

  subroutine grid_res(ss, dd, ff, uu, mm, face_type, lcross, uniform_dh)

    type(multifab), intent(in)    :: ff, ss
    type(multifab), intent(inout) :: dd, uu
    type(imultifab), intent(in)   :: mm
    integer, intent(in) :: face_type(:,:,:)
    logical, intent(in) :: lcross
    logical, intent(in) :: uniform_dh

    integer :: i, n
    real(kind=dp_t), pointer :: dp(:,:,:,:)
    real(kind=dp_t), pointer :: fp(:,:,:,:)
    real(kind=dp_t), pointer :: up(:,:,:,:)
    real(kind=dp_t), pointer :: sp(:,:,:,:)
    integer        , pointer :: mp(:,:,:,:)
    integer :: dm,nodal_ng
    type(bl_prof_timer), save :: bpt

    call build(bpt, "grid_res")

    nodal_ng = 0; if ( nodal_q(uu) ) nodal_ng = 1

    dm = get_dim(uu)

    call multifab_fill_boundary(uu, cross = lcross)

    do i = 1, nfabs(uu)
       dp  => dataptr(dd, i)
       fp  => dataptr(ff, i)
       up  => dataptr(uu, i)
       sp  => dataptr(ss, i)
       mp  => dataptr(mm, i)
       do n = 1, ncomp(uu)
          select case(dm)
          case (1)
             call grid_laplace_1d(sp(:,:,1,1), dp(:,1,1,n), fp(:,1,1,n), up(:,1,1,n), &
                                  mp(:,1,1,1),nghost(uu))
          case (2)
             call grid_laplace_2d(sp(:,:,:,1), dp(:,:,1,n), fp(:,:,1,n), up(:,:,1,n), &
                                  mp(:,:,1,1), nghost(uu), face_type(i,:,:))
          case (3)
             call grid_laplace_3d(sp(:,:,:,:), dp(:,:,:,n), fp(:,:,:,n), up(:,:,:,n), &
                                  mp(:,:,:,1), nghost(uu), face_type(i,:,:), uniform_dh)
          end select
       end do
    end do

    call destroy(bpt)

  end subroutine grid_res

  subroutine grid_laplace_1d(ss, dd, ff, uu, mm, ng)

    use impose_neumann_bcs_module

    integer           , intent(in   ) :: ng
    real (kind = dp_t), intent(in   ) :: ss(0:,:)
    real (kind = dp_t), intent(inout) :: dd(0:)
    real (kind = dp_t), intent(in   ) :: ff(0:)
    integer,            intent(in   ) :: mm(:)
    real (kind = dp_t), intent(inout) :: uu(1-ng:)

    integer :: i,nx,lo(1)
    nx = size(ss,dim=2)-1

    lo = 1
    call impose_neumann_bcs_1d(uu,mm,lo,ng)

    i = 1
    dd(i) = HALF*ff(i) - (ss(0,i)*uu(i) + ss(1,i)*(uu(i+1)-uu(i)))

    do i = 2,nx
      dd(i) = ff(i) - &
              (ss(0,i)*uu(i) + ss(1,i) * uu(i+1) + ss(2,i) * uu(i-1))
    end do

    i = nx+1
    dd(i) = HALF*ff(i) - (ss(2,i)*(uu(i-1)-uu(i)))

  end subroutine grid_laplace_1d

  subroutine grid_laplace_2d(ss, dd, ff, uu, mm, ng, face_type)

    use bc_module
    use impose_neumann_bcs_module

    integer           , intent(in   ) :: ng
    real (kind = dp_t), intent(in   ) :: ss(0:,:,:)
    real (kind = dp_t), intent(inout) :: dd(0:,0:)
    real (kind = dp_t), intent(in   ) :: ff(0:,0:)
    real (kind = dp_t), intent(inout) :: uu(1-ng:,1-ng:)
    integer,            intent(in   ) :: mm(:,:)
    integer,            intent(in ) :: face_type(:,:)
    integer :: i,j,nx,ny,lo(2)
    integer :: istart,iend,jstart,jend

    nx = size(ss,dim=2)-1
    ny = size(ss,dim=3)-1

    lo = 1
    call impose_neumann_bcs_2d(uu,mm,lo,ng)

    if (face_type(1,1) == BC_NEU) then
      istart = 1
    else
      istart = 2
    end if
    if (face_type(1,2) == BC_NEU) then
      iend = nx+1
    else
      iend = nx
    end if
    if (face_type(2,1) == BC_NEU) then
      jstart = 1
    else
      jstart = 2
    end if
    if (face_type(2,2) == BC_NEU) then
      jend = ny+1
    else
      jend = ny
    end if

    if (size(ss,dim=1) .eq. 9) then
!     Corners
      i = 1
      j = 1
      dd(i,j) = ss(8,i,j)*(uu(i+1,j+1)+HALF*uu(i,j+1)+HALF*uu(i+1,j)-TWO*uu(i,j))
      dd(i,j) = FOURTH*ff(i,j) - dd(i,j)

      i = 1
      j = ny+1
      dd(i,j) = ss(3,i,j)*(uu(i+1,j-1)+HALF*uu(i,j-1)+HALF*uu(i+1,j)-TWO*uu(i,j))
      dd(i,j) = FOURTH*ff(i,j) - dd(i,j)

      i = nx+1
      j = 1
      dd(i,j) = ss(6,i,j)*(uu(i-1,j+1)+HALF*uu(i,j+1)+HALF*uu(i-1,j)-TWO*uu(i,j))
      dd(i,j) = FOURTH*ff(i,j) - dd(i,j)

      i = nx+1
      j = ny+1
      dd(i,j) = ss(1,i,j)*(uu(i-1,j-1)+HALF*uu(i,j-1)+HALF*uu(i-1,j)-TWO*uu(i,j))
      dd(i,j) = FOURTH*ff(i,j) - dd(i,j)

!     Lo-x edge
      i = 1
      do j = jstart,jend
         dd(i,j) = ss(3,i,j)*(uu(i+1,j-1)+HALF*uu(i,j-1)+HALF*uu(i+1,j)-TWO*uu(i,j)) &
                  +ss(8,i,j)*(uu(i+1,j+1)+HALF*uu(i,j+1)+HALF*uu(i+1,j)-TWO*uu(i,j))
         dd(i,j) = HALF*ff(i,j) - dd(i,j)
      end do

!     Hi-x edge
      i = nx+1
      do j = jstart,jend
         dd(i,j) = ss(1,i,j)*(uu(i-1,j-1)+HALF*uu(i,j-1)+HALF*uu(i-1,j)-TWO*uu(i,j)) &
                  +ss(6,i,j)*(uu(i-1,j+1)+HALF*uu(i,j+1)+HALF*uu(i-1,j)-TWO*uu(i,j))
         dd(i,j) = HALF*ff(i,j) - dd(i,j)
      end do

!     Lo-y edge
      j = 1
      do i = istart,iend
         dd(i,j) = ss(6,i,j)*(uu(i-1,j+1)+HALF*uu(i,j+1)+HALF*uu(i-1,j)-TWO*uu(i,j)) &
                  +ss(8,i,j)*(uu(i+1,j+1)+HALF*uu(i,j+1)+HALF*uu(i+1,j)-TWO*uu(i,j))
         dd(i,j) = HALF*ff(i,j) - dd(i,j)
      end do

!     Hi-y edge
      j = ny+1
      do i = istart,iend
         dd(i,j) = ss(1,i,j)*(uu(i-1,j-1)+HALF*uu(i,j-1)+HALF*uu(i-1,j)-TWO*uu(i,j)) &
                  +ss(3,i,j)*(uu(i+1,j-1)+HALF*uu(i,j-1)+HALF*uu(i+1,j)-TWO*uu(i,j))
         dd(i,j) = HALF*ff(i,j) - dd(i,j)
      end do

!     Interior
      do j = jstart,jend
      do i = istart,iend
         dd(i,j) = ss(0,i,j)*uu(i,j) + ss(1,i,j) * uu(i-1,j-1) &
                                     + ss(2,i,j) * uu(i  ,j-1) &
                                     + ss(3,i,j) * uu(i+1,j-1) &
                                     + ss(4,i,j) * uu(i-1,j  ) &
                                     + ss(5,i,j) * uu(i+1,j  ) &
                                     + ss(6,i,j) * uu(i-1,j+1) &
                                     + ss(7,i,j) * uu(i  ,j+1) &
                                     + ss(8,i,j) * uu(i+1,j+1)
         dd(i,j) = ff(i,j) - dd(i,j)
      end do
      end do

    else if (size(ss,dim=1) .eq. 5) then

!     Corners
      i = 1
      j = 1
      dd(i,j) = ss(1,i,j)*(uu(i+1,j)-uu(i,j)) + ss(3,i,j)*(uu(i,j+1)-uu(i,j))
      dd(i,j) = FOURTH*ff(i,j) - dd(i,j)

      i = 1
      j = ny+1
      dd(i,j) = ss(1,i,j)*(uu(i+1,j)-uu(i,j)) + ss(4,i,j)*(uu(i,j-1)-uu(i,j))
      dd(i,j) = FOURTH*ff(i,j) - dd(i,j)

      i = nx+1
      j = 1
      dd(i,j) = ss(2,i,j)*(uu(i-1,j)-uu(i,j)) + ss(3,i,j)*(uu(i,j+1)-uu(i,j))
      dd(i,j) = FOURTH*ff(i,j) - dd(i,j)

      i = nx+1
      j = ny+1
      dd(i,j) = ss(2,i,j)*(uu(i-1,j)-uu(i,j)) + ss(4,i,j)*(uu(i,j-1)-uu(i,j))
      dd(i,j) = FOURTH*ff(i,j) - dd(i,j)

!     Lo-x edge
      i = 1
      do j = jstart,jend
         dd(i,j) =  ss(1,i,j)*(uu(i+1,j)-uu(i,j)) &
                   +ss(3,i,j)*(uu(i,j+1)-uu(i,j)) &
                   +ss(4,i,j)*(uu(i,j-1)-uu(i,j))
         dd(i,j) = HALF*ff(i,j) - dd(i,j)
      end do

!     Hi-x edge
      i = nx+1
      do j = jstart,jend
         dd(i,j) = ss(2,i,j)*(uu(i-1,j)-uu(i,j)) &
                  +ss(3,i,j)*(uu(i,j+1)-uu(i,j)) &
                  +ss(4,i,j)*(uu(i,j-1)-uu(i,j))
         dd(i,j) = HALF*ff(i,j) - dd(i,j)
      end do

!     Lo-y edge
      j = 1
      do i = istart,iend
         dd(i,j) = ss(3,i,j)*(uu(i,j+1)-uu(i,j)) &
                  +ss(1,i,j)*(uu(i+1,j)-uu(i,j)) &
                  +ss(2,i,j)*(uu(i-1,j)-uu(i,j))
         dd(i,j) = HALF*ff(i,j) - dd(i,j)
      end do

!     Hi-y edge
      j = ny+1
      do i = istart,iend
         dd(i,j) = ss(4,i,j)*(uu(i,j-1)-uu(i,j)) &
                  +ss(1,i,j)*(uu(i+1,j)-uu(i,j)) &
                  +ss(2,i,j)*(uu(i-1,j)-uu(i,j))
         dd(i,j) = HALF*ff(i,j) - dd(i,j)
      end do

!     Interior
      do j = jstart,jend
      do i = istart,iend
         dd(i,j) = ss(0,i,j)*uu(i,j) + ss(1,i,j) * uu(i+1,j  ) &
                                     + ss(2,i,j) * uu(i-1,j  ) &
                                     + ss(3,i,j) * uu(i  ,j+1) &
                                     + ss(4,i,j) * uu(i  ,j-1)
         dd(i,j) = ff(i,j) - dd(i,j)
      end do
      end do

    end if

  end subroutine grid_laplace_2d

  subroutine grid_laplace_3d(ss, dd, ff, uu, mm, ng, face_type, uniform_dh)
    use bc_module
    use impose_neumann_bcs_module
    integer, intent(in) :: ng
    real (kind = dp_t), intent(in   ) :: ff(0:,0:,0:)
    real (kind = dp_t), intent(inout) :: uu(1-ng:,1-ng:,1-ng:)
    real (kind = dp_t), intent(inout) :: dd(0:,0:,0:)
    real (kind = dp_t), intent(in   ) :: ss(0:,:,:,:)
    integer,            intent(in   ) :: mm(:,:,:)
    integer,            intent(in   ) :: face_type(:,:)
    logical,            intent(in)    :: uniform_dh

    integer :: i, j, k, lo(3)
    integer :: istart, iend, jstart, jend, kstart, kend
    integer :: nx, ny, nz

    nx = size(ss,dim=2)-1
    ny = size(ss,dim=3)-1
    nz = size(ss,dim=4)-1

    lo = 1
    call impose_neumann_bcs_3d(uu,mm,lo,ng)

    if (face_type(1,1) == BC_NEU) then
      istart = 1
    else
      istart = 2
    end if
    if (face_type(1,2) == BC_NEU) then
      iend = nx+1
    else
      iend = nx
    end if
    if (face_type(2,1) == BC_NEU) then
      jstart = 1
    else
      jstart = 2
    end if
    if (face_type(2,2) == BC_NEU) then
      jend = ny+1
    else
      jend = ny
    end if
    if (face_type(3,1) == BC_NEU) then
      kstart = 1
    else
      kstart = 2
    end if
    if (face_type(3,2) == BC_NEU) then
      kend = nz+1
    else
      kend = nz
    end if

    if (size(ss,dim=1) .eq. 27 .or. size(ss,dim=1) .eq. 21) then
       !
       !     Corners
       !
      i = 1
      j = 1
      k = 1
      dd(i,j,k) = ss(20,i,j,k)*(uu(i+1,j+1,k+1) + uu(i+1,j+1,k  ) &
                               +uu(i+1,j  ,k+1) + uu(i  ,j+1,k+1) &
                         - FOUR*uu(i  ,j  ,k) )
      dd(i,j,k) = EIGHTH*ff(i,j,k) - dd(i,j,k)

      i = 1
      j = ny+1
      k = 1
      dd(i,j,k) = ss(15,i,j,k)*(uu(i+1,j-1,k+1) + uu(i+1,j-1,k  ) &
                               +uu(i+1,j  ,k+1) + uu(i  ,j-1,k+1) &
                         - FOUR*uu(i  ,j  ,k) )
      dd(i,j,k) = EIGHTH*ff(i,j,k) - dd(i,j,k)

      i = nx+1
      j = 1
      k = 1
      dd(i,j,k) = ss(18,i,j,k)*(uu(i-1,j+1,k+1) + uu(i-1,j+1,k  ) &
                               +uu(i-1,j  ,k+1) + uu(i  ,j+1,k+1) &
                         - FOUR*uu(i  ,j  ,k) )
      dd(i,j,k) = EIGHTH*ff(i,j,k) - dd(i,j,k)

      i = nx+1
      j = ny+1
      k = 1
      dd(i,j,k) = ss(13,i,j,k)*(uu(i-1,j-1,k+1) + uu(i-1,j-1,k  ) &
                               +uu(i-1,j  ,k+1) + uu(i  ,j-1,k+1) &
                         - FOUR*uu(i  ,j  ,k) )
      dd(i,j,k) = EIGHTH*ff(i,j,k) - dd(i,j,k)

      i = 1
      j = 1
      k = nz+1
      dd(i,j,k) = ss( 8,i,j,k)*(uu(i+1,j+1,k-1) + uu(i+1,j+1,k  ) &
                               +uu(i+1,j  ,k-1) + uu(i  ,j+1,k-1) &
                         - FOUR*uu(i  ,j  ,k) )
      dd(i,j,k) = EIGHTH*ff(i,j,k) - dd(i,j,k)

      i = 1
      j = ny+1
      k = nz+1
      dd(i,j,k) = ss( 3,i,j,k)*(uu(i+1,j-1,k-1) + uu(i+1,j-1,k  ) &
                               +uu(i+1,j  ,k-1) + uu(i  ,j-1,k-1) &
                         - FOUR*uu(i  ,j  ,k) )
      dd(i,j,k) = EIGHTH*ff(i,j,k) - dd(i,j,k)

      i = nx+1
      j = 1
      k = nz+1
      dd(i,j,k) = ss( 6,i,j,k)*(uu(i-1,j+1,k-1) + uu(i-1,j+1,k  ) &
                               +uu(i-1,j  ,k-1) + uu(i  ,j+1,k-1) &
                         - FOUR*uu(i  ,j  ,k) )
      dd(i,j,k) = EIGHTH*ff(i,j,k) - dd(i,j,k)

      i = nx+1
      j = ny+1
      k = nz+1
      dd(i,j,k) = ss( 1,i,j,k)*(uu(i-1,j-1,k-1) + uu(i-1,j-1,k  ) &
                               +uu(i-1,j  ,k-1) + uu(i  ,j-1,k-1) &
                         - FOUR*uu(i  ,j  ,k) )
      dd(i,j,k) = EIGHTH*ff(i,j,k) - dd(i,j,k)

!     Lo-x / Lo-y edge
      i = 1
      j = 1
      do k = kstart,kend
         dd(i,j,k) = ss(20,i,j,k)*(uu(i+1,j+1,k+1) + uu(i+1,j+1,k  ) &
                                  +uu(i+1,j  ,k+1) + uu(i  ,j+1,k+1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss( 8,i,j,k)*(uu(i+1,j+1,k-1) + uu(i+1,j+1,k  ) &
                                  +uu(i+1,j  ,k-1) + uu(i  ,j+1,k-1) &
                            - FOUR*uu(i  ,j  ,k) )
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Hi-x / Lo-y edge
      !
      i = nx+1
      j = 1
      do k = kstart,kend
         dd(i,j,k) = ss(18,i,j,k)*(uu(i-1,j+1,k+1) + uu(i-1,j+1,k  ) &
                                   +uu(i-1,j  ,k+1) + uu(i  ,j+1,k+1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss( 6,i,j,k)*(uu(i-1,j+1,k-1) + uu(i-1,j+1,k  ) &
                                  +uu(i-1,j  ,k-1) + uu(i  ,j+1,k-1) &
                            - FOUR*uu(i  ,j  ,k) )
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Lo-x / Hi-y edge
      !
      i = 1
      j = ny+1
      do k = kstart,kend
         dd(i,j,k) = ss(15,i,j,k)*(uu(i+1,j-1,k+1) + uu(i+1,j-1,k  ) &
                                  +uu(i+1,j  ,k+1) + uu(i  ,j-1,k+1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss( 3,i,j,k)*(uu(i+1,j-1,k-1) + uu(i+1,j-1,k  ) &
                                  +uu(i+1,j  ,k-1) + uu(i  ,j-1,k-1) &
                            - FOUR*uu(i  ,j  ,k) )
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Hi-x / Hi-y edge
      !
      i = nx+1
      j = ny+1
      do k = kstart,kend
         dd(i,j,k) = ss(13,i,j,k)*(uu(i-1,j-1,k+1) + uu(i-1,j-1,k  ) &
                                  +uu(i-1,j  ,k+1) + uu(i  ,j-1,k+1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss( 1,i,j,k)*(uu(i-1,j-1,k-1) + uu(i-1,j-1,k  ) &
                                  +uu(i-1,j  ,k-1) + uu(i  ,j-1,k-1) &
                            - FOUR*uu(i  ,j  ,k) )
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Lo-x / Lo-z edge
      !
      i = 1
      k = 1
      do j = jstart,jend
         dd(i,j,k) = ss(20,i,j,k)*(uu(i+1,j+1,k+1) + uu(i+1,j+1,k  ) &
                                  +uu(i+1,j  ,k+1) + uu(i  ,j+1,k+1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss(15,i,j,k)*(uu(i+1,j-1,k+1) + uu(i+1,j-1,k  ) &
                                  +uu(i+1,j  ,k+1) + uu(i  ,j-1,k+1) &
                            - FOUR*uu(i  ,j  ,k) )
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Hi-x / Lo-z edge
      !
      i = nx+1
      k = 1
      do j = jstart,jend
         dd(i,j,k) = ss(18,i,j,k)*(uu(i-1,j+1,k+1) + uu(i-1,j+1,k  ) &
                                  +uu(i-1,j  ,k+1) + uu(i  ,j+1,k+1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss(13,i,j,k)*(uu(i-1,j-1,k+1) + uu(i-1,j  ,k+1) &
                                  +uu(i-1,j-1,k  ) + uu(i  ,j-1,k+1) &
                            - FOUR*uu(i  ,j  ,k) )
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Lo-x / Hi-z edge
      !
      i = 1
      k = nz+1
      do j = jstart,jend
         dd(i,j,k) = ss( 8,i,j,k)*(uu(i+1,j+1,k-1) + uu(i+1,j  ,k-1) &
                                  +uu(i+1,j+1,k  ) + uu(i  ,j+1,k-1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss( 3,i,j,k)*(uu(i+1,j-1,k-1) + uu(i+1,j-1,k  ) &
                                  +uu(i+1,j  ,k-1) + uu(i  ,j-1,k-1) &
                            - FOUR*uu(i  ,j  ,k) )
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Hi-x / Hi-z edge
      !
      i = nx+1
      k = nz+1
      do j = jstart,jend
         dd(i,j,k) = ss( 6,i,j,k)*(uu(i-1,j+1,k-1) + uu(i-1,j  ,k-1) &
                                  +uu(i-1,j+1,k  ) + uu(i  ,j+1,k-1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss( 1,i,j,k)*(uu(i-1,j-1,k-1) + uu(i-1,j-1,k  ) &
                                +uu(i-1,j  ,k-1) + uu(i  ,j-1,k-1) &
                            - FOUR*uu(i  ,j  ,k) )
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Lo-y / Lo-z edge
      !
      j = 1
      k = 1
      do i = istart,iend
         dd(i,j,k) = ss(20,i,j,k)*(uu(i+1,j+1,k+1) + uu(i+1,j+1,k  ) &
                                  +uu(i+1,j  ,k+1) + uu(i  ,j+1,k+1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss(18,i,j,k)*(uu(i-1,j+1,k+1) + uu(i-1,j+1,k  ) &
                                  +uu(i  ,j+1,k+1) + uu(i-1,j  ,k+1) &
                            - FOUR*uu(i  ,j  ,k) )
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Hi-y / Lo-z edge
      !
      j = ny+1
      k = 1
      do i = istart,iend
         dd(i,j,k) = ss(15,i,j,k)*(uu(i+1,j-1,k+1) + uu(i+1,j-1,k  ) &
                                  +uu(i  ,j-1,k+1) + uu(i+1,j  ,k+1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss(13,i,j,k)*(uu(i-1,j-1,k+1) + uu(i  ,j-1,k+1) &
                                  +uu(i-1,j-1,k  ) + uu(i-1,j  ,k+1) &
                            - FOUR*uu(i  ,j  ,k) )
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Lo-y / Hi-z edge
      !
      j = 1
      k = nz+1
      do i = istart,iend
         dd(i,j,k) = ss( 8,i,j,k)*(uu(i+1,j+1,k-1) + uu(i+1,j  ,k-1) &
                                  +uu(i+1,j+1,k  ) + uu(i  ,j+1,k-1) &
                          - FOUR*uu(i  ,j  ,k) ) &
                      +ss( 6,i,j,k)*(uu(i-1,j+1,k-1) + uu(i-1,j+1,k  ) &
                                  +uu(i  ,j+1,k-1) + uu(i-1,j  ,k-1) &
                            - FOUR*uu(i  ,j  ,k) )
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Hi-y / Hi-z edge
      !
      j = ny+1
      k = nz+1
      do i = istart,iend
         dd(i,j,k) = ss( 3,i,j,k)*(uu(i+1,j-1,k-1) + uu(i  ,j-1,k-1) &
                                  +uu(i+1,j-1,k  ) + uu(i+1,j  ,k-1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss( 1,i,j,k)*(uu(i-1,j-1,k-1) + uu(i-1,j-1,k  ) &
                                  +uu(i-1,j  ,k-1) + uu(i  ,j-1,k-1) &
                            - FOUR*uu(i  ,j  ,k) )
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Lo-x face
      !
      i = 1
      !$OMP PARALLEL DO PRIVATE(j,k)
      do k = kstart,kend
      do j = jstart,jend
         dd(i,j,k) = ss(20,i,j,k)*(uu(i+1,j+1,k+1) + uu(i+1,j+1,k  ) &
                                  +uu(i+1,j  ,k+1) + uu(i  ,j+1,k+1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss( 8,i,j,k)*(uu(i+1,j+1,k-1) + uu(i+1,j+1,k  ) &
                                  +uu(i+1,j  ,k-1) + uu(i  ,j+1,k-1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss(15,i,j,k)*(uu(i+1,j-1,k+1) + uu(i+1,j-1,k  ) &
                                  +uu(i+1,j  ,k+1) + uu(i  ,j-1,k+1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss( 3,i,j,k)*(uu(i+1,j-1,k-1) + uu(i+1,j-1,k  ) &
                                  +uu(i+1,j  ,k-1) + uu(i  ,j-1,k-1) &
                            - FOUR*uu(i  ,j  ,k) )
         dd(i,j,k) = HALF*ff(i,j,k) - dd(i,j,k)
      end do
      end do
      !$OMP END PARALLEL DO
      !
      !     Hi-x face
      !
      i = nx+1
      !$OMP PARALLEL DO PRIVATE(j,k)
      do k = kstart,kend
      do j = jstart,jend
         dd(i,j,k) = ss(18,i,j,k)*(uu(i-1,j+1,k+1) + uu(i-1,j+1,k  ) &
                                  +uu(i-1,j  ,k+1) + uu(i  ,j+1,k+1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss(13,i,j,k)*(uu(i-1,j-1,k+1) + uu(i-1,j  ,k+1) &
                                  +uu(i-1,j-1,k  ) + uu(i  ,j-1,k+1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss( 6,i,j,k)*(uu(i-1,j+1,k-1) + uu(i-1,j  ,k-1) &
                                  +uu(i-1,j+1,k  ) + uu(i  ,j+1,k-1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss( 1,i,j,k)*(uu(i-1,j-1,k-1) + uu(i-1,j-1,k  ) &
                                  +uu(i-1,j  ,k-1) + uu(i  ,j-1,k-1) &
                          - FOUR*uu(i  ,j  ,k) )
         dd(i,j,k) = HALF*ff(i,j,k) - dd(i,j,k)
      end do
      end do
      !$OMP END PARALLEL DO
      !
      !     Lo-y face
      !
      j = 1
      !$OMP PARALLEL DO PRIVATE(i,k)
      do k = kstart,kend
      do i = istart,iend
         dd(i,j,k) = ss(20,i,j,k)*(uu(i+1,j+1,k+1) + uu(i+1,j+1,k  ) &
                                  +uu(i+1,j  ,k+1) + uu(i  ,j+1,k+1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss( 8,i,j,k)*(uu(i+1,j+1,k-1) + uu(i+1,j+1,k  ) &
                                  +uu(i+1,j  ,k-1) + uu(i  ,j+1,k-1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss(18,i,j,k)*(uu(i-1,j+1,k+1) + uu(i-1,j+1,k  ) &
                                  +uu(i-1,j  ,k+1) + uu(i  ,j+1,k+1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss( 6,i,j,k)*(uu(i-1,j+1,k-1) + uu(i-1,j+1,k  ) &
                                  +uu(i-1,j  ,k-1) + uu(i  ,j+1,k-1) &
                            - FOUR*uu(i  ,j  ,k) )
         dd(i,j,k) = HALF*ff(i,j,k) - dd(i,j,k)
      end do
      end do
      !$OMP END PARALLEL DO
      !
      !     Hi-y face
      !
      j = ny+1
      !$OMP PARALLEL DO PRIVATE(i,k)
      do k = kstart,kend
      do i = istart,iend
         dd(i,j,k) = ss( 3,i,j,k)*(uu(i+1,j-1,k-1) + uu(i  ,j-1,k-1) &
                                  +uu(i+1,j-1,k  ) + uu(i+1,j  ,k-1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss( 1,i,j,k)*(uu(i-1,j-1,k-1) + uu(i-1,j-1,k  ) &
                                  +uu(i-1,j  ,k-1) + uu(i  ,j-1,k-1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss(15,i,j,k)*(uu(i+1,j-1,k+1) + uu(i+1,j-1,k  ) &
                                  +uu(i  ,j-1,k+1) + uu(i+1,j  ,k+1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss(13,i,j,k)*(uu(i-1,j-1,k+1) + uu(i  ,j-1,k+1) &
                                  +uu(i-1,j-1,k  ) + uu(i-1,j  ,k+1) &
                            - FOUR*uu(i  ,j  ,k) )
         dd(i,j,k) = HALF*ff(i,j,k) - dd(i,j,k)
      end do
      end do
      !$OMP END PARALLEL DO
      !
      !     Lo-z face
      !
      k = 1
      !$OMP PARALLEL DO PRIVATE(i,j)
      do j = jstart,jend
      do i = istart,iend
         dd(i,j,k) = ss(15,i,j,k)*(uu(i+1,j-1,k+1) + uu(i+1,j-1,k  ) &
                                  +uu(i  ,j-1,k+1) + uu(i+1,j  ,k+1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss(13,i,j,k)*(uu(i-1,j-1,k+1) + uu(i  ,j-1,k+1) &
                                  +uu(i-1,j-1,k  ) + uu(i-1,j  ,k+1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss(20,i,j,k)*(uu(i+1,j+1,k+1) + uu(i+1,j+1,k  ) &
                                  +uu(i+1,j  ,k+1) + uu(i  ,j+1,k+1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss(18,i,j,k)*(uu(i-1,j+1,k+1) + uu(i-1,j+1,k  ) &
                                  +uu(i  ,j+1,k+1) + uu(i-1,j  ,k+1) &
                            - FOUR*uu(i  ,j  ,k) )
         dd(i,j,k) = HALF*ff(i,j,k) - dd(i,j,k)
      end do
      end do
      !$OMP END PARALLEL DO
      !
      !     Hi-z face
      !
      k = nz+1
      !$OMP PARALLEL DO PRIVATE(i,j)
      do j = jstart,jend
      do i = istart,iend
         dd(i,j,k) = ss( 3,i,j,k)*(uu(i+1,j-1,k-1) + uu(i  ,j-1,k-1) &
                                  +uu(i+1,j-1,k  ) + uu(i+1,j  ,k-1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss( 1,i,j,k)*(uu(i-1,j-1,k-1) + uu(i-1,j-1,k  ) &
                                  +uu(i-1,j  ,k-1) + uu(i  ,j-1,k-1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss( 8,i,j,k)*(uu(i+1,j+1,k-1) + uu(i+1,j  ,k-1) &
                                  +uu(i+1,j+1,k  ) + uu(i  ,j+1,k-1) &
                            - FOUR*uu(i  ,j  ,k) ) &
                    +ss( 6,i,j,k)*(uu(i-1,j+1,k-1) + uu(i-1,j+1,k  ) &
                                  +uu(i  ,j+1,k-1) + uu(i-1,j  ,k-1) &
                            - FOUR*uu(i  ,j  ,k) )
         dd(i,j,k) = HALF*ff(i,j,k) - dd(i,j,k)
      end do
      end do
      !$OMP END PARALLEL DO
      !
      !     Interior
      !
      !$OMP PARALLEL DO PRIVATE(i,j,k)
      do k = kstart,kend
      do j = jstart,jend
      do i = istart,iend

        dd(i,j,k) = ss(0,i,j,k)*uu(i,j,k) &
            + ss( 1,i,j,k) * uu(i-1,j-1,k-1) + ss( 2,i,j,k) * uu(i  ,j-1,k-1) &
            + ss( 3,i,j,k) * uu(i+1,j-1,k-1) + ss( 4,i,j,k) * uu(i-1,j  ,k-1) &
            + ss( 5,i,j,k) * uu(i+1,j  ,k-1) + ss( 6,i,j,k) * uu(i-1,j+1,k-1) &
            + ss( 7,i,j,k) * uu(i  ,j+1,k-1) + ss( 8,i,j,k) * uu(i+1,j+1,k-1) &
            + ss( 9,i,j,k) * uu(i-1,j-1,k  ) + ss(10,i,j,k) * uu(i+1,j-1,k  ) &
            + ss(11,i,j,k) * uu(i-1,j+1,k  ) + ss(12,i,j,k) * uu(i+1,j+1,k  ) &
            + ss(13,i,j,k) * uu(i-1,j-1,k+1) + ss(14,i,j,k) * uu(i  ,j-1,k+1) &
            + ss(15,i,j,k) * uu(i+1,j-1,k+1) + ss(16,i,j,k) * uu(i-1,j  ,k+1) &
            + ss(17,i,j,k) * uu(i+1,j  ,k+1) + ss(18,i,j,k) * uu(i-1,j+1,k+1) &
            + ss(19,i,j,k) * uu(i  ,j+1,k+1) + ss(20,i,j,k) * uu(i+1,j+1,k+1)

        if ((size(ss,dim=1) .eq. 27) .and. (.not. uniform_dh) ) then
            !
            ! Add faces (only non-zero for non-uniform dx)
            !
            dd(i,j,k) = dd(i,j,k) + &
                  ss(21,i,j,k) * uu(i-1,j  ,k  ) + ss(22,i,j,k) * uu(i+1,j  ,k  ) &
                  + ss(23,i,j,k) * uu(i  ,j-1,k  ) + ss(24,i,j,k) * uu(i  ,j+1,k  ) &
                  + ss(25,i,j,k) * uu(i  ,j  ,k-1) + ss(26,i,j,k) * uu(i  ,j  ,k+1)
         end if

         dd(i,j,k) = ff(i,j,k) - dd(i,j,k)
      end do
      end do
      end do
      !$OMP END PARALLEL DO

    else if (size(ss,dim=1) .eq. 7) then
       !
       !     Corners
       !
      i = 1
      j = 1
      k = 1
      dd(i,j,k) =  ss(1,i,j,k)*(uu(i+1,j,k)-uu(i,j,k)) &
                 + ss(3,i,j,k)*(uu(i,j+1,k)-uu(i,j,k)) &
                 + ss(5,i,j,k)*(uu(i,j,k+1)-uu(i,j,k))
      dd(i,j,k) = EIGHTH*ff(i,j,k) - dd(i,j,k)

      i = 1
      j = ny+1
      k = 1
      dd(i,j,k) =  ss(1,i,j,k)*(uu(i+1,j,k)-uu(i,j,k)) &
                 + ss(4,i,j,k)*(uu(i,j-1,k)-uu(i,j,k)) &
                 + ss(5,i,j,k)*(uu(i,j,k+1)-uu(i,j,k))
      dd(i,j,k) = EIGHTH*ff(i,j,k) - dd(i,j,k)

      i = nx+1
      j = 1
      k = 1
      dd(i,j,k) =  ss(2,i,j,k)*(uu(i-1,j,k)-uu(i,j,k)) &
                 + ss(3,i,j,k)*(uu(i,j+1,k)-uu(i,j,k)) &
                 + ss(5,i,j,k)*(uu(i,j,k+1)-uu(i,j,k))
      dd(i,j,k) = EIGHTH*ff(i,j,k) - dd(i,j,k)

      i = nx+1
      j = ny+1
      k = 1
      dd(i,j,k) =  ss(2,i,j,k)*(uu(i-1,j,k)-uu(i,j,k)) &
                 + ss(4,i,j,k)*(uu(i,j-1,k)-uu(i,j,k)) &
                 + ss(5,i,j,k)*(uu(i,j,k+1)-uu(i,j,k))
      dd(i,j,k) = EIGHTH*ff(i,j,k) - dd(i,j,k)

      i = 1
      j = 1
      k = nz+1
      dd(i,j,k) =  ss(1,i,j,k)*(uu(i+1,j,k)-uu(i,j,k)) &
                 + ss(3,i,j,k)*(uu(i,j+1,k)-uu(i,j,k)) &
                 + ss(6,i,j,k)*(uu(i,j,k-1)-uu(i,j,k))
      dd(i,j,k) = EIGHTH*ff(i,j,k) - dd(i,j,k)

      i = 1
      j = ny+1
      k = nz+1
      dd(i,j,k) =  ss(1,i,j,k)*(uu(i+1,j,k)-uu(i,j,k)) &
                 + ss(4,i,j,k)*(uu(i,j-1,k)-uu(i,j,k)) &
                 + ss(6,i,j,k)*(uu(i,j,k-1)-uu(i,j,k))
      dd(i,j,k) = EIGHTH*ff(i,j,k) - dd(i,j,k)

      i = nx+1
      j = 1
      k = nz+1
      dd(i,j,k) =  ss(2,i,j,k)*(uu(i-1,j,k)-uu(i,j,k)) &
                 + ss(3,i,j,k)*(uu(i,j+1,k)-uu(i,j,k)) &
                 + ss(6,i,j,k)*(uu(i,j,k-1)-uu(i,j,k))
      dd(i,j,k) = EIGHTH*ff(i,j,k) - dd(i,j,k)

      i = nx+1
      j = ny+1
      k = nz+1
      dd(i,j,k) =  ss(2,i,j,k)*(uu(i-1,j,k)-uu(i,j,k)) &
                 + ss(4,i,j,k)*(uu(i,j-1,k)-uu(i,j,k)) &
                 + ss(6,i,j,k)*(uu(i,j,k-1)-uu(i,j,k))
      dd(i,j,k) = EIGHTH*ff(i,j,k) - dd(i,j,k)
      !
      !     Lo-x / Lo-y edge
      !
      i = 1
      j = 1
      do k = kstart,kend
         dd(i,j,k) =      ss(1,i,j,k)*(uu(i+1,j,k)-uu(i,j,k)) &
                        + ss(3,i,j,k)*(uu(i,j+1,k)-uu(i,j,k)) &
                        + ss(5,i,j,k)*(uu(i,j,k+1)-uu(i,j,k)) &
                        + ss(6,i,j,k)*(uu(i,j,k-1)-uu(i,j,k))
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Hi-x / Lo-y edge
      !
      i = nx+1
      j = 1
      do k = kstart,kend
         dd(i,j,k) =      ss(2,i,j,k)*(uu(i-1,j,k)-uu(i,j,k)) &
                        + ss(3,i,j,k)*(uu(i,j+1,k)-uu(i,j,k)) &
                        + ss(5,i,j,k)*(uu(i,j,k+1)-uu(i,j,k)) &
                        + ss(6,i,j,k)*(uu(i,j,k-1)-uu(i,j,k))
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Lo-x / Hi-y edge
      !
      i = 1
      j = ny+1
      do k = kstart,kend
         dd(i,j,k) =      ss(1,i,j,k)*(uu(i+1,j,k)-uu(i,j,k)) &
                        + ss(4,i,j,k)*(uu(i,j-1,k)-uu(i,j,k)) &
                        + ss(5,i,j,k)*(uu(i,j,k+1)-uu(i,j,k)) &
                        + ss(6,i,j,k)*(uu(i,j,k-1)-uu(i,j,k))
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Hi-x / Hi-y edge
      !
      i = nx+1
      j = ny+1
      do k = kstart,kend
         dd(i,j,k) =      ss(2,i,j,k)*(uu(i-1,j,k)-uu(i,j,k)) &
                        + ss(4,i,j,k)*(uu(i,j-1,k)-uu(i,j,k)) &
                        + ss(5,i,j,k)*(uu(i,j,k+1)-uu(i,j,k)) &
                        + ss(6,i,j,k)*(uu(i,j,k-1)-uu(i,j,k))
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Lo-x / Lo-z edge
      !
      i = 1
      k = 1
      do j = jstart,jend
         dd(i,j,k) =      ss(1,i,j,k)*(uu(i+1,j,k)-uu(i,j,k)) &
                        + ss(3,i,j,k)*(uu(i,j+1,k)-uu(i,j,k)) &
                        + ss(4,i,j,k)*(uu(i,j-1,k)-uu(i,j,k)) &
                        + ss(5,i,j,k)*(uu(i,j,k+1)-uu(i,j,k))
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Hi-x / Lo-z edge
      !
      i = nx+1
      k = 1
      do j = jstart,jend
         dd(i,j,k) =      ss(2,i,j,k)*(uu(i-1,j,k)-uu(i,j,k)) &
                        + ss(3,i,j,k)*(uu(i,j+1,k)-uu(i,j,k)) &
                        + ss(4,i,j,k)*(uu(i,j-1,k)-uu(i,j,k)) &
                        + ss(5,i,j,k)*(uu(i,j,k+1)-uu(i,j,k))
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Lo-x / Hi-z edge
      !
      i = 1
      k = nz+1
      do j = jstart,jend
         dd(i,j,k) =      ss(1,i,j,k)*(uu(i+1,j,k)-uu(i,j,k)) &
                        + ss(3,i,j,k)*(uu(i,j+1,k)-uu(i,j,k)) &
                        + ss(4,i,j,k)*(uu(i,j-1,k)-uu(i,j,k)) &
                        + ss(6,i,j,k)*(uu(i,j,k-1)-uu(i,j,k))
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Hi-x / Hi-z edge
      !
      i = nx+1
      k = nz+1
      do j = jstart,jend
         dd(i,j,k) =      ss(2,i,j,k)*(uu(i-1,j,k)-uu(i,j,k)) &
                        + ss(3,i,j,k)*(uu(i,j+1,k)-uu(i,j,k)) &
                        + ss(4,i,j,k)*(uu(i,j-1,k)-uu(i,j,k)) &
                        + ss(6,i,j,k)*(uu(i,j,k-1)-uu(i,j,k))
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Lo-y / Lo-z edge
      !
      j = 1
      k = 1
      do i = istart,iend
         dd(i,j,k) =      ss(1,i,j,k)*(uu(i+1,j,k)-uu(i,j,k)) &
                        + ss(2,i,j,k)*(uu(i-1,j,k)-uu(i,j,k)) &
                        + ss(3,i,j,k)*(uu(i,j+1,k)-uu(i,j,k)) &
                        + ss(5,i,j,k)*(uu(i,j,k+1)-uu(i,j,k))
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Hi-y / Lo-z edge
      !
      j = ny+1
      k = 1
      do i = istart,iend
         dd(i,j,k) =      ss(1,i,j,k)*(uu(i+1,j,k)-uu(i,j,k)) &
                        + ss(2,i,j,k)*(uu(i-1,j,k)-uu(i,j,k)) &
                        + ss(4,i,j,k)*(uu(i,j-1,k)-uu(i,j,k)) &
                        + ss(5,i,j,k)*(uu(i,j,k+1)-uu(i,j,k))
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Lo-y / Hi-z edge
      !
      j = 1
      k = nz+1
      do i = istart,iend
         dd(i,j,k) =      ss(1,i,j,k)*(uu(i+1,j,k)-uu(i,j,k)) &
                        + ss(2,i,j,k)*(uu(i-1,j,k)-uu(i,j,k)) &
                        + ss(3,i,j,k)*(uu(i,j+1,k)-uu(i,j,k)) &
                        + ss(6,i,j,k)*(uu(i,j,k-1)-uu(i,j,k))
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Hi-y / Hi-z edge
      !
      j = ny+1
      k = nz+1
      do i = istart,iend
         dd(i,j,k) =      ss(1,i,j,k)*(uu(i+1,j,k)-uu(i,j,k)) &
                        + ss(2,i,j,k)*(uu(i-1,j,k)-uu(i,j,k)) &
                        + ss(4,i,j,k)*(uu(i,j-1,k)-uu(i,j,k)) &
                        + ss(6,i,j,k)*(uu(i,j,k-1)-uu(i,j,k))
         dd(i,j,k) = FOURTH*ff(i,j,k) - dd(i,j,k)
      end do
      !
      !     Lo-x face
      !
      i = 1
      !$OMP PARALLEL DO PRIVATE(j,k)
      do k = kstart,kend
      do j = jstart,jend
         dd(i,j,k) =      ss(1,i,j,k)*(uu(i+1,j,k)-uu(i,j,k)) &
                        + ss(3,i,j,k)*(uu(i,j+1,k)-uu(i,j,k)) &
                        + ss(4,i,j,k)*(uu(i,j-1,k)-uu(i,j,k)) &
                        + ss(5,i,j,k)*(uu(i,j,k+1)-uu(i,j,k)) &
                        + ss(6,i,j,k)*(uu(i,j,k-1)-uu(i,j,k))
         dd(i,j,k) = HALF*ff(i,j,k) - dd(i,j,k)
      end do
      end do
      !$OMP END PARALLEL DO
      !
      !     Hi-x face
      !
      i = nx+1
      !$OMP PARALLEL DO PRIVATE(j,k)
      do k = kstart,kend
      do j = jstart,jend
         dd(i,j,k) =      ss(2,i,j,k)*(uu(i-1,j,k)-uu(i,j,k)) &
                        + ss(3,i,j,k)*(uu(i,j+1,k)-uu(i,j,k)) &
                        + ss(4,i,j,k)*(uu(i,j-1,k)-uu(i,j,k)) &
                        + ss(5,i,j,k)*(uu(i,j,k+1)-uu(i,j,k)) &
                        + ss(6,i,j,k)*(uu(i,j,k-1)-uu(i,j,k))
         dd(i,j,k) = HALF*ff(i,j,k) - dd(i,j,k)
      end do
      end do
      !$OMP END PARALLEL DO
      !
      !     Lo-y face
      !
      j = 1
      !$OMP PARALLEL DO PRIVATE(i,k)
      do k = kstart,kend
      do i = istart,iend
         dd(i,j,k) =      ss(3,i,j,k)*(uu(i,j+1,k)-uu(i,j,k)) &
                        + ss(1,i,j,k)*(uu(i+1,j,k)-uu(i,j,k)) &
                        + ss(2,i,j,k)*(uu(i-1,j,k)-uu(i,j,k)) &
                        + ss(5,i,j,k)*(uu(i,j,k+1)-uu(i,j,k)) &
                        + ss(6,i,j,k)*(uu(i,j,k-1)-uu(i,j,k))
         dd(i,j,k) = HALF*ff(i,j,k) - dd(i,j,k)
      end do
      end do
      !$OMP END PARALLEL DO
      !
      !     Hi-y face
      !
      j = ny+1
      !$OMP PARALLEL DO PRIVATE(i,k)
      do k = kstart,kend
      do i = istart,iend
         dd(i,j,k) =      ss(4,i,j,k)*(uu(i,j-1,k)-uu(i,j,k)) &
                        + ss(1,i,j,k)*(uu(i+1,j,k)-uu(i,j,k)) &
                        + ss(2,i,j,k)*(uu(i-1,j,k)-uu(i,j,k)) &
                        + ss(5,i,j,k)*(uu(i,j,k+1)-uu(i,j,k)) &
                        + ss(6,i,j,k)*(uu(i,j,k-1)-uu(i,j,k))
         dd(i,j,k) = HALF*ff(i,j,k) - dd(i,j,k)
      end do
      end do
      !$OMP END PARALLEL DO
      !
      !     Lo-z face
      !
      k = 1
      !$OMP PARALLEL DO PRIVATE(i,j)
      do j = jstart,jend
      do i = istart,iend
         dd(i,j,k) =      ss(5,i,j,k)*(uu(i,j,k+1)-uu(i,j,k)) &
                        + ss(1,i,j,k)*(uu(i+1,j,k)-uu(i,j,k)) &
                        + ss(2,i,j,k)*(uu(i-1,j,k)-uu(i,j,k)) &
                        + ss(3,i,j,k)*(uu(i,j+1,k)-uu(i,j,k)) &
                        + ss(4,i,j,k)*(uu(i,j-1,k)-uu(i,j,k))
         dd(i,j,k) = HALF*ff(i,j,k) - dd(i,j,k)
      end do
      end do
      !$OMP END PARALLEL DO
      !
      !     Hi-z face
      !
      k = nz+1
      !$OMP PARALLEL DO PRIVATE(i,j)
      do j = jstart,jend
      do i = istart,iend
         dd(i,j,k) =      ss(6,i,j,k)*(uu(i,j,k-1)-uu(i,j,k)) &
                        + ss(1,i,j,k)*(uu(i+1,j,k)-uu(i,j,k)) &
                        + ss(2,i,j,k)*(uu(i-1,j,k)-uu(i,j,k)) &
                        + ss(3,i,j,k)*(uu(i,j+1,k)-uu(i,j,k)) &
                        + ss(4,i,j,k)*(uu(i,j-1,k)-uu(i,j,k))
         dd(i,j,k) = HALF*ff(i,j,k) - dd(i,j,k)
      end do
      end do
      !$OMP END PARALLEL DO
      !
      !     Interior
      !
      !$OMP PARALLEL DO PRIVATE(i,j,k)
      do k = kstart,kend
      do j = jstart,jend
      do i = istart,iend
          dd(i,j,k) = ss(0,i,j,k)*uu(i,j,k) &
            + ss(1,i,j,k) * uu(i+1,j  ,k  ) + ss(2,i,j,k) * uu(i-1,j  ,k  ) &
            + ss(3,i,j,k) * uu(i  ,j+1,k  ) + ss(4,i,j,k) * uu(i  ,j-1,k  ) &
            + ss(5,i,j,k) * uu(i  ,j  ,k+1) + ss(6,i,j,k) * uu(i  ,j  ,k-1)
          dd(i,j,k) = ff(i,j,k) - dd(i,j,k)
      end do
      end do
      end do
      !$OMP END PARALLEL DO

    end if

  end subroutine grid_laplace_3d

end module ml_nd_module

max_stars_repo_name: memmett/BoxLib

max_issues_repo_name: memmett/BoxLib

max_forks_repo_name: memmett/BoxLib


text:
!Do Metropolis-Hastings and slice sampling algorithms
!Also directional gridding method for fast-slow parameters
!Provisional implementation of Wang-Landau and Multicanonical methods
!August 2006

module MonteCarlo
 use ParamDef
 use CalcLike
 use Random
 use propose
 use Readmodel

 implicit none

 integer :: indep_sample = 0
  !number of iterations between dumping full model data. If zero then skip.

 integer :: directional_grid_steps = 20
  !for sampling_method = sampling_slowgrid, number of steps per grid

 real, parameter    :: eps_1 = 1.00001
 real :: MaxLike

!Multicanonical/W-L sampling parameters
  real, parameter :: mc_logspace =  6. !1/ steps to use in log-likelihood * number of parameter
  integer :: mc_mini = 0, mc_maxi = 0
  integer :: mc_update_steps = 7000, mc_update_burn = 50
  integer :: mc_steps_inc = 1000
  real, dimension(:,:), allocatable :: mc_like_counts
  real, dimension(:), allocatable :: mc_lnweights
!Wang-Landau parameter
  real :: WL_f = 1, WL_min = 0.1, WL_minW
  real :: WL_maxL=8. !WL_maxL is log like from best to use, in units of number of parameters
  integer :: WL_update_steps = 4000
  real :: WL_flat_tol = 1.4 !factor by which smallest histogram can be smaller than mean
contains

 function UpdateParamsLike(Params, fast, dist, i, freeNew, Lik) result(NewLik)
  !For slice sampling movement and likelihood
  Type(ParamSet) tmp, Params
  integer, intent(in) :: i
  real, intent(in), optional :: Lik
  real NewLik
  logical, intent(in) :: fast, freeNew
  real, intent(in) :: dist

  tmp = Params

  call UpdateParamsDirection(tmp, fast, dist, i)

  NewLik = GetLogLike(tmp)
  if (freeNew) then
!    call AcceptReject(.true., tmp%Info, Params%Info)
  else
     if (NewLik <= Lik*eps_1) then
       !accept move
!       call AcceptReject(.false., tmp%Info, Params%Info)
       Params = tmp
     else
!       call AcceptReject(.true., tmp%Info, Params%Info)
     endif
  endif

 end function UpdateParamsLike

 subroutine SliceUpdate(Params, fast, i, Like)
 !Do slice sampling, see http://www.cs.toronto.edu/~radford/slice-aos.abstract.html
 !update parameter i, input Like updated to new value.
 !Use linear stepping for the moment
  Type(ParamSet) Params
  integer, intent(in) :: i
  logical, intent(in) :: fast
  real, intent(inout) :: Like
  real offset,  P
  real L, R
  real LikL, LikR, Range, LikT, w
  integer fevals

   w = propose_scale
   if (slice_randomsize) w =w * Gaussian1()
   Like = Like + randexp1()  !New vertical position (likelihood)
   fevals = 0

   offset = ranmar()

   L = -offset
   R = 1- offset

   !step out
   if (slice_stepout) then
       do
        LikL = UpdateParamsLike(Params, fast, w*L, i, .true.)
        fevals = fevals + 1
        if (LikL*eps_1 > Like) exit
        L = L  - 1
       end do

       do
        LikR = UpdateParamsLike(Params, fast, w*R, i,.true.)
        fevals = fevals + 1
        if (LikR*eps_1 > Like) exit
        R = R + 1
       end do
   end if

   !stepping in
   do
    Range =  R - L
    P  =  ranmar() * Range + L
    LikT = UpdateParamsLike(Params, fast, w*P, i, .false., Like)
    fevals = fevals + 1
    if (LikT < Like*eps_1) exit
    if (P > 0) then
        R = P
    else
        L = P
    end if
   end do
   Like = LikT
   if (.not. fast) slow_proposals = slow_proposals + fevals

 end subroutine SliceUpdate

 subroutine SliceSampleSlowParam(CurParams, CurLike)
  Type(ParamSet) CurParams
  real CurLike
  integer, save:: loopix = 0

   if (mod(loopix,num_slow)==0) then
        if (.not. allocated(Rot_slow)) allocate(Rot_slow(num_slow,num_slow))
        call RotMatrix(Rot_slow,num_slow)
        loopix = 0
   end if
   loopix = loopix + 1
   call SliceUpdate(CurParams, .false., loopix,CurLike)

 end subroutine SliceSampleSlowParam

 subroutine SliceSampleFastParams(CurParams, CurLike)
  Type(ParamSet) CurParams
  real CurLike
  integer j

   if (.not. allocated(Rot_fast)) allocate(Rot_fast(num_fast,num_fast))
   call RotMatrix(Rot_fast,num_fast)

   do j = 1, num_fast * max(1,oversample_fast)
      call SliceUpdate(CurParams, .true., mod(j-1,num_fast)+1,CurLike)
   end do

 end subroutine SliceSampleFastParams

 subroutine SampleSlowGrid(CurParams, CurLike, mult, MaxLike, num, num_accept)
  !sampling_method = method_slowgrid
  !Make proposals on a grid in each slow direction simultaneously with changes to
  !fast parameters. Allows fast paraters to 'adapt' to slow parameter values, thereby
  !increasing slow acceptance rate
  Type(ParamSet) Trial, CurParams
  real CurLike, MaxLike
  integer num, num_accept, mult
  real w
  integer r, r_min, r_max, last_r
  Type(ParamSet), dimension(:), allocatable :: grid
  integer i, accpt
  real Like
  integer, save:: loopix = 0

   if (mod(loopix,num_slow)==0) then
        if (.not. allocated(Rot_slow)) allocate(Rot_slow(num_slow,num_slow))
        call RotMatrix(Rot_slow,num_slow)
        loopix = 0
   end if
   loopix = loopix + 1

   w = max(1e-5,propose_r(num_slow)) * propose_scale / 2

   r = 0
   r_min = 0
   r_max = 0
   accpt = 0

   allocate(grid(-directional_grid_steps:directional_grid_steps))
   grid(0) = CurParams

   do i = 1, directional_grid_steps*oversample_fast

    last_r = r

    if (ranmar() < 0.5) then
     r =r  + 1
    else
     r =r  - 1
    end if

    Trial = Grid(last_r)
    call UpdateParamsDirection(Trial, .false., (r-last_r)*w, loopix)

    if (r > r_max  .or. r < r_min) then
     grid(r) = Trial
     slow_proposals = slow_proposals + 1
     r_min = min(r,r_min)
     r_max = max(r,r_max)
    else
     grid(r)%P(fast_params_used) = Trial%P(fast_params_used)
    end if

    call GetProposalProjFast(grid(r), grid(r), propose_scale)

    Like = GetLogLike(grid(r))

    if (Feedback > 1) write (*,*) r, 'Likelihood: ', Like, 'Current Like:', CurLike

    if ((Like /= logZero) .and. (CurLike > Like .or. randexp1() > Like - CurLike)) then
    !Accept

      if (num_accept> burn_in) then
          output_lines = output_lines +1
         call WriteParams(grid(last_r),real(mult),CurLike)
      end if

      CurLike = Like
      accpt = accpt + 1
      if (Like < MaxLike) MaxLike = Like
      mult=1

    else
    !Reject
      r = last_r
      mult=mult+1
    end if

    call AddMPIParams(grid(r)%P,CurLike)

   num = num + 1

   end do

   if (Feedback > 1) write(*,*) 'Done slow grid, direction', loopix
   if (Feedback > 1) write(*,*) 'grid steps moved:',abs(r), ' acc rate =',&
            real(accpt)/(directional_grid_steps*oversample_fast)

   !Have to be careful freeing up here.. (logZero results don't allocate new)
   do last_r = r_max, r+1, -1
!      call AcceptReject(.true., grid(last_r)%Info,grid(last_r-1)%Info)
   end do
   do last_r = r_min, r-1
!      call AcceptReject(.true., grid(last_r)%Info,grid(last_r+1)%Info)
   end do

   CurParams = grid(r)
   num_accept = num_accept + accpt

   deallocate(grid)

 end subroutine SampleSlowGrid

!Multicanonical/Wang-Landau

 subroutine MC_AddLike(L)
  real, intent(in) :: L
  integer like_ix
  integer, dimension(:,:), allocatable :: tmp
  integer sc

   sc =0

   like_ix = int(L*mc_logspace/num_params_used)
   if (.not. allocated(mc_like_counts)) then
     mc_mini=like_ix-100
     mc_maxi=like_ix+100
     allocate(mc_like_counts(mc_mini:mc_maxi,2))
     mc_like_counts=sc
   end if
   if (like_ix > mc_maxi) then
     allocate(tmp(mc_mini:like_ix+100,2))
     tmp=sc
     tmp(mc_mini:mc_maxi,:) = mc_like_counts
     deallocate(mc_like_counts)
     mc_maxi = like_ix + 100
     allocate(mc_like_counts(mc_mini:mc_maxi,2))
     mc_like_counts = tmp
     deallocate(tmp)
   end if
   if (like_ix < mc_mini) then
     allocate(tmp(mc_mini-100:mc_maxi,2))
     tmp=sc
     tmp(mc_mini:mc_maxi,:) = mc_like_counts
     deallocate(mc_like_counts)
     mc_mini = like_ix - 100
     allocate(mc_like_counts(mc_mini:mc_maxi,2))
     mc_like_counts = tmp
     deallocate(tmp)
   end if

   mc_like_counts(like_ix,1)=mc_like_counts(like_ix,1)+ WL_f
   mc_like_counts(like_ix,2)=mc_like_counts(like_ix,2)+ 1

 end subroutine MC_AddLike

 subroutine MC_UpdateWeights
    real, dimension(:), allocatable :: tmp
    real win(-2:2)
    integer i

    allocate(tmp(mc_mini:mc_maxi))
    tmp=0
    tmp(LBOUND(mc_lnweights,DIM=1):UBOUND(mc_lnweights,DIM=1)) = mc_lnweights
    if (allocated(mc_lnweights)) deallocate(mc_lnweights)
    allocate(mc_lnweights(mc_mini:mc_maxi))
    mc_lnweights = tmp

    tmp = 0
    do i=-2, 2
     win(i) = exp(-real(i**2)*2)
    end do
    do i = mc_mini+2 , mc_maxi-2
       tmp(i) = sum(mc_like_counts(i-2:i+2,1)*win)
    end do
    !tmp=mc_like_counts
    tmp = max(tmp,1.)
    mc_lnweights = mc_lnweights + log(tmp/real(maxval(tmp)))
    deallocate(tmp)
    mc_like_counts = 0

 end  subroutine MC_UpdateWeights

 function MC_LnWeight(L)
  real, intent(in) :: L
  real MC_LnWeight
  integer like_ix

  if (.not. allocated(mc_lnweights)) then
       MC_LnWeight=0
  else
      like_ix = int(L*mc_logspace/num_params_used)
      if (like_ix > UBound(mc_lnweights,dim=1)) then
       MC_LnWeight = mc_lnweights(UBound(mc_lnweights,dim=1))
      elseif (like_ix < LBound(mc_lnweights,dim=1)) then
       MC_LnWeight= mc_lnweights(LBound(mc_lnweights,dim=1))
      else
       MC_LnWeight = mc_lnweights(like_ix)
      end if
  endif

 end  function MC_LnWeight

 function MC_WeightLike(L) result (Like)
  real, intent(in) :: L
  real Like

   Like = L + MC_LnWeight(L)

 end function MC_WeightLike

 function MC_Weight(L) result(W)
  real, intent(in) :: L
  real mx, W
  integer like_ix

  if (.not. allocated(mc_lnweights)) then
   W=1
  else
   like_ix =int(L*mc_logspace/num_params_used)

   mx = maxval(mc_lnweights)
   W =  exp((MC_LnWeight(L) - mx))
  end if

 end function MC_Weight

 subroutine WL_UpdateWeights
  integer i,n
  real amin,asum,L

   if (WL_F > WL_min) then

     amin = LogZero
     asum = 0
     n=0
     do i=mc_mini,mc_maxi
      L = i*num_params_used/mc_logspace
      if (mc_like_counts(i,1) /= 0 .and. L <= MaxLike + WL_maxL*num_params_used) then
           n=n+1
           amin = min(amin,mc_like_counts(i,2))
           asum = asum + mc_like_counts(i,2)
      end if
     end do
    WL_update_steps = WL_update_steps + 30*num_params_used

    if (amin > asum / (WL_flat_tol*n)) then
          WL_f = WL_f*0.5
          mc_like_counts(:,2)=0
          WL_update_steps = WL_update_steps + 300*num_params_used

          if (WL_F < WL_min) then
             if (Feedback > 0) write(*,*) 'Outputting Markovian samples'
             WL_f = 0
             WL_minW=LogZero
!            call CreateTxtFile('density.txt',1)
             do i=mc_mini,mc_maxi
              L = i*num_params_used/mc_logspace
              if (mc_like_counts(i,1) /= 0 .and. L <= MaxLike + WL_maxL*num_params_used) then
                     WL_minW = min(WL_minW,L - mc_like_counts(i,1))
  !                   write (1,*) L,mc_like_counts(i,1)
              end if
             end do
 !            close(1)
          end if
    end if
   end if
 end  subroutine WL_UpdateWeights

 function WL_WeightLike(L) result (Like)
  real, intent(in) :: L
  real Like
  integer like_ix

  if (L > MaxLike + WL_maxL*num_params_used) then
    Like = LogZero  !!!!
    return
  end if

  like_ix = int(L*mc_logspace/num_params_used)
  if (like_ix >= mc_maxi .or. like_ix <= mc_mini) then
   Like = 0
  else
   Like = mc_like_counts(like_ix,1)
  end if

 end function WL_WeightLike

function WL_Weight(L) result (W)
  real, intent(in) :: L
  real W
  integer like_ix

  if (L > MaxLike + WL_maxL*num_params_used) then
   W=0
   return
  end if

  like_ix = int(L*mc_logspace/num_params_used)
  if (like_ix >= mc_maxi .or. like_ix <= mc_mini) then
   W = 0
  else
   W  = exp( -L + mc_like_counts(like_ix,1) + WL_minW)
  end if

 end function WL_Weight

 subroutine MCMCsample(Params, samples_to_get)

   use Readmodel

   integer samples_to_get
   Type(ParamSet) Trial, Params, CurParams
   real Like, CurLike
   logical accpt
   real rmult
   integer mult,num_accept, last_num
   integer num_metropolis,  num_metropolis_accept
   real testlike, testCurLike
   integer mc_update_burn_num, mc_updates

   MaxLike = LogZero
   CurLike = StartLike
   testCurLike = StartLike
   CurParams = Params
   num_accept = 0
   last_num = 0
   mult= 1
   num_metropolis = 0
   num_metropolis_accept = 0
   mc_update_burn_num = mc_update_burn
   mc_updates = 0

   call read_powlg_obsdata()

   do while (num <= samples_to_get)
      if (mod(num,1000)==1) write(*,*) num,samples_to_get

   num = num + 1

   if (indep_sample /= 0) then
      if (num_accept>= burn_in .and. mod(num, indep_sample)==0) then
       if (CurLike /= LogZero) call WriteIndepSample(CurParams, CurLike)
      end if
   end if

   if (CurLike /= LogZero .and. sampling_method == sampling_slice) then
    !Slice sampling
          if (num_accept> burn_in) then
           output_lines = output_lines +1
           call WriteParams(CurParams, real(mult), CurLike)
          end if
          if (Feedback > 1) write (*,*) instance, 'Slicing, Current Like:', CurLike
          mult = 1
          if (num_slow /=0) call SliceSampleSlowParam(CurParams, CurLike)
          if (num_fast /=0) call SliceSampleFastParams(CurParams, CurLike)
          num_accept = num_accept + 1
          if (CurLike < MaxLike) MaxLike = CurLike

  elseif (sampling_method == sampling_slowgrid .and. CurLike /= LogZero .and. num_fast /= 0 .and. num_slow /=0) then

          if (Feedback > 1) write (*,*) instance, 'Directional gridding, Like: ', CurLike
          call SampleSlowGrid(CurParams, CurLike,mult, MaxLike, num, num_accept)

  else

   !Do metropolis, except for (optional) slice sampling on fast parameters
   !(slicing dynamically adjusts for the marginal distribution width at the
   ! expense of additional function evaluations)
    if (sampling_method == sampling_fastslice) then
        if (CurLike /=LogZero .and. num_fast /=0) then
            if (num_accept> burn_in) then
              output_lines = output_lines +1
              call WriteParams(CurParams,real(mult), CurLike)
            end if
            call SliceSampleFastParams(CurParams, CurLike)
           if (CurLike < MaxLike) MaxLike = CurLike
            mult = 1
        end if
        call GetProposalProjSlow(CurParams, Trial)
    else
        call GetProposal(CurParams, Trial)
    end if

    Like = GetLogLike(Trial)

    if (Feedback > 1) write (*,*) 'Likelihood: ', Like, 'Current Like:', CurLike

!!!    accpt = (Like /= logZero) .and. (CurLike > Like .or. randexp1() > Like - CurLike)
             !Include the min() so that compilers not doing optimal compilation don't complain

    if (Like /= logZero) then
       if (sampling_method == sampling_multicanonical) then
        testLike = MC_WeightLike(Like)
       elseif (sampling_method == sampling_wang_landau) then
        testLike = WL_WeightLike(Like)
       else
        testLike = Like
       end if
       accpt = (testLike /= LogZero) .and. (testCurLike > testLike .or. randexp1() > testLike - testCurLike)
       if (.not. accpt .and. sampling_method == sampling_wang_landau) testCurLike = testCurLike + WL_f
    else
      accpt = .false.
    end if

!   call AcceptReject(accpt, CurParams%Info, Trial%Info)
   num_metropolis = num_metropolis +1

   if (accpt) then

      if (num_accept> burn_in .and. CurLike/=LogZero ) then
         if (sampling_method == sampling_multicanonical) then
          if (num>mc_update_burn_num) then
           rmult = mult*MC_Weight(CurLike)
           else
           rmult = 0.
           end if
         else
          rmult = mult
         end if
         if (sampling_method == sampling_wang_landau) then
           if (WL_f >= WL_min) then
             rmult = 0.
           else
             rmult = mult*WL_Weight(CurLike)
           end if
         end if

         if (rmult /= 0.) then
          output_lines = output_lines +1
          call WriteParams(CurParams,rmult,CurLike)
         end if
      end if

      CurLike = Like
      testCurLike = testLike
      CurParams = Trial
      num_accept = num_accept + 1
      num_metropolis_accept = num_metropolis_accept + 1
      if (Like < MaxLike) MaxLike = Like
      mult=1

      if (Feedback > 1) write (*,*) num_metropolis, ' accepting. ratio:', real(num_metropolis_accept)/num_metropolis

      if (mod(num_metropolis_accept,50) ==0) then
          write (*,*) MPIRank, 'rat:',real(num_metropolis_accept)/num_metropolis, ' in ',num_metropolis, &
               ' (M) best: ',real(MaxLike)
      end  if
      if (logfile_unit /=0 .and. mod(num_metropolis_accept,50) ==0) then
           write (logfile_unit,*) 'rat:',real(num_metropolis_accept)/num_metropolis, ' in ',num_metropolis, &
               ', (M) best: ',real(MaxLike)
           write (logfile_unit,*) 'local acceptance ratio:', 50./(num_metropolis - last_num)

           call FlushFile(logfile_unit)
           last_num = num_metropolis
      end if
    else
       mult=mult+1
       if (Feedback > 1) write (*,*) num_metropolis,' rejecting. ratio:', real(num_metropolis_accept)/num_metropolis
    end if

    end if !not slicing

    if (CurLike /= logZero) then
      if (sampling_method /= sampling_slowgrid) call AddMPIParams(CurParams%P,CurLike)
      if (mod(num,100)==0) call CheckParamChange

      if (sampling_method == sampling_multicanonical .and. &
              mod(num,mc_update_steps + mc_updates*mc_steps_inc)==0) then
        call MC_UpdateWeights
        Like=LogZero
        testCurLike = StartLike
        mc_update_burn_num = num + mc_update_burn
        mc_updates = mc_updates + 1
      end if

      if (sampling_method == sampling_wang_landau .and. &
              num > WL_update_steps) then
           call WL_UpdateWeights
      end if

      if (sampling_method == sampling_multicanonical .and. num>mc_update_burn_num .or. &
           sampling_method == sampling_wang_landau) then
         call MC_AddLike(CurLike)
      end if

     else
      if (num > 1000) then
        call DoAbort('MCMC.f90: Couldn''t start after 1000 tries - check starting ranges')
      end if
    end if

   end do

   if (Feedback > 0) then
     write(*,*) 'Stopping as have ',samples_to_get ,' samples. '
     write(*,*) MPIRank, ' Slow proposals: ',slow_proposals
   end if

 end subroutine MCMCsample

  subroutine WriteParams(P, mult, like)
     use settings
     use ParamDef
     implicit none
    Type(ParamSet) P
    real, intent(in) :: mult, like
    character(LEN =60) fmt
!    Type(CMBParams) C
    real r10

    if (outfile_unit ==0) return

    fmt = trim(numcat('(2E16.7,',num_params))//'E16.7)'
    write (outfile_unit,fmt) mult,like, P%P

    if (flush_write) call FlushFile(outfile_unit)

  end  subroutine WriteParams

end module MonteCarlo

max_stars_repo_name: chiaki-hikage/glenclustfit

max_issues_repo_name: chiaki-hikage/glenclustfit

max_forks_repo_name: chiaki-hikage/glenclustfit


text:
      SUBROUTINE zrinfo (IFLTAB, CPATH, LFOUND, IDTYPE, CDTYPE, LDOUB,
     *                   LQUAL, IPRECIS, CRTAG, CLWDATE, CLWTIME,
     *                   CPNAME, IVERS, NDATA, NSPACE, ICOMPRES, LPASS)

C
      implicit none
C
C     Get record information
C     Written by Bill Charley at HEC, 2006.
C
C     Input:   IFLTAB  - DSS file table
C              CPATH   - Pathname
C     Output:  LFOUND  - Logical, true if exists
C              IDTYPE  - Integer, record type
C              CDTYPE  - Character, record type description (50 chr)
C              LDOUB   - Logical, double precision
C              LQUAL   - Logical, quality flags stored
C              IPRECIS - Integer, decimal precision of data
C              CRTAG   - Character, pathname tag (8 chr)
C              CLWDATE - Character, last written date (8 chr)
C              CLWTIME - Character, last written time (10 chr)
C              CPNAME  - Character, program name that wrote data (8)
C              IVERS   - Integer, version of data (number time written)
C              NDATA   - Integer, number of data
C              NSPACE  - Integer, space allocated for data
C              ICOMPRES- Integer, compress method (TS only, 0 for none)
C              LPASS   - Logical, password applied to this record
C
C
C
      INTEGER IFLTAB(*), IDTYPE, IPRECIS, IVERS, NDATA, NSPACE, ICOMPRES
      integer zdssVersion
      CHARACTER CPATH*(*)
      CHARACTER CDTYPE*(*), CRTAG*(*), CLWDATE*(*), CLWTIME*(*)
      CHARACTER CPNAME*(*)
      LOGICAL LFOUND, LDOUB, LQUAL, LPASS
      CHARACTER PATHNAME*393
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     FIX ME - MOVE PATHCHECK TO FUNCTION AND CALL ONLY
C     IF RECORD NOT FOUND!!!!
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     Adjust the time interval for the DSS version, if necessary
      pathname = cpath
      call ztsPathCheckInterval(ifltab, PATHNAME)
C
C
      IF (zdssVersion(IFLTAB).EQ.6) THEN
          CALL zrinfo6 (IFLTAB, PATHNAME, LFOUND, IDTYPE, CDTYPE, LDOUB,
     *                  LQUAL, IPRECIS, CRTAG, CLWDATE, CLWTIME,
     *                  CPNAME, IVERS, NDATA, NSPACE, ICOMPRES, LPASS)
      ELSE
          CALL zrinfo7 (IFLTAB, PATHNAME, LFOUND, IDTYPE, CDTYPE, LDOUB,
     *                  LQUAL, IPRECIS, CRTAG, CLWDATE, CLWTIME,
     *                  CPNAME, IVERS, NDATA, NSPACE, ICOMPRES, LPASS)
      ENDIF
C
      RETURN
      END

ext: f

max_stars_repo_name: HydrologicEngineeringCenter/heclib

max_issues_repo_name: HydrologicEngineeringCenter/heclib

max_forks_repo_name: HydrologicEngineeringCenter/heclib


text:
c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine  adi(rho_i, us, vs, ws, qs, square, rhs,
     $                forcing, u, nx, nxmax, ny, nz)

      integer nx, nxmax, ny, nz
      double precision rho_i  (  0:nxmax-1,0:ny-1,0:nz-1),
     $                 us     (  0:nxmax-1,0:ny-1,0:nz-1),
     $                 vs     (  0:nxmax-1,0:ny-1,0:nz-1),
     $                 ws     (  0:nxmax-1,0:ny-1,0:nz-1),
     $                 qs     (  0:nxmax-1,0:ny-1,0:nz-1),
     $                 square (  0:nxmax-1,0:ny-1,0:nz-1),
     $                 rhs    (5,0:nxmax-1,0:ny-1,0:nz-1),
     $                 forcing(5,0:nxmax-1,0:ny-1,0:nz-1),
     $                 u      (5,0:nxmax-1,0:ny-1,0:nz-1)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      call compute_rhs(rho_i, us, vs, ws, qs, square, rhs,
     $                 forcing, u, nx, nxmax, ny, nz)

      call x_solve(rho_i, qs, square, u, rhs, nx, nxmax, ny, nz)

      call y_solve(rho_i, qs, square, u, rhs, nx, nxmax, ny, nz)

      call z_solve(rho_i, qs, square, u, rhs, nx, nxmax, ny, nz)

      call add(u, rhs, nx, nxmax, ny, nz)

      return
      end

max_stars_repo_name: robert-mijakovic/readex-ptf

max_issues_repo_name: robert-mijakovic/readex-ptf

max_forks_repo_name: robert-mijakovic/readex-ptf


text:
! =============================================================================
!               Low-storage 4th order Runge-Kutta method
!            (see https://doi.org/10.5194/gmd-10-3145-2017)
! =============================================================================
module ls_rk4
    use options, only : parcel, time
    use parcel_container
    use parcel_bc
    use rk4_utils, only: get_dBdt, get_time_step
    use utils, only : write_step
    use parcel_interpl, only : par2grid, grid2par, grid2par_add
    use fields, only : velgradg, velog, vortg, vtend, tbuoyg
    use inversion_mod, only : vor2vel, vorticity_tendency
    use parcel_diagnostics, only : calculate_parcel_diagnostics
    use field_diagnostics, only : calculate_field_diagnostics
    use parameters, only : nx, nz
    use timer, only : start_timer, stop_timer, timings
    implicit none

    integer, parameter :: dp=kind(zero)           ! double precision

    integer :: rk4_timer

    double precision, allocatable, dimension(:, :) :: &
        delta_pos, &
        strain,    &   ! strain at parcel location
        delta_b,   &   ! B matrix integration
        delta_vor      ! vorticity integration

    double precision, parameter :: &
        ca1 = - 567301805773.0_dp/1357537059087.0_dp,  &
        ca2 = -2404267990393.0_dp/2016746695238.0_dp,  &
        ca3 = -3550918686646.0_dp/2091501179385.0_dp,  &
        ca4 = -1275806237668.0_dp/842570457699.0_dp,   &
        ca5 = 0.0,   &  !dummy value, not actually used
        cb1 =  1432997174477.0_dp/9575080441755.0_dp,  &
        cb2 =  5161836677717.0_dp/13612068292357.0_dp, &
        cb3 =  1720146321549.0_dp/2090206949498.0_dp,  &
        cb4 =  3134564353537.0_dp/4481467310338.0_dp,  &
        cb5 =  2277821191437.0_dp/14882151754819.0_dp

    contains

        ! Allocate memory of low-storage RK-4 temporaries.
        ! This memory is deallocated at the end of the simulation.
        ! @param[in] num is the size to allocate
        subroutine ls_rk4_alloc(num)
            integer, intent(in) :: num

            allocate(delta_pos(3, num))
            allocate(delta_vor(3, num))
            allocate(strain(5, num))
            allocate(delta_b(5, num))

        end subroutine ls_rk4_alloc

        ! Deallocate memory of temporaries
        subroutine ls_rk4_dealloc

            deallocate(delta_pos)
            deallocate(delta_vor)
            deallocate(strain)
            deallocate(delta_b)

        end subroutine ls_rk4_dealloc

        ! Advances the parcels by a single ls-RK-4 step. It calls a
        ! function to obtain the current time step based on the velocity
        ! strain and the buoyancy gradient.
        ! @param[in] t is the time
        ! Precondition: this routine assumes that the fields are
        ! up-to-date for the first sub-step
        subroutine ls_rk4_step(t)
            double precision, intent(inout) :: t
            double precision                :: dt

            call par2grid((t > time%initial))

            ! need to be called in order to set initial time step;
            ! this is also needed for the first ls-rk4 substep
            call vor2vel(vortg, velog, velgradg)

            call vorticity_tendency(vortg, velog, tbuoyg, vtend)

            ! update the time step
            dt = get_time_step(t)

            call grid2par(delta_pos, delta_vor, strain)

            call calculate_parcel_diagnostics(delta_pos)

            call calculate_field_diagnostics

            call write_step(t)

            call ls_rk4_substep(ca1, cb1, dt, 1)

            call par2grid
            call ls_rk4_substep(ca2, cb2, dt, 2)

            call par2grid
            call ls_rk4_substep(ca3, cb3, dt, 3)

            call par2grid
            call ls_rk4_substep(ca4, cb4, dt, 4)

            call par2grid
            call ls_rk4_substep(ca5, cb5, dt, 5)

            call start_timer(rk4_timer)
            call apply_parcel_bc(parcels%position, parcels%B)
            call stop_timer(rk4_timer)

            ! we need to subtract 14 calls since we start and stop
            ! the timer multiple times which increments n_calls
            timings(rk4_timer)%n_calls =  timings(rk4_timer)%n_calls - 14

            t = t + dt
        end subroutine ls_rk4_step

        ! Do a ls-RK-4 substep.
        ! @param[in] ca is a numerical coefficient
        ! @param[in] cb is a numerical coefficient
        ! @param[in] dt is the time step
        ! @param[in] step is the number of the substep (1 to 5)
        subroutine ls_rk4_substep(ca, cb, dt, step)
            double precision, intent(in) :: ca
            double precision, intent(in) :: cb
            double precision, intent(in) :: dt
            integer,          intent(in) :: step
            integer                      :: n

            if (step == 1) then
                call start_timer(rk4_timer)

                !$omp parallel do default(shared) private(n)
                do n = 1, n_parcels
                    delta_b(:, n) = get_dBdt(parcels%B(:, n), strain(:, n), &
                                             parcels%vorticity(:, n), parcels%volume(n))
                enddo
                !$omp end parallel do

                call stop_timer(rk4_timer)
            else
                call vor2vel(vortg, velog, velgradg)

                call vorticity_tendency(vortg, velog, tbuoyg, vtend)

                call grid2par_add(delta_pos, delta_vor, strain)

                call start_timer(rk4_timer)

                !$omp parallel do default(shared) private(n)
                do n = 1, n_parcels
                    delta_b(:, n) = delta_b(:, n) &
                                  + get_dBdt(parcels%B(:, n), strain(:, n), &
                                             parcels%vorticity(:, n), parcels%volume(n))
                enddo
                !$omp end parallel do

                call stop_timer(rk4_timer)
            endif

            call start_timer(rk4_timer)

            !$omp parallel do default(shared) private(n)
            do n = 1, n_parcels
                parcels%position(:, n) = parcels%position(:, n) &
                                       + cb * dt * delta_pos(:, n)

                parcels%vorticity(:, n) = parcels%vorticity(:, n) + cb * dt * delta_vor(:, n)
                parcels%B(:, n) = parcels%B(:, n) + cb * dt * delta_b(:, n)
            enddo
            !$omp end parallel do

            call stop_timer(rk4_timer)

            if (step == 5) then
               return
            endif

            call start_timer(rk4_timer)

            !$omp parallel do default(shared) private(n)
            do n = 1, n_parcels
                delta_pos(:, n) = ca * delta_pos(:, n)
                delta_vor(:, n) = ca * delta_vor(:, n)
                delta_b(:, n) = ca * delta_b(:, n)
            enddo
            !$omp end parallel do

            call stop_timer(rk4_timer)

        end subroutine ls_rk4_substep

end module ls_rk4

max_stars_repo_name: matt-frey/epic

max_issues_repo_name: matt-frey/epic

max_forks_repo_name: matt-frey/epic


! { dg-do compile }
! { dg-options "-fdump-tree-original -fcoarray=lib" }
!
! Argument passing was wrong
!

program atomic
  use iso_fortran_env
  implicit none

  integer :: me
  integer(atomic_int_kind) :: atom[*]
  me = this_image()
  call atomic_define(atom[1],0)
  sync all
  call ATOMIC_ADD (atom[1], me)
  if(me == 1) call atomic_ref(me,atom[1])
  sync all
  write(*,*) me
end program

! { dg-final { scan-tree-dump-times "value.. = 0;" 1 "original" } }
! { dg-final { scan-tree-dump-times "_gfortran_caf_atomic_define \\(caf_token.0, 0, 1, &value.., 0B, 1, 4\\);" 1 "original" } }
! { dg-final { scan-tree-dump-times "_gfortran_caf_atomic_op \\(1, caf_token.0, 0, 1, &me, 0B, 0B, 1, 4\\);" 1 "original" } }
! { dg-final { scan-tree-dump-times "_gfortran_caf_atomic_ref \\(caf_token.0, 0, 1, &me, 0B, 1, 4\\);" 1 "original" } }