psmile_mg_next_level_3d_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_3d_real
!
! !INTERFACE:

      subroutine psmile_mg_next_level_3d_real (            &
                           grid_id, lev, nlev,             &
                           chmin1, chmin2, chmin3,         &
                           chmax1, chmax2, chmax3,         &
                           midp1,  midp2,  midp3,          &
                           levdim,                         &
                           found, loc, range,              &
                           coords1, coords2, coords3,      &
                           shape, control,                 &
                           ijkinc, ijkcoa, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_mg_next_level_3d_real

      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_3d)

!     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), 
                                                 0:levdim(3))
      Real, Intent (In)               :: chmin2 (0:levdim(1), 
                                                 0:levdim(2), 
                                                 0:levdim(3))
      Real, Intent (In)               :: chmin3 (0:levdim(1), 
                                                 0:levdim(2), 
                                                 0:levdim(3))

!     Minimum of grid coordinates per grid cell

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

!     Maximum of grid coordinates per grid cell

      Real, Intent (In)               ::  midp1 (0:levdim(1), 
                                                 0:levdim(2), 
                                                 0:levdim(3))
      Real, Intent (In)               ::  midp2 (0:levdim(1), 
                                                 0:levdim(2), 
                                                 0:levdim(3))
      Real, Intent (In)               ::  midp3 (0:levdim(1), 
                                                 0:levdim(2), 
                                                 0:levdim(3))

!     Mid point of the cell

      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,k.
!     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_3d,               
                                                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.

      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))
      Real, Intent (In)               :: coords3 (shape(1,1):shape(2,1), 
                                                  shape(1,2):shape(2,2), 
                                                  shape(1,3):shape(2,3))

!     Coordinates to be searched

      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

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

!     Increment

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

!     Coarsening
!
! !OUTPUT PARAMETERS:
!
      integer, Intent (Out)           :: ierror

!     Returns the error code of PSMILe_MG_coars_level;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !DEFINED PARAMETERS:
!
      Integer, parameter              :: ndtry = 64
      Integer, parameter              :: ntry1 = 8
      Integer, parameter              :: ntry3 = 27
      Integer, parameter              :: ntry2 = ndtry - ntry1
      Real, parameter                 :: remax = 1.0e20
!
! !LOCAL VARIABLES
!
!     ... for levels
!
      Integer                         :: levc
!
!     ... for locations searched
!
      Integer                         :: i, j, k
      Integer                         :: ibegl, n, nprev
      Integer                         :: ifound, nmin(1)
      Integer                         :: ijkf(ndim_3d), ijkold(ndim_3d), newijk(ndim_3d)
      Integer                         :: ijkctl (ndim_3d, ndtry), ijkdst (ndim_3d, ndtry)
      Integer                         :: ijkctl3 (ndim_3d, ndtry)
      Real                            :: dist (ndtry)
      Real                            :: xyz (ndim_3d)
!
#ifdef DEBUG
!     ... for locations searched
!
      Integer                      :: nfnd (0:4)
#endif /* DEBUG */
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_mg_next_level_3d_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_3d_real.F90 2325 2010-04-21 15:00:07Z valcke $
!  $Author: valcke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_mg_next_level_3d_real.F90 2325 2010-04-21 15:00:07Z valcke $'
!
!----------------------------------------------------------------------
!
!  Relative control indices for sub-cells of coarse grid cells
!
      data ((ijkctl (i, n), i=1,ndim_3d), n = 1,ntry1) &
                             / 0, 0, 0,   1, 0, 0,    &
                               0, 1, 0,   1, 1, 0,    &
                               0, 0, 1,   1, 0, 1,    &
                               0, 1, 1,   1, 1, 1 /
