psmile_get_faces_gauss2_dble.F90

Go to the documentation of this file.

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

      subroutine psmile_get_faces_gauss2_dble (search, extra_search, &
                      corners1, corners2, corners3,                  &
                      corner_shape, nbr_corners, grid_valid_shape,   &
                      neighbors_3d, nloc, num_neigh,                 &
                      faces, nreq, ierror)
!
! !USES:
!
      use PRISM
!
      use PSMILe, dummy_interface => PSMILe_Get_faces_gauss2_dble
#ifdef DEBUG_TRACE
      use psmile_debug_trace
#endif

      Implicit none
!
! !INPUT PARAMETERS:
!
      Type (Enddef_search),     Intent (In) :: search

!     Info's on coordinates to be searched
!
      Type (Extra_search_info), Intent (In) :: extra_search
!
!     Number of locations where
!       (*) required mask values were not true
!
      Integer,            Intent (In) :: corner_shape (2, ndim_3d)
      
!     Specifies the dimension of arrays "corners1", "corners2" and "corners3"

      Integer,            Intent (In) :: nbr_corners

!     Number of corners in control volume 
!     Must be 8 for a grid of type "PRISM_Reglonlatvrt".

      Integer,            Intent (In) :: grid_valid_shape (2, ndim_3d)
      
!     Specifies the valid block shape for the "ndim_3d"-dimensional block

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

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

      Double Precision,   Intent (In) :: 
         corners3 ( corner_shape(1,3):corner_shape(2,3), 2) 

!
!     Coordinates of corners of cell/control volume
!     Note: This is a grid of type PRISM_Gaussreduced_Regvrt, i.e. 
!           2 corner coordinates per regular direction
!
      Integer,            Intent (In)    :: nloc
!
!     Total  Number of locations to be transferred
!
!     Last dimension of neighbors array "neighbors_3d" and
!     number of neighbors to be searched.

      Integer,            Intent (In)    :: num_neigh
!
!     Last dimension of neighbors array "neighbors_3d" and
!     number of neighbors to be searched.

      Integer,            Intent (In)    :: neighbors_3d (ndim_3d, nloc, 
                                                          num_neigh)

!     Indices of neighbor locations (interpolation bases)
!
      Integer,            Intent (In) :: nreq
!
!     Total number of extra points to be searched.
!     Dimension of "faces" and "nfaces".
!
! !OUTPUT PARAMETERS:
!
      Double Precision,  Intent (Out) :: faces (2, ndim_3d, nreq)

!     Number of block faces which cell "i" is coinciding
!
      Integer,           Intent (Out) :: ierror

!     Returns the error code of PSMILe_Get_faces_gauss2_dble;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !LOCAL VARIABLES
!
!     ... loop variables
!
      Integer                         :: i, j, k, n
!
!     ... for partitions
!
      Integer                         :: ipart, nprev
!
      Integer, Pointer                :: indices_req (:)
      Integer, Pointer                :: required (:)
      Integer, Pointer                :: len_req (:)
!
!     ... locations controlled
!
      Integer                         :: code
!
      Double Precision                :: box (2, ndim_3d)
!
!     ... for error handling
!
!     Integer, Parameter              :: nerrp = 3
!     Integer                         :: ierrp (nerrp)
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_Get_faces_gauss2_dble" computes
!
! Subroutine "PSMILe_Get_faces_gauss2_dble" is designed for grids of
! type "PRISM_Gaussreduced_Regvrt"
!
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 02.02.05    H. Ritzdorf  created
! 24.02.11    M. Hanke     changed neighbors_3d from 2d representation to 1d
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_get_faces_gauss2_dble.F90 2988 2011-02-24 14:36:19Z hanke $
!  $Author: hanke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_get_faces_gauss2_dble.F90 2988 2011-02-24 14:36:19Z hanke $'
!
!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
      print 9990, trim(ch_id)
!     print *, 'corner_shape    ', corner_shape
!     print *, 'grid_valid_shape', grid_valid_shape

      call psmile_flushstd
#endif /* VERBOSE */
!
#ifdef PRISM_ASSERTION
      if (grid_valid_shape (1,2) /= 1 .or. &
          grid_valid_shape (2,2) /= 1) then
         call psmile_assert (__FILE__, __LINE__, &
              "Internal definition of GaussReduced grid not correct")
      endif
!
      if (grid_valid_shape (1,1) /= corner_shape(1,1) .or. &
          grid_valid_shape (2,1) /= corner_shape(2,1) .or. &
          grid_valid_shape (1,2) /= corner_shape(1,2) .or. &
          grid_valid_shape (2,2) /= corner_shape(2,2)) then
!
         print *, 'corner_shape    ', corner_shape
         print *, 'grid_valid_shape', grid_valid_shape
!
         call psmile_assert (__FILE__, __LINE__, &
              "corner_shape /= grid_valid_shape; don't know to address")
      endif
#endif
!
      ierror = 0
!
      len_req => extra_search%len_req
!
!===> 
!
!
      faces (1, 1, :) = 1.0
      faces (2, 1, :) = 0.0
