psmile_mg_method_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_mg_method_gauss2_dble
!
! !INTERFACE:
      subroutine psmile_mg_method_gauss2_dble (method_id, search_range,     &
                                               tgt_shape, tgt_coords_x,     &
                                               tgt_coords_y, loc_fnd_shape, &
                                               found, loc, virtual_cell,    &
                                               ierror)
!
! !USES:
!
      use PRISM_constants
      use psmile, dummy => psmile_mg_method_gauss2_dble
      use psmile_geometry, only : point_is_in_quadrangle
      use psmile_grid, only : common_grid_range, psmile_transform_cell_cyclic
      use psmile_grid_reduced_gauss
#ifdef DEBUG_TRACE
      use psmile_debug_trace
#endif

      implicit none
!
! !INPUT PARAMETERS:
!
      integer, intent (in)            :: method_id
!
      integer, intent (in)            :: search_range (2, ndim_3d)

!
      integer, intent (in)            :: loc_fnd_shape (2, ndim_3d)

!     Dimension of loc and found

      integer, intent (in)            :: tgt_shape (2, ndim_3d)

!     Dimension of coordinate arrays "coord1", "coord2"
!     which contain the coordinates to be searched.
!
      double precision, intent (in)   :: tgt_coords_x (tgt_shape(1,1):tgt_shape(2,1),  
                                                       tgt_shape(1,2):tgt_shape(2,2),  
                                                       tgt_shape(1,3):tgt_shape(2,3)), 
                                         tgt_coords_y (tgt_shape(1,1):tgt_shape(2,1),  
                                                       tgt_shape(1,2):tgt_shape(2,2),  
                                                       tgt_shape(1,3):tgt_shape(2,3))

!     Coordinates to be searched
!
! !INPUT/OUTPUT PARAMETERS:
!
      integer, intent (inout)         :: found (loc_fnd_shape(1,1):loc_fnd_shape(2,1), 
                                                loc_fnd_shape(1,2):loc_fnd_shape(2,2), 
                                                loc_fnd_shape(1,3):loc_fnd_shape(2,3))

!     Finest level number on which a grid cell was found for point i,j,k.
!     Input :
!     Level number < -1: Point was not found and
!                        and last level number was (-found(i,j,k))
!     Level number = -(nlev+1): Never found
!     Output :
!     found(i,j,k) = +1: Point (i,j,k) is located on a point
!                        of the method grid.
!     found(i,j,k) = -1: Point (i,j,k) is located in a cell
!                        of the method grid.
!     found(,j,k) = -(nlev+2): Not located in the method grid.

      integer, intent (inout)         :: loc   (loc_fnd_shape(1,1):loc_fnd_shape(2,1), 
                                                loc_fnd_shape(1,2):loc_fnd_shape(2,2), 
                                                loc_fnd_shape(1,3):loc_fnd_shape(2,3))
!
!     Indices of the grid cell in which the point was found.
!     The indices are relative to "shape".
!     Note: The indices are indices in the Gauss Reduced Grid.
!
! !OUTPUT PARAMETERS:
!
      integer, intent (out)           :: virtual_cell (loc_fnd_shape(1,1):loc_fnd_shape(2,1), 
                                                       loc_fnd_shape(1,2):loc_fnd_shape(2,2), 
                                                       loc_fnd_shape(1,3):loc_fnd_shape(2,3))

!     Is the searched point (i,j,k) located in a virtual cell ?
!     virtual_cell (i,j,k) returns a code for the location in
!     the virtual cell (outside the local method grid)
!     virtual_cell (i,j,k) = 0 : Not located in a virtual cell
!     virtual_cell (i,j,k) > 0 :     Located in a virtual cell
!                          = 1 : (-1,-1) direction
!                          = 2 : ( 0,-1) direction
!                          = 3 : ( 1,-1) direction
!                          = 4 : (-1, 0) direction
!                          = 6 : ( 1, 0) direction
!                          = 7 : (-1, 1) direction
!                          = 8 : ( 0, 1) direction
!                          = 9 : ( 1, 1) direction

      integer, intent (out)           :: ierror

!     Returns the error code of psmile_mg_method_gauss2_dble;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !DEFINED PARAMETERS:
!
!  val_coupler = Code for locations which should to be transferred
!                to the coupler process.
!
      integer, parameter              :: val_coupler = -1
