text:
!!  Routines that run the MGRIT algorithm
!
! This file is part of LIBPFASST.
!

!> Module of routines to run MGRIT
module pf_mod_MGRIT
  use pf_mod_interpolate
  use pf_mod_restrict
  use pf_mod_utils
  use pf_mod_timer
  use pf_mod_dtype
  use pf_mod_hooks
  use pf_mod_pfasst
  use pf_mod_comm
  use pf_mod_results
  use mpi

  implicit none

  type :: int_vector
     integer, allocatable :: val(:)
  end type int_vector

  type :: mgrit_level_data
     integer :: Nt
     integer :: Nt_glob
     real(pfdp) :: dt
     real(pfdp) :: T0_glob
     real(pfdp) :: Tfin_glob
     real(pfdp) :: t0
     real(pfdp) :: tfin
     real(pfdp) :: res_norm_loc(1)
     integer :: res_norm_index
     integer :: cycle_phase
     logical :: FCF_flag = .true.
     logical :: FAS_flag = .false.
     logical :: setup_start_coarse_flag = .true.
     type(int_vector), allocatable :: f_blocks(:)
     integer, allocatable :: c_pts(:)
     class(pf_encap_t), allocatable :: g(:)
     class(pf_encap_t), allocatable :: qc(:)
     class(pf_encap_t), allocatable :: qc_prev(:)
     class(pf_encap_t), allocatable :: qc_fas(:)
     class(pf_encap_t), allocatable :: r
     class(pf_encap_t), allocatable :: qc_boundary
     class(pf_encap_t), allocatable :: qc_fas_boundary
     class(pf_encap_t), allocatable :: q_temp
     class(pf_encap_t), allocatable :: Q0
     class(pf_encap_t), allocatable :: g_temp
     integer :: send_to_rank
     integer :: recv_from_rank
     integer :: rank_shifted
     integer :: coarsest_level
     logical :: f_pts_flag = .true.
     logical :: c_pts_flag = .true.
  end type