!
!  Relative control indices for neighbour cells
!
      data ((ijkctl (i, n), i=1,ndim_3d), n = ntry1+1,  ntry1+16) &
                  /-1,-1,-1,   0,-1,-1,   1,-1,-1,   2,-1,-1,    &
                   -1, 0,-1,   0, 0,-1,   1, 0,-1,   2, 0,-1,    &
                   -1, 1,-1,   0, 1,-1,   1, 1,-1,   2, 1,-1,    &
                   -1, 2,-1,   0, 2,-1,   1, 2,-1,   2, 2,-1/
      data ((ijkctl (i, n), i=1,ndim_3d), n = ntry1+17, ntry1+28) &
                  /-1,-1, 0,   0,-1, 0,   1,-1, 0,   2,-1, 0,    &
                   -1, 0, 0,                         2, 0, 0,    &
                   -1, 1, 0,                         2, 1, 0,    &
                   -1, 2, 0,   0, 2, 0,   1, 2, 0,   2, 2, 0/
      data ((ijkctl (i, n), i=1,ndim_3d), n = ntry1+29, ntry1+40) &
                  /-1,-1, 1,   0,-1, 1,   1,-1, 1,   2,-1, 1,    &
                   -1, 0, 1,                         2, 0, 1,    &
                   -1, 1, 1,                         2, 1, 1,    &
                   -1, 2, 1,   0, 2, 1,   1, 2, 1,   2, 2, 1/
      data ((ijkctl (i, n), i=1,ndim_3d), n = ntry1+41, ntry1+56) &
                  /-1,-1, 2,   0,-1, 2,   1,-1, 2,   2,-1, 2,    &
                   -1, 0, 2,   0, 0, 2,   1, 0, 2,   2, 0, 2,    &
                   -1, 1, 2,   0, 1, 2,   1, 1, 2,   2, 1, 2,    &
                   -1, 2, 2,   0, 2, 2,   1, 2, 2,   2, 2, 2/
!
!  Relative control indices for computed cells
!
      data ((ijkctl3 (i, n), i=1,ndim_3d), n = 1,ntry3) &
                 / -1,-1,-1,   0,-1,-1,   1,-1,-1,     &
                   -1, 0,-1,   0, 0,-1,   1, 0,-1,     &
                   -1, 1,-1,   0, 1,-1,   1, 1,-1,     &
                   -1,-1, 0,   0,-1, 0,   1,-1, 0,     &
                   -1, 0, 0,   0, 0, 0,   1, 0, 0,     &
                   -1, 1, 0,   0, 1, 0,   1, 1, 0,     &
                   -1,-1, 1,   0,-1, 1,   1,-1, 1,     &
                   -1, 0, 1,   0, 0, 1,   1, 0, 1,     &
                   -1, 1, 1,   0, 1, 1,   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
!
#ifdef DEBUG
      nfnd (0) = 0
      nfnd (1) = 0
      nfnd (2) = 0
      nfnd (3) = 0
      nfnd (4) = 0
#endif /* DEBUG */

         do k = control(1,3), control(2,3), ijkinc (3)
           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: Index list instead of while loop
!
!     TO TEST:
!         (i)   Compute minimal distance without checking chmin and chmax
!         (ii)  Check chmin and chmax for cell found
!         (iii) If outside, check all cells and mark cells by remax and
!               compute minimum again
!-----------------------------------------------------------------------
!
      nprev = 0

      ibegl = control (1, 1)
!
!===> ... Look for the next cell to be controlled
!
      do while (ibegl <= control(2,1))
!cdir vector
           do i = ibegl, control(2,1), ijkinc(1)
           if (found (i,j,k) == levc) go to 11
           end do
!
        if (i .lt. control(2,1)+ijkinc(1)) then
           i = control(2,1)
           if (found (i,j,k) == levc) go to 11
        endif
!
        exit
!
!===> ... Look in the fine parts of cell loc(i,j,k)
!
11         continue
!
           ijkf(:) = min (loc(:,i,j,k) * ijkcoa(:), levdim(:))
!
!cdir vector
              do n = 1, ntry1
              ijkdst (:, n) = max(0, min(ijkf(:) + ijkctl(:,n), levdim(:)))
              end do

!cdir vector
              do n = 1, ntry1
              if (coords1(i,j,k) >= chmin1(ijkdst(1,n),         &
                                           ijkdst(2,n),         &
                                           ijkdst(3,n)) .and.   &
                  coords2(i,j,k) >= chmin2(ijkdst(1,n),         &
                                           ijkdst(2,n),         &
                                           ijkdst(3,n)) .and.   &
                  coords3(i,j,k) >= chmin3(ijkdst(1,n),         &
                                           ijkdst(2,n),         &
                                           ijkdst(3,n)) .and.   &
                  coords1(i,j,k) <= chmax1(ijkdst(1,n),         &
                                           ijkdst(2,n),         &
                                           ijkdst(3,n)) .and.   &
                  coords2(i,j,k) <= chmax2(ijkdst(1,n),         &
                                           ijkdst(2,n),         &
                                           ijkdst(3,n)) .and.   &
                  coords3(i,j,k) <= chmax3(ijkdst(1,n),         &
                                           ijkdst(2,n),         &
                                           ijkdst(3,n))) then
                 dist (n) = (coords1(i,j,k) -                                   &
                               midp1(ijkdst(1,n), ijkdst(2,n), ijkdst(3,n)))**2 &
                          + (coords2(i,j,k) -                                   &
                               midp2(ijkdst(1,n), ijkdst(2,n), ijkdst(3,n)))**2 &
                          + (coords3(i,j,k) -                                   &
                               midp3(ijkdst(1,n), ijkdst(2,n), ijkdst(3,n)))**2

              else

                 dist (n) = remax

              endif
              end do
