psmile_mg_srch_nneigh_reg_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_MG_srch_nneigh_reg_dble
!
subroutine psmile_mg_srch_nneigh_reg_dble (grid_id, nn_srch,    &
                arrays, search_mode, nref_3d, grid_valid_shape, &
                neighbors_3d, nloc, num_neigh,                  &
                sin_search, cos_search, z_search,               &
                dist_dble, dim1, indices, jbeg, jend,           &
                mask_array, mask_shape, mask_available,         &
                tol, ierror)
!
! !USES:
!
  use PRISM_constants
!
  use PSMILe, dummy_interface => PSMILe_MG_srch_nneigh_reg_dble

  implicit none
!
! !INPUT PARAMETERS:
!
      Integer, Intent (In)            :: grid_id
!
!     Info on the component in which the donor cells
!     should be searched.
!
      Type (Extra_search_nn), Intent (In) :: nn_srch
!
!     Info about the size of the coarse grid problem
!
      Type (Extra_search_dble)        :: arrays
!
!     Memory for metric information
!
      Integer, Intent (In)            :: nref_3d
!
!     Something like the maximum number of neighbor to be searched for
!
      Integer, Intent (In)            :: search_mode
!
!     Specifies the search mode for nearest neigbours with
!        search_mode = 3 : Full 3d-search
!        search_mode = 2 : Search in 2d-hyperplane (1st 2nd direction lon, lat)
!        search_mode = 1 : Search in 1d-plane
!
      Integer, Intent (In)            :: grid_valid_shape (2, ndim_3d)
!
!     Specifies the valid block shape for the "ndim_3d"-dimensional block

      Integer, Intent (In)            :: nloc
!
!     Second dimension of neighbors array "neighbors_3d" and
!     Total number of locations to be transferred
!
      Integer, Intent (In)            :: num_neigh
!
!     Number of neighbors to be searched.
!     Last dimension of neighbors array "neighbors_3d" and
!     number of neighbors to be searched.
!
      Integer,          Intent (In)   :: jbeg, jend
!
!     Shape of input arrays
!
      Integer, Parameter              :: lat = 2
      Double Precision, Intent (In)   :: sin_search (jbeg:jend, lat)
!
!     Sin values of the points for which the nearest neighbours
!     should be searched.
!
      Double Precision, Intent (In)   :: cos_search (jbeg:jend, lat)
!
!     Cos values of the points for which the nearest neighbours
!     should be searched.
!
      Double Precision, Intent (In)   :: z_search (jbeg:jend)
!
!     Z values of the points for which the nearest neighbours
!     should be searched.
!
      Integer, Intent (In)            :: dim1 (2)
!
!     Dimensions of "dist_dble" in first direction 
!
      Integer, Intent (In)            :: indices (:)
!
!     Indices of the extra points in entire list of all points to be
!     searched
!
      Integer, Intent (In)            :: mask_shape (2, ndim_3d)
!
!    Dimension of (method) mask array
!
      Logical, Intent (In)            :: 
         mask_array (mask_shape (1,1):mask_shape (2,1), 
                     mask_shape (1,2):mask_shape (2,2), 
                     mask_shape (1,3):mask_shape (2,3))
!
!    Mask of the method
!
      Logical, Intent (In)             :: mask_available
!
!    Is mask specified in array "mask_array" ?
!    mask_available = .false. : Mask is not available
!                     .true.  : Mask is     available
!
      Double Precision, Intent (In)    :: tol
!
!    Some tolerance may be used later to evaluate distances
!
! !INPUT/OUTPUT PARAMETERS:
!
      Double Precision, Intent (InOut) :: dist_dble (dim1(1):dim1(2), num_neigh)
!
!    Initial distances and final distances
!
      Integer,           Intent (Out) :: neighbors_3d (ndim_3d, nloc, num_neigh)