contains

  subroutine mgrit_initialize(pf, mg_ld, T0_glob, Tfin_glob, Nt_start, coarsen_factor, FAS_flag, FCF_flag, start_coarse_flag)
     type(pf_pfasst_t), intent(inout) :: pf
     type(mgrit_level_data), allocatable, target, intent(inout) :: mg_ld(:)
     integer, intent(in) :: Nt_start, coarsen_factor
     logical, intent(in) :: FAS_flag, FCF_flag, start_coarse_flag
     real(pfdp) :: T0_glob, Tfin_glob, x, y
     integer :: level_index, nlevels, N, i, kk, p, Nt, p_rank, k
     type(pf_level_t) :: pf_lev, pf_f_lev, pf_c_lev
     type(mgrit_level_data), pointer :: mg_lev, mg_f_lev, mg_c_lev
     class(pf_encap_t), allocatable :: qtemp
     integer :: count_levels
     logical :: start_dropping_procs
     integer :: my_rank, rank_send, rank_recv, Nt_temp
     integer :: Nt_coarsest, coarsest_level

     if (start_coarse_flag .eqv. .false.) then
        x = real(pf%comm%nproc, pfdp)
        y = log(x) / log(real(coarsen_factor, pfdp))
        if (ceiling(y) .ne. floor(y)) then
           print *,'Error in MGRIT setup: number of procs must be a power of the coarsening factor'
           stop
        end if
        x = real(Nt_start * pf%comm%nproc, pfdp)
        y = log(x) / log(real(coarsen_factor, pfdp))
        !print *,y,ceiling(y),floor(y)
        if (ceiling(y) .ne. floor(y)) then
           print *,'Error in MGRIT setup: number of finest-grid points must be a power of the coarsening factor'
           stop
        end if
     end if

     Nt_coarsest = max(4, coarsen_factor)

     nlevels = pf%nlevels
     n = coarsen_factor

     allocate(mg_ld(nlevels))
     mg_ld(nlevels)%setup_start_coarse_flag = start_coarse_flag

     do level_index = 1,nlevels
        mg_lev => mg_ld(level_index)
        mg_lev%res_norm_loc(1) = 0.0_pfdp
        mg_lev%T0_glob = T0_glob
        mg_lev%Tfin_glob = Tfin_glob
     end do
     if (mg_ld(nlevels)%setup_start_coarse_flag .eqv. .true.) then
        coarsest_level = 1
        mg_ld(nlevels)%coarsest_level = 1
        N = Nt_start;
        do level_index = 1,nlevels
           mg_lev => mg_ld(level_index)
           mg_lev%Nt = N
           mg_lev%Nt_glob = N * pf%comm%nproc
           mg_lev%dt = (Tfin_glob - T0_glob)/real(mg_lev%Nt_glob, pfdp)
           mg_lev%t0 = T0_glob + real(N, pfdp) * mg_lev%dt * real(pf%rank, pfdp) + mg_lev%dt
           mg_lev%tfin = mg_lev%t0 + N * mg_lev%dt - mg_lev%dt
           pf%state%t0 = mg_lev%t0
           N = N * coarsen_factor

           if (level_index .eq. nlevels) then
              pf%state%t0 = mg_lev%tfin - mg_lev%dt
              pf%state%dt = mg_lev%dt
              pf%state%step = (pf%rank+1) * mg_lev%Nt - 1
           end if

           mg_lev%rank_shifted = pf%rank
           mg_lev%send_to_rank = pf%rank+1
           mg_lev%recv_from_rank = pf%rank-1
        end do
     else
        Nt_temp = Nt_start * pf%comm%nproc
        level_index = nlevels
        mg_ld(level_index)%Nt_glob = Nt_temp
        do while (level_index .gt. 1)
           Nt_temp = Nt_temp / coarsen_factor
           if (Nt_temp < Nt_coarsest) then
              exit
           end if
           level_index = level_index - 1
           mg_ld(level_index)%Nt_glob = Nt_temp
        end do
        coarsest_level = level_index
        mg_ld(nlevels)%coarsest_level = coarsest_level
        mg_ld(nlevels)%Nt = Nt_start
        mg_ld(nlevels)%rank_shifted = pf%rank
        mg_ld(nlevels)%send_to_rank = pf%rank+1
        mg_ld(nlevels)%recv_from_rank = pf%rank-1
        start_dropping_procs = .false.
        i = 0
        do level_index = (nlevels-1),coarsest_level,-1
           mg_lev => mg_ld(level_index)
           mg_f_lev => mg_ld(level_index+1)
           mg_lev%send_to_rank = -1
           mg_lev%recv_from_rank = -1
           mg_lev%rank_shifted = -1
           p = pf%rank+1
           if (mg_f_lev%Nt .eq. 1) then
              if (start_dropping_procs .eqv. .false.) then
                 start_dropping_procs = .true.
              end if
              if (mod(p, coarsen_factor**(i+1)) .eq. 0) then
                 mg_lev%Nt = 1
                 mg_lev%rank_shifted = (mg_f_lev%rank_shifted+1)/ coarsen_factor - 1
                 mg_lev%send_to_rank = (mg_lev%rank_shifted+2) * coarsen_factor**(i+1) - 1
                 mg_lev%recv_from_rank = mg_lev%rank_shifted * coarsen_factor**(i+1) - 1
                 mg_f_lev%f_pts_flag = .false.
                 mg_f_lev%c_pts_flag = .true.
              else
                 mg_f_lev%f_pts_flag = .true.
                 mg_f_lev%c_pts_flag = .false.
                 mg_lev%c_pts_flag = .false.
                 mg_lev%f_pts_flag = .false.
                 mg_lev%Nt = 0
              end if
           else if (mg_f_lev%Nt .gt. 1) then
              mg_lev%Nt = mg_f_lev%Nt / coarsen_factor
              mg_lev%rank_shifted = pf%rank
              mg_lev%send_to_rank = pf%rank+1
              mg_lev%recv_from_rank = pf%rank-1
           else
              mg_lev%Nt = 0
              mg_lev%c_pts_flag = .false.
              mg_lev%f_pts_flag = .false.
           end if
           if (start_dropping_procs .eqv. .true.) then
              i = i + 1;
           end if
        end do

        do level_index = coarsest_level,nlevels
           mg_lev => mg_ld(level_index)
           mg_lev%dt = (Tfin_glob - T0_glob)/mg_lev%Nt_glob
           mg_lev%t0 = T0_glob + mg_lev%dt * pf%rank * mg_lev%Nt + mg_lev%dt
           mg_lev%tfin = mg_lev%t0 + mg_lev%Nt * mg_lev%dt
           pf%state%t0 = mg_lev%t0
        end do

        !do level_index = nlevels,coarsest_level,-1
        !   mg_lev => mg_ld(level_index)
        !   print *,level_index,pf%rank,mg_lev%send_to_rank,mg_lev%recv_from_rank
        !end do
     end if

     do level_index = nlevels,coarsest_level+1,-1
        mg_lev => mg_ld(level_index)
        mg_c_lev => mg_ld(level_index-1)

        if (mg_ld(nlevels)%setup_start_coarse_flag .eqv. .true.) then
           call FC_Setup(mg_c_lev%Nt, mg_c_lev%Nt, coarsen_factor, mg_lev%f_blocks, mg_lev%c_pts)
        else
           call FC_Setup(mg_c_lev%Nt, mg_lev%Nt, coarsen_factor, mg_lev%f_blocks, mg_lev%c_pts)
        end if
     end do

     do level_index = nlevels,coarsest_level,-1
        mg_lev => mg_ld(level_index)
        pf_lev = pf%levels(level_index)
        !if ((level_index .lt. nlevels) .and. (level_index .gt. 1) .and. (mg_lev%Nt .gt. 0)) then
        if ((level_index .lt. nlevels) .and. (level_index .gt. 1)) then
           call pf_lev%ulevel%factory%create_array(mg_lev%g, max(1, mg_lev%Nt), level_index, pf_lev%lev_shape)
           do i = 1,mg_lev%Nt
              call mg_lev%g(i)%setval(0.0_pfdp)
           end do
        end if
        if (FAS_flag .eqv. .true.) then
           mg_ld(level_index)%FAS_flag = .true.
           !if ((level_index .lt. nlevels) .and. (mg_lev%Nt .gt. 0)) then
           if (level_index .lt. nlevels) then
              call pf_lev%ulevel%factory%create_array(mg_lev%qc_fas, max(1, mg_lev%Nt), level_index, pf_lev%lev_shape)
              do i = 1,mg_lev%Nt
                 call mg_lev%qc_fas(i)%setval(0.0_pfdp)
              end do
           end if
        end if
        if (FCF_flag .eqv. .false.) then
           mg_ld(level_index)%FCF_flag = .false.
        end if
        if (level_index .gt. coarsest_level) then
           mg_c_lev => mg_ld(level_index-1)
           !if ((mg_c_lev%Nt .eq. 0) .and. (mg_lev%Nt .eq. 1)) then
           !   call pf_lev%ulevel%factory%create_array(mg_lev%qc, 1, level_index, pf_lev%lev_shape)
           !   call pf_lev%ulevel%factory%create_array(mg_lev%qc_prev, 1, level_index, pf_lev%lev_shape)
           !   call mg_lev%qc(1)%setval(0.0_pfdp)
           !   call mg_lev%qc_prev(1)%setval(0.0_pfdp)
           !else if (mg_c_lev%Nt .gt. 0) then
              call pf_lev%ulevel%factory%create_array(mg_lev%qc, max(1, mg_c_lev%Nt), level_index, pf_lev%lev_shape)
              call pf_lev%ulevel%factory%create_array(mg_lev%qc_prev, max(1, mg_c_lev%Nt), level_index, pf_lev%lev_shape)
              do i = 1,mg_c_lev%Nt
                 call mg_lev%qc(i)%setval(0.0_pfdp)
                 call mg_lev%qc_prev(i)%setval(0.0_pfdp)
              end do
           !end if
        end if
        call pf_lev%ulevel%factory%create_single(mg_lev%qc_boundary, level_index, pf_lev%lev_shape)
        call pf_lev%ulevel%factory%create_single(mg_lev%q_temp, level_index, pf_lev%lev_shape)
        call pf_lev%ulevel%factory%create_single(mg_lev%g_temp, level_index, pf_lev%lev_shape)
        call pf_lev%ulevel%factory%create_single(mg_lev%r, level_index, pf_lev%lev_shape)
        call pf_lev%ulevel%factory%create_single(mg_lev%Q0, level_index, pf_lev%lev_shape)
        if (FAS_flag .eqv. .true.) then
           call pf_lev%ulevel%factory%create_single(mg_lev%qc_fas_boundary, level_index, pf_lev%lev_shape)
        end if
     end do

     if (mg_ld(nlevels)%coarsest_level .eq. nlevels) then
        print *,'Error in MGRIT setup: number of levels must be greater than 1'
        stop
     end if
  end subroutine mgrit_initialize

  subroutine FC_Setup(Nc, Nf, coarsen_factor, f_blocks, c_pts)
     type(int_vector), allocatable, intent(inout) :: f_blocks(:)
     integer, allocatable, intent(inout) :: c_pts(:)
     integer, intent(in) :: Nc, Nf, coarsen_factor
     integer, allocatable :: f_pts(:)
     integer :: i, j, n
     integer :: j_start, j_end, jj
     integer :: f_size

     n = coarsen_factor

     if (Nc .gt. 0) then
        allocate(c_pts(Nc))
        j = min(Nf, n)
        do i = 1,Nc
           c_pts(i) = j;
           j = j + n;
        end do
     end if

     if (Nf .gt. 0) then
        if (Nc .gt. 0) then
           allocate(f_blocks(Nc))
           do i = 1,Nc
              j_start = (i-1)*n+1;
              j_end = i*n-1;
              f_pts = [(jj, jj = j_start,j_end)]
              f_size = size(f_pts)
              allocate(f_blocks(i)%val(f_size))
              j = 1
              do jj = j_start,j_end
                 f_blocks(i)%val(j) = jj
                 j = j + 1
              end do
           end do
        else
           allocate(f_blocks(1))
           allocate(f_blocks(1)%val(1))
           f_blocks(1)%val(1) = 1
        end if
     end if
  end subroutine FC_Setup

  !>  Do the MGRIT algorithm
  subroutine pf_MGRIT_run(pf, mg_ld, Q0, qend)
    type(pf_pfasst_t), target, intent(inout)   :: pf   !!  The complete PFASST structure
    type(mgrit_level_data), allocatable, target, intent(inout) :: mg_ld(:)
    class(pf_encap_t), intent(in   )           :: Q0   !!  The initial condition
    class(pf_encap_t), intent(inout), optional :: qend    !!  The computed solution at tend
    class(pf_encap_t), pointer :: qc(:)
    type(mgrit_level_data), pointer :: mg_lev
    type(pf_level_t), pointer :: pf_lev
    logical :: zero_rhs_flag

    !  Local variables
    integer :: nproc  !!  Total number of processors
    integer :: nsteps_loc  !!  local number of time steps
    real(pfdp) :: tend_loc !!  The final time of run
    integer :: ierror, iter, k, i, level_index
    integer :: Nt, Nt_glob
    integer :: qc_len
    real(pfdp) :: t0, dt, big_dt, t0_glob
    integer :: nlevels, coarsest_level

    k = 1
    nlevels = pf%nlevels
    coarsest_level = mg_ld(nlevels)%coarsest_level
    pf%state%pfblock = 1
    pf%state%nsteps = pf%comm%nproc
    pf%state%sweep = 1
    pf%state%pstatus = PF_STATUS_ITERATING

    !>  Allocate stuff for holding results
    call initialize_results(pf)

    if (pf%save_timings > 0) call pf_start_timer(pf, T_TOTAL)

    mg_ld(nlevels)%cycle_phase = 0
    level_index = coarsest_level
    call pf_start_timer(pf, T_SWEEP, level_index)
    call InitExactSolve(pf, mg_ld, Q0, level_index)
    zero_rhs_flag = .true.
    call IdealInterp(pf, mg_ld, 0, zero_rhs_flag)
    call pf_stop_timer(pf, T_SWEEP, level_index)
    !qc => mg_ld(nlevels)%qc
    !call qend%copy(qc(size(qc)))
    !return

    do level_index = coarsest_level,nlevels
       pf_lev => pf%levels(level_index)
       call pf_lev%qend%setval(0.0_pfdp)
       call pf_lev%q0%setval(0.0_pfdp)
    end do

    !>  Try to sync everyone
    call mpi_barrier(pf%comm%comm, ierror)

    !> Start timer

    call pf_start_timer(pf, T_TOTAL)

    do iter = 1, pf%niters
       pf%state%iter = iter

       do level_index = coarsest_level,nlevels
          mg_lev => mg_ld(level_index)
          pf_lev => pf%levels(level_index)
          qc_len = size(mg_lev%qc)
          if (level_index .gt. coarsest_level) then
             call pf_lev%delta_q0%copy(mg_lev%qc(qc_len))
          else
             call pf_lev%delta_q0%copy(pf_lev%qend)
          end if
       end do

       !  Do a v_cycle
       call pf_start_timer(pf, T_ITERATION)
       call pf_MGRIT_v_cycle(pf, mg_ld, iter)
       call pf_set_iter(pf, iter)
       call pf_stop_timer(pf, T_ITERATION)

       do level_index = coarsest_level,nlevels
          mg_lev => mg_ld(level_index)
          pf_lev => pf%levels(level_index)
          qc_len = size(mg_lev%qc)
          if (level_index .gt. coarsest_level) then
             call pf_lev%delta_q0%axpy(-1.0_pfdp, mg_lev%qc(qc_len))
          else
             call pf_lev%delta_q0%axpy(-1.0_pfdp, pf_lev%qend)
          end if
          pf_lev%max_delta_q0 = pf_lev%delta_q0%norm()
          call pf_set_delta_q0(pf, level_index, pf_lev%delta_q0%norm())
          call pf_set_resid(pf, level_index, mg_lev%res_norm_loc(1))
       end do

       qc => mg_ld(nlevels)%qc
       call pf%levels(nlevels)%qend%copy(qc(size(qc)))
       call call_hooks(pf, nlevels, PF_POST_ITERATION)
       if (pf%state%pstatus .eq. PF_STATUS_CONVERGED) then
          exit
       end if
    end do
    call pf_stop_timer(pf, T_TOTAL)

    call pf_dump_stats(pf)

    qc => mg_ld(nlevels)%qc
    call qend%copy(qc(size(qc)))

   ! qc => mg_ld(pf%nlevels)%qc

   ! if (pf%rank .eq. 0) then
   !    do i = 1,size(qc)
   !       print *,i,qc(i)%norm()
   !    end do
   ! end if
   ! call mpi_barrier(pf%comm%comm, ierror)
   ! call mpi_barrier(pf%comm%comm, ierror)
   ! call mpi_barrier(pf%comm%comm, ierror)
   ! if (pf%rank .eq. 1) then
   !    do i = 1,size(qc)
   !       print *,i,qc(i)%norm()
   !    end do
   ! end if

   ! if (pf%rank .eq. pf%comm%nproc-1) then
   !    call qend%eprint()
   !
   !    t0 = mg_ld(pf%nlevels)%T0_glob
   !    Nt = mg_ld(pf%nlevels)%Nt * pf%comm%nproc
   !    dt = mg_ld(pf%nlevels)%dt
   !    big_dt = dt * real(Nt, pfdp)
   !    call pf%levels(pf%nlevels)%q0%copy(q0)
   !    call pf%levels(pf%nlevels)%ulevel%stepper%do_n_steps(pf, pf%nlevels, t0, pf%levels(pf%nlevels)%q0, &
   !                               pf%levels(pf%nlevels)%qend, big_dt, Nt)
   !    call qend%copy(pf%levels(pf%nlevels)%qend)
   ! end if
  end subroutine pf_MGRIT_run

  !> Execute a MGRIT V-cycle (iteration)
  !!  It is assumed that we have two levels and two nodes here
  !!  When this is called the previous coarse integrator result should be stored in Q(1)
  !!  and the MGRIT iteration in qend (both on coarse level).  If this is called
  !!  directly after the predictor, these will be the same thing
  subroutine pf_MGRIT_v_cycle(pf, mg_ld, iteration)
    type(pf_pfasst_t), target, intent(inout) :: pf
    type(mgrit_level_data), allocatable, target, intent(inout) :: mg_ld(:)
    integer, intent(in) :: iteration
    type(pf_level_t), pointer :: pf_lev, pf_f_lev, pf_c_lev
    type(mgrit_level_data), pointer :: mg_lev, mg_f_lev, mg_c_lev
    logical :: zero_rhs_flag, zero_c_pts_flag, interp_flag
    integer :: qc_len
    integer :: i, j, k, n, ierr
    integer :: nlevels, level_index, level_index_f
    integer :: ierror
    real(pfdp) :: res_norm_glob(1)
    class(pf_encap_t), pointer :: qc(:)
    real(pfdp), pointer :: qend(:)
    real(pfdp) :: qexact, maxerr
    integer :: coarsest_level

    nlevels = pf%nlevels;
    coarsest_level = mg_ld(nlevels)%coarsest_level

    !> FCF-relaxation on finest grid
    mg_ld(nlevels)%cycle_phase = 1
    level_index = nlevels
    mg_lev => mg_ld(level_index)
    call pf_start_timer(pf, T_SWEEP, level_index)
    call FCF_Relax(pf, mg_ld, level_index, iteration)
    call pf_stop_timer(pf, T_SWEEP, level_index)

    if (pf%state%pstatus .eq. PF_STATUS_CONVERGED) then
        return
    end if

    do level_index = nlevels-1,coarsest_level+1,-1
        level_index_f = level_index+1
        mg_lev => mg_ld(level_index)
        mg_f_lev => mg_ld(level_index_f)

        !> Restrict residual
        mg_f_lev%res_norm_loc(1) = 0.0_pfdp
        if (mg_lev%Nt .gt. 0) then
           call Restrict(pf, mg_ld, level_index, level_index_f)
           do i = 1,size(mg_f_lev%qc)
              call mg_f_lev%qc(i)%copy(mg_f_lev%qc_prev(i))
           end do
        end if
        !if ((level_index_f .eq. nlevels) .and. (iteration .gt. 1)) then
        !   call mpi_allreduce(mg_f_lev%res_norm_loc(1), res_norm_glob(1), 1, myMPI_Datatype, MPI_MAX, pf%comm%comm, ierr)
        !   if (res_norm_glob(1) .lt. pf%abs_res_tol) then
        !      pf%state%pstatus = PF_STATUS_CONVERGED
        !      return
        !   end if
        !end if
        !> FCF-relaxation on intermediate grids
        call pf_start_timer(pf, T_SWEEP, level_index)
        call FCF_Relax(pf, mg_ld, level_index, iteration)
        call pf_stop_timer(pf, T_SWEEP, level_index)
    end do

    !> Coarsest grid solve
    mg_ld(nlevels)%cycle_phase = 2
    level_index = coarsest_level
    if (mg_lev%Nt .gt. 0) then
       call pf_start_timer(pf, T_SWEEP, level_index)
       call ExactSolve(pf, mg_ld, level_index)
       call pf_stop_timer(pf, T_SWEEP, level_index)
       !if ((level_index+1 .eq. nlevels) .and. (iteration .gt. 1)) then
       !   mg_f_lev => mg_ld(nlevels)
       !   call mpi_allreduce(mg_f_lev%res_norm_loc(1), res_norm_glob(1), 1, myMPI_Datatype, MPI_MAX, pf%comm%comm, ierr)
       !   if (res_norm_glob(1) .lt. pf%abs_res_tol) then
       !      pf%state%pstatus = PF_STATUS_CONVERGED
       !      return
       !   end if
       !end if
    end if
    zero_rhs_flag = .false.
    mg_ld(nlevels)%cycle_phase = 3
    call pf_start_timer(pf, T_SWEEP, level_index)
    call IdealInterp(pf, mg_ld, iteration, zero_rhs_flag)
    call pf_stop_timer(pf, T_SWEEP, level_index)
  end subroutine pf_MGRIT_v_cycle

  subroutine FCF_Relax(pf, mg_ld, level_index, iteration)
     type(pf_pfasst_t), target, intent(inout) :: pf
     type(mgrit_level_data), allocatable, target, intent(inout) :: mg_ld(:)
     integer, intent(in) :: iteration
     integer, intent(in) :: level_index
     type(pf_level_t), pointer :: pf_lev
     type(mgrit_level_data), pointer :: mg_lev, mg_c_lev
     logical :: zero_rhs_flag, zero_c_pts_flag, interp_flag, send_flag, recv_flag
     integer :: qc_len, i, nlevels
     integer :: ierror
     class(pf_encap_t), pointer :: qc(:)

     nlevels = pf%nlevels
     mg_lev => mg_ld(level_index)
     mg_c_lev => mg_ld(level_index-1)
     pf_lev => pf%levels(level_index)

     if (level_index .lt. nlevels) then
        zero_c_pts_flag = .true.
     else
        zero_c_pts_flag = .false.
     end if

     if (level_index .lt. nlevels) then
        zero_rhs_flag = .false.
     else
        zero_rhs_flag = .true.
     end if

     interp_flag = .false.

     qc_len = size(mg_lev%qc)
     do i = 1,qc_len
        call mg_lev%qc_prev(i)%copy(mg_lev%qc(i))
     end do

     if (zero_c_pts_flag .eqv. .false.) then
        !if ((pf%rank .lt. pf%comm%nproc-1) .and. (mg_lev%c_pts_flag .eqv. .true.)) then
           call pf_lev%qend%copy(mg_lev%qc(qc_len))
        !end if
        send_flag = mg_lev%c_pts_flag
        recv_flag = mg_lev%f_pts_flag
        call mgrit_send_recv(pf, mg_ld, level_index, send_flag, recv_flag)
     else
        call pf_lev%q0%setval(0.0_pfdp)
     end if

     if (mg_lev%f_pts_flag .eqv. .true.) then
        call F_Relax(pf, mg_ld, level_index, zero_rhs_flag, interp_flag, zero_c_pts_flag)
     end if
     if (mg_lev%Nt .eq. 1) then
        !if ((pf%rank .lt. pf%comm%nproc-1) .and. (mg_lev%f_pts_flag .eqv. .true.)) then
           call pf_lev%qend%copy(mg_lev%qc(1))
        !end if
        send_flag = mg_lev%f_pts_flag
        recv_flag = mg_lev%c_pts_flag
        call mgrit_send_recv(pf, mg_ld, level_index, send_flag, recv_flag)
        if (recv_flag .eqv. .true.) then
           call mg_lev%qc(1)%copy(mg_lev%qc_boundary)
        end if
     end if

     if (level_index .eq. nlevels) then
        call pf_start_timer(pf, T_RESIDUAL, level_index)
        call ResidualNorm_Cpoints_L2norm(pf, mg_ld, level_index, zero_rhs_flag)
        call pf_stop_timer(pf, T_RESIDUAL, level_index)
        if (pf%state%pstatus .eq. PF_STATUS_CONVERGED) then
           return
        end if
     end if

     if (mg_lev%FCF_flag .eqv. .true.) then
        zero_c_pts_flag = .false.
        if (mg_lev%c_pts_flag .eqv. .true.) then
           call C_Relax(pf, mg_ld, level_index, zero_rhs_flag, interp_flag)
           call pf_lev%qend%copy(mg_lev%qc(qc_len))
        end if
        send_flag = mg_lev%c_pts_flag
        recv_flag = mg_lev%f_pts_flag
        call mgrit_send_recv(pf, mg_ld, level_index, send_flag, recv_flag)
        if (mg_ld(nlevels)%FAS_flag .eqv. .true.) then
           call mg_lev%qc_fas_boundary%copy(mg_lev%qc_boundary)
        end if
        do i = 1,qc_len
           call mg_lev%qc_prev(i)%copy(mg_lev%qc(i))
        end do
        if (mg_lev%f_pts_flag .eqv. .true.) then
           call F_Relax(pf, mg_ld, level_index, zero_rhs_flag, interp_flag, zero_c_pts_flag)
        end if
        if (mg_lev%Nt .eq. 1) then
           call pf_lev%qend%copy(mg_lev%qc(1))
           send_flag = mg_lev%f_pts_flag
           recv_flag = mg_lev%c_pts_flag
           call mgrit_send_recv(pf, mg_ld, level_index, send_flag, recv_flag)
           if (recv_flag .eqv. .true.) then
              call mg_lev%qc(1)%copy(mg_lev%qc_boundary)
           end if
        end if
     end if
  end subroutine FCF_Relax

  subroutine IdealInterp(pf, mg_ld, iteration, zero_rhs_flag)
    type(pf_pfasst_t), target, intent(inout) :: pf
    type(mgrit_level_data), allocatable, target, intent(inout) :: mg_ld(:)
    integer, intent(in) :: iteration
    logical, intent(in) :: zero_rhs_flag
    type(pf_level_t), pointer :: pf_lev, pf_f_lev, pf_c_lev
    type(mgrit_level_data), pointer :: mg_lev, mg_f_lev, mg_c_lev
    logical :: zero_c_pts_flag, interp_flag, send_flag, recv_flag
    integer :: qc_len
    integer :: i, j, k, n, ierr
    integer :: nlevels, level_index, level_index_f, level_index_c
    class(pf_encap_t), pointer :: qc(:)
    integer :: coarsest_level

    nlevels = pf%nlevels;
    coarsest_level = mg_ld(nlevels)%coarsest_level
    zero_c_pts_flag = .false.
    interp_flag = .true.
    do level_index = coarsest_level+1,(nlevels-1)
       level_index_f = level_index+1
       mg_lev => mg_ld(level_index)
       pf_lev => pf%levels(level_index)
       mg_f_lev => mg_ld(level_index_f)
       pf_f_lev => pf%levels(level_index_f)
       qc_len = size(mg_lev%qc)
       do i = 1,qc_len
          call mg_lev%qc_prev(i)%copy(mg_lev%qc(i))
       end do
       if (mg_lev%c_pts_flag .eqv. .true.) then
          call pf_lev%qend%copy(mg_lev%qc(qc_len))
       end if
       send_flag = mg_lev%c_pts_flag
       recv_flag = mg_lev%f_pts_flag
       call mgrit_send_recv(pf, mg_ld, level_index, send_flag, recv_flag)
       if (mg_lev%f_pts_flag .eqv. .true.) then
          call F_Relax(pf, mg_ld, level_index, zero_rhs_flag, interp_flag, zero_c_pts_flag)
       else
          call mg_lev%q_temp%copy(mg_lev%qc_prev(1))
          if (mg_ld(nlevels)%FAS_flag .eqv. .true.) then
             call mg_lev%q_temp%axpy(-1.0_pfdp, mg_lev%qc_fas(1))
          end if
          call pf_f_lev%ulevel%interpolate(pf_f_lev, pf_lev, mg_f_lev%q_temp, mg_lev%q_temp, mg_lev%t0)
          call mg_f_lev%qc(1)%axpy(1.0_pfdp, mg_f_lev%q_temp)
       end if

       do i = 1,qc_len
          call mg_lev%qc(i)%copy(mg_lev%qc_prev(i))
       end do

    end do
  end subroutine IdealInterp

  subroutine F_Relax(pf, mg_ld, level_index, zero_rhs_flag, interp_flag, zero_c_pts_flag)
     type(pf_pfasst_t), target, intent(inout) :: pf
     type(mgrit_level_data), allocatable, target, intent(inout) :: mg_ld(:)
     integer, intent(in) :: level_index
     logical, intent(in) :: zero_rhs_flag, zero_c_pts_flag, interp_flag
     type(mgrit_level_data), pointer :: mg_f_lev, mg_lev
     type(pf_level_t), pointer :: pf_lev, pf_f_lev
     integer :: i, j, n, ii
     integer :: nlevels
     logical :: q_zero_flag
     integer :: level_index_f

     nlevels = pf%nlevels
     level_index_f = level_index + 1;
     mg_lev => mg_ld(level_index)
     pf_lev => pf%levels(level_index)
     if (interp_flag .eqv. .true.) then
        mg_f_lev => mg_ld(level_index_f)
        pf_f_lev => pf%levels(level_index_f)
     end if

     do i = 1,size(mg_lev%f_blocks)
        if ((i .eq. 1) .and. (mg_lev%rank_shifted .eq. 0)) then
           if ((level_index .eq. nlevels) .or. (mg_ld(nlevels)%cycle_phase .eq. 0)) then
              call pf_lev%q0%copy(mg_lev%Q0)
           else if (mg_ld(nlevels)%FAS_flag .eqv. .true.) then
              if (zero_c_pts_flag .eqv. .true.) then
                 call pf_lev%q0%setval(0.0_pfdp)
              else
                 call pf_lev%q0%copy(mg_lev%Q0)
              end if
           else
              call pf_lev%q0%setval(0.0_pfdp)
           end if
        else
           if (zero_c_pts_flag .eqv. .true.) then
              call pf_lev%q0%setval(0.0_pfdp)
           else if (i .gt. 1) then
              call pf_lev%q0%copy(mg_lev%qc_prev(i-1))
           end if
        end if

        do n = 1,size(mg_lev%f_blocks(i)%val)
           j = mg_lev%f_blocks(i)%val(n)
           !> Do a single step
           call PointRelax(pf, mg_ld, level_index, j, pf_lev%q0, pf_lev%qend)
           !> Add g to the result
           if (zero_rhs_flag .eqv. .false.) then
              call pf_lev%qend%axpy(1.0_pfdp, mg_lev%g(j))
           end if

           !> Interpolate to fine level
           if (interp_flag .eqv. .true.) then
              call pf_start_timer(pf, T_INTERPOLATE, level_index_f)
              call mg_lev%q_temp%copy(pf_lev%qend)
              if ((mg_ld(nlevels)%FAS_flag .eqv. .true.) .and. (mg_ld(nlevels)%cycle_phase .gt. 0)) then
                 call mg_lev%q_temp%axpy(-1.0_pfdp, mg_lev%qc_fas(j))
              end if
              call pf_f_lev%ulevel%interpolate(pf_f_lev, pf_lev, mg_f_lev%q_temp, mg_lev%q_temp, mg_f_lev%t0)
              call mg_f_lev%qc(j)%axpy(1.0_pfdp, mg_f_lev%q_temp)
              call pf_stop_timer(pf, T_INTERPOLATE, level_index_f)
           end if
           call pf_lev%q0%copy(pf_lev%qend)
        end do
        call mg_lev%qc(i)%copy(pf_lev%qend)
        if ((interp_flag .eqv. .true.) .and. (mg_lev%c_pts_flag .eqv. .true.)) then
           call pf_start_timer(pf, T_INTERPOLATE, level_index_f)
           call mg_lev%q_temp%copy(mg_lev%qc_prev(i))
           if ((mg_ld(nlevels)%FAS_flag .eqv. .true.) .and. (mg_ld(nlevels)%cycle_phase .gt. 0)) then
              call mg_lev%q_temp%axpy(-1.0_pfdp, mg_lev%qc_fas(j+1))
           end if
           call pf_f_lev%ulevel%interpolate(pf_f_lev, pf_lev, mg_f_lev%q_temp, mg_lev%q_temp, mg_lev%t0)
           call mg_f_lev%qc(j+1)%axpy(1.0_pfdp, mg_f_lev%q_temp)
           call pf_stop_timer(pf, T_INTERPOLATE, level_index_f)
        end if
     end do
  end subroutine F_Relax

  subroutine C_Relax(pf, mg_ld, level_index, zero_rhs_flag, interp_flag)
     type(pf_pfasst_t), target, intent(inout) :: pf
     type(mgrit_level_data), allocatable, target, intent(inout) :: mg_ld(:)
     integer, intent(in) :: level_index
     logical, intent(in) :: zero_rhs_flag, interp_flag
     integer :: i, j, k, n, ii
     integer :: nlevels
     logical :: q_zero_flag
     type(mgrit_level_data), pointer :: mg_lev, mg_finest_lev
     type(pf_level_t), pointer :: pf_lev

     nlevels = pf%nlevels

     mg_lev => mg_ld(level_index)

     pf_lev => pf%levels(level_index)

     do i = 1,size(mg_lev%c_pts)
        j = mg_lev%c_pts(i)
        call pf_lev%q0%copy(mg_lev%qc(i))
        call PointRelax(pf, mg_ld, level_index, j, pf_lev%q0, pf_lev%qend)
        call mg_lev%qc(i)%copy(pf_lev%qend)
        if (zero_rhs_flag .eqv. .false.) then
           call mg_lev%qc(i)%axpy(1.0_pfdp, mg_lev%g(j))
        end if
    end do
  end subroutine C_Relax

  subroutine mgrit_send_recv(pf, mg_ld, level_index, send_flag, recv_flag)
    type(pf_pfasst_t), target, intent(inout) :: pf
    type(mgrit_level_data), allocatable, target, intent(inout) :: mg_ld(:)
    integer, intent(in) :: level_index
    logical, intent(in) :: send_flag, recv_flag
    type(pf_level_t), pointer :: pf_lev
    integer :: ierror
    type(mgrit_level_data), pointer :: mg_lev

    if (pf%comm%nproc .gt. 1) then
       pf_lev => pf%levels(level_index)
       mg_lev => mg_ld(level_index)
       if (mod(mg_lev%rank_shifted, 2) .eq. 0) then
          if (mg_lev%rank_shifted .gt. 0) then
             if (recv_flag .eqv. .true.) then
                call mgrit_recv(pf, mg_ld, level_index, 1, .true.)
                call mg_lev%qc_boundary%copy(pf_lev%q0)
             end if
          end if
          if (pf%rank .lt. pf%comm%nproc-1) then
             if (send_flag .eqv. .true.) then
                call mgrit_send(pf, mg_ld, level_index, 1, .true.)
             end if
          end if
       else
          if (pf%rank .lt. pf%comm%nproc-1) then
             if (send_flag .eqv. .true.) then
                call mgrit_send(pf, mg_ld, level_index, 1, .true.)
             end if
          end if
          if (mg_lev%rank_shifted .gt. 0) then
             if (recv_flag .eqv. .true.) then
                call mgrit_recv(pf, mg_ld, level_index, 1, .true.)
                call mg_lev%qc_boundary%copy(pf_lev%q0)
             end if
          end if
       end if
    end if
  end subroutine mgrit_send_recv

  subroutine mgrit_send(pf, mg_ld, level_index, tag, blocking)
    type(pf_pfasst_t), target, intent(inout) :: pf
    type(mgrit_level_data), allocatable, target, intent(inout) :: mg_ld(:)
    integer, intent(in) :: level_index, tag
    logical, intent(in) :: blocking
    type(pf_level_t), pointer :: pf_lev
    integer :: ierror, dir, which, stride

    pf_lev => pf%levels(level_index)

    !call pf_lev%qend%pack(pf_lev%send)
    !call pf_start_timer(pf, T_SEND, level_index)
    !call pf%comm%send(pf, pf_lev, tag, blocking, ierror, mg_ld(level_index)%send_to_rank)
    !call pf_stop_timer(pf, T_SEND, level_index)

    dir = 1
    which = 1
    stride = mg_ld(level_index)%send_to_rank - pf%rank
    call pf_send(pf, pf_lev, tag, blocking, dir, which, stride)
  end subroutine mgrit_send

  subroutine mgrit_recv(pf, mg_ld, level_index, tag, blocking)
    type(pf_pfasst_t), target, intent(inout) :: pf
    type(mgrit_level_data), allocatable, target, intent(inout) :: mg_ld(:)
    integer, intent(in) :: level_index, tag
    logical, intent(in) :: blocking
    type(pf_level_t), pointer :: pf_lev
    integer :: ierror, dir, which, stride

    pf_lev => pf%levels(level_index)

    !call pf_start_timer(pf, T_RECEIVE, level_index)
    !call pf%comm%recv(pf, pf_lev, tag, blocking, ierror, mg_ld(level_index)%recv_from_rank)
    !call pf_stop_timer(pf, T_RECEIVE, level_index)
    !call pf_lev%q0%unpack(pf_lev%recv)

    dir = 1
    which = 1
    stride = pf%rank - mg_ld(level_index)%recv_from_rank
    call pf_recv(pf, pf_lev, tag, blocking, dir, which, stride)
  end subroutine mgrit_recv

  subroutine ExactSolve(pf, mg_ld, level_index)
     type(pf_pfasst_t), target, intent(inout) :: pf
     type(mgrit_level_data), allocatable, target, intent(inout) :: mg_ld(:)
     integer, intent(in) :: level_index
     integer :: level_index_f
     type(pf_level_t), pointer :: lev
     class(pf_encap_t), allocatable :: gi
     integer :: i, j, k, n, ii
     integer :: nlevels
     logical :: zero_rhs_flag, send_flag, recv_flag
     type(mgrit_level_data), pointer :: mg_lev, mg_f_lev, mg_finest_lev
     type(pf_level_t), pointer :: pf_lev, pf_f_lev

     mg_lev => mg_ld(level_index)
     nlevels = pf%nlevels
     level_index_f = level_index + 1

     mg_f_lev => mg_ld(level_index_f)
     mg_finest_lev => mg_ld(nlevels)

     pf_lev => pf%levels(level_index)
     pf_f_lev => pf%levels(level_index_f)

     if ((mg_ld(nlevels)%FAS_flag .eqv. .true.) .and. (mg_f_lev%Nt .eq. 1)) then
        call pf_f_lev%qend%copy(mg_f_lev%qc_boundary)
        send_flag = mg_f_lev%f_pts_flag
        recv_flag = mg_f_lev%c_pts_flag
        call mgrit_send_recv(pf, mg_ld, level_index_f, send_flag, recv_flag)
        if (recv_flag .eqv. .true.) then
           call mg_f_lev%qc_fas_boundary%copy(mg_f_lev%qc_boundary)
        end if
     end if

     if (mg_lev%Nt .gt. 0) then
        mg_f_lev%res_norm_loc(1) = 0.0_pfdp

        if ((mg_lev%rank_shifted .gt. 0) .and. (pf%comm%nproc .gt. 1)) then
           call mgrit_recv(pf, mg_ld, level_index, 3, .true.)
        else
           if (mg_ld(nlevels)%FAS_flag .eqv. .true.) then
              call pf_lev%q0%copy(mg_lev%Q0)
           else
              call pf_lev%q0%setval(0.0_pfdp);
           end if
        end if

        !call pf_lev%ulevel%factory%create_single(gi, level_index, pf_lev%lev_shape)
        do i = 1,mg_lev%Nt
           ii = mg_f_lev%c_pts(i)
           if (level_index_f .eq. nlevels) then
               zero_rhs_flag = .true.
           else
               zero_rhs_flag = .false.
           end if
           call InjectRestrictPoint(pf, mg_ld, mg_lev%g_temp, level_index, level_index_f, i, ii, zero_rhs_flag)
           call mg_f_lev%qc(i)%copy(mg_f_lev%qc_prev(i))

           call PointRelax(pf, mg_ld, level_index, i, pf_lev%q0, pf_lev%qend)
           call pf_lev%qend%axpy(1.0_pfdp, mg_lev%g_temp)
           call pf_lev%q0%copy(pf_lev%qend)

           call mg_lev%q_temp%copy(pf_lev%qend)
           if ((mg_ld(nlevels)%FAS_flag .eqv. .true.) .and. (mg_ld(nlevels)%cycle_phase .gt. 0)) then
              call mg_lev%q_temp%axpy(-1.0_pfdp, mg_lev%qc_fas(i))
           end if
           call pf_start_timer(pf, T_INTERPOLATE, level_index_f)
           call pf_f_lev%ulevel%interpolate(pf_f_lev, pf_lev, mg_f_lev%q_temp, mg_lev%q_temp, mg_f_lev%t0)
           call pf_stop_timer(pf, T_INTERPOLATE, level_index_f)
           call mg_f_lev%qc(i)%axpy(1.0_pfdp, mg_f_lev%q_temp)
        end do
        !call pf_lev%ulevel%factory%destroy_single(gi)

        !pf_f_lev%residual = mg_f_lev%res_norm_loc(1)
        !pf_lev%residual = 0.0_pfdp

        if ((pf%rank .lt. pf%comm%nproc-1) .and. (pf%comm%nproc .gt. 1)) then
           call mgrit_send(pf, mg_ld, level_index, 3, .true.)
        end if

     end if
  end subroutine ExactSolve

  subroutine InitExactSolve(pf, mg_ld, Q0, level_index)
     type(pf_pfasst_t), target, intent(inout) :: pf
     type(mgrit_level_data), allocatable, target, intent(inout) :: mg_ld(:)
     class(pf_encap_t), intent(in) :: Q0
     integer, intent(in) :: level_index
     integer :: level_index_f
     type(pf_level_t), pointer :: lev
     integer :: i, j, k, n, ii, l
     integer :: nlevels
     type(mgrit_level_data), pointer :: mg_lev, mg_f_lev, mg_finest_lev
     type(pf_level_t), pointer :: pf_lev, pf_f_lev

     mg_lev => mg_ld(level_index)
     nlevels = pf%nlevels
     level_index_f = level_index + 1
     mg_f_lev => mg_ld(level_index_f)
     mg_finest_lev => mg_ld(nlevels)

     pf_lev => pf%levels(level_index)
     pf_f_lev => pf%levels(level_index_f)

     if ((mg_lev%rank_shifted .gt. 0) .and. (pf%comm%nproc .gt. 1)) then
        call mgrit_recv(pf, mg_ld, level_index, 3, .true.)
     else
        call mg_finest_lev%Q0%copy(Q0)
        do l = nlevels,level_index+1,-1
           call pf%levels(l)%ulevel%restrict(pf%levels(l), pf%levels(l-1), mg_ld(l)%Q0, mg_ld(l-1)%Q0, mg_ld(l)%t0)
        end do
        call pf_lev%q0%copy(mg_lev%Q0)
     end if

     if (mg_lev%Nt .gt. 0) then
        do i = 1,mg_lev%Nt
           call PointRelax(pf, mg_ld, level_index, i, pf_lev%q0, pf_lev%qend)
           call pf_start_timer(pf, T_INTERPOLATE, level_index_f)
           call pf_lev%ulevel%interpolate(pf_f_lev, pf_lev, mg_f_lev%qc(i), pf_lev%qend, mg_f_lev%t0)
           call pf_stop_timer(pf, T_INTERPOLATE, level_index_f)
           call pf_lev%q0%copy(pf_lev%qend)
        end do

        if ((pf%rank .lt. pf%comm%nproc-1) .and. (pf%comm%nproc .gt. 1)) then
           call mgrit_send(pf, mg_ld, level_index, 3, .true.)
        end if
     end if
  end subroutine InitExactSolve

  subroutine Restrict(pf, mg_ld, level_index_c, level_index_f)
     type(pf_pfasst_t), target, intent(inout) :: pf
     type(mgrit_level_data), allocatable, target, intent(inout) :: mg_ld(:)
     integer, intent(in) :: level_index_c, level_index_f
     class(pf_encap_t), pointer :: gi
     integer :: i_c, i_f, k, n
     integer :: nlevels
     logical :: zero_rhs_flag, send_flag, recv_flag
     real(pfdp) :: r_norm
     type(mgrit_level_data), pointer :: mg_c_lev, mg_f_lev
     type(pf_level_t), pointer :: pf_lev, pf_f_lev

     nlevels = pf%nlevels

     mg_c_lev => mg_ld(level_index_c)
     mg_f_lev => mg_ld(level_index_f)

     pf_f_lev => pf%levels(level_index_f)

     mg_f_lev%res_norm_loc(1) = 0.0_pfdp
     mg_f_lev%res_norm_index = 0

     if ((mg_ld(nlevels)%FAS_flag .eqv. .true.) .and. (mg_f_lev%Nt .eq. 1)) then
        call pf_f_lev%qend%copy(mg_f_lev%qc_boundary)
        send_flag = mg_f_lev%f_pts_flag
        recv_flag = mg_f_lev%c_pts_flag
        call mgrit_send_recv(pf, mg_ld, level_index_f, send_flag, recv_flag)
        if (recv_flag .eqv. .true.) then
           call mg_f_lev%qc_fas_boundary%copy(mg_f_lev%qc_boundary)
        end if
     end if

     do i_c = 1,mg_c_lev%Nt
        i_f = mg_f_lev%c_pts(i_c)
        if (level_index_f .eq. nlevels) then
           zero_rhs_flag = .true.;
        else
           zero_rhs_flag = .false.;
        end if
        call InjectRestrictPoint(pf, mg_ld, mg_c_lev%g(i_c), level_index_c, level_index_f, i_c, i_f, zero_rhs_flag);
     end do

     !pf_f_lev%residual = mg_f_lev%res_norm_loc(1)
  end subroutine Restrict

  subroutine InjectRestrictPoint(pf, mg_ld, gci, level_index_c, level_index_f, i_c, i_f, zero_rhs_flag)
     type(pf_pfasst_t), target, intent(inout) :: pf
     type(mgrit_level_data), allocatable, target, intent(inout) :: mg_ld(:)
     class(pf_encap_t), intent(inout) :: gci
     integer, intent(in) :: level_index_c, level_index_f, i_c, i_f
     logical, intent(in) :: zero_rhs_flag
     type(mgrit_level_data), pointer :: mg_c_lev, mg_f_lev
     type(pf_level_t), pointer :: pf_f_lev, pf_c_lev
     integer :: nlevels

     mg_c_lev => mg_ld(level_index_c)
     mg_f_lev => mg_ld(level_index_f)

     pf_f_lev => pf%levels(level_index_f)
     pf_c_lev => pf%levels(level_index_c)
     nlevels = pf%nlevels

     !if (mg_f_lev%Nt .eq. 1) then
     !   call mg_f_lev%q_temp%copy(mg_f_lev%qc_boundary)
     !else
        call mg_f_lev%q_temp%copy(mg_f_lev%qc(i_c))
     !end if
     call PointRelax(pf, mg_ld, level_index_f, i_f, mg_f_lev%q_temp, mg_f_lev%r)
     call mg_f_lev%r%axpy(-1.0_pfdp, mg_f_lev%qc_prev(i_c))
     if (zero_rhs_flag .eqv. .false.) then
        call mg_f_lev%r%axpy(1.0_pfdp, mg_f_lev%g(i_f));
     end if
     if (i_c .eq. mg_c_lev%Nt) then
        call ResNorm(pf, mg_ld, level_index_f, mg_f_lev%r, i_c, 2)
     end if
     call pf_start_timer(pf, T_RESTRICT, level_index_f)
     call pf_f_lev%ulevel%restrict(pf_f_lev, pf_c_lev, mg_f_lev%r, mg_c_lev%r, mg_f_lev%t0)
     if (mg_ld(nlevels)%FAS_flag .eqv. .true.) then
        call pf_f_lev%ulevel%restrict(pf_f_lev, pf_c_lev, mg_f_lev%qc_prev(i_c), mg_c_lev%qc_fas(i_c), mg_f_lev%t0)
        call mg_c_lev%r%axpy(1.0_pfdp, mg_c_lev%qc_fas(i_c))
        if (i_c .eq. 1) then
           if (mg_c_lev%rank_shifted .eq. 0) then
              call mg_c_lev%q_temp%copy(mg_c_lev%Q0)
           else
              call mg_c_lev%q_temp%copy(mg_f_lev%qc_fas_boundary)
           end if
        else
           call mg_c_lev%q_temp%copy(mg_c_lev%qc_fas(i_c-1))
        end if
        call PointRelax(pf, mg_ld, level_index_c, i_c, mg_c_lev%q_temp, pf_c_lev%qend)
        call mg_c_lev%r%axpy(-1.0_pfdp, pf_c_lev%qend)
     end if
     call pf_stop_timer(pf, T_RESTRICT, level_index_f)
     call gci%copy(mg_c_lev%r)
  end subroutine InjectRestrictPoint

  subroutine ResidualNorm_Cpoints_L2norm(pf, mg_ld, level_index, zero_rhs_flag)
     type(pf_pfasst_t), target, intent(inout) :: pf
     type(mgrit_level_data), allocatable, target, intent(inout) :: mg_ld(:)
     integer, intent(in) :: level_index
     logical, intent(in) :: zero_rhs_flag
     type(mgrit_level_data), pointer :: mg_lev, mg_c_lev
     type(pf_level_t), pointer :: pf_lev, pf_c_lev
     integer :: level_index_c, nlevels, i, i_c, ierr
     real(pfdp) :: r_norm
     real(pfdp) :: res_norm_glob(1)

     nlevels = pf%nlevels
     level_index_c = level_index-1
     mg_lev => mg_ld(level_index)
     mg_c_lev => mg_ld(level_index_c)

     pf_lev => pf%levels(level_index)
     pf_c_lev => pf%levels(level_index_c)

     do i_c = 1,mg_c_lev%Nt
        i = mg_lev%c_pts(i_c)
        call mg_lev%q_temp%copy(mg_lev%qc(i_c))
        call PointRelax(pf, mg_ld, level_index, i, mg_lev%q_temp, mg_lev%r)
        call mg_lev%r%axpy(-1.0_pfdp, mg_lev%qc_prev(i_c))
        if (zero_rhs_flag .eqv. .false.) then
           call mg_lev%r%axpy(1.0_pfdp, mg_lev%g(i));
        end if
        mg_lev%res_norm_loc(1) = mg_lev%res_norm_loc(1) + (mg_lev%r%norm())**2
     end do

     call mpi_allreduce(mg_lev%res_norm_loc(1), res_norm_glob(1), 1, myMPI_Datatype, MPI_SUM, pf%comm%comm, ierr)
     res_norm_glob(1) = sqrt(res_norm_glob(1))
     pf_lev%residual = res_norm_glob(1)

     if (res_norm_glob(1) .lt. pf%abs_res_tol) then
        pf%state%pstatus = PF_STATUS_CONVERGED
     end if
  end subroutine ResidualNorm_Cpoints_L2norm

  subroutine PointRelax(pf, mg_ld, level_index, n, q0, qend)
     type(pf_pfasst_t), target, intent(inout) :: pf
     type(mgrit_level_data), allocatable, target, intent(inout) :: mg_ld(:)
     integer, intent(in) :: level_index, n
     class(pf_encap_t), intent(inout) :: q0, qend
     type(pf_level_t), pointer :: pf_lev
     type(mgrit_level_data), pointer :: mg_lev
     real(pfdp) :: t0n

     mg_lev => mg_ld(level_index)
     pf_lev => pf%levels(level_index)

     t0n = mg_lev%t0 + mg_lev%dt * real(n-1,pfdp)
     call pf_lev%ulevel%stepper%do_n_steps(pf, level_index, t0n, q0, qend, mg_lev%dt, 1)
  end subroutine PointRelax

  subroutine ResNorm(pf, mg_ld, level_index, r, i, norm_type)
     type(pf_pfasst_t), target, intent(inout) :: pf
     type(mgrit_level_data), allocatable, target, intent(inout) :: mg_ld(:)
     class(pf_encap_t), intent(in) :: r
     integer, intent(in) :: level_index, i
     integer, intent(in) :: norm_type
     type(mgrit_level_data), pointer :: mg_lev
     type(pf_level_t), pointer :: pf_lev
     integer :: nlevels
     real(pfdp) :: r_norm

     nlevels = pf%nlevels

     if (norm_type == 2) then
        r_norm = r%norm()
     else
        mg_lev => mg_ld(level_index)
        pf_lev => pf%levels(level_index)
        r_norm = r%norm()
        if (r_norm .gt. mg_lev%res_norm_loc(1)) then
           mg_lev%res_norm_loc(1) = r_norm
           mg_lev%res_norm_index = i
        end if
     end if
  end subroutine ResNorm