!
!===> ... Look for the minimum distance
!
              nmin = MINLOC (dist(1:ntry1))
#ifdef MINLOCFIX
              if (nmin(1).eq.0) nmin=1
#endif /* MINLOCFIX */
!
              if (dist(nmin(1)) .ne. remax) then
                 found (i,j,k)  = lev
                 loc (:,i,j,k) = ijkdst (:,nmin(1))
!
#ifdef DEBUG
                 nfnd (1) = nfnd (1) + 1
#endif /* DEBUG */
                 go to 95
!
              else
!
           if (levc == nlev) then
              found (i,j,k) = - found(i,j,k)
#ifdef DEBUG
              nfnd (0) = nfnd (0) + 1
#endif /* DEBUG */
              go to 95
           endif
        endif
!
!-----------------------------------------------------------------------
!          Look in the neigbourhood of the coarse grid cell
!          loc(i,j,k)
!-----------------------------------------------------------------------
!
!cdir vector
              do n = ntry1+1, ndtry
              ijkdst (:, n) = max (0, min(ijkf(:) + ijkctl(:,n), levdim(:)))
              end do
!
              do n = ntry1+1, ndtry
              if (coords1(i,j,k) >= chmin1(ijkdst(1,n),         &
                                           ijkdst(2,n),         &
                                           ijkdst(3,n)) .and.   &
                  coords2(i,j,k) >= chmin2(ijkdst(1,n),         &
                                           ijkdst(2,n),         &
                                           ijkdst(3,n)) .and.   &
                  coords3(i,j,k) >= chmin3(ijkdst(1,n),         &
                                           ijkdst(2,n),         &
                                           ijkdst(3,n)) .and.   &
                  coords1(i,j,k) <= chmax1(ijkdst(1,n),         &
                                           ijkdst(2,n),         &
                                           ijkdst(3,n)) .and.   &
                  coords2(i,j,k) <= chmax2(ijkdst(1,n),         &
                                           ijkdst(2,n),         &
                                           ijkdst(3,n)) .and.   &
                  coords3(i,j,k) <= chmax3(ijkdst(1,n),         &
                                           ijkdst(2,n),         &
                                           ijkdst(3,n))) then
                 dist (n) = (coords1(i,j,k) -                                   &
                               midp1(ijkdst(1,n), ijkdst(2,n), ijkdst(3,n)))**2 &
                          + (coords2(i,j,k) -                                   &
                               midp2(ijkdst(1,n), ijkdst(2,n), ijkdst(3,n)))**2 &
                          + (coords3(i,j,k) -                                   &
                               midp3(ijkdst(1,n), ijkdst(2,n), ijkdst(3,n)))**2

              else

                 dist (n) = remax

              endif
              end do
!
!===> ... Look for the minimum distance
!
           nmin = MINLOC (dist(ntry1+1:ndtry)) + ntry1
#ifdef MINLOCFIX
           if (nmin(1).eq.ntry1) nmin=ntry1+1
#endif /* MINLOCFIX */
!
           if (dist(nmin(1)) .ne. remax) then
              found (i,j,k)  = lev
              loc (:,i,j,k) = ijkdst (:,nmin(1))
!
#ifdef DEBUG
              nfnd (2) = nfnd (2) + 1
#endif /* DEBUG */
!
           else
!
              nprev = nprev + 1
           endif
!
!===> ... Start index of the next search
!
95         if (i == control(2,1)) exit
           ibegl = min (control(2,1), i+ijkinc(1))

         end do ! while
!
!===> ... Have points to be controlled o a coarser level ?
!
         if (nprev == 0) go to 20
