psmile_bbcells_virt_2d_dble.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_Bbcells_virt_2d_dble
!
! !INTERFACE:

      subroutine psmile_bbcells_virt_2d_dble (method_id,                &
                      coords1, coords2, coords_shape, coords_range,     &
                      corners1, corners2, corner_shape, nbr_corners,    &
                      chmin1_corner, chmin2_corner,                     &
                      chmax1_corner, chmax2_corner, levdim_corner,      &
                      chmin1,  chmin2, chmax1, chmax2, midp1, midp2,    &
                      levdim,  period, bmaski, bmaskj, ierror)
!
! !USES:
!
      use PRISM_constants
      use psmile_grid, only : psmile_transform_cell_cyclic
      use PSMILe, dummy_interface => PSMILe_Bbcells_virt_2d_dble

      Implicit none
!
! !INPUT PARAMETERS:
!
      Integer,           Intent (In)  :: method_id

!     Method id to access halo information

      Integer,           Intent (In)  :: coords_shape (2, ndim_2d)

!     Dimension of (method/grid) coordinate "coords1", ...

      Integer,           Intent (In)  :: coords_range (2, ndim_2d)

!     Definition of the subgrid for which the (bounding) cells should
!     be computed.
!     In general, this should be identical with "grid_valid_shape".

      Integer,           Intent (In)  :: corner_shape (2, ndim_2d)

!     Dimension of the array "corners1" containing the corner points
!     of the control volume (for a single coordinate).

      Integer,           Intent (In)  :: nbr_corners

!     Number of corners per control volume

      Double Precision,  Intent (In)  :: 
         corners1 (corner_shape(1,1):corner_shape(2,1), 
                   corner_shape(1,2):corner_shape(2,2), 
                   nbr_corners)
      Double Precision,  Intent (In)  :: 
         corners2 (corner_shape(1,1):corner_shape(2,1), 
                   corner_shape(1,2):corner_shape(2,2), 
                   nbr_corners)

!     Definition of the subgrid for which the (bounding)
!     cells should be computed.

      Integer, Intent (In)            :: levdim_corner (2, ndim_2d)

!     Dimensions of chmin, chmax and midp on level 1

      Double Precision,  Intent (In)  :: 
         chmin1_corner (levdim_corner(1,1):levdim_corner(2,1), 
                        levdim_corner(1,2):levdim_corner(2,2))
      Double Precision,  Intent (In)  :: 
         chmin2_corner (levdim_corner(1,1):levdim_corner(2,1), 
                        levdim_corner(1,2):levdim_corner(2,2))
      Double Precision,  Intent (In)  :: 
         chmax1_corner (levdim_corner(1,1):levdim_corner(2,1), 
                        levdim_corner(1,2):levdim_corner(2,2))
      Double Precision,  Intent (In)  :: 
         chmax2_corner (levdim_corner(1,1):levdim_corner(2,1), 
                        levdim_corner(1,2):levdim_corner(2,2))

!     Bounding box info of corner cells

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

!     Dimension of output arrays "chmin1", "chmax1" and "midp1".
!     For testing.
!
      Double Precision,  Intent (In)  :: period
!
!     Type (Split_Information), Intent (InOut) :: split_info
!
! !INPUT/OUTPUT PARAMETERS:
!
      Double Precision, Intent (InOut) :: 
         chmin1 (coords_range(1,1)-1:coords_range(2,1), 
                 coords_range(1,2)-1:coords_range(2,2))
      Double Precision,  Intent (InOut) :: 
         chmin2 (coords_range(1,1)-1:coords_range(2,1), 
                 coords_range(1,2)-1:coords_range(2,2))
!
!     Minimum value of bounding box of method cell
!        chmin1 (i,j) contains the minimal value of bounding box of
!                     method cell ((i,j), (i+1,j), (i,j+1), (i+1,j+1))
!                     of "coords1"
!        chmin2 (i,j) contains the minimal value of bounding box of
!                     method cell ((i,j), (i+1,j), (i,j+1), (i+1,j+1))
!                     of "coords2"

      Double Precision, Intent (InOut) :: 
         chmax1 (coords_range(1,1)-1:coords_range(2,1), 
                 coords_range(1,2)-1:coords_range(2,2))
      Double Precision, Intent (InOut) :: 
         chmax2 (coords_range(1,1)-1:coords_range(2,1), 
                 coords_range(1,2)-1:coords_range(2,2))
!
!     Maximum value of bounding box of method cell
!
      Double Precision, Intent (InOut) :: 
         midp1  (coords_range(1,1)-1:coords_range(2,1), 
                 coords_range(1,2)-1:coords_range(2,2))
      Double Precision, Intent (InOut) :: 
         midp2  (coords_range(1,1)-1:coords_range(2,1), 
                 coords_range(1,2)-1:coords_range(2,2))
!
!     Mid point of bounding box of method cell
!
! !OUTPUT PARAMETERS:
!
      Double Precision, Intent (InOut)   ::           
         coords1(coords_shape(1,1):coords_shape(2,1), 
                 coords_shape(1,2):coords_shape(2,2))
      Double Precision, Intent (InOut)   ::           
         coords2(coords_shape(1,1):coords_shape(2,1), 
                 coords_shape(1,2):coords_shape(2,2))

!        Input: Coordinates of the method
!
!        Output:
!        coords(coords_range(1,1)-1:coords_range(2,1)+1, coords_range(1,2)-1)
!        coords(coords_range(1,1)-1:coords_range(2,1)+1, coords_range(2,2)+1)
!        coords(coords_range(1,1)-1, coords_range(1,2)-1:coords_range(2,2)+1)
!        coords(coords_range(2,1)+1, coords_range(1,2)-1:coords_range(2,2)+1)
!
!        Coordinates of the method points over the boundary
!        which were received or inter/extrapolated.
!
      Integer,           Intent (Out) :: 
         bmaski (coords_range(1,1)-1:coords_range(2,1)+1, 2)
      Integer,           Intent (Out) :: 
         bmaskj (coords_range(1,2)-1:coords_range(2,2)+1, 2)
!
!     Mask for virtual method cells with
!     bmaski (i,lb) = Mask value for vitual method points
!        coords(coords_range(1,1)-1:coords_range(2,1)+1, coords_range(1,2)-1)
!                     and
!        coords(coords_range(1,1)-1:coords_range(2,1)+1, coords_range(2,2)+1)
!     bmaskj (j,lb) = Mask value for virtual method points
!        coords(coords_range(1,1)-1, coords_range(1,2)-1:coords_range(2,2)+1)
!                     and
!        coords(coords_range(2,1)+1, coords_range(1,2)-1:coords_range(2,2)+1)
!     bmask? (i,lb) =
!         neighbour_value        : Value was received from neighbouring block
!         neighbour_interpolated : Value was interpolated from
!                                  neighbouring method points
!         neighbour_extrapolated : Value was extrapolated from 
!                                  interior method points
!
      Integer,           Intent (Out) :: ierror

!     Returns the error code of PSMILe_Bbcells_virt_2d_dble
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !DEFINED PARAMETERS:
!
! r_nbr = 1 / Number of corners1 per grid cell.
!
      Integer,          Parameter  :: n_corners = 4
      Double Precision, Parameter  :: r_nbr = 0.25d0
!
      Integer,          Parameter  :: no_value = 0
      Integer,          Parameter  :: neighbour_value = 1
      Integer,          Parameter  :: neighbour_interpolated = 2
      Integer,          Parameter  :: neighbour_extrapolated = 3
!
! !LOCAL VARIABLES
!
      Integer                      :: i, i1, i2, iv, ibeg, iend, ifirst
      Integer                      :: j, j1, j2, jv, jbeg, jend, jfirst
      Integer                      :: lb, l, n
!
      Double Precision             :: dmn, dmx, fac, r
      Double Precision             :: d (ndim_2d)