end module pf_mod_MGRIT

max_stars_repo_name: schreiberx/LibPFASST

max_issues_repo_name: schreiberx/LibPFASST

max_forks_repo_name: schreiberx/LibPFASST


text:
module dynamicalLoop

  use constants, only: real12
  use definedtypes, only: settingparam, storedparam, kappagrid
  use greensroutines, only: greensfunc, allocate_3DGF, calculateGF
  use gf_fourier, only: greensfunc_initplan, greensfunc_killplan
  implicit none

contains

  subroutine loop(settings, stored, kgrid_coarse, kgrid_fine, GF_mom,&
       & hybridisation, ADOS, TDOS, ierr)
    type(settingparam), intent(in)  :: settings
    type(storedparam), intent(in)   :: stored
    type(kappagrid), intent(in)     :: kgrid_coarse(:, :, :)
    type(kappagrid), intent(in)     :: kgrid_fine(:, :, :)
    type(greensfunc), intent(inout) :: GF_mom(:, :, :)
    type(greensfunc), intent(inout) :: hybridisation(:, :, :)
    real(real12), intent(out)       :: ADOS(:)
    real(real12), intent(out)       :: TDOS(:)
    integer, intent(out)            :: ierr

    !routine variables
    type(greensfunc), allocatable  :: GF_real(:, :)
    type(greensfunc), allocatable  :: GF_typ(:, :)
    type(greensfunc), allocatable  :: new_hybrid(:, :, :)
    integer                        :: ix, iy, iz, iom
    integer                        :: np_coarse, np_fine
    integer                        :: iloop
    integer                        :: ierr1
    logical                        :: converged

    !allocate variables
    np_coarse = settings%ncell(1)*settings%ncell(2)*settings%ncell(3)
    np_fine = settings%nfpoints(1)*settings%nfpoints(2)*settings%nfpoints(3)

    allocate(GF_real(np_coarse, np_coarse))
    do ix = 1, np_coarse
       do iy = 1, np_coarse
          allocate(GF_real(ix, iy)%GF(settings%nomega))
       enddo
    enddo

    call allocate_3DGF(new_hybrid, settings%ncell(1), settings%ncell(2),&
         & settings%ncell(3), settings%nomega, ierr1)
    ! error handling
    call allocate_3DGF(GF_typ, settings%ncell(1), settings%ncell(2),&
         & settings%ncell(3), settings%nomega, ierr1)
    ! error handling

    !initialise fourier transform settings
    call greensfunc_initplan(GF_mom, ierr1)
    ! error handling

    self_consistent_loop:do iloop = 1, settings%max_loops

       call calculateGF(GF_mom, stored%omega_diff, kgrid_coarse,&
            & hybridisation, ierr1)
       ! error handling

       call momspace_to_realspace(GF_mom, GF_real, ierr1)
       ! error handling

       call disorder_effects(GF_real, GF_typ, ADOS, PDOS, ierr1)
       ! error handling

       call new_hyb(GF_mom, GF_typ, hybridisation, new_hybrid,&
            & settings%mix, ierr1)
       ! error handling

       converged = converge(new_hybrid, hybridisation, &
            & settings%conv_thresh)

       if (converged) exit

    end do self_consistent_loop

    call greensfunc_killplan(ierr1)
    ! error handling

    if (.not.converged) then
       ! not converged error handling
    end if

  end subroutine loop