!     Indices of neighbor locations (interpolation bases)
!     Index array for points that have been found
!
! !OUTPUT PARAMETERS:
!
     Integer, Intent (Out)           :: ierror

!    Returns the error code of PSMILE_Neigh_nearx_sub_reg_dble;
!            ierror = 0 : No error
!            ierror > 0 : Severe error!
!
! !DEFINED PARAMETERS:
!
!  dble_earth = Earth radius in meter
!
  Double Precision, Parameter     :: dble_earth = 6400000.0
!
  Double Precision, Parameter     :: acosp1 =  1.0
  Double Precision, Parameter     :: acosm1 = -1.0

  Integer, Parameter              :: lon = 1
  Integer, Parameter              :: maxlen = 64
!
! !LOCAL VARIABLES
!
  Double Precision                :: dist, maxdist, val, fac, fac2

  Double Precision                :: sin_corner_lat(8,3)
  Double Precision                :: cos_corner_lat(8,3)
  Double Precision                :: sin_corner_lon(8,3)
  Double Precision                :: cos_corner_lon(8,3)

  Double Precision                :: sin_midp_lat(3)
  Double Precision                :: cos_midp_lat(3)
  Double Precision                :: sin_midp_lon(3)
  Double Precision                :: cos_midp_lon(3)

  Double Precision                :: distrad(maxlen), radius (maxlen)

  Logical                         :: maskrad(maxlen)

  Integer                         :: ijk0 (ndim_3d)
  Integer                         :: ijk  (2,ndim_3d)
  Integer                         :: ii, jj, kk
  Integer                         :: inext, jnext, knext
  Integer                         :: i, j, k, n
  Integer                         :: ncount
  Integer                         :: iloc(1)

  Type (Grid), Pointer            :: grid_info
  Type (Enddef_mg_double), Pointer :: my_arrays
!
!     ... multi-grid control
!
      Integer                         :: lev, nlev
      Integer, Allocatable            :: icoarse (:, :)
!
!     ... number of points in i, j, k directions
!
      Integer                         :: leni, lenij, lenijk, lenj, lenk
!
!     ... loop indices in neighbors_3d and compact vector
!
      Integer                         :: ind, jind
!
!     ... for error handling
!
!     Integer, Parameter              :: nerrp = 2
!     Integer                         :: ierrp (nerrp)
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_MG_srch_nneigh_reg_dble" searches the next "num_neigh"
! nearest neighbours on the method-grid. num_neigh = 1 for the moment only.
!
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 05/02/01    R. Redler    created
!
!EOP
!----------------------------------------------------------------------
!
!  Initialization

#ifdef VERBOSE
  print 9990, trim(ch_id)

  call psmile_flushstd
#endif /* VERBOSE */

  ierror = 0

  grid_info => Grids(grid_id)
  ijk0 = grid_info%ijk0

  nlev = grid_info%nlev - 3

#ifdef PRISM_ASSERTION

! This should never happen

  if ( nn_srch%leni*nn_srch%lenj*nn_srch%lenk > maxlen ) then
     print *, "lenijk, maxlen", nn_srch%leni*nn_srch%lenj*nn_srch%lenk, maxlen
     call psmile_assert ( __FILE__, __LINE__, "lenijk > maxlen")
  endif
#endif

!
!
! compute coarsening factor to get
! from level "lev" onto the finest level "1"
!
  if (mask_available) then
     allocate (icoarse (ndim_3d, nlev), stat = ierror)
!
     icoarse (1:ndim_3d, 1) = 1
!
        do lev = 2, nlev
           do i = 1, ndim_3d
           if ( grid_info%mg_infos(lev  )%levdim(i) /= &
                grid_info%mg_infos(lev-1)%levdim(i) ) then
              icoarse(i, lev) = icoarse(i, lev-1) * 2
           else
              icoarse(i, lev) = icoarse(i, lev-1)
           end if
           end do
        end do
   endif

