psmile_mg_setup.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_setup
!
! !INTERFACE:

      subroutine psmile_mg_setup (grid_id, range, tol, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_MG_setup

      implicit none

!
! !INPUT PARAMETERS:
!
      Integer, Intent (In)               :: grid_id

!     Grid Id for which the multigrid sequence is generated

      Integer, Intent (In)               :: range (2, ndim_3d)

!     Range for which MG sequence should be setup

      Double precision, Intent (In)      :: tol

!     Tolerance for search
!
! !OUTPUT PARAMETERS:
!
      Integer, Intent (Out)              :: ierror
!
!     Returns the error code of PSMILe_MG_setup;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !LOCAL VARIABLES
!
!  ... for components
!
      integer                            :: comp_id

!  ... for grids

      integer                            :: lev, levf, ndlev
      integer                            :: levbeg
      integer, allocatable               :: levdim (:, :)
      integer                            :: icoarse (ndim_3d)
      Integer                            :: my_range (2, ndim_3d)

      integer                            :: j, len, n
!
      logical                            :: simplified_grid

      Type(Grid), Pointer                :: gp

!  ... for error parameters

      integer, parameter                 :: nerrp = 2
      integer                            :: ierrp (nerrp)
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_MG_setup" sets up the multigrid sequence for
! searching the neighbourhood information. The multigrid sequence is
! setup for the grid with grid id "grid_id" and restricted to subgrid
! "range".
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 03.06.20    H. Ritzdorf  created
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_mg_setup.F90 3016 2011-03-14 15:17:02Z hanke $
!  $Autor$
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_mg_setup.F90 3016 2011-03-14 15:17:02Z hanke $'
!
!----------------------------------------------------------------------

#ifdef VERBOSE
      print 9990, trim(ch_id), Grids(grid_id)%comp_id, range

      call psmile_flushstd
#endif /* VERBOSE */

#ifdef PRISM_ASSERTION
      if (Associated (Grids(grid_id)%mg_infos)) then
         call psmile_assert (__FILE__, __LINE__, &
                             "Multigrid info's already created for this grid")
      endif
#endif
!
!  Initialization
!
      ierror = 0
      gp => Grids(grid_id)
      my_range = range
!
      comp_id = gp%comp_id
!     gp%cyclic = .false.
!
! Define maximal number of multigrid levels
! ndlev = 1: Fertig ?
!

!
!  Special section for Gaussian reduced grids
!
!     The assumption is that each latitude band is described as a block on its own
!     and that there is only one block per latitude.
!
      simplified_grid = gp%grid_type == PRISM_Gaussreduced_regvrt
!
      if ( gp%grid_type == PRISM_Gaussreduced_regvrt ) then

         my_range(1:2,1:2) = gp%gcorner_pointer%corner_shape(1:2,1:2)
         my_range(1:2,3)   = range(1:2,3)

!        Add an extra level for psmile_mg_final_gauss2

         levbeg = 0

      else

!        Generate multigrid sequence for really defined grids

         levbeg = 1

         if (gp%grid_type == PRISM_Irrlonlat_regvrt) then
!           Full range is required for definition of virtual cells
            my_range (1:2, 1:ndim_2d) = gp%grid_shape (1:2, 1:ndim_2d)
         endif

      endif
!
!     Get number of levels to be generated

      len = MAXVAL (my_range(2,1:ndim_3d)-my_range(1,1:ndim_3d)+1)
!
      ndlev = log (dble(len)) / log (2.0d0)
      if (2**ndlev .lt. len) ndlev = ndlev + 1
      ndlev = ndlev + 1

#ifdef DEBUG
         PRINT *, TRIM(ch_id), ': psmile_mg_setup: we have ndlev', ndlev, 'to be generated'
         call psmile_flushstd
#endif /* DEBUG */
!
      Allocate (levdim(levbeg:ndim_3d, 1:ndlev), STAT = ierror)
      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = ndim_3d * (ndlev + (levbeg-1))
         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'levdim', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif
!
      Allocate (gp%mg_infos(levbeg:ndlev), STAT = ierror)
      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = ndlev + (levbeg-1)
         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'gp%mg_infos', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif
!
!  Initialize Pointers in structure Grids(grid_id)%mg_infos
!
      do lev = levbeg, ndlev
         Nullify (gp%mg_infos(lev)%real_arrays)
         Nullify (gp%mg_infos(lev)%double_arrays)
#if defined ( PRISM_QUAD_TYPE )
         Nullify (gp%mg_infos(lev)%quad_arrays)
#endif
      end do

      gp%nlev = ndlev
      gp%ijk0 (:) = my_range(1,:)
!
!-----------------------------------------------------------------------
!     Allocate and generate grid functions of first level "1"
!     using the corner volumes.
!     The coarse levels "2:ndlev" are generated using the data of
!     the previously generated multigrid level.
!-----------------------------------------------------------------------
!
!     Create first level
!
      lev = 1
!
      levdim (1:ndim_3d, lev) = my_range(2,1:ndim_3d)-my_range(1,1:ndim_3d)
      gp%mg_infos(lev)%levdim = levdim (1:ndim_3d, lev)
!
      call psmile_mg_first_level (grid_id, my_range, &
                                  gp%mg_infos(lev),  &
                                  tol, simplified_grid, ierror)
      if (ierror /= 0) return
!
!     For all coarser levels

         do lev = 2, ndlev
         levf = lev - 1

!     ... Compute dimensions of coarser level "lev"
!
         icoarse (:) = 2
!
            do j = 1, ndim_3d
            n = levdim(j,levf) + 1
!
            if (n .eq. 1) then
               icoarse (j) = 1
               levdim (j,lev) = levdim (j,levf)
            else if ((n/2) * 2 == n) then
               levdim (j,lev) = n/2 - 1
            else
!              add a dummy line
               levdim (j,lev) = n/2
            endif
            end do
!
!     ... Create level "lev" from finer level "levf" = "lev-1"
!
            gp%mg_infos(lev)%levdim = levdim (1:ndim_3d, lev)

            call psmile_mg_coars_level (grid_id,           &
                                        gp%mg_infos(levf), &
                                        gp%mg_infos(lev),  &
                                        icoarse, ierror)
            if (ierror /= 0) return
!
!     if (levdim (*, lev) = 0 break

         end do ! for all levels

#ifdef PRISM_ASSERTION
      if (levdim (1, ndlev) /= 0 .or. &
          levdim (2, ndlev) /= 0 .or. &
          levdim (3, ndlev) /= 0) then
         print *, 'ndlev', ndlev
         print *, 'levdim(:, 1)', levdim (:, 1)
         print *, 'levdim(:, n)', levdim (:, ndlev)
         call psmile_assert (__FILE__, __LINE__, &
                             "incorrect coarsest levdim")
      endif
#endif

      Deallocate (levdim)
!
!=======================================================================
!     Allocate MG info for final grid in Gaussreduced storage.
!     This final level is stored as level 0 and the data is
!     computed using the really defined corner volumes.
!
!     Note: Since the face info "gp%face" has to be used
!           in routine PSMILe_mg_final_gauss2_type,
!           the data is generated for the entire grid !
!
!     Note: my_range(:,3) is set to gp%grid_shape(1,3)
!           since this level is only used for lonlat part of
!           GaussReduced grid
!=======================================================================
!
      if ( gp%grid_type == PRISM_Gaussreduced_regvrt ) then
         lev = 0
!
         my_range      = gp%grid_shape
         my_range(2,3) = gp%grid_shape(1,3)

         gp%mg_infos(lev)%levdim = my_range(2,1:ndim_3d)-my_range(1,1:ndim_3d)
         simplified_grid = .false.

         call psmile_mg_first_level (grid_id, my_range,   &
                                     gp%mg_infos(lev), &
                                     tol, simplified_grid, ierror)
         if (ierror /= 0) return
      endif
!
!=======================================================================
!
!=======================================================================
!
      call psmile_mg_get_cyclic (grid_id, my_range, tol, ierror)
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id), ierror
      call psmile_flushstd
#endif /* VERBOSE */
!
#ifdef VERBOSE

9990 format (1x, a, ': psmile_mg_setup: comp_id', i3, ', range ', 6i6)
9980 format (1x, a, ': psmile_mg_setup: eof ierror =', i3)

#endif /* VERBOSE */

      end subroutine PSMILe_MG_setup