psmile_mg_first_level_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_first_level_dble
!
! !INTERFACE:

      subroutine psmile_mg_first_level_dble (grid_id, range, &
                      mg_info, tol, simplified_grid, ierror)
!
! !USES:
!
      use PRISM_constants
      use psmile_grid, only : common_grid_range
      use PSMILe, dummy_interface => PSMILe_MG_first_level_dble

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

!     Specifies the handle to the grid information.

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

!     Range for which MG sequence should be setup

      Real (PSMILe_float_kind), Intent (In) :: tol

!     Search tolerance
!
      Logical,                  Intent (In) :: simplified_grid
!
!     Should the simplified grids for grids of type
!     "PRISM_Gaussreduced_regvrt" be used ?
!
! !INPUT/OUTPUT PARAMETERS:
!
      Type (Enddef_mg),      Intent (InOut) :: mg_info

!     Structure returning the data on the first multigrid level
!     which is created by this routine.
!
! !OUTPUT PARAMETERS:
!
      Integer,               Intent (Out)   :: ierror

!     Returns the error code of PSMILe_MG_first_level_dble;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !DEFINED PARAMETERS:
!
      Integer, Parameter                  :: nc_reg = 2

!     Number of corners for regular directions

      Integer, Parameter                  :: nd_cyclic = 2

!     Number of possible cylic coordinates (lon lat)
!
! !LOCAL VARIABLES
!
      Type (Grid), Pointer                :: gp
      Type (Corner_Block), Pointer        :: corner_pointer
      Type (Enddef_mg_double), Pointer    :: arrays
!
      Integer                             :: grid_type
!
!  ... for dimensions
!
      integer                             :: index_c ! index in corner shape
      integer                             :: index_g ! index in grid   shape
      integer                             :: i, ii, jj, kk
      integer                             :: leni, lenij
      integer                             :: dim_i, dim_ij
      integer                             :: dim (ndim_3d)
!
!  ... for lonlat periods
!
      Double Precision                    :: len_cyclic (nd_cyclic)
      Double Precision                    :: bounds_lon (2)
!
!  ... for error parameters
!
      Integer, Parameter                  :: nerrp = 2
      Integer                             :: ierrp (nerrp)
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_MG_first_level_dble" creates the multi-grid
! information for the first level of the multigrid sequence.
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 03.06.05    H. Ritzdorf  created
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_mg_first_level_dble.F90 2687 2010-10-28 15:15:52Z coquart $
!  $Author: coquart $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_mg_first_level_dble.F90 2687 2010-10-28 15:15:52Z coquart $'
!
!----------------------------------------------------------------------

#ifdef VERBOSE
      print *, trim(ch_id), ': psmile_mg_first_level_dble: grid_id', grid_id

      call psmile_flushstd
#endif /* VERBOSE */
!
!  Initialization
!
      ierror = 0
      gp => Grids(grid_id)
!
      if ( gp%grid_type == PRISM_Gaussreduced_regvrt .and. &
           simplified_grid) then

!        use the simplified corner volumes

         corner_pointer => gp%gcorner_pointer
         grid_type = PRISM_Reglonlatvrt

      else

!        use the real corner volumes

         corner_pointer => gp%corner_pointer
         grid_type = gp%grid_type

      endif
!
!     lon lat range for cyclic indices is 360 or 180 degrees, respectively.
!
      len_cyclic (1:nd_cyclic) = common_grid_range(2,1:nd_cyclic) - &
                                 common_grid_range(1,1:nd_cyclic)
!
#ifdef PRISM_ASSERTION

     if (corner_pointer%corner_datatype /= MPI_DOUBLE_PRECISION) then
        call psmile_assert ( __FILE__, __LINE__, &
                            'Corner type is not MPI_DOUBLE_PRECISION')

     endif