!
! !LOCAL VARIABLES
!
      integer, pointer                :: local_1d_stencil_lookup(:,:)

      integer                         :: i, j, k, n, n_end

      integer                         :: not_found

      double precision, pointer       :: src_coords_x(:), 
                                         src_coords_y(:) ! coordinates of the method points
                                                         ! these can also contain data of
                                                         ! remote points

      integer, pointer                :: local_coords_shape(:,:)
      integer                         :: num_local_points, num_remote_points, num_global_lats

      integer, pointer                :: curr_stencil(:)
      integer                         :: curr_stencil_idx
      integer                         :: curr_3d_base_idx(ndim_3d)
      type(point_dble)                :: method_cell (4), point
      double precision                :: grid_extent
      integer                         :: grid_id

      integer, allocatable            :: loc_3d(:,:,:,:)

      integer, parameter              :: n_to_stencil_base_idx(9) = (/1,2,3, 5,6,7, 
                                                                      9,10,11/)
      integer, parameter              :: stencil_idx_to_virtual_cell(16) = (/1,2,3,-1, 
                                                                             4,5,6,-1, 
                                                                             7,8,9,-1, 
                                                                             -1,-1,-1,-1/)
!
#ifdef DEBUG_TRACE
!
!     ... for debugging
!
      logical                         :: debug_trace
      type(point_dble)                :: debug_point
      character (len=*), parameter    :: out_prefix = "ictl psmile_mg_method_gauss2_dble:"
#endif
!
! !DESCRIPTION:
!
! Subroutine "psmile_mg_method_gauss2_dble" gets a list of target points and the matching
! local source cell (grid cells). The routine searches for corresponding method cells. If
! a found method cell is not on the local process, it is indicated through the virtual_cell
! array
!
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 03.07.05    H. Ritzdorf  created
! 06.01.11    M. Hanke     rewritten
! 21.02.11    M. Hanke     introduced module psmile_grid_reduced_gauss
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_mg_method_gauss2_dble.F90 3011 2011-03-10 17:50:02Z hanke $
!  $Author: hanke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_mg_method_gauss2_dble.F90 3011 2011-03-10 17:50:02Z hanke $'

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

      call psmile_flushstd
#endif /* VERBOSE */
!
!  Initialization
!
      ierror = 0

      grid_id = methods(method_id)%grid_id
      local_1d_stencil_lookup => &
         grids(grid_id)%reduced_gauss_data%local_1d_stencil_lookup
      num_global_lats = &
         size (grids(grid_id)%reduced_gauss_data%nbr_points_per_lat)
      not_found = - (grids(grid_id)%nlev + 2)
      local_coords_shape => Methods(method_id)%coords_pointer%coords_shape

      ! get the number of remote point that have a local copy
      num_local_points  = local_coords_shape(2,1) - local_coords_shape(1,1)
      if (associated (grids(grid_id)%remote_index)) then
         num_remote_points = size (grids(grid_id)%remote_index)
      else
         num_remote_points = 0
      endif

      ! set the local source coordinate pointer
      ! if there are no remote points that have local copy
      if (num_remote_points <= 0) then

         src_coords_x => Methods(method_id)%coords_pointer%coords_dble(1)%vector
         src_coords_y => Methods(method_id)%coords_pointer%coords_dble(2)%vector
      else

         ! allocate memory of an array that contains the local and remote source coordinate  data
         allocate (src_coords_x(local_coords_shape(1,1):local_coords_shape(2,1)+num_remote_points), &
                   src_coords_y(local_coords_shape(1,1):local_coords_shape(2,1)+num_remote_points), &
                   stat = ierror)

         if ( ierror > 0 ) then
            call psmile_error (PRISM_Error_Alloc, 'src_coords_x and src_coords_y', (/ierror, &
                               2*(local_coords_shape(2,1)-local_coords_shape(1,1)+1+num_remote_points)/), &
                               2, __FILE__, __LINE__ )
            return
         endif

         src_coords_x(local_coords_shape(1,1):local_coords_shape(2,1)) = &
            Methods(method_id)%coords_pointer%coords_dble(1)%vector
         src_coords_y(local_coords_shape(1,1):local_coords_shape(2,1)) = &
            Methods(method_id)%coords_pointer%coords_dble(2)%vector

         src_coords_x(local_coords_shape(2,1)+1:local_coords_shape(2,1)+num_remote_points) = &
             Methods(method_id)%gauss2_dble(1)%vector
         src_coords_y(local_coords_shape(2,1)+1:local_coords_shape(2,1)+num_remote_points) = &
             Methods(method_id)%gauss2_dble(2)%vector
      endif ! (num_remote_points == 0)

      ! initialise with default value (no virtual cell)
      virtual_cell(:,:,:) = 0

      ! this should normally be 360
      grid_extent = common_grid_range(2,1) - common_grid_range(1,1)
      
      ! some optimisation
      allocate (loc_3d(ndim_3d, loc_fnd_shape(1,1):loc_fnd_shape(2,1), &
                                loc_fnd_shape(1,2):loc_fnd_shape(2,2), &
                                loc_fnd_shape(1,3):loc_fnd_shape(2,3)))

      loc_3d(:,:,:,:) = reshape (psmile_gauss_local_1d_to_3d(grid_id,                     &
                                                             reshape(loc, (/size(loc)/)), &
                                                             size(loc)),                  &
                                 (/ndim_3d, size(loc,1), size(loc,2), size(loc,3)/))
      ! for all tgt points
      do i = search_range(1, 1), search_range(2,1)
         do j = search_range(1, 2), search_range(2,2)
            do k = search_range(1, 3), search_range(2,3)
