psmile_gauss_setup.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_gauss_setup
!
! !INTERFACE:

      subroutine psmile_gauss_setup (grid_id, ierror)
!
! !USES:
!
      use prism_constants
      use psmile, dummy_interface => psmile_gauss_setup
      use psmile_grid_reduced_gauss

      implicit none
!
! !INPUT PARAMETERS:
!
      integer, intent (in)                :: grid_id
!
! !OUTPUT PARAMETERS:
!
      integer, intent (out)               :: ierror

!     Returns the error code of PSMILe_MG_first_level;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !LOCAL VARIABLES
!
      type(grid), pointer                 :: gp
!
!  ... for Gaussian reduced grids
!
      integer                             :: j, n, i
      integer                             :: leny, lenz
      type (block_info), pointer          :: local_block_info(:)
!
!  ... for error parameters
!
      integer, parameter                  :: nerrp = 2
      integer                             :: ierrp (nerrp)
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_gauss_setup" computes the auxiliary grid
! for search on reduced Gauss grids. In addition, it collect
! corner and point (method) data from neighbouring processes
! and builts the face and star arrays for reduced gauss grids,
! which are used to construct the stencils for linear and
! bicuibis interpolation.
!
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 05.12.02    R. Redler    created
! 09.01.11    M. Hanke     rewrote computation of aux. grid
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_gauss_setup.F90 3010 2011-03-10 13:26:49Z hanke $
!  $Author: hanke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_gauss_setup.F90 3010 2011-03-10 13:26:49Z hanke $'
!
!----------------------------------------------------------------------

#ifdef VERBOSE
      print *, trim(ch_id), ': psmile_gauss_setup: grid_id', grid_id
      call psmile_flushstd
#endif /* VERBOSE */
!
!  Initialization
!
      ierror = 0

      gp => Grids(grid_id)

! from mg_setup

      lenz = gp%grid_shape(2,3)-gp%grid_shape(1,3)+1

      if ( associated(gp%extent) ) then ! partitions are defined

         leny = size(gp%extent(:,1))

#ifdef DEBUG
         print *, trim(ch_id), ': psmile_gauss_setup: we have ', leny, ' partitions.'
#endif /* DEBUG */


      else ! we have to assume only one global grid partition

#ifdef DEBUG
         print *, trim(ch_id), ': psmile_gauss_setup: we assume to have only one global grid partition'
#endif /* DEBUG */

         leny = size(gp%nbr_points_per_lat)

         allocate(gp%extent(leny,2), gp%partition(leny,2), STAT = ierror)
         if ( ierror > 0 ) then
            ierrp (1) = ierror
            ierrp (2) = 2*leny
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'partition and extent', &
                 ierrp, 2, __FILE__, __LINE__ )
            return
         endif

         gp%extent(:,1) = gp%nbr_points_per_lat(:)
         gp%partition(1,1) = 0

         do n = 2, leny
            gp%partition(n,1) = sum(gp%nbr_points_per_lat(1:n-1))
         enddo

         gp%extent(:,2) = lenz
         gp%partition(:,2) = 0

      endif ! ( associated(gp%extent) )
      
      ! generate the auxiliary grid (used by the multigrid search instead of the reduced gauss grid)
      call psmile_gauss_gen_aux_grid (grid_id)
      ! generate a mapping between the auxiliary grid and the reduced gauss grid
      call psmile_gauss_gen_aux_grid_map (grid_id)

      Allocate (gp%gcorner_pointer, STAT = ierror)
      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = leny*2
         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'gcorner_pointer', &
              ierrp, 2, __FILE__, __LINE__ )
         return
      endif

      gp%gcorner_pointer%corner_shape = &
                        reshape ((/1,size (gp%reduced_gauss_data%aux_corners(1)%vector)/2,   &
                                   1,size (gp%reduced_gauss_data%aux_corners(2)%vector)/2,   &
                                   1,size (gp%reduced_gauss_data%aux_corners(3)%vector)/2/), &
                                 (/2,ndim_3d/))
      gp%gcorner_pointer%corner_datatype = gp%corner_pointer%corner_datatype
      gp%gcorner_pointer%corners_dble(1)%vector => gp%reduced_gauss_data%aux_corners(1)%vector
      gp%gcorner_pointer%corners_dble(2)%vector => gp%reduced_gauss_data%aux_corners(2)%vector
      gp%gcorner_pointer%corners_dble(3)%vector => gp%reduced_gauss_data%aux_corners(3)%vector
!
!===> Generate some additional data on the blocks of the local partition (used for optimisation)
!
      allocate (local_block_info(size (gp%partition, 1)))
      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = size (gp%partition, 1)
         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'local_block_info', &
              ierrp, 2, __FILE__, __LINE__ )
         return
      endif

      local_block_info(:)%global_1d_start_idx = gp%partition(:,1) + 1
      local_block_info(:)%global_1d_end_idx = gp%partition(:,1) + gp%extent(:,1)

      local_block_info(:)%local_1d_start_idx =                           &
         psmile_gauss_1d_global_to_local(grid_id, gp%partition(:,1) + 1, &
                                         size (gp%partition, 1),         &
                                         psmile_undef, .false.)

      local_block_info(:)%local_1d_end_idx =                          &
         psmile_gauss_1d_global_to_local(grid_id, gp%partition(:,1) + &
                                         gp%extent(:,1),              &
                                         size (gp%partition, 1),      &
                                         psmile_undef, .false.)

      do i = 1, size (gp%partition, 1)

         local_block_info(i)%global_3d_start_idx = &
            psmile_gauss_local_1d_to_3d(grid_id, local_block_info(i)%local_1d_start_idx)
            
         local_block_info(i)%global_3d_end_idx = &
            psmile_gauss_local_1d_to_3d(grid_id, local_block_info(i)%local_1d_end_idx)
      enddo ! i
      
      gp%reduced_gauss_data%local_block_info => local_block_info
      nullify (local_block_info)
!
!===> Define range of auxiliary grid (needed in tricu and trili interpolation)
!  nbr_lats = leny
!  points_per_lat = gp%extent
!
      Allocate (gp%g_irange(leny,2), STAT = ierror)
      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = leny * 2
         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'g_irange', &
              ierrp, 2, __FILE__, __LINE__ )
         return
      endif
!
      gp%g_irange(1,1) = 1
      gp%g_irange(1,2) = gp%extent(1,1)

!cdir vector
      do j = 2, leny
         gp%g_irange(j,1) = gp%g_irange(j-1,2) + 1
         gp%g_irange(j,2) = gp%g_irange(j-1,2) + gp%extent(j,1)
      enddo

      ! construct face and star array and collect corner information
      ! from neighbouring processes

      call psmile_gauss_get_neighbours ( grid_id, ierror )

      !
      !===> Free Memory
      !
      if ( associated(gp%get_list) ) then
         do n = 0, Comps(gp%comp_id)%size-1
            if (associated(gp%get_list(n)%vector)) &
                deallocate(gp%get_list(n)%vector)
         enddo
         deallocate (gp%get_list)
      endif

      if ( associated(gp%put_list) ) then
         do n = 0, Comps(gp%comp_id)%size-1
            if (associated(gp%put_list(n)%vector)) &
                deallocate(gp%put_list(n)%vector)
         enddo
         deallocate (gp%put_list)
      endif
!
!===> All done
!
#ifdef VERBOSE
      print *, trim(ch_id), ': psmile_gauss_setup eof: grid_id',&
                             grid_id, ', ierror =', ierror

      call psmile_flushstd
#endif /* VERBOSE */

      end subroutine psmile_gauss_setup