!    Internal control: Is the data really available

      if (.not. Associated(corner_pointer%corners_dble(1)%vector) ) then
         call psmile_assert ( __FILE__, __LINE__, &
                'Pointer corners_dble(1)%vector is not set')
      endif

      if (.not. Associated(corner_pointer%corners_dble(2)%vector) ) then
         call psmile_assert ( __FILE__, __LINE__, &
                'Pointer corners_dble(2)%vector is not set')
      endif

      if (.not. Associated(corner_pointer%corners_dble(3)%vector) ) then
         call psmile_assert ( __FILE__, __LINE__, &
                'Pointer corners_dble(3)%vector is not set')
      endif
#endif /* PRISM_ASSERTION */
!
!   Allocate pointer to double arrays
!
      Allocate (mg_info%double_arrays, STAT = ierror)
      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = 1
         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'mg_info%double_arrays', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif
!
      arrays => mg_info%double_arrays
!
! ====================================================================
!     Get dimensions of arrays to be allocated and
!     allocate the arrays.
!
!     "levdim" specifies the dimension of the grid level to be created.
!     The dimension is (0:levdim(i)) for each direction.
! ====================================================================

      select case ( grid_type )

! -----------------------------------------------------------------------
!         Regular in all directions
!         ===> Number of corners in all directions = 2
!              (cf. prism_set_corners_3d_double.F90)
! -----------------------------------------------------------------------
!
      case (PRISM_Reglonlatvrt)
!
         dim(1) = mg_info%levdim(1) + 1
         dim(2) = mg_info%levdim(2) + 1
         dim(3) = mg_info%levdim(3) + 1
!
! -----------------------------------------------------------------------
!       Irregular in lonlat   direction
!         Regular in vertical direction
!         ===> Number of corners in vertical direction = 2
!              (cf. prism_set_corners_3d_double.F90)
! -----------------------------------------------------------------------
!
      case (PRISM_Irrlonlat_regvrt)

         dim(1) = (mg_info%levdim(1)+1) * (mg_info%levdim(2)+1)
         dim(2) = dim(1)
         dim(3) = mg_info%levdim (3) + 1
!
! -----------------------------------------------------------------------
!       Gaussreduced in lonlat   direction
!         Regular in vertical direction
!         ===> Number of corners in vertical direction = 2
!              (cf. prism_set_corners_3d_double.F90)
! -----------------------------------------------------------------------
!
      case (PRISM_Gaussreduced_regvrt)

         dim(1) = mg_info%levdim(1) + 1
         dim(2) = dim (1)
         dim(3) = mg_info%levdim(3) + 1
!
! -----------------------------------------------------------------------
!       Irregular in lonlat   and vertical direction
! -----------------------------------------------------------------------
!
      case (PRISM_Irrlonlatvrt)

         dim(1) = (mg_info%levdim(1)+1) * (mg_info%levdim(2)+1) * &
                  (mg_info%levdim(3)+1)
         dim(2) = dim(1)
         dim(3) = dim(2)
!
! -----------------------------------------------------------------------
!     Error: unsupported grid type
! -----------------------------------------------------------------------
!
      case DEFAULT
!
          ierrp (1) = grid_type

          ierror = PRISM_Error_Internal

          call psmile_error ( ierror, 'unsupported grid generation type', &
                              ierrp, 1, __FILE__, __LINE__ )

      end select
!
#ifdef DEBUG
      print *, trim(ch_id), ': psmile_mg_first_level_dble: dim(1:3)', dim(1:3)
      call psmile_flushstd
#endif
!  ... Allocate arrays
!
         do i = 1, ndim_3d

         Allocate (arrays%chmin(i)%vector(dim(i)), &
                   arrays%chmax(i)%vector(dim(i)), &
                   arrays%midp (i)%vector(dim(i)), STAT = ierror)
         if ( ierror > 0 ) then
            ierrp (1) = ierror
            ierrp (2) = dim (i) * 3
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'mg_info%{chmin,chmax,midp}(i)%vector', &
                                ierrp, 2, __FILE__, __LINE__ )
            return
         endif

         end do