!
!-----------------------------------------------------------------------
!          Some points were not found.
!          Try to interpolate the fine point a from neighbouring points
!-----------------------------------------------------------------------
!
        ibegl = control (1, 1) + ijkinc(1)
!
        do while (ibegl <= control(2,1))
!cdir vector
              do i = ibegl, 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) go to 211
              end do
!
           if (i .lt. 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) go to 211
           endif
!
           exit
!
!===> ... Interpolate cell IJKF and
!         look in the neighbourhood of cell IJKF
!
211        continue

           ijkf (:) = min((loc (:, i-ijkinc(1),j,k) + loc (:, i+ijkinc(1),j,k))/2, &
                          levdim(:))
!
           ijkold (:) = min (loc(:, i,j,k) * ijkcoa(:), levdim(:))
!
           if (ijkold(1) == ijkf(1) .and. &
               ijkold(2) == ijkf(2) .and. &
               ijkold(3) == ijkf(3)) go to 260

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

              else

                 dist (n) = remax

              endif
              end do
!
!===> ... Look for the minimum distance
!
           nmin = MINLOC (dist(1:ntry3))
#ifdef MINLOCFIX
           if (nmin(1).eq.0) nmin=1
#endif /* MINLOCFIX */
!
           if (dist(nmin(1)) .ne. remax) then
              found (i,j,k)  = lev
              loc (:, i,j,k) = ijkdst (:, nmin(1))
!
              nprev = nprev - 1
!
#ifdef DEBUG
              nfnd (3) = nfnd (3) + 1
#endif /* DEBUG */
        endif
!
!===> ... Start index of the next search
!
260        if (i >= control(2,1)-ijkinc(1)) exit
           ibegl = min (control(2,1)-ijkinc(1), i+ijkinc(1))
!
         end do ! while
!
         if (nprev == 0) go to 20

!-----------------------------------------------------------------------
!     Some points were still not found.
!     Look again on the next coarser level whether the points
!     are located in other cells.
!-----------------------------------------------------------------------
!
      ibegl = control (1, 1)

      do while (ibegl <= control(2,1))
!
!===> ... Look for the next cell to be controlled
!
!cdir vector
            do i = ibegl, control(2,1), ijkinc(1)
            if (found (i,j,k) == levc) go to 311
            end do
!
         if (i .lt. control(2,1)+ijkinc(1)) then
            i = control(2,1)
            if (found (i,j,k) == levc) go to 311
         endif
!
         exit
!
!===> ... Look in the coarser parts of cell loc(i,j,k)
!
311      continue
!
         if (lev+1 < nlev) then
            xyz (1) = coords1 (i,j,k)
            xyz (2) = coords2 (i,j,k)
            xyz (3) = coords3 (i,j,k)
!
            call psmile_mg_prev_levels_3d_real (grid_id, lev, nlev, &
                              loc(1,i,j,k), xyz, ifound, newijk)
!
            if (ifound .gt. 0) then
               loc (:, i,j,k) = newijk (:)
!
               found (i,j,k) = lev
#ifdef DEBUG
               nfnd (4) = nfnd (4) + 1
#endif /* DEBUG */
            else
               found (i,j,k) = - found (i,j,k)
#ifdef DEBUG
               nfnd (0) = nfnd (0) + 1
#endif /* DEBUG */
            endif
         else
!           This is the second coarsest level.
!           The coarsest level contains only 1 cell.
!           All points which are not found are outside.

            found (i,j,k) = - found (i,j,k)
#ifdef DEBUG
            nfnd (0) = nfnd (0) + 1
#endif /* DEBUG */

         endif
!
!===> ...
!
         nprev = nprev - 1
         if (nprev == 0) go to 20
!
!===> ... Start index of the next search
!
         if (i == control(2,1)) exit

         ibegl = min (control(2,1), i+ijkinc(1))
      end do ! while
!
20    end do ! j
      end do ! k
!
!===> All done
!
#ifdef DEBUG
      print 9970, trim(ch_id), lev, nfnd, ijkinc
9970  format (1x, a, ': PSMILe_mg_next_level_3d_real: lev =', i3, &
              ': fnd =', 5i5, ', 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_3d_real: level', i3, &
                    ', control', 6i6)
9980 format (1x, a, ': PSMILe_mg_next_level_3d_real: eof, level', i3)

      end subroutine PSMILe_mg_next_level_3d_real