! jind = loop index in compact vector "sin_search", "cos_search",
!        "z_search" and "dist_dble"
!
! ind  = Corresponding index in "neighbors_3d"
!
  do jind = jbeg, jend

     ! reinitialise lengths and next indices

     inext = 0
     jnext = 0
     knext = 0
!
! Get number of points in i,j,k direction on finest level
!
     my_arrays => grid_info%mg_infos(nlev)%double_arrays
!
     leni   = nn_srch%leni
     lenj   = nn_srch%lenj
     lenk   = nn_srch%lenk
!
     lenij  = leni*lenj
     lenijk = lenij*lenk

     ind = indices (jind)

     maxdist = maxval (dist_dble (jind, 1:num_neigh) )

     ! distance from point to mid point of control volume on start level 
     ! Das ist zu teuer (lenijk !)

     fac = dble_earth + z_search (jind)
!
     n = 0
        do j = 1, lenj

!cdir vector
           do i = 1, leni
           n = n+1

           val = (  arrays%sin_cmidp_lat%vector(j) * sin_search(jind,lat)    &
                 +  arrays%cos_cmidp_lat%vector(j) * cos_search(jind,lat)    &
                 * (arrays%cos_cmidp_lon%vector(i) * cos_search(jind,lon)    &
                 +  arrays%sin_cmidp_lon%vector(i) * sin_search(jind,lon)) )

           val = min (acosp1, val)
           val = max (acosm1, val)

           dist = fac * acos (val)
           distrad (n) = dist * dist
           enddo
        enddo
!
!       Replicate distance for other k levels
!
        do k = 2, lenk
        distrad ((k-1)*lenij+1:k*lenij) = distrad (1:lenij)
        enddo

        if ( search_mode == 3 ) then
           n = 0
              do k = 1, lenk
              distrad (n+1:n+lenij) = distrad (n+1:n+lenij) + &
                   (my_arrays%midp(3)%vector(k)-z_search(jind))**2

              n = n + lenij
              enddo
        endif

        distrad (1:lenijk) = sqrt (distrad (1:lenijk))

     ! distance from point to sphere
     ! Note: This distance may be negative.
     
     distrad(1:lenijk) = distrad(1:lenijk) - arrays%radius(1:lenijk)

     ! ignore control volumes that are farther away than what we have already

     maskrad (1:lenijk) = distrad (1:lenijk) < maxdist

     iloc(:) = minloc (distrad(1:lenijk), maskrad(1:lenijk))
#ifdef MINLOCFIX
     if (iloc(1).eq.0) iloc(1)=1
#endif /* MINLOCFIX */

     if ( iloc(1) > lenijk ) iloc(1) = 0

     do lev = nlev, 3, -1

        ! sphere around the control volume for all but the start level

        if ( lev < nlev ) then

        n = 0

        do k = 1, lenk
           fac2 = dble_earth + my_arrays%midp(3)%vector(k)
           do j = 1, lenj
!cdir vector
              do i = 1, leni
                 n = n+1

                 val = (  sin_corner_lat(1,j) * sin_midp_lat(j)    &
                       +  cos_corner_lat(1,j) * cos_midp_lat(j)    &
                       * (cos_corner_lon(1,i) * cos_midp_lon(i)    &
                       +  sin_corner_lon(1,i) * sin_midp_lon(i)) )

                 val = min (acosp1, val)
                 val = max (acosm1, val)

                 dist = fac2 * acos (val)
                 radius (n) = dist * dist

              enddo ! i
           enddo ! j
        enddo ! k

        if ( search_mode == 3 ) then
           n = 0

           do k = 1, lenk
           radius(n+1:n+lenij) = radius(n+1:n+lenij) + &
              (my_arrays%midp(3)%vector(k) - my_arrays%chmin(3)%vector(k))**2
           n = n + lenij
           enddo
        endif