!
! ====================================================================
!  Now, the multi-grid informations has to be generated depending
!  on the grid type "grid_type". 
!
!  Waere es nicht sinnvoller, immer die 3d Routine zu verwenden
!  und die anderen Dimensionen auf 1 zu setzen.
! ====================================================================

      select case ( grid_type )

! -----------------------------------------------------------------------
!         Regular in all directions
!         ===> Number of corners in all directions = 2
!              (cf. prism_set_corners_3d_double.F90)
!
!         Note: This case is also used for
!               (*) grids of type "PRISM_Gaussreduced_regvrt"
!                   simplified corner volumes.
!                
! -----------------------------------------------------------------------
!
      case (PRISM_Reglonlatvrt)
!
            do i = 1, ndim_3d
            call psmile_mg_first_subgrid_1d_dble (   &
                corner_pointer%corners_dble(i)%vector,     &
                corner_pointer%corner_shape(1,i), &
                corner_pointer%corner_shape(2,i), &
                nc_reg,                           &
                range (1, i),                     &
                arrays%chmin(i)%vector, arrays%chmax(i)%vector, &
                arrays%midp(i)%vector, &
                mg_info%levdim (i), ierror)
            if (ierror > 0) return
            end do
!
!           ... Control multigrid requirements for lonlat directions
!
            do i = 1, nd_cyclic
            call psmile_mg_ctrl_subgrid_1d_dble (          &
                corner_pointer%corners_dble(i)%vector,     &
                corner_pointer%corner_shape(:,i), nc_reg,  &
                range (:, i),                              &
                arrays%chmin(i)%vector, arrays%chmax(i)%vector, &
                mg_info%levdim (i),                        &
                len_cyclic(i), grid_id, i, ierror)
            if (ierror > 0) return
            end do
!
! -----------------------------------------------------------------------
!       Irregular in lonlat   direction
!         Regular in vertical direction
!         ===> Number of corners in vertical direction = 2
!              (cf. prism_set_corners_3d_double.F90)
! -----------------------------------------------------------------------
!
      case (PRISM_Irrlonlat_regvrt)
!
         do i = 1, ndim_2d
         call psmile_mg_first_subgrid_2d_dble (   &
             corner_pointer%corners_dble(i)%vector,       &
             corner_pointer%corner_shape(1,1),   &
             corner_pointer%corner_shape(2,1),   &
             corner_pointer%corner_shape(1,2),   &
             corner_pointer%corner_shape(2,2),   &
             corner_pointer%nbr_corners/nc_reg,  &
             range (1, 1),                     &
             arrays%chmin(i)%vector, arrays%chmax(i)%vector, &
             arrays%midp(i)%vector, &
             mg_info%levdim (1), mg_info%levdim (2), ierror)
         if (ierror > 0) return
         end do
!
         call psmile_mg_first_subgrid_1d_dble (   &
             corner_pointer%corners_dble(3)%vector,     &
             corner_pointer%corner_shape(1,3), &
             corner_pointer%corner_shape(2,3), &
             nc_reg,                           &
             range (1, 3),                     &
             arrays%chmin(3)%vector, arrays%chmax(3)%vector, arrays%midp(3)%vector, &
             mg_info%levdim (3), ierror)
         if (ierror > 0) return
!
!           ... Control multigrid requirements for lonlat directions
!
            do i = 1, nd_cyclic
            call psmile_mg_ctrl_subgrid_2d_dble (          &
                corner_pointer%corners_dble(i)%vector,     &
                corner_pointer%corner_shape(:,1:2),        &
                corner_pointer%nbr_corners/nc_reg,         &
                range (:, 1:2),                            &
                arrays%chmin(i)%vector, arrays%chmax(i)%vector, &
                mg_info%levdim (1:2), &
                len_cyclic(i), grid_id, i, ierror)
            if (ierror > 0) return
            end do
