psmile_mg_next_level_2d_real.F90

Go to the documentation of this file.

!-----------------------------------------------------------------------
! Copyright 2006-2010, NEC Europe Ltd., London, UK.
! All rights reserved. Use is subject to OASIS4 license terms.
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: PSMILe_MG_Next_level_2d_real
!
! !INTERFACE:

      subroutine psmile_mg_next_level_2d_real (            &
                           grid_id, lev, nlev,             &
                           chmin1, chmin2,                 &
                           chmax1, chmax2,                 &
                           midp1,  midp2,                  &
                           levdim,                         &
                           found, loc, range,              &
                           coords1, coords2,               &
                           shape, control,                 &
                           ijkinc, ijkcoa, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_mg_next_level_2d_real
#ifdef DEBUG_TRACE
      use psmile_debug_trace
#endif

      implicit none
!
! !INPUT PARAMETERS:
!
      Integer, Intent (In)            :: grid_id

!     Info on the component in which the donor
!     should be searched.

      Integer, Intent (In)            :: lev

!     Level number of the next grid

      Integer, Intent (In)            :: nlev

!     Total number of levels

      Integer, Intent (In)            :: levdim (ndim_2d)

!     Dimension of chmin, chmax and midp

      Integer, Intent (In)            :: range (2, ndim_3d)

!     Dimension of loc and found

      Integer, Intent (In)            :: shape (2, ndim_3d)

!     Dimension of coordinate arrays

      Real, Intent (In)               :: chmin1 (0:levdim(1), 
                                                 0:levdim(2))
      Real, Intent (In)               :: chmin2 (0:levdim(1), 
                                                 0:levdim(2))

!     Minimum of grid coordinates per grid cell

      Real, Intent (In)               :: chmax1 (0:levdim(1), 
                                                 0:levdim(2))
      Real, Intent (In)               :: chmax2 (0:levdim(1), 
                                                 0:levdim(2))

!     Maximum of grid coordinates per grid cell

      Real, Intent (In)               ::  midp1 (0:levdim(1), 
                                                 0:levdim(2))
      Real, Intent (In)               ::  midp2 (0:levdim(1), 
                                                 0:levdim(2))
!     Mid point of the cell


      Real, Intent (In)               :: coords1 (shape(1,1):shape(2,1), 
                                                  shape(1,2):shape(2,2), 
                                                  shape(1,3):shape(2,3))
      Real, Intent (In)               :: coords2 (shape(1,1):shape(2,1), 
                                                  shape(1,2):shape(2,2), 
                                                  shape(1,3):shape(2,3))

!     Coordinates to be searched

      Integer, Intent (In)            :: ijkinc (ndim_3d)

!     Increment for each coordinate direction to be searched

      Integer, Intent (In)            :: ijkcoa (ndim_2d)

!     Coarsening of source grid
!
! !INPUT/OUTPUT PARAMETERS:
!
      Integer, Intent (InOut)         :: found (range(1,1):range(2,1), 
                                                range(1,2):range(2,2), 
                                                range(1,3):range(2,3))

!     Finest level number on which a grid cell was found for point i,j.
!     Level number < -1: Point was not found and
!                        and last level number was (-found(i,j,k))
!     Level number = (-nlev+1): Never found (input value)

      Integer, Intent (InOut)         :: loc   (ndim_2d,               
                                                range(1,1):range(2,1), 
                                                range(1,2):range(2,2), 
                                                range(1,3):range(2,3))

!     Indices of the grid cell in which the point was found.

      Integer, Intent (InOut)         :: control (2, ndim_3d)

!     Input  value: index range in "coords" to be searched
!     Output value: index range in "coords" for which a cell was found
!
! !OUTPUT PARAMETERS:
!
      integer, Intent (Out)           :: ierror

!     Returns the error code of PSMILe_mg_next_level_2d_real;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !DEFINED PARAMETERS:
!
      Integer, parameter              :: ndtry = 16
      Integer, parameter              :: ntry1 = 4
      Integer, parameter              :: ntry3 = 9
      Integer, parameter              :: ntry2 = ndtry - ntry1
!
      Integer, Parameter              :: temp_len = 256
!
      Real, parameter                 :: remax = 1.0e20
!
! !LOCAL VARIABLES
!
!     ... for levels

      Integer                         :: levc

!     ... for locations searched

      Integer                         :: i, j, k
      Integer                         :: ibegl, ibeglx, iendl, n, nprev
      Integer                         :: nmin (temp_len, 1)
      Integer                         :: ijkf(temp_len, ndim_2d)
      Integer                         :: ijkold(temp_len, ndim_2d)
      Integer                         :: newijk(ndim_2d, temp_len)
      Integer                         :: ijkctl (ndim_2d, ndtry)
      Integer                         :: ijkdst (temp_len, ndim_2d, ndtry)
      Integer                         :: ijkctl3 (ndim_2d, ndtry)
      Real                            :: dist (temp_len, ndtry)
      Real                            :: xyz (ndim_2d, temp_len)
!
!     ...
!
      Integer                         :: ii, nc, np, nc_part
      Integer                         :: list (temp_len*2)
      Integer                         :: part_list (temp_len)
      Integer                         :: loc2 (ndim_2d, temp_len)
!
#ifdef DEBUG
!     ... for locations searched
!
      Integer                         :: nfnd (0:4)
#endif /* DEBUG */
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_mg_next_level_2d_real" searches the donor cells
! on the next MG grid for the subgrid coords by the sending process.
!
!
! !REVISION HISTORY:
!   Date      Programmer   Description
! ----------  ----------   -----------
! 03.07.21    H. Ritzdorf  created
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_mg_next_level_2d_real.F90 2927 2011-01-28 14:04:12Z hanke $
!  $Author: hanke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_mg_next_level_2d_real.F90 2927 2011-01-28 14:04:12Z hanke $'
!
!----------------------------------------------------------------------
!
!  Relative control indices for sub-cells of coarse grid cells
!
      data ((ijkctl (i, n), i=1,ndim_2d), n = 1,ntry1) &
                             / 0, 0,   1, 0,    &
                               0, 1,   1, 1/
!
!  Relative control indices for neighbour cells
!
      data ((ijkctl (i, n), i=1,ndim_2d), n = ntry1+1, ntry1+12) &
                  /-1,-1,   0,-1,   1,-1,   2,-1,    &
                   -1, 0,                   2, 0,    &
                   -1, 1,                   2, 1,    &
                   -1, 2,   0, 2,   1, 2,   2, 2/
!
!  Relative control indices for computed cells
!
      data ((ijkctl3 (i, n), i=1,ndim_2d), n = 1,ntry3) &
                 / -1,-1,   0,-1,   1,-1,     &
                   -1, 0,   0, 0,   1, 0,     &
                   -1, 1,   0, 1,   1, 1/
!
!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
      print 9990, trim(ch_id), lev, control

      call psmile_flushstd
#endif /* VERBOSE */
!
#ifdef PRISM_ASSERTION
#endif
!
      ierror = 0
      levc = lev + 1
!
      iendl = control(2,1) + ijkinc(1)
!
#ifdef DEBUG
      nfnd (:) = 0
#endif /* DEBUG */

         do k = control(1,3), control (2,3), ijkinc (3)
loop_j:     do j = control(1,2), control (2,2), ijkinc (2)
!
!-----------------------------------------------------------------------
!     Look in the fine parts of the coarse grid cell loc(:,i,j,k)
!
!     TODO: List "uber 1-dimensionalen Vector (i,j,k)
!-----------------------------------------------------------------------
!
            nprev = 0

            ibegl = control (1, 1)
!
!===> ... Look for the next cells to be controlled
!         Store indices of cells in list "list"
!
loop_i:        do while (ibegl < iendl)
!
               nc = 0
               i = ibegl
!
                  do while (nc < temp_len .and. i < iendl)
                     ibegl = i
!cdir vector
                        do i = ibegl, min (ibegl+((temp_len*2-1)-nc)*ijkinc(1), &
                                           control(2,1)), ijkinc(1)
                        if (found(i, j, k) == levc) then
                           nc = nc + 1
                           list (nc) = i
                        endif
                        end do
!
                     if (nc < temp_len .and. &
                         i > control(2,1) .and. i < iendl) then
                        i = control(2,1)
                        if (found (i,j,k) == levc) then
                           nc = nc + 1
                           list (nc) = i
!
                           exit ! do while (nc < temp_len .and.
                        endif
                     endif
                  end do
!
               if (nc == 0) exit loop_i ! exit while loop (ibegl < ...)
!
               if (nc < temp_len) then
                  ibegl = iendl
               else if (nc > temp_len) then
                  nc = temp_len
                  ibegl = list (temp_len+1)
               else
                  ibegl = list (temp_len) + ijkinc(1)
               endif
!
!===> ... Look in the fine parts of cell loc(:,i,j,k)
!
               ijkf(1:nc, 1) = min (loc(1,list(1:nc),j,k) * ijkcoa(1), levdim(1))
               ijkf(1:nc, 2) = min (loc(2,list(1:nc),j,k) * ijkcoa(2), levdim(2))
!
                  do n = 1, ntry1
                     do ii = 1, ndim_2d
                     ijkdst (1:nc, ii, n) = max (min(ijkf(1:nc, ii)+ijkctl(ii,n), &
                                                     levdim(ii)), 0)
                     end do ! ii
                  end do ! n
!
!cdir vector
                  do n = 1, ntry1
                     do i = 1, nc
                     if (coords1(list(i),j,k) >= chmin1(ijkdst(i,1,n),         &
                                                        ijkdst(i,2,n)) .and.   &
                         coords2(list(i),j,k) >= chmin2(ijkdst(i,1,n),         &
                                                        ijkdst(i,2,n)) .and.   &
                         coords1(list(i),j,k) <= chmax1(ijkdst(i,1,n),         &
                                                        ijkdst(i,2,n)) .and.   &
                         coords2(list(i),j,k) <= chmax2(ijkdst(i,1,n),         &
                                                        ijkdst(i,2,n))) then
                        dist (i, n) = (coords1(list(i),j,k) -                  &
                                       midp1(ijkdst(i,1,n), ijkdst(i,2,n)))**2 &
                                    + (coords2(list(i),j,k) -                  &
                                       midp2(ijkdst(i,1,n), ijkdst(i,2,n)))**2

                     else

                        dist (i, n) = remax

                     endif
                     end do ! i
                  end do ! n
!
!===> ... Look for the minimum distance
!
#if 0
                  do i = 1, nc
                  nmin (i, :) = MINLOC (dist(i, 1:ntry1))
                  end do
#ifdef MINLOCFIX
                  do i = 1, nc
                  if (nmin(i, 1).eq.0) nmin (i, 1) = 1
                  end do
#endif /* MINLOCFIX */
#else
               nmin (1:nc, 1) = 1
!
                  do n = 2, ntry1
!cdir vector
                     do i = 1, nc
                     if (dist(i, n) < dist(i, nmin(i,1))) then
                        nmin (i,1) = n
                     endif
                     end do
                  end do
#endif 
!
!     ... Look for points which were found in a fine grid cell
!
               nc_part = 0
!cdir vector
                  do i = 1, nc
                  if (dist(i, nmin(i, 1)) /= remax) then
                     nc_part = nc_part + 1
                     part_list (nc_part) = i
                  endif
                  end do
!
               if (nc_part > 0) then
!cdir vector
                     do np = 1, nc_part
                     i = part_list (np)

                     found  (list(i),j,k) = lev
                     loc (:, list(i),j,k) = ijkdst (i, :, nmin(i, 1))
                     end do

#ifdef DEBUG
                  nfnd (1) = nfnd (1) + nc_part
#endif /* DEBUG */
!
                  if (nc == nc_part) cycle loop_i
!
!     ... Get list of points which require further search
!
                  nc_part = 0
!cdir vector
                     do i = 1, nc
                     if (dist(i, nmin(i, 1)) == remax) then
                        nc_part = nc_part + 1
                        part_list (nc_part) = i
                     endif
                     end do
!
!===> ...... Get data on remaining points out of previous list
!
                  nc = nc_part
!cdir vector
                     do i = 1, nc
! Hilfsvector einfuehren
                     list (i) = list (part_list (i))
!
                     ijkf (i,1) = ijkf (part_list(i), 1)
                     ijkf (i,2) = ijkf (part_list(i), 2)
                     enddo
!
               endif
!
!     ... Is the coarse grid the coarsest grid ?
!
               if (levc == lev) then
!cdir vector
                     do i = 1, nc
                     found (list(i),j,k) = - found(list(i),j,k)
                     end do
#ifdef DEBUG
                  nfnd (0) = nfnd (0) + nc
#endif /* DEBUG */

                  cycle loop_i
               endif
!
!-----------------------------------------------------------------------
!          Look in the neigbourhood of the coarse grid cell
!          loc(:,list(i),j,k)
!-----------------------------------------------------------------------
!
                 do n = ntry1+1, ndtry
!cdir vector
                    do ii = 1, ndim_2d
                    ijkdst (1:nc, ii, n) = max(0, min(ijkf(1:nc, ii)+ijkctl(ii,n), &
                                                      levdim(ii)))
                    end do
!
                    do i = 1, nc
                    if (coords1(list(i),j,k) >= chmin1(ijkdst(i, 1,n),         &
                                                       ijkdst(i, 2,n)) .and.   &
                        coords2(list(i),j,k) >= chmin2(ijkdst(i, 1,n),         &
                                                       ijkdst(i, 2,n)) .and.   &
                        coords1(list(i),j,k) <= chmax1(ijkdst(i, 1,n),         &
                                                       ijkdst(i, 2,n)) .and.   &
                        coords2(list(i),j,k) <= chmax2(ijkdst(i, 1,n),         &
                                                       ijkdst(i, 2,n))) then

                    dist (i, n) = (coords1(list(i),j,k) -                    &
                                   midp1(ijkdst(i, 1,n), ijkdst(i, 2,n)))**2 &
                                + (coords2(list(i),j,k) -                    &
                                   midp2(ijkdst(i, 1,n), ijkdst(i, 2,n)))**2

                 else

                    dist (i, n) = remax

                 endif
                 end do ! i
              end do ! n
!
!===> ... Look for the minimum distance
!
#if 0
              do i = 1, nc
              nmin(i,:) = MINLOC (dist(i, ntry1+1:ndtry))
              end do
!
#ifdef MINLOCFIX
              do i = 1, nc
              if (nmin(i, 1).eq.0) nmin (i, 1) = 1
              end do

#endif /* MINLOCFIX */
#else
               nmin (1:nc, 1) = ntry1+1
!
                  do n = ntry1+2, ndtry
!cdir vector
                     do i = 1, nc
                     if (dist(i, n) < dist(i, nmin(i,1))) then
                        nmin (i,1) = n
                     endif
                     end do
                  end do
#endif
!
!     ... Look for points which were found in a fine grid cell
!
            nc_part = 0
!cdir vector
               do i = 1, nc
               if (dist(i, nmin(i, 1)) /= remax) then
                  found  (list(i),j,k) = lev
                  loc (:, list(i),j,k) = ijkdst (i, :, nmin(i, 1))
!
                  nc_part = nc_part + 1
               endif
               end do
!
#ifdef DEBUG
            nfnd (2) = nfnd (2) + nc_part
#endif /* DEBUG */
!
            nprev = nprev + (nc - nc_part)
!
         end do loop_i ! while
!
!===> ... Have points to be controlled o a coarser level ?
!
         if (nprev == 0) cycle loop_j
!
!-----------------------------------------------------------------------
!          Some points were not found.
!          Try to interpolate the fine point from neighbouring points
!-----------------------------------------------------------------------
!
            ibeglx = control (1, 1) +  ijkinc(1)
!
!===> ... Look for the next cells to be controlled
!         Store indices of cells in list "list"
!
loop_2:        do while (ibeglx < iendl)
!
               nc = 0
               i = ibeglx
!
                  do while (nc < temp_len .and. i < control(2,1))
                     ibeglx = i
!cdir vector
                        do i = ibeglx, &
                               min (ibeglx+((temp_len*2-1)-nc)*ijkinc(1), &
                                    control(2,1)-ijkinc(1)), ijkinc(1)
                        if (found (i, j, k) == levc .and.        &
                            found (i-ijkinc(1),j,k) == lev .and. &
                            found (i+ijkinc(1),j,k) == lev) then
                           nc = nc + 1
                           list (nc) = i
                        endif
                        end do
!
                     if (nc < temp_len .and. &
                         i > control(2,1)-ijkinc(1) .and. i < control(2,1)) then
                        i = control(2,1) - ijkinc(1)
                        if (found (i, j, k) == levc .and.        &
                            found (i-ijkinc(1),j,k) == lev .and. &
                            found (i+ijkinc(1),j,k) == lev) then
                           nc = nc + 1
                           list (nc) = i
!
                           exit ! while (nc < temp_len .and. ...)
                        endif
                     endif
                  end do
!
               if (nc == 0) exit loop_2 ! exit while loop (ibeglx < ...)
!
               if (nc < temp_len) then
                  ibeglx = control(2,1) + ijkinc (1)
               else if (nc > temp_len) then
                  nc = temp_len
                  ibeglx = list (temp_len+1)
               else
                  ibeglx = list (temp_len) + ijkinc(1)
               endif
!
              do ii = 1, ndim_2d
              ijkf (1:nc, ii) = min ((loc (ii, list(1:nc)-ijkinc(1),j,k) +   &
                                      loc (ii, list(1:nc)+ijkinc(1),j,k))/2, &
                                     levdim(ii))
!
              ijkold (1:nc, ii) = min (loc(ii, list(1:nc),j,k) * ijkcoa(ii), &
                                       levdim(ii))
              end do
!
!===> ... Interpolate cell IJKF and
!         look in the neighbourhood of cell IJKF
!
#ifdef HUHU
              do i = 1, nc
!
!! rausnehmen aus der liste gar nicht erzeugen ? oder ignorieren
                 if (ijkold(i, 1) == ijkf(i, 1) .and. &
                     ijkold(i, 2) == ijkf(i, 2)) cycle
!
              end do
#endif

              do n = 1, ntry3
!cdir vector
                 do ii = 1, ndim_2d
                 ijkdst (1:nc, ii, n) = max(0, min(ijkf(1:nc, ii) + ijkctl3(ii,n), &
                                                   levdim(ii)))
                 end do
!
!===> ... Compute distances to the cell midpoints
!
!cdir vector
                 do i = 1, nc
                 if (coords1(list(i),j,k) >= chmin1(ijkdst(i,1,n),         &
                                                    ijkdst(i,2,n)) .and.   &
                     coords2(list(i),j,k) >= chmin2(ijkdst(i,1,n),         &
                                                    ijkdst(i,2,n)) .and.   &
                     coords1(list(i),j,k) <= chmax1(ijkdst(i,1,n),         &
                                                    ijkdst(i,2,n)) .and.   &
                     coords2(list(i),j,k) <= chmax2(ijkdst(i,1,n),         &
                                                    ijkdst(i,2,n))) then
                    dist (i, n) = (coords1(list(i),j,k) - midp1(ijkdst(i,1,n),     &
                                                                ijkdst(i,2,n)))**2 &
                                + (coords2(list(i),j,k) - midp2(ijkdst(i,1,n),     &
                                                                ijkdst(i,2,n)))**2

                 else

                    dist (i, n) = remax

                 endif
                 end do
              end do ! n
!
!===> ... Look for the minimum distance
!
#if 0
              do i = 1, nc
              nmin (i, :) = MINLOC (dist(i, 1:ntry3))
              end do
!
#ifdef MINLOCFIX
              do i = 1, nc
              if (nmin(i, 1).eq.0) nmin (i, 1) = 1
              end do
#endif /* MINLOCFIX */
#else
              nmin (1:nc, 1) = 1
!
                 do n = 2, ntry3
!cdir vector
                    do i = 1, nc
                    if (dist(i, n) < dist(i, nmin(i,1))) then
                       nmin (i,1) = n
                    endif
                    end do
                 end do
#endif
!
              nc_part = 0
!cdir vector
                 do i = 1, nc
                 if (dist(i, nmin(i, 1)) /= remax) then
                    found ( list(i),j,k) = lev
                    loc (:, list(i),j,k) = ijkdst (i, :, nmin(i, 1))
   
                    nc_part = nc_part + 1
                 endif
                 end do
!
                 nprev = nprev - nc_part
#ifdef DEBUG
              nfnd (3) = nfnd (3) + nc_part
#endif /* DEBUG */
!
              end do loop_2 ! while
!
           if (nprev == 0) cycle loop_j

!-----------------------------------------------------------------------
!     Some points were still not found.
!     Look again on the next coarser level whether the points
!     are located in other cells.
!-----------------------------------------------------------------------
!
           if (lev+1 < nlev) then
!
         ibeglx = control (1, 1)
!
!===> ... Look for the next cells to be controlled
!         Store indices of cells in list "list"
!
loop_3:        do while (ibeglx < iendl)
!
               nc = 0
               i = ibeglx
!
                  do while (nc < temp_len .and. i < iendl)
                     ibeglx = i
!cdir vector
                        do i = ibeglx, min (ibeglx+((temp_len*2-1)-nc)*ijkinc(1), &
                                           control(2,1)), ijkinc(1)
                        if (found(i, j, k) == levc) then
                           nc = nc + 1
                           list (nc) = i
                        endif
                        end do
!
                     if (nc < temp_len .and. &
                         i > control(2,1) .and. i < iendl) then
                        i = control(2,1)
                        if (found (i,j,k) == levc) then
                           nc = nc + 1
                           list (nc) = i
                        endif
                     endif
                  end do
!
               if (nc == 0) exit loop_3 ! exit while loop (ibeglx < ...)
!
               if (nc < temp_len) then
                  ibeglx = control(2,1) + ijkinc (1)
               else if (nc > temp_len) then
                  nc = temp_len
                  ibeglx = list (temp_len+1)
               else
                  ibeglx = list (temp_len) + ijkinc(1)
               endif
!
!===> ... Look in the coarser parts of cell loc(:,list(i),j,k)
!     !!! TODO: collect xyz and reduce number of subroutine calls
!               to PSMILE_mg_prev_levels_2d_real
! 3d loop mit nprev ?
!cdir vector
                  do i = 1, nc
                  xyz (1, i) = coords1 (list(i),j,k)
                  xyz (2, i) = coords2 (list(i),j,k)
!
                  loc2 (1, i) = loc(1,list(i),j,k)
                  loc2 (2, i) = loc(2,list(i),j,k)
                  end do
!
               call psmile_mg_prev_levels_2d_real (grid_id, lev, nlev, &
                                 loc2, xyz, part_list, &
                                 newijk, nc)
!
               do i = 1, nc
               if (part_list(i) > 0) then
!
!                 Point (list(i), j, k) was found in cell newijk (*, i)
!
                  loc (:, list(i),j,k) = newijk (:, i)
!
                  found (list(i),j,k) = lev
#ifdef DEBUG
                  nfnd (4) = nfnd (4) + 1
#endif /* DEBUG */
               else
                  found (list(i),j,k) = - found (list(i),j,k)
!                 print *, 'not found a', list(i),j,k, coords1 (list(i),j,k), coords2 (list(i),j,k)
#ifdef DEBUG
                  nfnd (0) = nfnd (0) + 1
#endif /* DEBUG */
               endif
               end do
!
            end do loop_3 ! while
            else
!              This is the second coarsest level.
!              The coarsest level contains only 1 cell.
!              All points which are not found are outside.

!cdir vector
               do i = control (1,1), control (2,1)
               if (found(i, j, k) == levc) then
                  found(i, j, k) = -levc
!                 print *, 'not found b', i,j,k
#ifdef DEBUG
                  nfnd (0) = nfnd (0) + 1
#endif /* DEBUG */
               endif
               end do
            endif ! if (lev+1 < nlev)
!
         end do loop_j ! j
      end do ! k
!
!===> All done
!
#ifdef DEBUG_TRACE
      if (range(1,1) <= ictl_ind(1) .and. range(2,1) >= ictl_ind(1) .and. &
          range(1,2) <= ictl_ind(2) .and. range(2,2) >= ictl_ind(2)) then
         i = ictl_ind(1)
         j = ictl_ind(2)
         k = control(1,3)

         if ( found (i,j,k) == lev .or. found (i,j,k) == -levc) then
            print 8990, ictl_ind (1:ndim_2d), &
                     found (i,j,k), loc(:, i,j,k)
!
            if (found (i,j,k) == lev) then
               print *, 'chmin ', &
                        chmin1 (loc(1,i,j,k), loc(2,i,j,k)), &
                        chmin2 (loc(1,i,j,k), loc(2,i,j,k))
               print *, 'coords', coords1 (i,j,k), coords2 (i,j,k)
               print *, 'chmax ', &
                        chmax1 (loc(1,i,j,k), loc(2,i,j,k)), &
                        chmax2 (loc(1,i,j,k), loc(2,i,j,k))
            endif

8990  format (1x, '### mg_next ictl_ind', 2i5, '; found ', i3, &
              :, 1x, ', loc', 2i6)
         endif
      endif
#endif
!
#ifdef DEBUG
      print 9970, trim(ch_id), lev, nfnd, ijkinc
9970  format (1x, a, ': PSMILe_mg_next_level_2d_real: lev =', i3, &
              ': not fnd', i5, ', fnd =', i6, 3i5, ', ijkinc ', 3i4)
#endif /* DEBUG */

#ifdef VERBOSE
      print 9980, trim(ch_id), lev

      call psmile_flushstd
#endif /* VERBOSE */
!
!  Formats:
!
9990 format (1x, a, ': PSMILe_mg_next_level_2d_real: level', i3, &
                    ', control', 6i6)
9980 format (1x, a, ': PSMILe_mg_next_level_2d_real: eof, level', i3)

      end subroutine PSMILe_mg_next_level_2d_real