#ifdef DEBUG_TRACE
               if (i == ictl_ind(1) .and. j == ictl_ind(2)) then
                  print *, out_prefix, ' ictl_ind', &
                     ictl_ind (1:ndim_2d), ' old found = ', &
                     found (ictl_ind(1),ictl_ind(2),search_range(1,3)), &
                     ', old loc =', loc   (ictl_ind(1),ictl_ind(2),search_range(1,3))
                  debug_trace = .true.
               else
                  debug_trace = .false.
               endif
#endif
               ! if we have a valid location for the given point
               if (found (i,j,k) == 1) then

                  point%x = tgt_coords_x (i,j,k)
                  point%y = tgt_coords_y (i,j,k)

                  curr_stencil => local_1d_stencil_lookup(:, loc (i,j,k))
                  curr_3d_base_idx(:) = loc_3d (:,i,j,k)

#ifdef DEBUG_TRACE
                  if (debug_trace) then
                     print '(2(1x,a),2f8.2)', out_prefix, 'point to be searched: ', point
                     debug_point = point
                  endif
#endif

                  ! skip the southern points if the base point is already in the most southern latitude
                  if (curr_3d_base_idx(2) == 1) then
                     n = 4
                  else
                     n = 1
                  endif

                  ! skip the northern points if we are 1 row below the most northern row
                  if (curr_3d_base_idx(2) == num_global_lats - 1) then
                     n_end = 6
                  else
                     n_end = 9
                  endif

                  ! for all points in the current stencil of the src point associated with the current tgt point
                  do n = n, n_end

                     curr_stencil_idx = n_to_stencil_base_idx(n)

                     ! generate method cell
                     method_cell = get_method_cell (curr_stencil_idx, curr_stencil, &
                                                    curr_3d_base_idx(2))

                     ! move the point into the x range of the cell
                     do while (point%x < minval (method_cell(:)%x))
                        point%x = point%x + grid_extent
                     enddo
                     do while (point%x > maxval (method_cell(:)%x))
                        point%x = point%x - grid_extent
                     enddo

#ifdef DEBUG_TRACE
                     if (debug_trace) then
                        if (debug_point%x /= point%x) then
                           print '(2(1x,a),/,2f8.2)', out_prefix, &
                                 "point coordinates were adjusted to method cell range:", point
                           debug_point = point
                        endif
                     endif
#endif

                     ! check point in cell
                     if (point_is_in_quadrangle (point, method_cell)) then

                        found (i,j,k) = val_coupler
                        exit
                     endif
                  enddo ! n

                  ! if a suitable method cell was found
                  if (found (i,j,k) == val_coupler) then
#ifdef DEBUG_TRACE
                     if (debug_trace) then
                        print *, out_prefix, ' found matching cell'
                     endif
#endif
                     ! if the point is located in the local grid
                     if (curr_stencil(curr_stencil_idx) <= local_coords_shape(2,1)) then

                        loc (i,j,k) = curr_stencil(curr_stencil_idx)
#ifdef DEBUG_TRACE
                        if (debug_trace) then
                           print *, out_prefix, ' new loc =', loc (i,j,k)
                        endif
#endif
                     else

                        virtual_cell (i,j,k) = &
                           stencil_idx_to_virtual_cell(curr_stencil_idx)
#ifdef DEBUG_TRACE
                        if (debug_trace) then
                           print *, out_prefix, ' cell is virtual cell at stencil idx', &
                                    curr_stencil_idx, 'virt. cell idx', virtual_cell (i,j,k)
                        endif