!
! -----------------------------------------------------------------------
!    GaussReduced in lonlat   direction
!         Regular in vertical direction
!         ===> Number of corners in all directions = 2
!              (cf. prism_set_corners_3d_double.F90)
!
!    The lonlat direction is handled as 1--dimensional vector with
!    the same dimension in all lon and lat direction.
!                
! -----------------------------------------------------------------------
!
      case (PRISM_Gaussreduced_regvrt)
!
            do i = 1, ndim_2d
            call psmile_mg_first_subgrid_1d_dble (     &
                corner_pointer%corners_dble(i)%vector, &
                corner_pointer%corner_shape(1,1),      &
                corner_pointer%corner_shape(2,1),      &
                nc_reg,                                &
                range (1, 1),                          &
                arrays%chmin(i)%vector, arrays%chmax(i)%vector, &
                arrays%midp(i)%vector,                 &
                mg_info%levdim (1), ierror)
            if (ierror > 0) return
            end do
!
         call psmile_mg_first_subgrid_1d_dble (       &
               corner_pointer%corners_dble(3)%vector, &
               corner_pointer%corner_shape(1,3),      &
               corner_pointer%corner_shape(2,3),      &
               nc_reg, range (1, 3),                  &
               arrays%chmin(3)%vector, arrays%chmax(3)%vector, &
               arrays%midp(3)%vector,                 &
               mg_info%levdim (3), ierror)
          if (ierror > 0) return
!
!           ... Control multigrid requirements for lonlat directions
!
            do i = 1, nd_cyclic
            call psmile_mg_ctrl_subgrid_1d_dble (          &
                corner_pointer%corners_dble(i)%vector,     &
                corner_pointer%corner_shape(:,1), nc_reg,  &
                range (:, 1),                              &
                arrays%chmin(i)%vector, arrays%chmax(i)%vector, &
                mg_info%levdim (i),                        &
                len_cyclic(i), grid_id, i, ierror)
            if (ierror > 0) return
            end do
!
! -----------------------------------------------------------------------
!       Irregular in lonlat   and vertical direction
! -----------------------------------------------------------------------
!
      case (PRISM_Irrlonlatvrt)

            do i = 1, ndim_3d
            call psmile_mg_first_subgrid_3d_dble ( &
                corner_pointer%corners_dble(i)%vector,       &
                corner_pointer%corner_shape(1,1),   &
                corner_pointer%corner_shape(2,1),   &
                corner_pointer%corner_shape(1,2),   &
                corner_pointer%corner_shape(2,2),   &
                corner_pointer%corner_shape(1,3),   &
                corner_pointer%corner_shape(2,3),   &
                corner_pointer%nbr_corners,         &
                range,                              &
                arrays%chmin(i)%vector, arrays%chmax(i)%vector, &
                arrays%midp(i)%vector, &
                mg_info%levdim (1), mg_info%levdim (2), mg_info%levdim (3), ierror)
            if (ierror > 0) return
            end do
!
!           ... Control multigrid requirements for lonlat directions
!
            do i = 1, nd_cyclic
            call psmile_mg_ctrl_subgrid_3d_dble (      &
                corner_pointer%corners_dble(i)%vector, &
                corner_pointer%corner_shape,           &
                corner_pointer%nbr_corners,            &
                range,                                 &
                arrays%chmin(i)%vector,                &
                arrays%chmax(i)%vector,                &
                mg_info%levdim, &
                len_cyclic(i), grid_id, i, ierror)
            if (ierror > 0) return
            end do
!
! -----------------------------------------------------------------------
!     Error: unsupported grid type
! -----------------------------------------------------------------------
!
      case DEFAULT
!
          ierrp (1) = grid_type

          ierror = PRISM_Error_Internal

          call psmile_error ( ierror, 'unsupported grid generation type', &
                              ierrp, 1, __FILE__, __LINE__ )

      end select