!
        radius (1:lenijk) = sqrt (radius (1:lenijk))

           ! distance from point to mid point of control volume
           !  for all but the start level

        n = 0

           do j = 1, lenj
              do i = 1, leni
                 n = n+1

                 val = (  sin_midp_lat(j) * sin_search(jind,lat)    &
                       +  cos_midp_lat(j) * cos_search(jind,lat)    &
                       * (cos_midp_lon(i) * cos_search(jind,lon)    &
                       +  sin_midp_lon(i) * sin_search(jind,lon)) )

              val = min (acosp1, val)
              val = max (acosm1, val)

              dist = fac * acos (val)
              distrad (n) = dist * dist
           enddo
           enddo

           do k = 2, lenk
           distrad ((k-1)*lenij+1:k*lenij) = distrad (1:lenij)
           enddo

           if ( search_mode == 3 ) then
              n = 0
              do k = 1, lenk
                 distrad (n+1:n+lenij) = distrad (n+1:n+lenij) + &
                      (my_arrays%midp(3)%vector(k)-z_search(jind))**2
                 n = n + lenij
              end do
           endif

           distrad (1:lenijk) = sqrt (distrad (1:lenijk))

           ! distance from point to sphere

           distrad (1:lenijk) = distrad (1:lenijk) - radius(1:lenijk)

           ! ignore control volumes that are farther away than what we have already

           maskrad (1:lenijk) = distrad (1:lenijk) < maxdist

           iloc(:) = minloc(distrad(1:lenijk),maskrad(1:lenijk))
#ifdef MINLOCFIX
           if (iloc(1).eq.0) iloc(1)=1
#endif /* MINLOCFIX */

           if ( iloc(1) > lenijk ) iloc(1) = 0

        endif ! (lev < nlev)
!
!-------------------------------------------------------------------------
!    Ignore points maked out
!-------------------------------------------------------------------------
!
        if ( mask_available .and. iloc(1) > 0 ) then

           ! check whether we have valid points in this volume

           do i = 1, lenijk

              kk = (iloc(1)                - 1) / lenij + 1
              jj = (iloc(1) - (kk-1)*lenij - 1) / leni  + 1
              ii =  iloc(1) - (kk-1)*lenij - (jj-1)*leni

              ii = ii - 1 + inext
              jj = jj - 1 + jnext
              kk = kk - 1 + knext

              ijk(1,1) = ijk0(1) + ii*icoarse(1,lev)
              ijk(1,2) = ijk0(2) + jj*icoarse(2,lev)
              ijk(1,3) = ijk0(3) + kk*icoarse(3,lev)

              ijk(2,1) = min (grid_valid_shape(2,1), ijk(1,1)+icoarse(1,lev)-1)
              ijk(2,2) = min (grid_valid_shape(2,2), ijk(1,2)+icoarse(2,lev)-1)
              ijk(2,3) = min (grid_valid_shape(2,3), ijk(1,3)+icoarse(3,lev)-1)

              ncount = count( mask_array(ijk(1,1):ijk(2,1), &
                                         ijk(1,2):ijk(2,2), &
                                         ijk(1,3):ijk(2,3)) )
#ifdef DEBUGX
              print *, '     masks:                     ', maskrad(1:lenijk)
              print *, '     ii, jj, kk:                ', ii, jj, kk
              print *, '     lev, i, iloc(1), ncount : ', lev, i, iloc(1), ncount
              print *, '     => ijk                     ', ijk
#endif

              if ( ncount > 0 ) exit

              ! reject control volume without any valid points

              maskrad (iloc(1)) = .false.

              iloc(:) = minloc (distrad(1:lenijk), maskrad(1:lenijk))
#ifdef MINLOCFIX
              if (iloc(1).eq.0) iloc(1)=1
#endif /* MINLOCFIX */

              if ( iloc(1) > lenijk ) then
                 iloc(1) = 0
                 exit
              endif