!
!  For all part "ipart"
!
!  indices_req = Indices of the extra points in entire list of all points
!                to be searched (i.e. for all parts "ipart")
!  required    = Code for the extra points which are additionally required.

      nprev = 0
      do ipart = 1, search%npart

         if (len_req (ipart) > 0) then
            indices_req => extra_search%indices_req(ipart)%vector
            required    => extra_search%required(ipart)%vector
!
            code = 1
!
            do j = 1, len_req(ipart)

               i = indices_req (j)

               do n = 1, num_neigh
#ifdef DEBUG_TRACE
                  if ( indices_req (j) == ictl ) then
                     print *, ' j, i, required         ', j, indices_req(j), required(j)
                     print *, ' neighs                 ', neighbors_3d (1,i,n)
                     print *, ' required -code         ', required (j), code, IAND (required (j), code)
                     print *, ' grid_valid_shape (2,1) ', grid_valid_shape (2,1)
                     code = code * 2
                  endif
#endif
                  if ( neighbors_3d (1,i,n) >= 1          .and. &
                       neighbors_3d (1,i,n) <= grid_valid_shape (2,1) ) then

#ifdef PRISM_ASSERTION
                     if (neighbors_3d (1,i,n) < 1 .or. &
                         neighbors_3d (1,i,n) > grid_valid_shape(2,1)) then
                        print *, "neighbors_3d (1,i,n)", neighbors_3d (1,i,n)
                        call psmile_assert (__FILE__, __LINE__, &
                           "Index out of range")
                     endif
#endif
                     if (neighbors_3d (3,i,n) < grid_valid_shape(1,3)) then
                        k = grid_valid_shape (1,3)
                     else
                        k = min (grid_valid_shape(2,3), neighbors_3d (3,i,n))
                     endif

!               HR: This was already calculated. This is the finest multi-grid.

                     box (1, 1) = minval(corners1 (neighbors_3d (1,i,n), :))
                     box (2, 1) = maxval(corners1 (neighbors_3d (1,i,n), :))

                     box (1, 2) = minval(corners2 (neighbors_3d (1,i,n), :))
                     box (2, 2) = maxval(corners2 (neighbors_3d (1,i,n), :))

                     box (1, 3) = minval(corners3 (k, :))
                     box (2, 3) = maxval(corners3 (k, :))

                     if (faces (1, 1, nprev+j) > faces (2, 1, nprev+j)) then

!                       ... first entry

                        faces (:, :, nprev+j) = box (:, :)
#ifdef DEBUG_TRACE
                        if ( indices_req (j) == ictl ) then
                           print *, ' face first entry '
                           print *, faces (:, :, nprev+j)
                        endif
#endif
                     else

                        faces (1, 1, nprev+j) = min (faces (1,1,nprev+j), box (1,1))
                        faces (2, 1, nprev+j) = max (faces (2,1,nprev+j), box (2,1))
                        faces (1, 2, nprev+j) = min (faces (1,2,nprev+j), box (1,2))
                        faces (2, 2, nprev+j) = max (faces (2,2,nprev+j), box (2,2))
                        faces (1, 3, nprev+j) = min (faces (1,3,nprev+j), box (1,3))
                        faces (2, 3, nprev+j) = max (faces (2,3,nprev+j), box (2,3))
#ifdef DEBUG_TRACE
                        if ( indices_req (j) == ictl ) then
                           print *, ' face following entries ', n
                           print *, faces (:, :, nprev+j)
                        endif
#endif
                     endif

                  endif ! neighbours_3d

               end do ! n

            end do ! j

         endif !( len_req (ipart) > 0)
         !
         nprev = nprev + len_req (ipart)
         !
      end do ! ipart
!
#ifdef PRISM_ASSERTION
      do j = 1, nreq
         if (faces (1,1,j) > faces (2,1,j)) exit
      end do
!
      if (j < nreq) then
         nprev = 0
         do ipart = 1, search%npart
            if (nprev + len_req(ipart) >= j) exit
            nprev = nprev + len_req(ipart)
         end do 
         
         indices_req => extra_search%indices_req(ipart)%vector
         i = indices_req (j)
         print *, 'i, j, nreq', i, j, nreq, grid_valid_shape
         do n = 1, num_neigh
            print *, 'i, n, ng', i, n, &
                              neighbors_3d (1,j,n), &
                              neighbors_3d (2,j,n), &
                              neighbors_3d (3,j,n)
         end do
         print *, ' x faces ', faces (:,1,j)
         print *, ' y faces ', faces (:,2,j)
         call psmile_assert (__FILE__, __LINE__, &
              "neighbour info is not consistent")
      endif
#endif
!
! Hier muss ich noch die faces auf -180:180, .... mappen
! oder 
!!
!!
!
!===> Muss ich hier nicht die Faces mit der cyclic behandeln
!     Oder mache ich das wieder dort, wo die Punkte gesucht werden.
!
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id), ierror

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

#ifdef VERBOSE

9990 format (1x, a, ': psmile_get_faces_gauss2_dble:')
9980 format (1x, a, ': psmile_get_faces_gauss2_dble: eof ierror =', i3)

#endif /* VERBOSE */

#ifdef DEBUG
#endif

      end subroutine PSMILe_Get_faces_gauss2_dble