! -----------------------------------------------------------------------
!        ... extent cell over the pole onto the pole point
!     HUHU: Why no pole cells for PRISM_Reglonlatvrt ?
! -----------------------------------------------------------------------

      if ( gp%grid_type == PRISM_Irrlonlat_regvrt .or. &
           gp%grid_type == PRISM_Irrlonlatvrt) then

         bounds_lon(1) = -180.0
         bounds_lon(2) =  180.0

         if ( associated ( corner_pointer%pole_array ) ) then

            leni  = corner_pointer%corner_shape(2,1) - &
                    corner_pointer%corner_shape(1,1) + 1
            lenij = leni * & 
                   (corner_pointer%corner_shape(2,2) - &
                    corner_pointer%corner_shape(1,2) + 1)

            dim_i  = (mg_info%levdim(1)+1)
            dim_ij = dim_i * (mg_info%levdim(2)+1)

            do i = 1, size(corner_pointer%pole_array)

               index_c = corner_pointer%pole_array(i)

!              Transfer 1-dimensional index into 3-dimensional index 
!              within the corner shape

               kk =  (index_c-1) / lenij
               jj = ((index_c-1) -kk*lenij) / leni
               ii =  (index_c-1) -kk*lenij -jj*leni
!
               ii = corner_pointer%corner_shape(1,1) + ii
               jj = corner_pointer%corner_shape(1,2) + jj
               kk = corner_pointer%corner_shape(1,3) + kk

!              Transfer 3-dimensional index into 1-dimensional index 
!              in chmin, chmax and midp

               index_g = (kk-gp%ijk0(3)) * dim_ij &
                       + (jj-gp%ijk0(2)) * dim_i  &
                       +  ii-gp%ijk0(1)  + 1
!
!              Extent chmin and chmax only when the pole is included in the subsection
!
               if ( index_g < dim(1) .and. index_g < dim(2) ) then
#ifdef DEBUG
                  print *, trim(ch_id), ': psmile_mg_first_level_dble: chmax extended to pole'
                  print *, trim(ch_id), '   for index in valid_shape:'
                  print *, trim(ch_id), '   i ', ii, ' j ', jj, ' k ', kk
                  call psmile_flushstd
#endif /* DEBUG */

                  if ( corner_pointer%corners_dble(2)%vector(index_c) > 0.0 ) then
                     ! HUHU: TODO: Improve setting of lon extension; edge on pole; more than 1 cell
                     !       There seems to be a problem in nearest neighbour search
                     !       setting to -360:360 changes the result of nearest neighbour search
                     !       in example cicle_1 using 1 proc
                     !
                     arrays%chmin(1)%vector(index_g) =  min (bounds_lon(1), &
                          arrays%chmin(1)%vector(index_g))
                     arrays%chmax(1)%vector(index_g) =  max (bounds_lon(2), &
                          arrays%chmax(1)%vector(index_g))
                     !                 arrays%midp (1)%vector(index_g) =     0.0

                     arrays%chmax(2)%vector(index_g) =  90.0
                     arrays%midp (2)%vector(index_g) =  90.0

                  else if ( corner_pointer%corners_dble(2)%vector(index_c) < 0.0 ) then
                     arrays%chmin(1)%vector(index_g) =  min (bounds_lon(1), &
                          arrays%chmin(1)%vector(index_g))
                     arrays%chmax(1)%vector(index_g) =  max (bounds_lon(2), &
                          arrays%chmax(1)%vector(index_g))
                     !                 arrays%midp (1)%vector(index_g) =     0.0

                     arrays%chmin(2)%vector(index_g) = -90.0
                     arrays%midp(2)%vector(index_g)  = -90.0
                  endif

               endif

            enddo

         endif

      endif
!
!===> All done
!
#ifdef VERBOSE
      print *, trim(ch_id), ': psmile_mg_first_level_dble eof: grid_id',&
                             grid_id, ', ierror =', ierror

      call psmile_flushstd
#endif /* VERBOSE */
!
      end subroutine PSMILe_MG_first_level_dble