!             if ( iloc(1) == 0 ) exit
              
           enddo

           if ( ncount == 1 ) then
              do k = ijk(1,3),ijk(2,3)
                 do j = ijk(1,2),ijk(2,2)
                    do i = ijk(1,1),ijk(2,1)
                       if ( mask_array(i,j,k) ) then
                          neighbors_3d (1, ind, 1) = i
                          neighbors_3d (2, ind, 1) = j
                          neighbors_3d (3, ind, 1) = k
                       endif
                    enddo
                 enddo
              enddo

! ist das korreckt ?
           if ( lev < nlev ) then
              dist_dble (jind, 1) = distrad (iloc(1)) + radius(iloc(1))
           else
              dist_dble (jind, 1) = distrad (iloc(1)) + arrays%radius(iloc(1))
           endif 

              iloc(1) = 0 ! do not enter level 3

              exit ! lev

           endif

        endif

        if ( iloc(1) == 0 ) exit
!
!----------------------------------------------------------------------
!       Prepare next level
!----------------------------------------------------------------------
!
        if ( lev > 3 ) then 

           ! prepare next level

           inext = ii * (icoarse(1,lev)/icoarse(1,lev-1))
           jnext = jj * (icoarse(2,lev)/icoarse(2,lev-1))
           knext = kk * (icoarse(3,lev)/icoarse(3,lev-1))

           my_arrays => grid_info%mg_infos(lev-1)%double_arrays

           leni   = min(inext+2,grid_info%mg_infos(lev-1)%levdim(1)) - inext + 1
           lenj   = min(jnext+2,grid_info%mg_infos(lev-1)%levdim(2)) - jnext + 1
           lenk   = min(knext+2,grid_info%mg_infos(lev-1)%levdim(3)) - knext + 1
           lenij  = leni*lenj
           lenijk = lenij*lenk

#ifdef PRISM_ASSERTION
!HR ich halte die Assertion fuer falsch
           if ( lenijk > maxlen ) then
              print *, "lenijk, maxlen", lenijk, maxlen
              call psmile_assert ( __FILE__, __LINE__, "lenijk > maxlen")
           endif
#endif
        sin_corner_lon(1,1:leni) = sin(my_arrays%chmin(1)%vector(inext+1:inext+leni)*dble_deg2rad)
        cos_corner_lon(1,1:leni) = cos(my_arrays%chmin(1)%vector(inext+1:inext+leni)*dble_deg2rad)
        sin_corner_lat(1,1:lenj) = sin(my_arrays%chmin(2)%vector(jnext+1:jnext+lenj)*dble_deg2rad)
        cos_corner_lat(1,1:lenj) = cos(my_arrays%chmin(2)%vector(jnext+1:jnext+lenj)*dble_deg2rad)

        sin_midp_lon(1:leni) = sin(my_arrays%midp(1)%vector(inext+1:inext+leni)*dble_deg2rad)
        cos_midp_lon(1:leni) = cos(my_arrays%midp(1)%vector(inext+1:inext+leni)*dble_deg2rad)
        sin_midp_lat(1:lenj) = sin(my_arrays%midp(2)%vector(jnext+1:jnext+lenj)*dble_deg2rad)
        cos_midp_lat(1:lenj) = cos(my_arrays%midp(2)%vector(jnext+1:jnext+lenj)*dble_deg2rad)

        endif

     enddo ! lev

!=========================================================================
! go over all points of the final selected control volume on level 3
! and find the nearest points.
!=========================================================================
     
     if ( iloc(1) > 0 ) then

        my_arrays => grid_info%mg_infos(1)%double_arrays

        inext  = ijk(1,1) - ijk0(1)
        jnext  = ijk(1,2) - ijk0(2)
        knext  = ijk(1,3) - ijk0(3)

        leni   = ijk(2,1) - ijk(1,1) + 1
        lenj   = ijk(2,2) - ijk(1,2) + 1
        lenk   = ijk(2,3) - ijk(1,3) + 1

        lenij  = leni*lenj
        lenijk = lenij*lenk

