psmile_mg_final_prev_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_final_prev_2d_real
!
! !INTERFACE:

      subroutine psmile_mg_final_prev_2d_real (grid_id, nlev,        &
                      lstijk, xyz, found, newijk, nc,                &
                      corners1, corners2, corner_shape, nbr_corners, &
                      tol, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_MG_final_prev_2d_real

      Implicit none
!
! !INPUT PARAMETERS:
!
      Integer, Intent (In)            :: grid_id
!
!     Grid Id for which the donor cells should be searched.
!
      Integer, Intent (In)            :: nlev
!
!     Total number of levels
!
      Integer, Intent (In)            :: nc
!
!     Number of points to be searched.
!
      Integer, Intent (In)            :: lstijk (ndim_2d, nc)
!
!     Last indices of cell on finest level "lev == 1".
!     !!! This routine assumes that all parts on the corrsponding
!     coarse grid cells on level "LEV+1" were already controlled and
!     point "xyz" was not found.
!     Note: "lstijk" is relative to "ijk0"
!
      Real, Intent (In)               :: xyz (ndim_2d, nc)
!
!     Coordinates of the points to be searched
!
      Integer, Intent (In)            :: corner_shape (2, ndim_2d)

!     Dimension of corner coordinate arrays

      Integer, Intent (In)            :: nbr_corners

!     Number of corners
!     Note: This routine assumes nbr_corners == 4.

      Real,   Intent (In)             ::                 
         corners1 ( corner_shape(1,1):corner_shape(2,1), 
                    corner_shape(1,2):corner_shape(2,2), nbr_corners)

      Real,   Intent (In)             ::                 
         corners2 ( corner_shape(1,1):corner_shape(2,1), 
                    corner_shape(1,2):corner_shape(2,2), nbr_corners)

!     Corner coordinates
!
      Real, Intent (In)               :: tol
!
!     Absolute tolerance for search of "identical" points
!     TOL >= 0.0
!
! !OUTPUT PARAMETERS:
!
      Integer, Intent (Out)           :: found (nc)
!
!     Was the point xyz found ?
!     found = 0: Not found
!     found > 0: Part of the cell where point XYZ was found.
!
      Integer, Intent (Out)           :: newijk (ndim_2d, nc)
!
!     Indices on level LEV if found > 0 is returned.
!
      Integer, Intent (Out)           :: ierror

!     Returns the error code of PSMILE_mg_final_prev_2d_real;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !DEFINED PARAMETERS:
!
!  NTET = Number of triangles used to determine whether
!         point is within the grid cell
!  EPS  = Tolerance for determinate of deformed grid cells
!
      Integer, Parameter              :: ntet = 3
!
      Real, Parameter                 :: eps = 1.0d-9
      Real, Parameter                 :: rmxeps = 1.0d30
      Real, Parameter                 :: dzero = 0.0
!
      Integer, Parameter              :: lev  = 1
      Integer, Parameter              :: levc = 2
      Logical, Parameter              :: wide = .true.
      Logical, Parameter              :: small = .false.
!
! !LOCAL VARIABLES
!
      Integer                         :: i, n

!     ... for Grid

      Type (Enddef_mg), Pointer       :: mg_infos (:)
      Integer, Pointer                :: ijk0 (:)

!     ... for levels

      Integer                         :: levd, levdbg, levu, levubg
      Integer                         :: nlev1
      Integer                         :: nold (nc, nlev-1)
      Integer                         :: ijk (ndim_2d, nc, nlev-1)
      Integer                         :: ignore (ndim_2d, nlev-1)
!
!     ... for tolerances
!
      Real                            :: tol1
!
!     ... for locations controlled
!
      Integer                         :: j, m, nt
!
!     ... Logical reference point and edges (constant)
! 
      Integer, Save                   :: ijkedg (2, ntet)
!
!     ... Reference point and edges and distances
! 
      Integer                         :: ind (ndim_2d)
      Real                            :: xyzpnt (ntet, ndim_2d), 
                                         xyzref (ntet, ndim_2d)
      Real                            :: xyzedg (ntet, ndim_2d, 2)
      Real                            :: d1, d2
      Real                            :: det (ntet)
      Real                            :: dnorm
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_MG_final_prev_2d_real" searches for the donor cell
! of points "xyz (*, 1:nc)". The points was last found on final level
! in donor cell "lstijk (*, 1:nc)". 
! If the i-th point is found in another cell of the real grid,
! "found(i)" > 0 is returned and
! "newijk (*, i)" returns the donor cell index on final level.
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 12.03.07    H. Ritzdorf  created
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_mg_final_prev_2d_real.F90,v 1.1.2.2 2008/12/08 18:36:56 ritzdorf Exp $
!  $Author: ritzdorf $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_mg_final_prev_2d_real.F90,v 1.1.2.2 2008/12/08 18:36:56 ritzdorf Exp $'
!
!----------------------------------------------------------------------
!
!  Relative indices of triangles
!
!     The cell is split into 3 triangles
!
      data ((ijkedg (j, n), j = 1, 2), n = 1, ntet) &
           / 2, 4,  3, 1,  4, 2/
!             1      2      3

!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
      print 9990, trim(ch_id), nc

      call psmile_flushstd
#endif /* VERBOSE */
!
      ierror = 0
!
      mg_infos => Grids(grid_id)%mg_infos
      ijk0     => Grids(grid_id)%ijk0

      nlev1 = nlev - 1
      tol1 = tol + 1.0
!
#ifdef PRISM_ASSERTION
      if (levc >= nlev) then
         call psmile_assert (__FILE__, __LINE__, 'incorrect level')
      endif
!
      if (tol < 0.0) then
         call psmile_assert (__FILE__, __LINE__, &
                             "tol must be >= 0.0")
      endif
!
!cdir vector
         do i = 1, nc
         if (lstijk(1,i)+ijk0(1) < corner_shape(1,1) .or. &
             lstijk(1,i)+ijk0(1) > corner_shape(2,1) .or. &
             lstijk(2,i)+ijk0(2) < corner_shape(1,2) .or. &
             lstijk(2,i)+ijk0(2) > corner_shape(2,2) ) then
            exit
         endif
         end do
!
      if (i <= nc) then
         print *, 'i, lstijk', i, lstijk (:, i), ijk0
         print *, 'corner_shape', corner_shape
         call psmile_assert (__FILE__, __LINE__, &
                             'wrong index')
      endif
#endif
!
         do levd = levc, nlev1
         nold (:, levd) = 0 ! not clear to which part "lstijk" belongs
         end do
!
!===> ... Compute indices on the coarser levels
!         for coarsening factor 2 !
!
         do i = 1, ndim_2d
         if (mg_infos(lev )%levdim(i) /= &
             mg_infos(levc)%levdim(i) ) then
!           Note: Calling routines have controlled cells
!                (lstijk()-1:lstijk()+1)
!           Thus index on level levc cannot lie on a coarse grid line
            ijk (:, :, levc) = lstijk (:, :) / 2
         else
            ijk (:, :, levc) = lstijk (:, :)
         endif
         end do
!
         do levd = levc+1, nlev1
            do i = 1, ndim_2d
            if (mg_infos(levd  )%levdim(i) /= &
                mg_infos(levd-1)%levdim(i) ) then
!
!              make sure that ijk lies on coarse grid line
!
               ijk (i, :, levd) = (ijk (i, :, levd-1) / 4) * 2
            else
               ijk (i, :, levd) = ijk (i, :, levd-1)
            endif
            end do
         end do
!
!---------------------------------------------------------------------
!     Go down in the level hierarchie
!
!     TODO: transfer nc loop into PSMILe_mg_control_cell_2d_real
!           Use mask for points to be controlled and/or
!           gather/scatter of points to be controlled.
!
!     Note: The "wide" control may cause some additional work
!           (repeated control of fine grid cells) 
!           at 2**i boundaries.
!           Otherwise, the "wide" control may find the correct
!           cell faster instead of looping till to the coarsest level.
!           It's not totally clear, whether "small" control would
!           be better should.
!           Should performance control done with example
!           "simple-orca-atmland" since it generates many controls.
!           Possibly, adaptive setting of "wide" or "small" control.
!---------------------------------------------------------------------
!
loop_n: do n = 1, nc

      levdbg = levc

!     Define cells to be ignored
      ignore (:, lev)  = lstijk (:, n)
!
         do i = 1, ndim_2d
            do levd = levc, nlev1
            if (mg_infos(levd-1)%levdim(i) > &
                mg_infos(levd  )%levdim(i)) then
               ignore (i, levd) = ignore (i, levd-1) / 2
            else
               ignore (i, levd) = ignore (i, levd-1)
            endif
            end do
         end do
!
100   continue
!
         do levd = levdbg, nlev1
!
         call psmile_mg_control_cell_2d_real (                   &
                   mg_infos(levd)%real_arrays%chmin(1)%vector, &
                   mg_infos(levd)%real_arrays%chmin(2)%vector, &
                   mg_infos(levd)%real_arrays%chmax(1)%vector, &
                   mg_infos(levd)%real_arrays%chmax(2)%vector, &
                   mg_infos(levd)%real_arrays%midp (1)%vector, &
                   mg_infos(levd)%real_arrays%midp (2)%vector, &
                   mg_infos(levd)%levdim,                        &
                   ijk(1, n, levd), xyz(:,n), nold(n, levd),     &
                   ignore (1, levd), wide, found (n), newijk (:, n))
!
         if (found (n) > 0) then
            go to 1990
         endif
!
         end do
!
!===> Not found
!
      cycle loop_n
!
!---------------------------------------------------------------------
!     Up again; startpoint is newijk on level "levd"
!---------------------------------------------------------------------
!
1990  nold (n, levd) = found (n)
!
      levdbg = levd
      levubg = levd - 1
!
      nold (n, levubg) = 0
!
         do i = 1, ndim_2d
         if (mg_infos(levubg  )%levdim(i) > &
             mg_infos(levubg+1)%levdim(i)) then
            ijk (i, n, levubg) = min (newijk(i,n)*2, &
                                      mg_infos(levubg)%levdim(i))
         else
            ijk (i, n, levubg) = newijk (i,n)
         endif
         enddo
!
1995  continue
!
         do levu = levubg, lev, -1
!
2000     continue
!
         call psmile_mg_control_cell_2d_real (                   &
                   mg_infos(levu)%real_arrays%chmin(1)%vector, &
                   mg_infos(levu)%real_arrays%chmin(2)%vector, &
                   mg_infos(levu)%real_arrays%chmax(1)%vector, &
                   mg_infos(levu)%real_arrays%chmax(2)%vector, &
                   mg_infos(levu)%real_arrays%midp (1)%vector, &
                   mg_infos(levu)%real_arrays%midp (2)%vector, &
                   mg_infos(levu)%levdim,                        &
                   ijk(1, n, levu), xyz(:,n), nold(n, levu),     &
                   ignore (1, levu), small, found(n), newijk (:,n))
!
         if (found(n) == 0) then
            if (levu < levd-1) then
!
!===> ... go back to the previous level; and try again
!
               levubg = levu + 1
               go to 1995
            else
!
!===> ... go down
!
               go to 100
            endif
!
         else if (levu > lev) then
!           Not the finest level
            nold (n, levu) = found (n)
!
!           Initialize search on next finer level "levu-1"
!
            nold (n, levu-1) = 0
!
               do i = 1, ndim_2d
               if (mg_infos(levu-1)%levdim(i) > &
                   mg_infos(levu  )%levdim(i)) then
                  ijk (i, n, levu-1) = min (newijk(i,n)*2, &
                                            mg_infos(levu-1)%levdim(i))
               else
                  ijk (i, n, levu-1) = newijk (i,n)
               endif
               enddo
!
         else ! levu == lev
!           Finest level
!-----------------------------------------------------------------------
!           Point "n" was found in box newijk (*, n)
!           Is the point contained in the real corner cell ?
!-----------------------------------------------------------------------
!
!===>       Compute reference points, corners and edges of transformed
!           triangles.
!
               ind (:) = newijk (:, n) + ijk0 (1:ndim_2d)
!
               do nt = 1, ntet
               xyzref (nt, 1) = corners1 (ind(1), ind(2), nt)
               xyzref (nt, 2) = corners2 (ind(1), ind(2), nt)
!
               xyzpnt (nt, 1) = xyz (1, n) - xyzref (nt, 1)
               xyzpnt (nt, 2) = xyz (2, n) - xyzref (nt, 2)
!
                  do m = 1, 2
                  xyzedg (nt, 1, m) = corners1 (ind(1), ind(2), &
                                                ijkedg(m, nt)) &
                                    - xyzref (nt, 1)

                  xyzedg (nt, 2, m) = corners2 (ind(1), ind(2), &
                                                ijkedg(m, nt)) &
                                    - xyzref (nt, 2)
                  end do
               end do ! nt
!
!===> Solve bilinear system; compute spat product
!
!   (xyzedg(*,1), xyzedg(*,2)) (x) = xyzpnt
!
            det(:) = xyzedg(:, 1,1)*xyzedg(:, 2,2) &
                   - xyzedg(:, 1,2)*xyzedg(:, 2,1)
!
#define ALLOW_DEGRADED_CELL
#ifdef ALLOW_DEGRADED_CELL
!
! ??? ist das koscher ?? wird d1(nt), d2(nt), d3(nt) == 0
!
               do nt = 1, ntet
               dnorm = abs(xyzedg(nt,1,1)) + abs(xyzedg(nt,2,1)) + &
                       abs(xyzedg(nt,1,2)) + abs(xyzedg(nt,2,2))

               if (abs(det(nt)) < eps*dnorm*dnorm .or. dnorm == dzero) then
                  det (nt) = sign (rmxeps, det(nt))
               else
                  det (nt) = 1.0 / det (nt)
               endif
               end do ! nt

#else 

loop_n:        do nt = 1, ntet
               dnorm = abs(xyzedg(nt,1,1)) + abs(xyzedg(nt,2,1)) + &
                       abs(xyzedg(nt,1,2)) + abs(xyzedg(nt,2,2))

               if (abs(det(i, nt)) < eps*dnorm*dnorm .or. &
                   dnorm(i,nt) == dzero) then
                  det (i, nt) = dnorm * 1e-25
! eigentlich exit loop_n und print_degraded cell
! ist bei vektorisierung verschwunden
! sollte man in part_list suchen
               endif
               end do ! nt
!
#ifdef PRINT_DEGRADED_CELL
            if (nt < ntet) then
!
!===> cell is degraded
!
               print 9960, trim(ch_id), grid_id, newijk(:, n), &
                        (corners1(ind(1), ind(2), nt),         &
                         corners2(ind(1), ind(2), nt), nt = 1, nbr_corners)
!
                  do nt = 1, ntet
                  if (abs(det(nt)) < eps) then
                     print 9950, nt, det (nt), dnorm, &
                       xyzedg (nt,1, 1), xyzedg (nt,1, 2), &
                       xyzedg (nt,2, 1), xyzedg (nt,2, 2)
                  endif
                  end do
!
               nicht gefunden fnd (i) = .false. versucht, durch
               ganz kleine Determinante zu loesen
            endif
#endif /* PRINT_DEGRADED_CELL */
!
            det (:) = 1.0 / det (:)
!
#endif /* ALLOW_DEGRADED_CELL */
!
!===>          Is the point contained in a triangle ?
!
                do nt = 1, ntet
                d1 = (xyzpnt(nt,1)*xyzedg(nt,2,2)- &
                      xyzpnt(nt,2)*xyzedg(nt,1,2)) * det (nt)
                d2 = (xyzpnt(nt,2)*xyzedg(nt,1,1)- &
                      xyzpnt(nt,1)*xyzedg(nt,2,1)) * det (nt)
!
                if (min(d1, d2) >= -tol .and. &
                    max(d1,dzero) + max(d2,dzero) <= tol1) exit
                end do ! ntet
!
            if (nt > ntet) then
!                 Point was not found in the real cell.
!                 Search again on same level.
!
               nold (n, levu) = found (n)
!
               go to 2000
            endif
         endif
         end do ! levu
!
         end do loop_n
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id)

      call psmile_flushstd
#endif /* VERBOSE */
!
      return
!
!  Formats:
!

9990 format (1x, a, ': psmile_mg_final_prev_2d_real: nc', i5)
9980 format (1x, a, ': psmile_mg_final_prev_2d_real: eof')

#ifdef PRINT_DEGRADED_CELL
9960 format (1x, a, ': psmile_MG_final_prev_2d_real: grid id =', i5, &
            /1x,  ' Grid cell (', i4, ',', i4,                    &
                 ') is degraded! Cell coordinates are:', 1p,      &
            4 (/1x, '(', 1 (e17.9, ','), e17.9, ')'))
9950 format (1x, i1, '-th spawn; det =', 1p, e17.9, ' (', e17.9, &
             ') of', 2 (/1x, 2e17.9, ','))
#endif

      end subroutine psmile_mg_final_prev_2d_real