#endif
                     endif ! (loc (i,j,k) > local_coords_shape(2,1))
                  else
                     found (i,j,k) = not_found
                  endif ! (found (i,j,k) == val_coupler)
               endif ! (found (i,j,k) == 1)
            enddo ! k = search_range(1, 3), search_range(2,3)
         enddo ! j = search_range(1, 2), search_range(2,2)
      enddo ! i = search_range(1, 1), search_range(2,1)

      ! free memory if necessary
      deallocate (loc_3d)
      if (num_remote_points > 0) deallocate (src_coords_x, src_coords_y)

#ifdef VERBOSE
      print 9980, trim(ch_id)

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

      contains

!-----------------------------------------------------------------------------------

         ! this function generates a method cell based on the point index in the
         ! given stencil array of the lower left corner of the cell
         function get_method_cell(idx, stencil, g_latitude_idx)

            integer, intent (in) :: idx, stencil(16), g_latitude_idx

            type (point_dble)    :: get_method_cell(4)

            type (point_dble)    :: method_cell(4)
            integer, parameter   :: neighbours(3,16) = 
                                       reshape ((/ 2, 6, 5,  3, 7, 6,  4, 8, 7, -1,-1,-1,   
                                                   6,10, 9,  7,11,10,  8,12,11, -1,-1,-1,   
                                                  10,14,13, 11,15,14, 12,16,15, -1,-1,-1,   
                                                  -1,-1,-1, -1,-1,-1, -1,-1,-1, -1,-1,-1/), 
                                                (/3,16/)) ! (see definition of the stencil in
                                                          !  psmile_grid_reduced_gauss)
            integer, parameter   :: latitude_off(16) = (/-1,-1,-1, 0, 
                                                          0, 0, 0, 0, 
                                                          1, 1, 1, 0, 
                                                          0, 0, 0, 0/)

#ifdef DEBUG_TRACE
            if (debug_trace) then
               print *, out_prefix, ' stencil     =', stencil
               print *, out_prefix, ' stencil idx =', idx
               print *, out_prefix, ' generating method cell with indices:', &
                        stencil(idx), stencil(neighbours(:,idx))
            endif
#endif
            ! get method point of current point
            method_cell(1)%x = src_coords_x(stencil(idx))

            ! get method point of right neighbour
            ! (see definition of stencil in psmile_grid_reduced_gauss)
            method_cell(2)%x = src_coords_x(stencil(neighbours(1,idx)))

            ! get method point of right-upper neighbour
            method_cell(3)%x = src_coords_x(stencil(neighbours(2,idx)))

            ! get method point of upper neighbour
            method_cell(4)%x = src_coords_x(stencil(neighbours(3,idx)))

            ! method cell is extended to the pole
            if (g_latitude_idx + latitude_off(idx) == 1) then

               method_cell(1:2)%y = dsign (90.d0, src_coords_y(stencil(idx)))

#ifdef DEBUG_TRACE
               if (debug_trace) then
                  print *, out_prefix, ' extending method cell to pole'
               endif
#endif
            else

               method_cell(1)%y = src_coords_y(stencil(idx))
               method_cell(2)%y = src_coords_y(stencil(neighbours(1,idx)))
            endif

            ! method cell is extended to the pole
            if (g_latitude_idx + latitude_off(idx) == num_global_lats - 1) then

               method_cell(3:4)%y = dsign (90.d0, src_coords_y(stencil(idx)))

#ifdef DEBUG_TRACE
               if (debug_trace) then
                  print *, out_prefix, ' extending method cell to pole'
               endif
#endif
            else

               method_cell(3)%y = src_coords_y(stencil(neighbours(2,idx)))
               method_cell(4)%y = src_coords_y(stencil(neighbours(3,idx)))
            endif

#ifdef DEBUG_TRACE
            if (debug_trace) then
               print "(2(1x,a),4(/,2f8.2))", out_prefix, 'generated method cell:', &
                                             method_cell
            endif
#endif

            ! apply cyclic coordinate handling
            call psmile_transform_cell_cyclic(method_cell(:)%x, grid_extent, ierror)

#ifdef DEBUG_TRACE
            if (debug_trace) then
               print "(2(1x,a),4(/,2f8.2))", &
                     out_prefix, 'after cyclic transformation:', method_cell
            endif
#endif

            get_method_cell = method_cell

         end function get_method_cell

      end subroutine psmile_mg_method_gauss2_dble