!
      Integer, Allocatable         :: liste(:)
!
!     Auxillary for halo regions
!
      Integer                      :: grid_id
      Integer                      :: face
      Integer                      :: n_seg
      Integer                      :: ishift
      Integer                      :: nbr_halo_segments
      Integer                      :: halo_shape(2,ndim_3d)
      Type (Halo_Block), Pointer   :: halo_pointer
!
!     Auxiliary fields for cyclic method cells
!
      Double Precision             :: cell (n_corners)
!
      Double Precision             :: period2
!
! Error parameters
!
      Integer, Parameter           :: nerrp = 2
      Integer                      :: ierrp (nerrp)
!
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_Bbcells_virt_2d_dble" computes the bounding box
! for each virtual cell of the 2-dimensional block. The virtual cells
! of the block are the cells
! (coords_range(1,1)-1, :), (coords_range(2,1), :),
! (:, coords_range(1,2)-1), (:, coords_range(2,2))
!
! chmin1 (i, j) returns the minimal value of bounding box of virtual
!               method cell ((i,j), (i+1,j), (i,j+1), (i+1,j+1))
!
! The first box
! chmin1 (coords_range(1,1)-1,j) returns the minimal value of the bounding box
!                         of virtual cell
!       (coords1(coords_range(1,1), j), coords1(coords_range(1,1), j+1),
!        bounding box of corner control volume of (coords_range(1,1),j)
!
! The last box
! chmin1 (coords_range(2,1),  j) returns the minimal value of the bounding box
!                         of virtual cell
!       (coords1(coords_range(2,1), j), coords1(coords_range(2,1), j+1),
!        bounding box of corner control volume of (coords_range(2,1),j)
!
! Corresponding values are returned for other virtual cells at block
! border.
!
! In addition, the virtual method points (received from neighbouring
! block or inter/extrapolated) are returned together with a mask
! which contains a code how the virtual boundary point was created.
!
! Subroutine "PSMILe_Bbcells_virt_2d_dble" is designed for grids of
! type "PRISM_Irrlonlat_regvrt".
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 05.05.09    H. Ritzdorf  created
! 13.10.10    M. Hanke     changed handling of cyclic coordinates
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_bbcells_virt_2d_dble.F90 2694 2010-10-29 14:58:50Z hanke $
!  $Author: hanke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_bbcells_virt_2d_dble.F90 2694 2010-10-29 14:58:50Z hanke $'
!
!----------------------------------------------------------------------
!
#ifdef VERBOSE
      print 9990, trim(ch_id)

      call psmile_flushstd
#endif /* VERBOSE */
!
#ifdef PRISM_ASSERTION
      if (coords_range(1,1) < corner_shape(1,1) .or. &
          coords_range(2,1) > corner_shape(2,1) .or. &
          coords_range(1,2) < corner_shape(1,2) .or. &
          coords_range(2,2) > corner_shape(2,2)) then
         print *, trim(ch_id), "PSMILe_bbcells_virt_2d_dble: range", &
                  coords_range
         print *, trim(ch_id), "PSMILe_bbcells_virt_2d_dble: corner_shape", &
                  corner_shape

         call psmile_assert (__FILE__, __LINE__, &
              "range must be in corner_shape")
      endif
!
      if (minval(levdim(:)-(coords_range(2,:)-coords_range(1,:)+2)) < 0) then
         print *, trim(ch_id), "PSMILe_bbcells_virt_2d_dble: levdim", levdim
         print *, trim(ch_id), "PSMILe_bbcells_virt_2d_dble: range ", &
                  coords_range

         call psmile_assert (__FILE__, __LINE__, &
              "levdim < coords_range(2,:)-coords_range(1,:)+2")
      endif

      if (coords_range(1,1)-1 < coords_shape (1,1) .or. &
          coords_range(2,1)+1 > coords_shape (2,1) .or. &
          coords_range(1,2)-1 < coords_shape (1,2) .or. &
          coords_range(2,2)+1 > coords_shape (2,2)) then

         print *, trim(ch_id), "PSMILe_bbcells_virt_2d_dble: coords_shape", &
                               coords_shape
         print *, trim(ch_id), "PSMILe_bbcells_virt_2d_dble: coords_range", &
                               coords_range

         call psmile_assert (__FILE__, __LINE__, &
              "coords_shape is not sufficent")
      endif

      if (coords_range(1,1) < levdim_corner (1,1) .or. &
          coords_range(2,1) > levdim_corner (2,1) .or. &
          coords_range(1,2) < levdim_corner (1,2) .or. &
          coords_range(2,2) > levdim_corner (2,2)) then

         print *, trim(ch_id), "PSMILe_bbcells_virt_2d_dble: levdim_corner", &
                               levdim_corner
         print *, trim(ch_id), "PSMILe_bbcells_virt_2d_dble: coords_range", &
                               coords_range

         call psmile_assert (__FILE__, __LINE__, &
              "levdim_corner is not sufficent")
      endif
#endif
!
!  Initialization
!
      ierror = 0
!
!  The corner cells (chmin?_corner, corners?, ...) are used in the
!      range of (ibeg:iend, jbeg:jend).
!  The method points (coords1, coords2) are used in the same range
!      such as the corner cells.
!  The method cells (chmin?, chmax?, ...) are used in the
!      range of (ibeg-1:iend, jbeg-1:jend)
!      where the cell (ibeg-1, :), (iend, :), (:, jbeg-1), (:. jend)
!      are virtual method cells.
!  The boundary (method) points are corners of the virtual method
!      cells which are outside of the interior of corner cells and
!      are used in the range of (ibeg-1:iend+1, *) and
!      (jbeg-1:jend+1, *) where boundary? (:, ij, *) contains
!      the coordinates of lower left corner of method cell 
!      chmin? (i, j), ...
!     
      ibeg = coords_range (1,1)
      iend = coords_range (2,1)
!
      jbeg = coords_range (1,2)
      jend = coords_range (2,2)

#ifdef DEBUG
      PRINT*, 'ibeg, iend :',ibeg,iend
      PRINT*, 'jbeg, jend :',jbeg,jend
      CALL psmile_flushstd
#endif
!
      period2 = period * 0.5d0
!
      bmaski (:, :) = no_value
      bmaskj (:, :) = no_value
!
      face = PSMILe_Undef
!=======================================================================
!     Apply halo information if available
!
!     (ibeg:iend,   jbeg:jend)  is the cell range
!     (ibeg:iend+1, jbeg:jend+1) is the range for boundary points
!=======================================================================
!
      grid_id = Methods(method_id)%grid_id
      nbr_halo_segments = Grids(grid_id)%nbr_halo_segments
!
!
!             face 3
!          -------------
!  face 4  |           |  face 2
!          -------------
!             face 1
!
      do n_seg = 1, nbr_halo_segments

         halo_pointer => Methods(method_id)%halo_pointer(n_seg)
         halo_shape = halo_pointer%halo_shape
#ifdef DEBUG
         PRINT*, 'number of halo segments :',nbr_halo_segments
         PRINT*, 'Halo_shape :',halo_shape
         call psmile_flushstd
#endif

! UNSCHOEN !!
         if ( halo_shape(1,1) == PSMILe_Undef ) cycle

         if ( halo_shape(1,1) == halo_shape(2,1) ) then
            if ( halo_shape(1,1) == ibeg-1 ) then
               face = 4
            else if ( halo_shape(1,1) == iend+1 ) then
               face = 2
            endif
         else if ( halo_shape(1,2) == halo_shape(2,2) ) then
            if ( halo_shape(1,2) == jbeg-1 ) then
               face = 1
            else if ( halo_shape(1,2) == jend+1 ) then
               face = 3
            endif
         endif

#ifdef PRISM_ASSERTION
         IF (face == PSMILe_Undef) THEN

             CALL psmile_assert (__FILE__, __LINE__, &
                "face is undefined ! ")
         ENDIF
#endif

#ifdef DEBUG
         PRINT *, TRIM(ch_id), ' - applying halo for face ', face
         CALL psmile_flushstd
#endif /* DEBUG */

         select case (face)
         case (1)

            ishift = 1-halo_shape(1,1)

!cdir vector
!DEC$ IVDEP
               do i = max(ibeg-1,halo_shape(1,1)), min(iend+1,halo_shape(2,1))
               bmaski  (i,1)  = neighbour_value
               coords1 (i,jbeg-1) = halo_pointer%halo_dble(1)%vector(i+ishift)
               coords2 (i,jbeg-1) = halo_pointer%halo_dble(2)%vector(i+ishift)
               end do

#ifdef DEBUG
         print *, trim(ch_id), 'face ', face, 'ibeg:iend, jbeg-1', &
                  max(ibeg-1,halo_shape(1,1)), &
                  min(iend+1,halo_shape(2,1)), jbeg-1
#endif /* DEBUG */

         case (2)

            ishift = 1-halo_shape(1,2)
#ifdef DEBUG
         print *, trim(ch_id), 'face ', face, 'iend+1, jbeg:jend', &
                  iend+1, max(jbeg-1,halo_shape(1,2)), &
                          min(jend+1,halo_shape(2,2))
#endif /* DEBUG */

!cdir vector
!DEC$ IVDEP
               do j = max(jbeg-1,halo_shape(1,2)), min(jend+1,halo_shape(2,2))
               bmaskj  (j,2) = neighbour_value
               coords1 (iend+1,j) = halo_pointer%halo_dble(1)%vector(j+ishift)
               coords2 (iend+1,j) = halo_pointer%halo_dble(2)%vector(j+ishift)
               end do

         case (3)

            ishift = 1-halo_shape(1,1)
#ifdef DEBUG
         print *, trim(ch_id), 'face ', face, 'ibeg:iend, jend+1', &
                  max(ibeg-1,halo_shape(1,1)), &
                  min(iend+1,halo_shape(2,1)), jend+1
#endif /* DEBUG */

!cdir vector
!DEC$ IVDEP
               do i = max(ibeg-1,halo_shape(1,1)), min(iend+1,halo_shape(2,1))
               bmaski  (i,2) = neighbour_value
               coords1 (i,jend+1) = halo_pointer%halo_dble(1)%vector(i+ishift)
               coords2 (i,jend+1) = halo_pointer%halo_dble(2)%vector(i+ishift)
               end do

         case (4)

            ishift = 1-halo_shape(1,2)
#ifdef DEBUG
         print *, trim(ch_id), 'face ', face, 'ibeg-1, jbeg:jend', &
                  ibeg-1, max(jbeg-1,halo_shape(1,2)), &
                          min(jend+1,halo_shape(2,2))
#endif /* DEBUG */

!cdir vector
!DEC$ IVDEP
               do j = max(jbeg-1,  halo_shape(1,2)), min(jend+1,halo_shape(2,2))
               bmaskj  (j,1) = neighbour_value
               coords1 (ibeg-1,j) = halo_pointer%halo_dble(1)%vector(j+ishift)
               coords2 (ibeg-1,j) = halo_pointer%halo_dble(2)%vector(j+ishift)
               end do

         end select

      enddo ! n_seg
!
! Are values in the corners already available ?
! Currently not available.
!
      if (bmaski (ibeg-1, 1) == neighbour_value) then
          bmaskj (jbeg-1, 1) = neighbour_value
      else if (bmaskj (jbeg-1, 1) == neighbour_value) then
          bmaski (ibeg-1, 1) = neighbour_value
      endif
!
      if (bmaski (iend+1, 1) == neighbour_value) then
          bmaskj (jbeg-1, 2) = neighbour_value
      else if (bmaskj (jbeg-1, 2) == neighbour_value) then
          bmaski (iend+1, 1) = neighbour_value
      endif
!
      if (bmaski (ibeg-1, 2) == neighbour_value) then
          bmaskj (jend+1, 1) = neighbour_value
      else if (bmaskj (jend+1, 1) == neighbour_value) then
          bmaski (ibeg-1, 2) = neighbour_value
      endif
!
      if (bmaski (iend+1, 2) == neighbour_value) then
          bmaskj (jend+1, 2) = neighbour_value
      else if (bmaskj (jend+1, 2) == neighbour_value) then
          bmaski (iend+1, 2) = neighbour_value
      endif

!-----------------------------------------------------------------------
!     Interpolate missing virtual method cell points
!     (*) which are not the corner points (i.e. from (ibeg:iend)) and
!     (*) which are neighboured by 2 neighbouring values
!     for faces 1 and 3 (j = jbeg and j = jend)
!-----------------------------------------------------------------------
!
      do lb = 1, 2
      if ( .not. ANY (bmaski(ibeg:iend, lb) == no_value) ) cycle

      if (lb == 1) then 
!        Face 1 (j = jbeg)
!        j  = Index of the first interior method points/corner cell
!        jv = Index of the virtual method point
         j  = jbeg
         jv = j - 1
      else 
!        Face 2 (j = jend);
!        j = Index of last interior method point/corner cell
!        jv = Index of the virtual method point
         j  = jend
         jv = j + 1
      endif

!cdir vector
!DEC$ IVDEP
         do i = ibeg, iend
         if (bmaski(i, lb) == no_value .and. &
            (bmaski(i-1, lb) == neighbour_value .and. &
             bmaski(i+1, lb) == neighbour_value)) then
!        Interpolate missing point by neighboured method points
! bei der Interpolation koennte man noch den methodenpunkt (i,j)
! und (i-1,j) und (i+1, j) einfliessen lassen
!
            bmaski  (i,lb) = neighbour_interpolated
            coords1 (i,jv) = (coords1 (i+1,jv) + coords1 (i-1,jv)) * 0.5
            coords2 (i,jv) = (coords2 (i+1,jv) + coords2 (i-1,jv)) * 0.5
         endif
         end do
!
      end do ! lb
!
!-----------------------------------------------------------------------
!     Interpolate missing virtual method cell points
!     (*) which are not the corner points (i.e. from (jbeg:jend)) and
!     (*) which are neighboured by 2 neighbouring values
!     for faces 2 and 4 (i = ibeg and i = iend)
!-----------------------------------------------------------------------
!
         do lb = 1, 2
         if (.not. ANY (bmaskj(jbeg:jend, lb) == no_value)) cycle

         if (lb == 1) then 
!           Face 3 (i = ibeg)
!           i  = Index of the first interior method points/corner cell
!           iv = Index of the virtual method point
            i  = ibeg
            iv = i - 1
         else 
!           Face 4 (i = iend)
!           i  = Index of the last interior method points/corner cell
!           iv = Index of the virtual method point
            i  = iend
            iv = i + 1
         endif

!cdir vector
!DEC$ IVDEP
         do j = jbeg, jend
         if (bmaskj(j, lb) == no_value .and. &
            (bmaskj(j-1, lb) == neighbour_value .and. &
             bmaskj(j+1, lb) == neighbour_value)) then
!        Interpolate missing point by neighboured method points
! bei der Interpolation koennte man noch den methodenpunkt (i,j)
! und (i,j-1) und (i, j+1) einfliessen lassen
!
            bmaskj  (j,lb) = neighbour_interpolated
            coords1 (iv,j) = (coords1 (iv,j+1) + coords1 (iv,j-1)) * 0.5
            coords2 (iv,j) = (coords2 (iv,j+1) + coords2 (iv,j-1)) * 0.5
         endif
         end do
!
      end do ! lb
!
!=======================================================================
!     Special case: 1-dimensional
!=======================================================================
!
!     if (coords_range(1,2) == coords_range (2,2)) then
      if (jbeg == jend) then
!
         j = jbeg
         if (ibeg < iend) then
!
!  Try to interpolate all points
!  Simple linear logical interpolation for this exotic case
!  Implement something better if required
!
           do lb = 1, 2
           if (lb == 1) then
              jv = j - 1
           else
              jv = j + 1
           endif
!
           n = 0
!
!cdir vector
               do i1 = ibeg-1, iend+1
               if (bmaski(i1, lb) == neighbour_value) exit
               end do
!
            if (i1 < iend+1) then
               ifirst = i1
               i2 = i1 + 1
               n = 1
!
               do while (i2 < iend+1)
                  do i2 = i1+1, iend+1
                  if (bmaski(i2, lb) == neighbour_value) exit
                  end do
!
               if (i2 > iend+1) exit

               n =  n + 1
!
               if (i2 == i1 + 1) then 
                  i1 = i2
                  cycle
               endif
!
               d(1) = (coords1 (i2, jv) - coords1 (i1, jv)) / (i2-i1)
               d(2) = (coords2 (i2, jv) - coords2 (i1, jv)) / (i2-i1)
!
!cdir vector
!DEC$ IVDEP
                  do i = i1+1, i2-1
                  bmaski(i, lb) = neighbour_interpolated
                  coords1 (i,jv) = coords1 (i1,jv) + d(1)*(i-i1)
                  coords2 (i,jv) = coords2 (i1,jv) + d(2)*(i-i1)
                  end do
!
               i1 = i2
               end do ! while
!
            if (n > 1) then
               if (ifirst > ibeg-1) then
!
                  d(1) = coords1 (ifirst+1, jv) - coords1 (ifirst, jv)
                  d(2) = coords2 (ifirst+1, jv) - coords2 (ifirst, jv)
!cdir vector
!DEC$ IVDEP
                     do i = ibeg-1, ifirst-1
                     bmaski(i, lb) = neighbour_extrapolated
                     coords1 (i,jv) = coords1 (ifirst,jv) - d(1)*(ifirst-i)
                     coords2 (i,jv) = coords2 (ifirst,jv) - d(2)*(ifirst-i)
                     end do
               endif
!
               if (i2 < iend+1) then
!
                  d(1) = coords1 (i2, jv) - coords1 (i2-1, jv)
                  d(2) = coords2 (i2, jv) - coords2 (i2-1, jv)
!cdir vector
!DEC$ IVDEP
                     do i = i2+1, iend+1
                     bmaski(i, lb) = neighbour_extrapolated
                     coords1 (i,jv) = coords1 (i2,jv) + d(1)*(i-iend)
                     coords2 (i,jv) = coords2 (i2,jv) + d(2)*(i-iend)
                     end do
               endif
            endif
         endif ! if (i1 < iend+1)
!
         if (n <= 1) then

!        Extrapolate method point to exterior in rectangular direction
!cdir vector
!DEC$ IVDEP
            do i = ibeg+1, iend-1
            if (bmaski(i, lb) == no_value) then
               d(1) = (coords1 (i+1, j) - coords1 (i-1, j)) * 0.5
               d(2) = (coords2 (i+1, j) - coords2 (i-1, j)) * 0.5

               coords1 (i,jv) = coords1 (i, j) - d(2)*(j-jv)
               coords2 (i,jv) = coords2 (i, j) + d(1)
            endif
            end do
!
!        Extrapolate method point to exterior in rectangular direction
!
            i = ibeg
!
            d(1) = coords1 (i+1, j) - coords1 (i, j)
            d(2) = coords2 (i+1, j) - coords2 (i, j)
!
            if (bmaski(i, lb) == no_value) then
               coords1 (i,jv) = coords1 (i, j) - d(2)*(j-jv)
               coords2 (i,jv) = coords2 (i, j) + d(1)
            endif

            if (bmaskj (j, 1) == no_value) then
               coords1 (ibeg-1,j) = coords1 (i, j) - d(1)
               coords2 (ibeg-1,j) = coords2 (i, j) - d(2)
               bmaskj (j,1) = neighbour_extrapolated
            endif
!
!        Extrapolate method point to exterior in rectangular direction

            i = iend
            d(1) = coords1 (i, j) - coords1 (i-1, j)
            d(2) = coords2 (i, j) - coords2 (i-1, j)
!
            if (bmaski(i, lb) == no_value) then
               coords1 (i,jv) = coords1 (i, j) - d(2)*(j-jv)
               coords2 (i,jv) = coords2 (i, j) + d(1)
            endif

            if (bmaskj (j, 2) == no_value) then
               coords1 (iend+1, j) = coords1 (i, j) + d(1)
               coords2 (iend+1, j) = coords2 (i, j) + d(2)
               bmaskj (j,2) = neighbour_extrapolated
            endif
!
! corners (extrapolate virtual line)
!
            if (bmaski(ibeg-1, lb) == no_value) then
               d(1) = coords1 (ibeg+1, jv) - coords1 (ibeg, jv)
               d(2) = coords2 (ibeg+1, jv) - coords2 (ibeg, jv)
!
               coords1 (ibeg-1,jv) = coords1 (ibeg, jv) - d(1)
               coords2 (ibeg-1,jv) = coords2 (ibeg, jv) - d(2)
            endif

            if (bmaski(iend+1, lb) == no_value) then
               d(1) = coords1 (iend, jv) - coords1 (iend-1, jv)
               d(2) = coords2 (iend, jv) - coords2 (iend-1, jv)
!
               coords1 (iend+1,jv) = coords1 (iend, jv) + d(1)
               coords2 (iend+1,jv) = coords2 (iend, jv) + d(2)
            endif

!cdir vector
            do i = ibeg-1, iend+1
            if (bmaski(i, lb) == no_value) &
                bmaski(i, lb) = neighbour_extrapolated
            end do
         endif
!
         end do ! lb

! ??? Was ist mit den corners ??? TODO

         else
!           ibeg == iend, jbeg == jend

            i = ibeg

            d(1) = chmax1_corner (i, j) - chmin1_corner (i, j)
            d(2) = chmax2_corner (i, j) - chmin2_corner (i, j)

            if (bmaski(i, 1) == no_value) then
               coords1 (i, j-1) = coords1 (i, j)
               coords2 (i, j-1) = coords2 (i, j) - d(2)
               bmaski (i, 1) = neighbour_extrapolated
            endif

            if (bmaski(i, 2) == no_value) then
               coords1 (i, j+1) = coords1 (i, j)
               coords2 (i, j+1) = coords2 (i, j) + d(2)
               bmaski (i, 2) = neighbour_extrapolated
            endif

            if (bmaskj(j, 1) == no_value) then
               coords1 (i-1, j) = coords1 (i, j) - d(1)
               coords2 (i-1, j) = coords2 (i, j)
               bmaskj (j, 1) = neighbour_extrapolated
            endif

            if (bmaskj(j, 2) == no_value) then
               coords1 (i+1, j) = coords1 (i, j) + d(1)
               coords2 (i+1, j) = coords2 (i, j)
               bmaskj (j, 2) = neighbour_extrapolated
            endif

         endif
!
!     else if (coords_range(1,1) == coords_range(2,1)) then
      else if (ibeg == iend) then
!
         i = ibeg
!
!  Try to interpolate all points
!  Simple linear logical interpolation for this exotic case
!  Implement something better if required
!
           do lb = 1, 2
           if (lb == 1) then
              iv = i - 1
           else
              iv = i + 1
           endif
!
           n = 0
!
!cdir vector
               do j1 = jbeg-1, jend+1
               if (bmaskj(j1, lb) == neighbour_value) exit
               end do
!
            if (j1 < jend+1) then
               jfirst = j1
               j2 = j1 + 1
               n = 1
!
               do while (j2 < jend+1)
                  do j2 = j1+1, jend+1
                  if (bmaskj(j2, lb) == neighbour_value) exit
                  end do
!
               if (j2 > jend+1) exit

               n =  n + 1
!
               if (j2 == j1 + 1) then 
                  j1 = j2
                  cycle
               endif
!
               d(1) = (coords1 (iv, j2) - coords1 (iv, j1)) / (j2-j1)
               d(2) = (coords2 (iv, j2) - coords2 (iv, j1)) / (j2-j1)
!
!cdir vector
!DEC$ IVDEP
                  do j = j1+1, j2-1
                  bmaskj(j, lb) = neighbour_interpolated
                  coords1 (iv,j) = coords1 (iv,j1) + d(1)*(j-j1)
                  coords2 (iv,j) = coords2 (iv,j1) + d(2)*(j-j1)
                  end do
!
               j1 = j2
               end do ! while
!
            if (n > 1) then
               if (jfirst > jbeg-1) then
!
                  d(1) = coords1 (iv, jfirst+1) - coords1 (iv, jfirst)
                  d(2) = coords2 (iv, jfirst+1) - coords2 (iv, jfirst)
!cdir vector
!DEC$ IVDEP
                     do j = jbeg-1, jfirst-1
                     bmaskj(j, lb) = neighbour_extrapolated
                     coords1 (iv,j) = coords1 (iv,jfirst) - d(1)*(jfirst-j)
                     coords2 (iv,j) = coords2 (iv,jfirst) - d(2)*(jfirst-j)
                     end do
               endif
!
               if (j2 < jend+1) then
!
                  d(1) = coords1 (iv, j2) - coords1 (iv, j2-1)
                  d(2) = coords2 (iv, j2) - coords2 (iv, j2-1)
!cdir vector
!DEC$ IVDEP
                     do j = j2+1, jend+1
                     bmaskj(j, lb) = neighbour_extrapolated
                     coords1 (iv,j) = coords1 (iv,j2) + d(1)*(j-jend)
                     coords2 (iv,j) = coords2 (iv,j2) + d(2)*(j-jend)
                     end do
               endif
            endif
         endif ! if (j1 < jend+1)
!
         if (n <= 1) then

!        Extrapolate method point to exterior in rectangular direction
!cdir vector
!DEC$ IVDEP
            do j = jbeg+1, jend-1
            if (bmaskj(j, lb) == no_value) then
               d(1) = (coords1 (i, j+1) - coords1 (i, j-1)) * 0.5
               d(2) = (coords2 (i, j+1) - coords2 (i, j-1)) * 0.5

               coords1 (iv,j) = coords1 (i, j) - d(2)*(i-iv)
               coords2 (iv,j) = coords2 (i, j) + d(1)
            endif
            end do
!
!        Extrapolate method point to exterior in rectangular direction
!
            j = jbeg
!
            d(1) = coords1 (i, j+1) - coords1 (i, j)
            d(2) = coords2 (i, j+1) - coords2 (i, j)
!
            if (bmaskj(j, lb) == no_value) then
               coords1 (iv,j) = coords1 (i, j) - d(2)*(i-iv)
               coords2 (iv,j) = coords2 (i, j) + d(1)
            endif

            if (bmaski (i, 1) == no_value) then
               coords1 (i,jbeg-1) = coords1 (i, j) - d(1)
               coords2 (i,jbeg-1) = coords2 (i, j) - d(2)
               bmaski (i,1) = neighbour_extrapolated
            endif
!
!        Extrapolate method point to exterior in rectangular direction

            j = jend
            d(1) = coords1 (i, j) - coords1 (i, j-1)
            d(2) = coords2 (i, j) - coords2 (i, j-1)
!
            if (bmaskj(j, lb) == no_value) then
               coords1 (iv,j) = coords1 (i, j) - d(2)*(i-iv)
               coords2 (iv,j) = coords2 (i, j) + d(1)
            endif

            if (bmaski (i, 2) == no_value) then
               coords1 (i, jend+1) = coords1 (i, j) + d(1)
               coords2 (i, jend+1) = coords2 (i, j) + d(2)
               bmaski (i,2) = neighbour_extrapolated
            endif
!
! corners (extrapolate virtual line)
!
            if (bmaskj(jbeg-1, lb) == no_value) then
               d(1) = coords1 (iv, jbeg+1) - coords1 (iv, jbeg)
               d(2) = coords2 (iv, jbeg+1) - coords2 (iv, jbeg)
!
               coords1 (iv, jbeg-1) = coords1 (ibeg, jv) - d(1)
               coords2 (iv, jbeg-1) = coords2 (ibeg, jv) - d(2)
            endif

            if (bmaskj(jend+1, lb) == no_value) then
               d(1) = coords1 (iv, jend) - coords1 (iv, jend-1)
               d(2) = coords2 (iv, jend) - coords2 (iv, jend-1)
!
               coords1 (iv, jend+1) = coords1 (iv, jend) + d(1)
               coords2 (iv, jend+1) = coords2 (iv, jend) + d(2)
            endif

!cdir vector
            do j = jbeg-1, jend+1
            if (bmaskj(j, lb) == no_value) &
                bmaskj(j, lb) = neighbour_extrapolated
            end do
         endif
!
         end do ! lb
!
      else
!
!=======================================================================
!     Regular 2d array
!=======================================================================
!
!-----------------------------------------------------------------------
!     Interpolate missing virtual method cell points which are not
!     the corner points (i.e. from (ibeg:iend))
!     for faces 1 and 3 (j = jbeg and j = jend)
!-----------------------------------------------------------------------
!
         do lb = 1, 2
         n = Count (bmaski(ibeg:iend, lb) /= neighbour_value)
         if (n == 0) cycle

         Allocate (liste (n), STAT = ierror)
         if (ierror > 0) then
            ierrp (1) = ierror
            ierrp (2) = n

            ierror = PRISM_Error_Alloc

            call psmile_error ( ierror, 'liste', &
                                ierrp, 2, __FILE__, __LINE__ )
            return
         endif

         if (lb == 1) then 
!           Face 1 (j = jbeg)
!           j  = Index of the first interior method points/corner cell
!           j1 = Index of the next interior method point
!           jv = Index of the virtual method point
            j  = jbeg
            j1 = j + 1
            jv = j - 1
         else 
!           Face 2 (j = jend);
!           j = Index of last interior method point/corner cell
!           j1 = Index of the next interior method point
!           jv = Index of the virtual method point
            j  = jend
            j1 = j - 1
            jv = j + 1
         endif

!cdir vector
!DEC$ IVDEP
         do i = ibeg, iend
         if (bmaski(i, lb) == no_value) then
!        Extrapolate from interior method point to exterior
            bmaski  (i,lb) = neighbour_extrapolated
            coords1 (i,jv) = coords1 (i, j) - &
                            (coords1 (i, j1) - coords1 (i, j))
            coords2 (i,jv) = coords2 (i, j) - &
                            (coords2 (i, j1) - coords2 (i, j))
         endif
         end do
!
!   Get number and indices of inter/extrapolated points
!   which are contained in the cell box.
!
!   HR: Handle extrapolated and interpolated values in a
!       different way ?
!
         n = 0
         do i = ibeg, iend
         if (bmaski(i, lb) /= neighbour_value .and. &
            (coords1(i,jv) >= chmin1_corner (i,j) .and. &
             coords1(i,jv) <= chmax1_corner (i,j) .and. &
             coords2(i,jv) >= chmin2_corner (i,j) .and. &
             coords2(i,jv) <= chmax2_corner (i,j)) ) then
            n = n + 1
            liste (n) = i
         endif
         end do
!
!   ... and move these points out of the box
!
         do l = 1, n
         i = liste(n)
         d(1) = coords1(i,jv) - coords1 (i,j)
         d(2) = coords2(i,jv) - coords2 (i,j)
!
         if (d(1) /= 0.0d0) then
            r = max (chmax1_corner (i,j) - coords1 (i,j), &
                     coords1 (i,j) - chmin1_corner (i,j))
! dangerous
            fac = r / d(1);
         else
            fac = 1.0d0
         endif
!
         if (d(2) /= 0.0d0) then
            r = max (chmax2_corner (i,j) - coords2 (i,j), &
                     coords2 (i,j) - chmin2_corner (i,j))
! dangerous
            fac = max (fac, r/d(2))
         endif

         fac = fac * 1.0078125d0 ! Move it really outside
         coords1(i,jv) = coords1 (i,j) + fac*d(1)
         coords2(i,jv) = coords2 (i,j) + fac*d(2)
         end do ! l
!
         Deallocate (liste)
!
      end do ! lb
!
!-----------------------------------------------------------------------
!     Interpolate missing virtual method cell points which are not
!     the corner points (i.e. from (jbeg:jend))
!     for faces 2 and 4 (i = ibeg and i = iend)
!-----------------------------------------------------------------------
!
         do lb = 1, 2
         n = Count (bmaskj(jbeg:jend, lb) /= neighbour_value)
         if (n == 0) cycle

         Allocate (liste (n), STAT = ierror)
         if (ierror > 0) then
            ierrp (1) = ierror
            ierrp (2) = n

            ierror = PRISM_Error_Alloc

            call psmile_error ( ierror, 'liste', &
                                ierrp, 2, __FILE__, __LINE__ )
            return
         endif

         if (lb == 1) then 
!           Face 3 (i = ibeg)
!           i  = Index of the first interior method points/corner cell
!           i1 = Index of the next interior method point
!           iv = Index of the virtual method point
            i  = ibeg
            i1 = i + 1
            iv = i - 1
         else 
!           Face 4 (i = iend)
!           i  = Index of the last interior method points/corner cell
!           i1 = Index of the next interior method point
!           iv = Index of the virtual method point
            i  = iend
            i1 = i - 1
            iv = i + 1
         endif

!cdir vector
!DEC$ IVDEP
         do j = jbeg, jend
         if (bmaskj(j, lb) == no_value) then
!        Extrapolate from interior method point to exterior
            bmaskj  (j,lb) = neighbour_extrapolated
            coords1 (iv,j) = coords1 (i,  j) - &
                            (coords1 (i1, j) - coords1 (i, j))
            coords2 (iv,j) = coords2 (i,  j) - &
                            (coords2 (i1, j) - coords2 (i, j))
         endif
         end do
!
!   Get number and indices of inter/extrapolated points
!   which are contained in the cell box.
!
!   HR: Handle extrapolated and interpolated values in a
!       different way ?
!
         n = 0
         do j = jbeg, jend
         if (bmaskj(j, lb) /= neighbour_value .and. &
            (coords1(iv,j) >= chmin1_corner (i,j) .and. &
             coords1(iv,j) <= chmax1_corner (i,j) .and. &
             coords2(iv,j) >= chmin2_corner (i,j) .and. &
             coords2(iv,j) <= chmax2_corner (i,j)) ) then
            n = n + 1
            liste (n) = j
         endif
         end do
!
!   ... and move these points out of the box
!
         do l = 1, n
         j = liste(n)
         d(1) = coords1(iv,j) - coords1 (i,j)
         d(2) = coords2(iv,j) - coords2 (i,j)
!
         if (d(1) /= 0.0d0) then
            r = max (chmax1_corner (i,j) - coords1 (i,j), &
                     coords1 (i,j) - chmin1_corner (i,j))
! dangerous
            fac = r / d(1);
         else
            fac = 1.0d0
         endif
!
         if (d(2) /= 0.0d0) then
            r = max (chmax2_corner (i,j) - coords2 (i,j), &
                     coords2 (i,j) - chmin2_corner (i,j))
! dangerous
            fac = max (fac, r/d(2))
         endif

         fac = fac * 1.0078125d0 ! Move it really outside
         coords1(iv,j) = coords1 (i,j) + fac*d(1)
         coords2(iv,j) = coords2 (i,j) + fac*d(2)
         end do ! l
!
         Deallocate (liste)
         end do ! lb
!
!-----------------------------------------------------------------------
!  .. for corners of Faces
!     Extrapolate from the interior in diagonal direction
!
!     The virtual corner points to be determined are
!     (ibeg-1, jbeg-1), (ibeg-1, jend+1),
!     (iend-1, jbeg-1), (iend-1, jend+1)
!-----------------------------------------------------------------------

! Extrapolate in diagonal direction
! Corner of Faces 1 and 3

      if (bmaski(ibeg-1, 1) == no_value) then
!        ASSERT (bmaskj (jbeg-1,1) = no_value)
         bmaski (ibeg-1,1) = neighbour_extrapolated
         bmaskj (jbeg-1,1) = neighbour_extrapolated

         d(1) = coords1 (ibeg+1, jbeg+1) - coords1 (ibeg, jbeg)
         d(2) = coords2 (ibeg+1, jbeg+1) - coords2 (ibeg, jbeg)
         coords1 (ibeg-1,jbeg-1) = coords1 (ibeg, jbeg) - d(1)
         coords2 (ibeg-1,jbeg-1) = coords2 (ibeg, jbeg) - d(2)

         if (coords1(ibeg-1,jbeg-1) >= chmin1_corner (ibeg,jbeg) .and. &
             coords1(ibeg-1,jbeg-1) <= chmax1_corner (ibeg,jbeg) .and. &
             coords2(ibeg-1,jbeg-1) >= chmin2_corner (ibeg,jbeg) .and. &
             coords2(ibeg-1,jbeg-1) <= chmax2_corner (ibeg,jbeg)) then
!
            if (d(1) /= 0.0d0) then
               r = max (chmax1_corner (ibeg,jbeg) - coords1 (ibeg,jbeg), &
                        coords1 (ibeg,jbeg) - chmin1_corner (ibeg,jbeg))
! dangerous
               fac = r / d(1);
            else
               fac = 1.0d0
            endif
!
            if (d(2) /= 0.0d0) then
               r = max (chmax2_corner (ibeg,jbeg) - coords2 (ibeg,jbeg), &
                        coords2 (ibeg,jbeg) - chmin2_corner (ibeg,jbeg))
! dangerous
               fac = max (fac, r/d(2))
            endif

            fac = fac * 1.0078125d0 ! Move it really outside
            coords1(ibeg-1,jbeg-1) = coords1 (ibeg,jbeg) + fac*d(1)
            coords2(ibeg-1,jbeg-1) = coords2 (ibeg,jbeg) + fac*d(2)
         endif
      endif

! Corner of Faces 1 and 4

      if (bmaski(iend+1, 1) == no_value) then
!        ASSERT (bmaskj (jbeg-1,2) = no_value)
         bmaski (iend+1,1) = neighbour_extrapolated
         bmaskj (jbeg-1,2) = neighbour_extrapolated

         d(1) = coords1 (iend-1, jbeg+1) - coords1 (iend, jbeg)
         d(2) = coords2 (iend-1, jbeg+1) - coords2 (iend, jbeg)
         coords1 (iend+1,jbeg-1) = coords1 (iend, jbeg) - d(1)
         coords2 (iend+1,jbeg-1) = coords2 (iend, jbeg) - d(2)

         if (coords1(iend+1,jbeg-1) >= chmin1_corner (iend,jbeg) .and. &
             coords1(iend+1,jbeg-1) <= chmax1_corner (iend,jbeg) .and. &
             coords2(iend+1,jbeg-1) >= chmin2_corner (iend,jbeg) .and. &
             coords2(iend+1,jbeg-1) <= chmax2_corner (iend,jbeg)) then
!
            if (d(1) /= 0.0d0) then
               r = max (chmax1_corner (iend,jbeg) - coords1 (iend,jbeg), &
                        coords1 (iend,jbeg) - chmin1_corner (iend,jbeg))
! dangerous
               fac = r / d(1);
            else
               fac = 1.0d0
            endif
!
            if (d(2) /= 0.0d0) then
               r = max (chmax2_corner (iend,jbeg) - coords2 (iend,jbeg), &
                        coords2 (iend,jbeg) - chmin2_corner (iend,jbeg))
! dangerous
               fac = max (fac, r/d(2))
            endif

            fac = fac * 1.0078125d0 ! Move it really outside
            coords1(iend+1,jbeg-1) = coords1 (iend,jbeg) + fac*d(1)
            coords2(iend+1,jbeg-1) = coords2 (iend,jbeg) + fac*d(2)
         endif
      endif

! Corner of Faces 2 and 3

      if (bmaski(ibeg-1, 2) == no_value) then
!        ASSERT (bmaskj (jend+1,1) = no_value)
         bmaski (ibeg-1,2) = neighbour_extrapolated
         bmaskj (jend+1,1) = neighbour_extrapolated

         d(1) = coords1 (ibeg+1, jend-1) - coords1 (ibeg, jend)
         d(2) = coords2 (ibeg+1, jend-1) - coords2 (ibeg, jend)
         coords1 (ibeg-1,jend+1) = coords1 (ibeg, jend) - d(1)
         coords2 (ibeg-1,jend+1) = coords2 (ibeg, jend) - d(2)

         if (coords1(ibeg-1,jend+1) >= chmin1_corner (ibeg,jend) .and. &
             coords1(ibeg-1,jend+1) <= chmax1_corner (ibeg,jend) .and. &
             coords2(ibeg-1,jend+1) >= chmin2_corner (ibeg,jend) .and. &
             coords2(ibeg-1,jend+1) <= chmax2_corner (ibeg,jend)) then
!
            if (d(1) /= 0.0d0) then
               r = max (chmax1_corner (ibeg,jend) - coords1 (ibeg,jend), &
                        coords1 (ibeg,jend) - chmin1_corner (ibeg,jend))
! dangerous
               fac = r / d(1);
            else
               fac = 1.0d0
            endif
!
            if (d(2) /= 0.0d0) then
               r = max (chmax2_corner (ibeg,jend) - coords2 (ibeg,jend), &
                        coords2 (ibeg,jend) - chmin2_corner (ibeg,jend))
! dangerous
               fac = max (fac, r/d(2))
            endif

            fac = fac * 1.0078125d0 ! Move it really outside
            coords1(ibeg-1,jend+1) = coords1 (ibeg,jend) + fac*d(1)
            coords2(ibeg-1,jend+1) = coords2 (ibeg,jend) + fac*d(2)
         endif
      endif

! Corner of Faces 2 and 4

      if (bmaski(iend+1, 2) == no_value) then
!        ASSERT (bmaskj (jend+1,2) = no_value)
         bmaski (iend+1,2) = neighbour_extrapolated
         bmaskj (jend+1,2) = neighbour_extrapolated

         d(1) = coords1 (iend-1, jend-1) - coords1 (iend, jend)
         d(2) = coords2 (iend-1, jend-1) - coords2 (iend, jend)
         coords1 (iend+1,jend+1) = coords1 (iend, jend) - d(1)
         coords2 (iend+1,jend+1) = coords2 (iend, jend) - d(2)

         if (coords1(iend+1,jend+1) >= chmin1_corner (iend,jend) .and. &
             coords1(iend+1,jend+1) <= chmax1_corner (iend,jend) .and. &
             coords2(iend+1,jend+1) >= chmin2_corner (iend,jend) .and. &
             coords2(iend+1,jend+1) <= chmax2_corner (iend,jend)) then
!
            if (d(1) /= 0.0d0) then
               r = max (chmax1_corner (iend,jend) - coords1 (iend,jend), &
                        coords1 (iend,jend) - chmin1_corner (iend,jend))
! dangerous
               fac = r / d(1);
            else
               fac = 1.0d0
            endif
!
            if (d(2) /= 0.0d0) then
               r = max (chmax2_corner (iend,jend) - coords2 (iend,jend), &
                        coords2 (iend,jend) - chmin2_corner (iend,jend))
! dangerous
               fac = max (fac, r/d(2))
            endif

            fac = fac * 1.0078125d0 ! Move it really outside
            coords1(iend+1,jend+1) = coords1 (iend,jend) + fac*d(1)
            coords2(iend+1,jend+1) = coords2 (iend,jend) + fac*d(2)
         endif
      endif !  if (bmaski(iend+1, 2) == no_value)
!
      endif !  if (jbeg == jend)
!
!-----------------------------------------------------------------------
!  ...... Set bounding box of virtual method cells
!     for faces 1 and 3 (j = jbeg and j = jend)
!
!     The virtual corner cells are located in
!     (ibeg-1, jbeg-1), (ibeg-1, jend),
!     (iend,   jbeg-1), (iend,   jend)
!-----------------------------------------------------------------------
!
      do lb = 1, 2
      if (lb == 1) then 
!        Face 1 (j = jbeg)
!        j  = Index of the first interior method point/corner cell
!        j1 = Index of the virtual method point
!        jv = Index of the virtual cell
         j  = jbeg
         j1 = j - 1
         jv = j - 1
      else 
!        Face 2 (j = jend)
!        j  = Index of the last interior method point/corner cell
!        j1 = Index of the virtual method point
!        jv = Index of the virtual cell
         j  = jend
         j1 = j + 1
         jv = j
      endif
!
!     Note: This loop sets also the corner values
!           chmin(ibeg-1, jv) and chmin1(iend, jv)
!cdir vector
!DEC$ IVDEP
      do i = ibeg-1, iend
         chmin1 (i, jv) = min (coords1(i,j),  coords1(i+1,j), &
                               coords1(i,j1), coords1(i+1,j1))
         chmax1 (i, jv) = max (coords1(i,j),  coords1(i+1,j), &
                               coords1(i,j1), coords1(i+1,j1))
      enddo

!     Look for cyclic cells

!cdir vector
!DEC$ IVDEP
      do i = ibeg-1, iend
         ! if the bounding box is too big, we most probably need to apply cyclic coordinate handling
         if (chmax1(i,jv)-chmin1(i,jv) > period2) then

            ! get the cell
            cell(1) = coords1(i  ,j ); cell (2) = coords1(i+1,j )
            cell(3) = coords1(i  ,j1); cell (4) = coords1(i+1,j1)

            ! apply cyclic coordinate handling
            call psmile_transform_cell_cyclic(cell, period, ierror)

            ! get new bounding box
            chmin1 (i, jv) = minval (cell(:))
            chmax1 (i, jv) = maxval (cell(:))
         endif
      enddo
!
!        ... limit virtual cell volume by corner volumes
!        ... HR: Soll man das wirklich machen ?
!                Fuer die globale Suche waeren die echten Zellen besser.
!                Sonst nimmt mg_method vielleicht eine andere Zelle ?
!
#define LIMIT
#ifdef LIMIT
!cdir vector
!DEC$ IVDEP
      do i = ibeg, iend-1
         dmn = min (chmin1_corner(i,j), chmin1_corner (i+1,j))
         chmin1 (i, jv) = max (chmin1(i, jv), dmn)

         dmx = max (chmax1_corner(i,j), chmax1_corner (i+1,j))
         chmax1 (i, jv) = min (chmax1(i, jv), dmx)
      enddo

!      Corners of Face 3 (ibeg) and Face 4 (iend)

       chmin1 (ibeg-1, jv) = max (chmin1 (ibeg-1, jv), &
                                  chmin1_corner (ibeg, j))
       chmax1 (ibeg-1, jv) = min (chmax1 (ibeg-1, jv), &
                                  chmax1_corner (ibeg, j))

       chmin1 (iend,   jv) = max (chmin1 (iend, jv), &
                                  chmin1_corner (iend, j))
       chmax1 (iend,   jv) = min (chmax1 (iend, jv), &
                                  chmax1_corner (iend, j))
#endif
!
!cdir vector
!DEC$ IVDEP
          do i = ibeg-1, iend
          chmin2 (i, jv) = min (coords2(i,j),  coords2(i+1,j), &
                                coords2(i,j1), coords2(i+1,j1))
          chmax2 (i, jv) = max (coords2(i,j),  coords2(i+1,j), &
                                coords2(i,j1), coords2(i+1,j1))
          enddo
!
#ifdef LIMIT
!cdir vector
!DEC$ IVDEP
          do i = ibeg, iend-1
          dmn = min (chmin2_corner(i,j), chmin2_corner (i+1,j))
          chmin2 (i, jv) = max (chmin2(i, jv), dmn)

          dmx = max (chmax2_corner(i,j), chmax2_corner (i+1,j))
          chmax2 (i, jv) = min (chmax2(i, jv), dmx)
          enddo

!      Corners of Face 3 (ibeg) and Face 4 (iend)

      chmin2 (ibeg-1, jv) = max (chmin2 (ibeg-1, jv), &
                                 chmin2_corner (ibeg, j))
      chmax2 (ibeg-1, jv) = min (chmax2 (ibeg-1, jv), &
                                 chmax2_corner (ibeg, j))

      chmin2 (iend,   jv) = max (chmin2 (iend, jv), &
                                 chmin2_corner (iend, j))
      chmax2 (iend,   jv) = min (chmax2 (iend, jv), &
                                 chmax2_corner (iend, j))
#endif

      end do ! lb
!
!-----------------------------------------------------------------------
!  ...... Set bounding box of virtual method cells
!     for faces 3 and 4 (i = ibeg and i = iend-1)
!-----------------------------------------------------------------------
!
      do lb = 1, 2

      if (lb == 1) then 
!        Face 3 (i = ibeg)
!        i  = Index of the first interior method point/corner cell
!        j1 = Index of the virtual method point
!        iv = Index of the virtual cell
         i  = ibeg
         i1 = i - 1
         iv = i - 1
      else 
!        Face 4 (i = iend)
!        i  = Index of the last interior method point/corner cell
!        j1 = Index of the virtual method point
!        iv = Index of the virtual cell
         i  = iend
         i1 = i + 1
         iv = iend
      endif
!
!     Note: chmin1(iv, jbeg-1) and chmin1(iv, jend)
!           were already set in the loop above.
!
!cdir vector
!DEC$ IVDEP
      do j = jbeg, jend-1
         chmin1 (iv, j) = min (coords1(i, j), coords1(i, j+1), &
                               coords1(i1,j), coords1(i1,j+1))
         chmax1 (iv, j) = max (coords1(i, j), coords1(i, j+1), &
                               coords1(i1,j), coords1(i1,j+1))
      enddo

!     Look for cyclic cells

!cdir vector
!DEC$ IVDEP
      do j = jbeg, jend-1
         ! if the bounding box is too big, we most probably need to apply cyclic coordinate handling
         if (chmax1(iv,j)-chmin1(iv,j) > period2) then

            ! get the cell
            cell(1) = coords1(i  ,j  ); cell (2) = coords1(i1,j  )
            cell(3) = coords1(i  ,j+1); cell (4) = coords1(i1,j+1)

            ! apply cyclic coordinate handling
            call psmile_transform_cell_cyclic(cell, period, ierror)

            ! get new bounding box
            chmin1 (i, jv) = minval (cell(:))
            chmax1 (i, jv) = maxval (cell(:))
         endif
      enddo
!
!        ... limit virtual cell volume by corner volumes
!
#ifdef LIMIT
!cdir vector
!DEC$ IVDEP
      do j = jbeg, jend-1
         dmn = min (chmin1_corner(i,j), chmin1_corner (i,j+1))
         chmin1 (iv, j) = max (chmin1(iv, j), dmn)

         dmx = max (chmax1_corner(i,j), chmax1_corner (i,j+1))
         chmax1 (iv, j) = min (chmax1(iv, j), dmx)
      enddo
#endif
!
!cdir vector
!DEC$ IVDEP
         do j = jbeg, jend-1
            chmin2 (iv, j) = min (coords2(i,j),  coords2 (i,j+1), &
                                  coords2(i1,j), coords2(i1,j+1))
            chmax2 (iv, j) = max (coords2(i,j),  coords2 (i,j+1), &
                                  coords2(i1,j), coords2(i1,j+1))
         enddo
!
#ifdef LIMIT
!cdir vector
!DEC$ IVDEP
          do j = jbeg, jend-1
          dmn = min (chmin2_corner(i,j), chmin2_corner (i,j+1))
          chmin2 (iv, j) = max (chmin2(iv, j), dmn)

          dmx = max (chmax2_corner(i,j), chmax2_corner (i,j+1))
          chmax2 (iv, j) = min (chmax2(iv, j), dmx)
          enddo
#endif

      end do ! lb
!
!-----------------------------------------------------------------------
!  ... Set mid point of virtual method cells
!-----------------------------------------------------------------------
!
         j = jbeg-1
!cdir vector
!DEC$ IVDEP
            do i = ibeg-1, iend
            midp1  (i, j) = (chmin1(i,j) + chmax1 (i,j)) * 0.5
            midp2  (i, j) = (chmin2(i,j) + chmax2 (i,j)) * 0.5
            enddo
!
!        ... Move extrapolated points away
!            to give other points higher preference

!cdir vector
!DEC$ IVDEP
            do i = ibeg-1, iend
            if (bmaski(i,  1) == neighbour_extrapolated .or. &
                bmaski(i+1,1) == neighbour_extrapolated) then
               midp1  (i, j) = midp1 (i, j) + 720.0
            endif
            enddo
!
         j = jend
!cdir vector
!DEC$ IVDEP
            do i = ibeg-1, iend
            midp1  (i, j) = (chmin1(i,j) + chmax1 (i,j)) * 0.5
            midp2  (i, j) = (chmin2(i,j) + chmax2 (i,j)) * 0.5
            enddo

!cdir vector
!DEC$ IVDEP
            do i = ibeg-1, iend
            if (bmaski(i,  2) == neighbour_extrapolated .or. &
                bmaski(i+1,2) == neighbour_extrapolated) then
               midp1  (i, j) = midp1 (i, j) + 720.0
            endif
            enddo
!
         i = ibeg-1
!cdir vector
!DEC$ IVDEP
            do j = jbeg, jend-1
            midp1  (i, j) = (chmin1(i,j) + chmax1 (i,j)) * 0.5
            midp2  (i, j) = (chmin2(i,j) + chmax2 (i,j)) * 0.5
            enddo

!cdir vector
!DEC$ IVDEP
            do j = jbeg-1, jend
            if (bmaskj(j,  1) == neighbour_extrapolated .or. &
                bmaskj(j+1,1) == neighbour_extrapolated) then
               midp1  (i, j) = midp1 (i, j) + 720.0
            endif
            enddo
!
         i = iend
!cdir vector
!DEC$ IVDEP
            do j = jbeg, jend-1
            midp1  (i, j) = (chmin1(i,j) + chmax1 (i,j)) * 0.5
            midp2  (i, j) = (chmin2(i,j) + chmax2 (i,j)) * 0.5
            enddo

!cdir vector
!DEC$ IVDEP
            do j = jbeg-1, jend
            if (bmaskj(j,  2) == neighbour_extrapolated .or. &
                bmaskj(j+1,2) == neighbour_extrapolated) then
               midp1  (i, j) = midp1 (i, j) + 720.0
            endif
            enddo
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id), ierror

      call psmile_flushstd
#endif /* VERBOSE */
!
!  Formats:
!
9990 format (1x, a, ': psmile_bbcells_virt_2d_dble')
9980 format (1x, a, ': psmile_bbcells_virt_2d_dble: eof ierror =', i3)
!
      end subroutine PSMILe_Bbcells_virt_2d_dble