#ifdef PRISM_ASSERTION
!HR ich halte die Assertion fuer falsch
        if ( lenijk > maxlen ) then
           print *, "lenijk, maxlen", lenijk, maxlen
           call psmile_assert ( __FILE__, __LINE__, "lenijk > maxlen")
        endif
#endif

     arrays%sin_fmidp_lon%vector(1:leni) = sin(my_arrays%midp(1)%vector(inext+1:inext+leni)*dble_deg2rad)
     arrays%cos_fmidp_lon%vector(1:leni) = cos(my_arrays%midp(1)%vector(inext+1:inext+leni)*dble_deg2rad)
     arrays%sin_fmidp_lat%vector(1:lenj) = sin(my_arrays%midp(2)%vector(jnext+1:jnext+lenj)*dble_deg2rad)
     arrays%cos_fmidp_lat%vector(1:lenj) = cos(my_arrays%midp(2)%vector(jnext+1:jnext+lenj)*dble_deg2rad)

     n = 0

        do j = 1, lenj
!cdir vector
           do i = 1, leni

              n = n+1

           ! distance from point to potential neighbor

              val = (  arrays%sin_fmidp_lat%vector(j) * sin_search(jind,lat)    &
                    +  arrays%cos_fmidp_lat%vector(j) * cos_search(jind,lat)    &
                    * (arrays%cos_fmidp_lon%vector(i) * cos_search(jind,lon)    &
                    +  arrays%sin_fmidp_lon%vector(i) * sin_search(jind,lon)) )

           val = min (acosp1, val)
           val = max (acosm1, val)

           dist = fac * acos (val)

           distrad(n) = dist*dist
        end do
        enddo

     do k = 2, lenk
        distrad ((k-1)*lenij+1:k*lenij) = distrad (1:lenij)
     enddo

        if ( search_mode == 3 ) then
           n = 0
           do k = 1, lenk
              distrad (n+1:n+lenij) = distrad (n+1:n+lenij) + &
                                   (my_arrays%midp(3)%vector(k)-z_search(jind))**2
              n = n + lenij
           end do
        endif

        ! ignore invalid points
        if ( mask_available ) then
        n = 0
           do k = 1, lenk
           do j = 1, lenj
!cdir vector
              do i = 1, leni
              n = n + 1

                 maskrad (n) = mask_array(ijk(1,1)+i-1, &
                                          ijk(1,2)+j-1, &
                                          ijk(1,3)+k-1) 
              end do
           end do
        end do

           i = min(num_neigh, count(maskrad(1:lenijk)))
        else
           maskrad (1:lenijk) = .true.

           i = min(num_neigh, lenijk)
        endif

        do n = 1, i

           iloc(:) = minloc (distrad(1:lenijk), maskrad(1:lenijk))
#ifdef MINLOCFIX
           if (iloc(1).eq.0) iloc(1)=1
#endif /* MINLOCFIX */

           kk = (iloc(1)                - 1) / lenij + 1
           jj = (iloc(1) - (kk-1)*lenij - 1) / leni  + 1
           ii =  iloc(1) - (kk-1)*lenij - (jj-1) * leni

           neighbors_3d (1, ind, n) = ijk(1,1)+ii-1
           neighbors_3d (2, ind, n) = ijk(1,2)+jj-1
           neighbors_3d (3, ind, n) = ijk(1,3)+kk-1

! ist das korrekt ?
           dist_dble (jind, n) = sqrt ( distrad(iloc(1)) )

           maskrad (iloc(1)) = .false.

           end do

     endif
  end do ! jind
 
  !
  !===> All done
  !
  if (mask_available) then
     Deallocate (icoarse)
  endif

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

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

end subroutine PSMILe_MG_srch_nneigh_reg_dble