end module dynamicalLoop

max_stars_repo_name: IOThomas/Anderson-Phonon

max_issues_repo_name: IOThomas/Anderson-Phonon

max_forks_repo_name: IOThomas/Anderson-Phonon


text:
! MinimalDX version 0.1.4 (https://www.github.com/dmey/minimal-dx).
! Copyright 2018-2020 D. Meyer and R. Raustad. Licensed under MIT.

program make_test_data_cooling
  !+ Generates test data for testing the cooling coils.

  ! Using Energy Plus psychrometric routines
  use Psychrometrics, only: PsyWFnTdbRhPb, PsyRhoAirFnPbTdbW

  implicit none

  !+ @note
  !+ The number one ('1') next to the variable description denotes that the variable is dimensionless - e.g.   COP `[1]`

  real :: OutdoorTDryBulb
      !+ Outdoor dry bulb air temperature `[°C]`
  real :: OutdoorRelHum
      !+ Outdoor relative humidity `[1]`
  real :: OutdoorHumRatio
      !+ Outdoor air humidity ratio `[kgH₂O kgAIR⁻¹]`
  real :: OutdoorPressure
      !+ Outdoor barometric pressure `[Pa]`
  real :: InletTDryBulb
      !+ Indoor (inlet) dry bulb air temperature `[°C]`
  real :: InletRelHum
      !+ Indoor (inlet) relative humidity `[1]`
  real :: InletHumRatio
      !+ Indoor (inlet) air humidity ratio `[kgH₂O kgAIR⁻¹]`
  real :: RatedCOP
      !+ Rated Coefficient Of Performance (COP) `[1]`
  real :: RatedTotCap
      !+ Rated (total) system capacity `[W]`
  real :: SensibleCoolingLoad
      !+ Building sensible load to be met `[W]`
  real :: RatedAirMassFlowRate
      !+ Rated air mass flow rate `[kg s⁻¹]
  integer ::  IdxOutdoorTDryBulb
      !+ Index for OutdoorTDryBulb `[1]`
  integer ::  IdxOutdoorRelHum
      !+ Index for OutdoorRelHum `[1]`
  integer ::  IdxOutdoorPressure
      !+ Index for OutdoorPressure `[1]`
  integer ::  IdxInletTDryBulb
      !+ Index for InletTDryBulb `[1]`
  integer ::  IdxInletRelHum
      !+ Index for InletRelHum `[1]`
  integer ::  IdxRatedCOP
      !+ Index for RatedCOP `[1]`
  integer ::  IdxRatedTotCap
      !+ Index for RatedTotCap `[1]`
  integer ::  IdxSensibleCoolingLoad
      !+ Index for SensibleCoolingLoad `[1]`
  integer ::  IterationCounter
      !+ Iteration counter `[1]`
  integer ::  num_args
      !+ Number of command line arguments `[1]`
  character(255) ::  DataSetPath
      !+ Name of dataset including path `[1]`

  IterationCounter = 0

  ! Check that the command line argument is given and that it is exactly 1
  num_args = command_argument_count()
  if ( num_args == 1 ) then
    call get_command_argument(1, DataSetPath)
  else
    error stop "You must give the filename path of the data set you want to save as a command line argument"
  end if

  open (unit=10, file=trim(adjustl(DataSetPath)), action="write", status="replace")
  write (10, '(A220)') "OutdoorTDryBulb|°C OutdoorHumRatio|kg/kg OutdoorPressure|Pa &
                        &InletTDryBulb|°C InletHumRatio|kg/kg RatedCOP|1 &
                        &RatedTotCap|W SensibleCoolingLoad|W RatedAirMassFlowRate|m3/(kg*s) &
                        &OutdoorRelHum|[0-1] InletRelHum|[0-1]"

  do IdxOutdoorTDryBulb = 20, 40, 2
      do IdxOutdoorRelHum = 0, 100, 5
          do IdxOutdoorPressure = 60000, 120000, 10000
              ! The internal temperature is always greater than the outdoor temperature
              do IdxInletTDryBulb = 20, 40, 2
                  do IdxInletRelHum = 0, 100, 5
                      do IdxRatedCOP = 4, 4, 4
                          do IdxRatedTotCap = 4000, 4000, 4000
                              ! SensibleCoolingLoad scaled between 0% and 100% of RatedTotCap with steps of 2000 W
                              do IdxSensibleCoolingLoad   = IdxRatedTotCap * 0, IdxRatedTotCap * 1, IdxRatedTotCap * 0.5

                                  OutdoorTDryBulb     = DBLE(IdxOutdoorTDryBulb)
                                  OutdoorRelHum       = DBLE(IdxOutdoorRelHum) / 100
                                  OutdoorPressure     = DBLE(IdxOutdoorPressure)
                                  InletTDryBulb       = DBLE(IdxInletTDryBulb)
                                  InletRelHum         = DBLE(IdxInletRelHum) / 100
                                  RatedCOP            = DBLE(IdxRatedCOP)
                                  RatedTotCap         = DBLE(IdxRatedTotCap)
                                  SensibleCoolingLoad = DBLE(IdxSensibleCoolingLoad)

                                  OutdoorHumRatio = PsyWFnTdbRhPb(OutdoorTDryBulb, OutdoorRelHum, OutdoorPressure)
                                  InletHumRatio = PsyWFnTdbRhPb(InletTDryBulb, InletRelHum, OutdoorPressure)
                                  RatedAirMassFlowRate = 0.00005 * RatedTotCap * &
                                                          PsyRhoAirFnPbTdbW(OutdoorPressure, InletTDryBulb, InletHumRatio)

                                  write (10, '(11(ES19.10E2))') OutdoorTDryBulb, OutdoorHumRatio, OutdoorPressure,        &
                                                                  InletTDryBulb, InletHumRatio, RatedCOP,                   &
                                                                  RatedTotCap, SensibleCoolingLoad, RatedAirMassFlowRate,   &
                                                                  OutdoorRelHum, InletRelHum
                                  IterationCounter = IterationCounter + 1
                              end do
                          end do
                      end do
                  end do
              end do
          end do
      end do
  end do

  print *, 'Total number of iterations: ', IterationCounter

end program make_test_data_cooling

max_stars_repo_name: dmey/minimal-dx

max_issues_repo_name: dmey/minimal-dx

max_forks_repo_name: dmey/minimal-dx


text:
!----------------------Subroutine--------------------------------------!
!
      SUBROUTINE CRNTNB
!
!-------------------------Disclaimer-----------------------------------!
!
!     This material was prepared as an account of work sponsored by
!     an agency of the United States Government. Neither the
!     United States Government nor the United States Department of
!     Energy, nor Battelle, nor any of their employees, makes any
!     warranty, express or implied, or assumes any legal liability or
!     responsibility for the accuracy, completeness, or usefulness
!     of any information, apparatus, product, software or process
!     disclosed, or represents that its use would not infringe
!     privately owned rights.
!
!----------------------Acknowledgement---------------------------------!
!
!     This software and its documentation were produced with Government
!     support under Contract Number DE-AC06-76RLO-1830 awarded by the
!     United Department of Energy. The Government retains a paid-up
!     non-exclusive, irrevocable worldwide license to reproduce,
!     prepare derivative works, perform publicly and display publicly
!     by or for the Government, including the right to distribute to
!     other Government contractors.
!
!---------------------Copyright Notices--------------------------------!
!
!            Copyright Battelle Memorial Institute, 1996
!                    All Rights Reserved.
!
!----------------------Description-------------------------------------!
!
!     Compute the local maximum Courant numbers.
!
!----------------------Authors-----------------------------------------!
!
!     Written by MD White, Battelle, PNL, October, 1995.
!     Last Modified by MD White, Battelle, PNL, October 10, 1995.
!     $Id: crntnb.F 1080 2017-03-14 16:22:02Z d3c002 $
!

!----------------------Fortran 90 Modules------------------------------!
!
      USE GLB_PAR
      USE TRNSPT
      USE SOLTN
      USE HYST
      USE GRID
      USE FLUXP
      USE FLUXN
      USE FDVP
      USE CONST
!

!----------------------Implicit Double Precision-----------------------!
!
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
!
!----------------------Parameter Statements----------------------------!
!

!
!----------------------Common Blocks-----------------------------------!
!

!
!----------------------Executable Lines--------------------------------!
!
      ISUB_LOG = ISUB_LOG+1
      SUB_LOG(ISUB_LOG) = '/CRNTNB'
      IF( INDEX(SVN_ID(27)(1:1),'$').EQ.0 ) SVN_ID(27) =
     & '$Id: crntnb.F 1080 2017-03-14 16:22:02Z d3c002 $'
!
!---  Aqueous Phase Courant Number  ---
!
      IF( IEQW.GT.0 ) THEN
        DO 100 N = 1,NFLD
          IF( IXP(N).EQ.0 ) GOTO 100
          NPX = NSX(N)
          NPY = NSY(N)
          NPZ = NSZ(N)
          NQX = NSX(N)+1
          IF( INBS(4,N).NE.0 ) NQX = INBS(4,N)
          NQY = NSY(N)+IFLD
          IF( INBS(5,N).NE.0 ) NQY = INBS(5,N)
          NQZ = NSZ(N)+IJFLD
          IF( INBS(6,N).NE.0 ) NQZ = INBS(6,N)
          I = ID(N)
          J = JD(N)
          K = KD(N)
          CRNTL(N) = MAX( ABS(UL(1,NPX)*DT/DXGP(NPX)),
     &      ABS(UL(1,NQX)*DT/DXGP(NQX)),
     &      ABS(VL(1,NPY)*DT/(DYGP(NPY)*RP(I))),
     &      ABS(VL(1,NQY)*DT/(DYGP(NQY)*RP(I))),
     &      ABS(WL(1,NPZ)*DT/DZGP(NPZ)),
     &      ABS(WL(1,NQZ)*DT/DZGP(NQZ)) )
          IF( SL(2,N)*PORD(2,N).GT.EPSL ) THEN
            CRNTL(N) = CRNTL(N)/(SL(2,N)*PORD(2,N))
          ELSE
            CRNTL(N) = 0.D+0
          ENDIF
  100   CONTINUE
      ENDIF
!
!---  Gas Phase Courant Number  ---
!
      IF( IEQA.GT.0 ) THEN
        DO 200 N = 1,NFLD
          IF( IXP(N).EQ.0 ) GOTO 200
          NPX = NSX(N)
          NPY = NSY(N)
          NPZ = NSZ(N)
          NQX = NSX(N)+1
          IF( INBS(4,N).NE.0 ) NQX = INBS(4,N)
          NQY = NSY(N)+IFLD
          IF( INBS(5,N).NE.0 ) NQY = INBS(5,N)
          NQZ = NSZ(N)+IJFLD
          IF( INBS(6,N).NE.0 ) NQZ = INBS(6,N)
          I = ID(N)
          J = JD(N)
          K = KD(N)
          CRNTG(N) = MAX( ABS(UG(1,NPX)*DT/DXGP(NPX)),
     &      ABS(UG(1,NQX)*DT/DXGP(NQX)),
     &      ABS(VG(1,NPY)*DT/(DYGP(NPY)*RP(I))),
     &      ABS(VG(1,NQY)*DT/(DYGP(NQY)*RP(I))),
     &      ABS(WG(1,NPZ)*DT/DZGP(NPZ)),
     &      ABS(WG(1,NQZ)*DT/DZGP(NQZ)) )
          IF( (SG(2,N)-SGT(2,N))*PORD(2,N).GT.EPSL ) THEN
            CRNTG(N) = CRNTG(N)/((SG(2,N)-SGT(2,N))*PORD(2,N))
          ELSE
            CRNTG(N) = 0.D+0
          ENDIF
  200   CONTINUE
      ENDIF
!
!---  NAPL Courant Number  ---
!
      IF( IEQO.GT.0 ) THEN
        DO 300 N = 1,NFLD
          IF( IXP(N).EQ.0 ) GOTO 300
          NPX = NSX(N)
          NPY = NSY(N)
          NPZ = NSZ(N)
          NQX = NSX(N)+1
          IF( INBS(4,N).NE.0 ) NQX = INBS(4,N)
          NQY = NSY(N)+IFLD
          IF( INBS(5,N).NE.0 ) NQY = INBS(5,N)
          NQZ = NSZ(N)+IJFLD
          IF( INBS(6,N).NE.0 ) NQZ = INBS(6,N)
          I = ID(N)
          J = JD(N)
          K = KD(N)
          CRNTN(N) = MAX( ABS(UN(1,NPX)*DT/DXGP(NPX)),
     &      ABS(UN(1,NQX)*DT/DXGP(NQX)),
     &      ABS(VN(1,NPY)*DT/(DYGP(NPY)*RP(I))),
     &      ABS(VN(1,NQY)*DT/(DYGP(NQY)*RP(I))),
     &      ABS(WN(1,NPZ)*DT/DZGP(NPZ)),
     &      ABS(WN(1,NQZ)*DT/DZGP(NQZ)) )
          IF( (SN(2,N)-SNT(2,N))*PORD(2,N).GT.EPSL ) THEN
            CRNTN(N) = CRNTN(N)/((SN(2,N)-SNT(2,N))*PORD(2,N))
          ELSE
            CRNTN(N) = 0.D+0
          ENDIF
  300   CONTINUE
      ENDIF
      ISUB_LOG = ISUB_LOG-1
!
!---  End of CRNTNB group  ---
!
      RETURN
      END

!----------------------Subroutine--------------------------------------!
!
      SUBROUTINE CRN_LIM( NSL )
!
!-------------------------Disclaimer-----------------------------------!
!
!     This material was prepared as an account of work sponsored by
!     an agency of the United States Government. Neither the
!     United States Government nor the United States Department of
!     Energy, nor Battelle, nor any of their employees, makes any
!     warranty, express or implied, or assumes any legal liability or
!     responsibility for the accuracy, completeness, or usefulness
!     of any information, apparatus, product, software or process
!     disclosed, or represents that its use would not infringe
!     privately owned rights.
!
!----------------------Acknowledgement---------------------------------!
!
!     This software and its documentation were produced with Government
!     support under Contract Number DE-AC06-76RLO-1830 awarded by the
!     United Department of Energy. The Government retains a paid-up
!     non-exclusive, irrevocable worldwide license to reproduce,
!     prepare derivative works, perform publicly and display publicly
!     by or for the Government, including the right to distribute to
!     other Government contractors.
!
!---------------------Copyright Notices--------------------------------!
!
!            Copyright Battelle Memorial Institute, 1996
!                    All Rights Reserved.
!
!----------------------Description-------------------------------------!
!
!     Compute a time step that globally satisfies the Courant
!     number limit.
!
!----------------------Authors-----------------------------------------!
!
!     Written by MD White, Battelle, PNNL, October, 2000.
!     Last Modified by MD White, Battelle, PNNL, October 17, 2000.
!     $Id: crntnb.F 1080 2017-03-14 16:22:02Z d3c002 $

!

!----------------------Fortran 90 Modules------------------------------!
!
      USE GLB_PAR
      USE TRNSPT
      USE SOLTN
      USE GRID
      USE FDVP
      USE CONST
!

!----------------------Implicit Double Precision-----------------------!
!
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
!
!----------------------Parameter Statements----------------------------!
!

!
!----------------------Common Blocks-----------------------------------!
!

!
!----------------------Executable Lines--------------------------------!
!
      ISUB_LOG = ISUB_LOG+1
      SUB_LOG(ISUB_LOG) = '/CRN_LIM'
      IF( INDEX(SVN_ID(27)(1:1),'$').EQ.0 ) SVN_ID(27) =
     & '$Id: crntnb.F 1080 2017-03-14 16:22:02Z d3c002 $'
!
!---  Compute the maximum Courant number ---
!
      CRNLX = 0.D+0
      CRNNX = 0.D+0
      CRNGX = 0.D+0
      IF( ISLC(17).EQ.1 ) THEN
        IF( IEQW.GT.0 ) THEN
          DO 100 N = 1,NFLD
            IF( IXP(N).EQ.0 ) GOTO 100

            CRNLX = MAX( CRNTL(N),CRNLX )
  100     CONTINUE
        ENDIF
        IF( IEQA.GT.0 ) THEN
          DO 110 N = 1,NFLD
            IF( IXP(N).EQ.0 ) GOTO 110

            CRNGX = MAX( CRNTG(N),CRNGX )
  110     CONTINUE
        ENDIF
        IF( IEQO.GT.0 ) THEN
          DO 120 N = 1,NFLD
            IF( IXP(N).EQ.0 ) GOTO 120

            CRNNX = MAX( CRNTN(N),CRNNX )
  120     CONTINUE
        ENDIF
      ELSEIF( ISLC(17).EQ.2 ) THEN
        IF( IEQW.GT.0 ) THEN
          DO 210 N = 1,NFLD
            IF( IXP(N).EQ.0 ) GOTO 210
            IF( 1.D+0-SL(2,N).LT.EPSL ) GOTO 210
            CLX = C(N,NSL)*YL(N,NSL)/(SL(2,N)*PORD(2,N)+EPSL)
            IF( CLX.GT.CCL_CRN(NSL) ) CRNLX = MAX( CRNTL(N),CRNLX )
  210     CONTINUE
        ENDIF
      ELSEIF( ISLC(17).EQ.3 ) THEN
        IF( IEQW.GT.0 ) THEN
          DO 300 N = 1,NFLD
            IF( IXP(N).EQ.0 ) GOTO 300

            CRNLX = MAX( CRNTL(N),CRNLX )
  300     CONTINUE
        ENDIF
      ENDIF
      CRNMX = MAX( CRNLX,CRNGX,CRNNX )
      DT_CRN = DT
      DTI_CRN = DTI
      TM_CRN = TM
      TM = TM-DT
      IF( CRNMX.GT.CRNTMXT ) THEN
        N_CRN(NSL) = INT(CRNMX/CRNTMXT) + 1
        REALX = REAL(N_CRN(NSL))
        DT = DT/REALX
        DTI = 1.D+0/(DT+EPSL)
      ENDIF
      ISUB_LOG = ISUB_LOG-1
!
!---  End of CRN_LIM group  ---
!
      RETURN
      END

!----------------------Subroutine--------------------------------------!
!
      SUBROUTINE CRN_G
!
!-------------------------Disclaimer-----------------------------------!
!
!     This material was prepared as an account of work sponsored by
!     an agency of the United States Government. Neither the
!     United States Government nor the United States Department of
!     Energy, nor Battelle, nor any of their employees, makes any
!     warranty, express or implied, or assumes any legal liability or
!     responsibility for the accuracy, completeness, or usefulness
!     of any information, apparatus, product, software or process
!     disclosed, or represents that its use would not infringe
!     privately owned rights.
!
!----------------------Acknowledgement---------------------------------!
!
!     This software and its documentation were produced with Government
!     support under Contract Number DE-AC06-76RLO-1830 awarded by the
!     United Department of Energy. The Government retains a paid-up
!     non-exclusive, irrevocable worldwide license to reproduce,
!     prepare derivative works, perform publicly and display publicly
!     by or for the Government, including the right to distribute to
!     other Government contractors.
!
!---------------------Copyright Notices--------------------------------!
!
!            Copyright Battelle Memorial Institute, 1996
!                    All Rights Reserved.
!
!----------------------Description-------------------------------------!
!
!     Compute the maximum gas Courant number.
!
!----------------------Authors-----------------------------------------!
!
!     Written by MD White, PNNL, 01 Feburary 2007
!     Last Modified by MD White, PNNL, 01 Feburary 2007
!     $Id: crntnb.F 1080 2017-03-14 16:22:02Z d3c002 $
!

!----------------------Fortran 90 Modules------------------------------!
!
      USE GLB_PAR
      USE SOURC
      USE SOLTN
      USE OUTPU
      USE HYST
      USE GRID
      USE FLUXP
      USE FILES
      USE FDVP
      USE CONST
!

!----------------------Implicit Double Precision-----------------------!
!
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
!
!----------------------Parameter Statements----------------------------!
!

!
!----------------------Common Blocks-----------------------------------!
!

!
!----------------------Executable Lines--------------------------------!
!
      ISUB_LOG = ISUB_LOG+1
      SUB_LOG(ISUB_LOG) = '/CRN_G'
      IF( INDEX(SVN_ID(27)(1:1),'$').EQ.0 ) SVN_ID(27) =
     & '$Id: crntnb.F 1080 2017-03-14 16:22:02Z d3c002 $'
!
!---  Gas Courant number  ---
!
      CRNGMX = 0.D+0
      DO 100 N = 1,NFLD
        IF( IXP(N).EQ.0 ) GOTO 100
        NPX = NSX(N)
        NPY = NSY(N)
        NPZ = NSZ(N)
        NQX = NSX(N)+1
        IF( INBS(4,N).NE.0 ) NQX = INBS(4,N)
        NQY = NSY(N)+IFLD
        IF( INBS(5,N).NE.0 ) NQY = INBS(5,N)
        NQZ = NSZ(N)+IJFLD
        IF( INBS(6,N).NE.0 ) NQZ = INBS(6,N)
        I = ID(N)
        J = JD(N)
        K = KD(N)
        IF( (SG(2,N)-SGT(2,N))*PORD(2,N).GT.EPSL ) THEN
          CRNGX = MAX( ABS(UG(1,NPX)*DT/DXGP(NPX)),
     &    ABS(UG(1,NQX)*DT/DXGP(NQX)),
     &    ABS(VG(1,NPY)*DT/(DYGP(NPY)*RP(I))),
     &    ABS(VG(1,NQY)*DT/(DYGP(NQY)*RP(I))),
     &    ABS(WG(1,NPZ)*DT/DZGP(NPZ)),
     &    ABS(WG(1,NQZ)*DT/DZGP(NQZ)) )
     &    /((SG(2,N)-SGT(2,N))*PORD(2,N))
          IF( CRNGX.GT.CRNGMX ) THEN
            CRNGMX = CRNGX
            NX = N
          ENDIF
        ENDIF
  100 CONTINUE
!
!---  Reduce time step  ---
!
      IF( CRNGMX.GT.CRNTMXC ) THEN
        ICNV = 1
        DTX = DT
        DT = 9.D-1*DT*CRNTMXC/CRNGMX
        TM = TM - DTX + DT
        DTO = DT
        DTI = 1.D+0/DT
        VAR = DT
        VARX = DTX
        IF( UNTM.NE.'null' ) THEN
          INDX = 1
          IUNS = 1
          CALL RDUNIT(UNTM,VAR,INDX)
          IUNS = 1
          CALL RDUNIT(UNTM,VARX,INDX)
          NCH = INDEX( UNTM,'  ')-1
        ENDIF
        WRITE(ISC,'(A,I6)') 'Execution Note: Gas Courant ' //
     &    'Limit Exceeded: Node: ',NX
        WRITE(IWR,'(A,I6)') 'Execution Note: Gas Courant ' //
     &    'Limit Exceeded: Node: ',NX
        WRITE(ISC,'(A,1PE11.4,1X,2A,1PE11.4,1X,A)')
     &    '  Time step reduced from ',VARX,
     &    UNTM(1:NCH),' to ',VAR,UNTM(1:NCH)
        WRITE(IWR,'(A,1PE11.4,1X,2A,1PE11.4,1X,A)')
     &    '  Time step reduced from ',VARX,
     &    UNTM(1:NCH),' to ',VAR,UNTM(1:NCH)
        DO 200 N = 1,NFLD
          PL(2,N) = PL(1,N)
          PN(2,N) = PN(1,N)
          XMLO(2,N) = XMLO(1,N)
          SN(2,N) = SN(1,N)
          NPHAZ(2,N) = NPHAZ(1,N)
  200   CONTINUE

      ENDIF
      ISUB_LOG = ISUB_LOG-1
!
!---  End of CRN_G group  ---
!
      RETURN
      END

!----------------------Subroutine--------------------------------------!
!
      SUBROUTINE CRN_L
!
!-------------------------Disclaimer-----------------------------------!
!
!     This material was prepared as an account of work sponsored by
!     an agency of the United States Government. Neither the
!     United States Government nor the United States Department of
!     Energy, nor Battelle, nor any of their employees, makes any
!     warranty, express or implied, or assumes any legal liability or
!     responsibility for the accuracy, completeness, or usefulness
!     of any information, apparatus, product, software or process
!     disclosed, or represents that its use would not infringe
!     privately owned rights.
!
!----------------------Acknowledgement---------------------------------!
!
!     This software and its documentation were produced with Government
!     support under Contract Number DE-AC06-76RLO-1830 awarded by the
!     United Department of Energy. The Government retains a paid-up
!     non-exclusive, irrevocable worldwide license to reproduce,
!     prepare derivative works, perform publicly and display publicly
!     by or for the Government, including the right to distribute to
!     other Government contractors.
!
!---------------------Copyright Notices--------------------------------!
!
!            Copyright Battelle Memorial Institute, 1996
!                    All Rights Reserved.
!
!----------------------Description-------------------------------------!
!
!     Compute the maximum aqueous Courant number.
!
!----------------------Authors-----------------------------------------!
!
!     Written by MD White, PNNL, 01 Feburary 2007
!     Last Modified by MD White, PNNL, 01 Feburary 2007
!     $Id: crntnb.F 1080 2017-03-14 16:22:02Z d3c002 $
!

!----------------------Fortran 90 Modules------------------------------!
!
      USE GLB_PAR
      USE SOURC
      USE SOLTN
      USE OUTPU
      USE HYST
      USE GRID
      USE FLUXP
      USE FILES
      USE FDVP
      USE CONST
!

!----------------------Implicit Double Precision-----------------------!
!
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
!
!----------------------Parameter Statements----------------------------!
!

!
!----------------------Common Blocks-----------------------------------!
!

!
!----------------------Executable Lines--------------------------------!
!
      ISUB_LOG = ISUB_LOG+1
      SUB_LOG(ISUB_LOG) = '/CRN_L'
      IF( INDEX(SVN_ID(27)(1:1),'$').EQ.0 ) SVN_ID(27) =
     & '$Id: crntnb.F 1080 2017-03-14 16:22:02Z d3c002 $'
!
!---  Aqueous Courant number  ---
!
      CRNLMX = 0.D+0
      DO 100 N = 1,NFLD
        IF( IXP(N).EQ.0 ) GOTO 100
        NPX = NSX(N)
        NPY = NSY(N)
        NPZ = NSZ(N)
        NQX = NSX(N)+1
        IF( INBS(4,N).NE.0 ) NQX = INBS(4,N)
        NQY = NSY(N)+IFLD
        IF( INBS(5,N).NE.0 ) NQY = INBS(5,N)
        NQZ = NSZ(N)+IJFLD
        IF( INBS(6,N).NE.0 ) NQZ = INBS(6,N)
        I = ID(N)
        J = JD(N)
        K = KD(N)
        IF( SL(2,N)*PORD(2,N).GT.EPSL ) THEN
          CRNLX = MAX( ABS(UL(1,NPX)*DT/DXGP(NPX)),
     &    ABS(UL(1,NQX)*DT/DXGP(NQX)),
     &    ABS(VL(1,NPY)*DT/(DYGP(NPY)*RP(I))),
     &    ABS(VL(1,NQY)*DT/(DYGP(NQY)*RP(I))),
     &    ABS(WL(1,NPZ)*DT/DZGP(NPZ)),
     &    ABS(WL(1,NQZ)*DT/DZGP(NQZ)) )
     &    /(SL(2,N)*PORD(2,N))
          IF( CRNLX.GT.CRNLMX ) THEN
            CRNLMX = CRNLX
            NX = N
          ENDIF
        ENDIF
  100 CONTINUE
!
!---  Reduce time step  ---
!
      IF( CRNLMX.GT.CRNTMXC ) THEN
        ICNV = 1
        DTX = DT
        DT = 9.D-1*DT*CRNTMXC/CRNLMX
        TM = TM - DTX + DT
        DTO = DT
        DTI = 1.D+0/DT
        VAR = DT
        VARX = DTX
        IF( UNTM.NE.'null' ) THEN
          INDX = 1
          IUNS = 1
          CALL RDUNIT(UNTM,VAR,INDX)
          IUNS = 1
          CALL RDUNIT(UNTM,VARX,INDX)
          NCH = INDEX( UNTM,'  ')-1
        ENDIF
        WRITE(ISC,'(A,I6)') 'Execution Note: Aqueous Courant ' //
     &    'Limit Exceeded: Node: ',NX
        WRITE(IWR,'(A,I6)') 'Execution Note: Aqueous Courant ' //
     &    'Limit Exceeded: Node: ',NX
        WRITE(ISC,'(A,1PE11.4,1X,2A,1PE11.4,1X,A)')
     &    '  Time step reduced from ',VARX,
     &    UNTM(1:NCH),' to ',VAR,UNTM(1:NCH)
        WRITE(IWR,'(A,1PE11.4,1X,2A,1PE11.4,1X,A)')
     &    '  Time step reduced from ',VARX,
     &    UNTM(1:NCH),' to ',VAR,UNTM(1:NCH)
        DO 200 N = 1,NFLD
          PL(2,N) = PL(1,N)
          PN(2,N) = PN(1,N)
          XMLO(2,N) = XMLO(1,N)
          SN(2,N) = SN(1,N)
          NPHAZ(2,N) = NPHAZ(1,N)
  200   CONTINUE

      ENDIF
      ISUB_LOG = ISUB_LOG-1
!
!---  End of CRN_L group  ---
!
      RETURN
      END

!----------------------Subroutine--------------------------------------!
!
      SUBROUTINE CRN_N
!
!-------------------------Disclaimer-----------------------------------!
!
!     This material was prepared as an account of work sponsored by
!     an agency of the United States Government. Neither the
!     United States Government nor the United States Department of
!     Energy, nor Battelle, nor any of their employees, makes any
!     warranty, express or implied, or assumes any legal liability or
!     responsibility for the accuracy, completeness, or usefulness
!     of any information, apparatus, product, software or process
!     disclosed, or represents that its use would not infringe
!     privately owned rights.
!
!----------------------Acknowledgement---------------------------------!
!
!     This software and its documentation were produced with Government
!     support under Contract Number DE-AC06-76RLO-1830 awarded by the
!     United Department of Energy. The Government retains a paid-up
!     non-exclusive, irrevocable worldwide license to reproduce,
!     prepare derivative works, perform publicly and display publicly
!     by or for the Government, including the right to distribute to
!     other Government contractors.
!
!---------------------Copyright Notices--------------------------------!
!
!            Copyright Battelle Memorial Institute, 1996
!                    All Rights Reserved.
!
!----------------------Description-------------------------------------!
!
!     Compute the maximum NAPL Courant number.
!
!----------------------Authors-----------------------------------------!
!
!     Written by MD White, PNNL, 01 Feburary 2007
!     Last Modified by MD White, PNNL, 01 Feburary 2007
!     $Id: crntnb.F 1080 2017-03-14 16:22:02Z d3c002 $
!

!----------------------Fortran 90 Modules------------------------------!
!
      USE GLB_PAR
      USE SOURC
      USE SOLTN
      USE OUTPU
      USE HYST
      USE GRID
      USE FLUXN
      USE FILES
      USE FDVP
      USE CONST
!

!----------------------Implicit Double Precision-----------------------!
!
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
!
!----------------------Parameter Statements----------------------------!
!

!
!----------------------Common Blocks-----------------------------------!
!

!
!----------------------Executable Lines--------------------------------!
!
      ISUB_LOG = ISUB_LOG+1
      SUB_LOG(ISUB_LOG) = '/CRN_N'
      IF( INDEX(SVN_ID(27)(1:1),'$').EQ.0 ) SVN_ID(27) =
     & '$Id: crntnb.F 1080 2017-03-14 16:22:02Z d3c002 $'
!
!---  NAPL Courant number  ---
!
      CRNNMX = 0.D+0
      DO 100 N = 1,NFLD
        IF( IXP(N).EQ.0 ) GOTO 100
        NPX = NSX(N)
        NPY = NSY(N)
        NPZ = NSZ(N)
        NQX = NSX(N)+1
        IF( INBS(4,N).NE.0 ) NQX = INBS(4,N)
        NQY = NSY(N)+IFLD
        IF( INBS(5,N).NE.0 ) NQY = INBS(5,N)
        NQZ = NSZ(N)+IJFLD
        IF( INBS(6,N).NE.0 ) NQZ = INBS(6,N)
        I = ID(N)
        J = JD(N)
        K = KD(N)
        IF( (SN(2,N)-SNT(2,N))*PORD(2,N).GT.EPSL ) THEN
          CRNNX = MAX( ABS(UN(1,NPX)*DT/DXGP(NPX)),
     &    ABS(UN(1,NQX)*DT/DXGP(NQX)),
     &    ABS(VN(1,NPY)*DT/(DYGP(NPY)*RP(I))),
     &    ABS(VN(1,NQY)*DT/(DYGP(NQY)*RP(I))),
     &    ABS(WN(1,NPZ)*DT/DZGP(NPZ)),
     &    ABS(WN(1,NQZ)*DT/DZGP(NQZ)) )
     &    /((SN(2,N)-SNT(2,N))*PORD(2,N))
          IF( CRNNX.GT.CRNNMX ) THEN
            CRNNMX = CRNNX
            NX = N
          ENDIF
        ENDIF
  100 CONTINUE
!
!---  Reduce time step  ---
!
      IF( CRNNMX.GT.CRNTMXC ) THEN
        ICNV = 1
        DTX = DT
        DT = 9.D-1*DT*CRNTMXC/CRNNMX
        TM = TM - DTX + DT
        DTO = DT
        DTI = 1.D+0/DT
        VAR = DT
        VARX = DTX
        IF( UNTM.NE.'null' ) THEN
          INDX = 1
          IUNS = 1
          CALL RDUNIT(UNTM,VAR,INDX)
          IUNS = 1
          CALL RDUNIT(UNTM,VARX,INDX)
          NCH = INDEX( UNTM,'  ')-1
        ENDIF
        WRITE(ISC,'(A,I6)') 'Execution Note: NAPL Courant ' //
     &    'Limit Exceeded: Node: ',NX
        WRITE(IWR,'(A,I6)') 'Execution Note: NAPL Courant ' //
     &    'Limit Exceeded: Node: ',NX
        WRITE(ISC,'(A,1PE11.4,1X,2A,1PE11.4,1X,A)')
     &    '  Time step reduced from ',VARX,
     &    UNTM(1:NCH),' to ',VAR,UNTM(1:NCH)
        WRITE(IWR,'(A,1PE11.4,1X,2A,1PE11.4,1X,A)')
     &    '  Time step reduced from ',VARX,
     &    UNTM(1:NCH),' to ',VAR,UNTM(1:NCH)
        DO 200 N = 1,NFLD
          PL(2,N) = PL(1,N)
          PN(2,N) = PN(1,N)
          XMLO(2,N) = XMLO(1,N)
          SN(2,N) = SN(1,N)
          NPHAZ(2,N) = NPHAZ(1,N)
  200   CONTINUE

      ENDIF
      ISUB_LOG = ISUB_LOG-1
!
!---  End of CRN_N group  ---
!
      RETURN
      END

ext: f

max_stars_repo_name: pnnl/STOMP-W

max_issues_repo_name: pnnl/stomp-w

max_forks_repo_name: pnnl/stomp-w


text:
program ch1002
    ! Reading Real Data
    implicit none

    integer, parameter :: n = 10
    real, dimension(1:n) :: h, w, bmi
    integer :: i

    open(100,file='ch1001.out',status='old')
    do i = 1, n
        read(100,'(1x,f5.2,2x,f4.1)') h(i), w(i)
    end do
    close(100)

    bmi = w/(h*h)
    do i = 1, n
        write(*,'(1x, f4.1)') bmi(i)
    end do
end program

max_stars_repo_name: hechengda/fortran

max_issues_repo_name: hechengda/fortran

max_forks_repo_name: hechengda/fortran


text:
      SUBROUTINE G0DIVL(IAXIS)
C
C          ------------------------------------------------
C          ROUTINE NO. ( 236)   VERSION (A7.5)    11:FEB:85
C          ------------------------------------------------
C
C          THIS FINDS THE POSITIONS OF THE INTERVALS ON THE GIVEN AXIS
C          (FOR LOG. SCALING) WHICH LIE WITHIN THE CURRENT WINDOW AREA.
C
C
C          <IAXIS>  GIVES THE REQUIRED AXIS BY ITS MODULUS:
C                   = 1, THE X-AXIS IS TAKEN, OR
C                   = 2, THE Y-AXIS IS TAKEN.
C
C
C          THE FOLLOWING ARGUMENTS ARE SUPPLIED THROUGH COMMON:
C
C          <X1WND0>    THE COORDINATES
C          <X2WND0>    OF
C          <Y1WND0>    THE
C          <Y2WND0>    WINDOW RECTANGLE.
C
C          <AXPOSX> THE POSITION OF THE X-AXIS ALONG Y.
C          <AXPOSY> THE POSITION OF THE Y-AXIS ALONG X.
C          <NOTATE> IF ZERO, NO ANNOTATION IS REQUIRED,
C                   IF NON-ZERO, ANNOT. FORMAT MUST BE CALCULATED.
C
C
C          THE FOLLOWING ARGUMENTS ARE RETURNED THROUGH COMMON:
C          (ONLY THE ARGS. RELEVANT TO THE GIVEN AXIS ARE CHANGED):
C
C          <KTYPEX> IS THE X-AXIS TYPE (= 2 FOR LOG. AXIS).
C          <KTYPEY> IS THE Y-AXIS TYPE (= 2 FOR LOG. AXIS).
C          <DIVLX>  GIVES THE X-AXIS SIGN.
C          <DIVLY>  GIVES THE Y-AXIS SIGN.
C          <NSKIPX> IS THE X-AXIS START-POINT EXPONENT.
C          <NSKIPY> IS THE Y-AXIS START-POINT EXPONENT.
C          <NTIKLX> IS THE X-AXIS START-POINT MANTISSA.
C          <NTIKLY> IS THE Y-AXIS START-POINT MANTISSA.
C          <NTIKHX> IS THE NO. OF X-AXIS TICK MARKS.
C          <NTIKHY> IS THE NO. OF Y-AXIS TICK MARKS.
C          <NDECSX> IS THE NO. OF DECADE MARKS ON THE X-AXIS.
C          <NDECSY> IS THE NO. OF DECADE MARKS ON THE Y-AXIS.
C          <NCHRSX> IS THE NO. OF CHARS. IN X-AXIS ANNOTATION.
C          <NCHRSY> IS THE NO. OF CHARS. IN Y-AXIS ANNOTATION.
C          <KAXIS>  IS SET BY <IAXIS> FOR SUBSEQUENT USE.
C
C
      LOGICAL ERRON
C
      COMMON /T0AARG/ KAXIS
      COMMON /T0ADIX/ DIVLX,NTIKLX,NTIKHX
      COMMON /T0ADIY/ DIVLY,NTIKLY,NTIKHY
      COMMON /T0ANOX/ KANNX,NCHRSX,NAFTPX
      COMMON /T0ANOY/ KANNY,NCHRSY,NAFTPY
      COMMON /T0ASKX/ NSKIPX,NDECSX
      COMMON /T0ASKY/ NSKIPY,NDECSY
      COMMON /T0ATYP/ KTYPEX,KTYPEY
      COMMON /T0WNDO/ X1WND0,X2WND0,Y1WND0,Y2WND0
      COMMON /T3ERRS/ ERRON,NUMERR
C
C
C          THE AXIS TYPE IS SET, THE SIGN OF THE AXIS REGION IS
C          FOUND, AND THE END POINTS ARE APPROPRIATELY SET.
C
      KAXIS= IAXIS
C
      ENDMIN= AMIN1(X1WND0,X2WND0)
      ENDMAX= AMAX1(X1WND0,X2WND0)
      IF (IABS(KAXIS).EQ.1) GO TO 1
      ENDMIN= AMIN1(Y1WND0,Y2WND0)
      ENDMAX= AMAX1(Y1WND0,Y2WND0)
C
    1 AXSIGN= +1.0
      IF (ENDMIN.GT.0.0) GO TO 2
      AXSIGN= -1.0
      SWOP=   ABS(ENDMIN)
      ENDMIN= ABS(ENDMAX)
      ENDMAX= SWOP
C
C          THE START-POINT EXPONENT AND MANTISSA ARE FOUND.
C          THE SECOND SECTION ALLOWS FOR ROUNDING ERRORS.
C
    2 EXP= ALOG10(ENDMIN)
      IF (EXP.LT.0.0) EXP= EXP-1.0
      IDEC1= EXP
      TIKPT= ENDMIN*(10.0**(-IDEC1))
      ITIK1= TIKPT+0.99999
C
      IF (ITIK1.LT.10) GO TO 3
      ITIK1= 1
      IDEC1= IDEC1+1
      GO TO 4
    3 IF (ITIK1.GE.1) GO TO 4
      ITIK1= 9
      IDEC1= IDEC1-1
C
C          THE SAME IS THEN DONE FOR THE END-POINT VALUES.
C
    4 EXP= ALOG10(ENDMAX)
      IF (EXP.LT.0.0) EXP= EXP-1.0
      IDEC2= EXP
      TIKPT= ENDMAX*(10.0**(-IDEC2))
      ITIK2= TIKPT+0.00001
C
      IF (ITIK2.LT.10) GO TO 5
      ITIK2= 1
      IDEC2= IDEC2+1
      GO TO 6
    5 IF (ITIK2.GE.1) GO TO 6
      ITIK2= 9
      IDEC2= IDEC2-1
C
C          THE NO. OF TICK AND DECADE MARKS CAN THEN BE CALCULATED.
C
    6 NDECS= IDEC2-IDEC1
      IF (NDECS.GT.20) GO TO 901
      NTICKS= 9*NDECS+(ITIK2-ITIK1)+1
      IF (ITIK1.EQ.1) NDECS= NDECS+1
C
C          THE NO. OF CHARS. IN THE ANNOTATION FORMAT IS FOUND
C          WHEN THE VALUES WILL BE PRINTED IN NON-EXPONENT FORM.
C
      NCHARS= 2
      IF (NTICKS.GT.10)   GO TO 7
      NCHARS= IDEC2+2
      IF (IDEC1.GE.0)     GO TO 7
      NCHARS= 3-IDEC1
      IF (IDEC1.EQ.IDEC2) GO TO 7
      NCHARS= NCHARS+1
C
C          ONCE THE VALUES HAVE BEEN FOUND, THEY ARE PLACED
C          INTO THE APPROPRIATE VARIABLES IN THE COMMON BLOCK.
C
    7 IF (IABS(KAXIS).EQ.2) GO TO 8
      KTYPEX= 2
      NTIKLX= ITIK1
      NSKIPX= IDEC1
      NTIKHX= NTICKS
      NDECSX= NDECS
      DIVLX= AXSIGN
      NCHRSX= NCHARS
      RETURN
C
    8 KTYPEY= 2
      NTIKLY= ITIK1
      NSKIPY= IDEC1
      NTIKHY= NTICKS
      NDECSY= NDECS
      DIVLY= AXSIGN
      NCHRSY= NCHARS
      RETURN
C
C          IF MORE THAN 20 CYCLES IS NECESSARY, AN ERROR
C          MESSAGE IS PRINTED AND THE MARKING IS SUPPRESSED.
C
  901 NUMERR= 32
      IF (ERRON) CALL G0ERMS
      IF (IABS(KAXIS).EQ.1) NTIKHX= 0
      IF (IABS(KAXIS).EQ.2) NTIKHY= 0
      RETURN
      END

max_stars_repo_name: ZedThree/GHOST

max_issues_repo_name: ZedThree/GHOST

max_forks_repo_name: ZedThree/GHOST