psmile_gridless_func_dble.F90

Go to the documentation of this file.

!-----------------------------------------------------------------------
! $COPYRIGHT_CERFACS
! $COPYRIGHT_SGI
! $COPYRIGHT_CCRL
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: psmile_gridless_func_dble
!
! !INTERFACE:

   subroutine psmile_gridless_func_dble ( field_id, userdef_id, il_side, data_array, ierror )
!
! !USES:
!
      use PSMILe

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

!  handle in Userdefs global structures   field_id of the geographical function

      Integer, Intent (In)  :: userdef_id

!  handle in Userdefs global structures   userdef_id associated to field_id

      Integer, Intent (In)  :: il_side

!  flag = 0 for the "source" side,  = 1 for the "target" side

      Double Precision, Intent (InOut)  :: data_array(*)

!  The data (geographical function) itself

!  OUTPUT PARAMETER:

      Integer, Intent (Out) :: ierror

!     error return code

!  LOCAL VARIABLES:

      Integer                   :: i
      Integer                   :: imaxg, il_celldim
      Integer                   :: iloc, ilink
      Integer                   :: indi, indj, indk
      Integer                   :: dim1, dim2, dim3
      Integer                   :: ijklin, ijkdat, ndim
      Integer                   :: il_nbdl, ibd


      Type(Userdef), POINTER    :: ug
      Type(GridFunction), POINTER    :: fp


!   user-defined interpolation : set pointer to Userdefs structures

      ierror = 0
      ug => Userdefs(userdef_id)
      imaxg      = ug%ig_nb_ppp
      il_celldim = ug%ig_celldim
      print *, ' imaxg = ',imaxg

!   Fill relevant values from the user-defined links
!   data_array must be used as a one-dimension array.
    
      fp => Fields(field_id)
      ndim = Grids(Methods(fp%method_id)%grid_id)%n_dim   ! dimension of geographical function
      print *, ' ndim = ',ndim
!
!  Verifier (prism assertion) que il_celldim = ndim
!
         dim1 = fp%var_shape(2,1) - fp%var_shape(1,1) + 1
      print *, ' dim1 = ', dim1
         dim2 = fp%var_shape(2,2) - fp%var_shape(1,2) + 1
      print *, ' dim2 = ', dim2
      if ( ndim == 3 ) then
         dim3 = fp%var_shape(2,3) - fp%var_shape(1,3) + 1
      print *, ' dim3 = ', dim3
      endif

!  Number of fields ( = 1, if no bundle,  > 1 if bundles )
      il_nbdl = ug%ig_nbr_fields
!
!  Side-dependent code : 0 (source) data_array is input
!                        1 (target) data_array is output

      if ( il_side == 0 ) then
!   Source side : ug%dble_gridless is output (gridless function)
!!!        if ( ug%lg_nolink ) then
!     Case of a fake link (nothing to send in fact)
!!!           ug%dble_gridless(1,1,1,1) = 0.0 
!!!        else
!     Normal case : user_defined links with non-zero weigh
          if ( ndim == 3 ) then
!   Specific to Three-dimensional grid
            do iloc = 1, imaxg
               indi  = ug%iga_igl(iloc,1) - fp%var_shape(1,1) + 1
               indj  = ug%iga_igl(iloc,2) - fp%var_shape(1,2) + 1 
               indk  = ug%iga_igl(iloc,3) - fp%var_shape(1,3) + 1
               ilink = ug%iga_igl(iloc,4)
               ijklin = indi + dim1*(indj-1) + dim1*dim2*(indk-1)
#ifdef DEBUG
               print *, 'indi-j-k  ilink = ',indi, indj, indk, ilink
               print *, 'ijkdat  et data = ',ijklin, data_array(ijklin)
#endif
               do ibd = 1, il_nbdl
                  ijkdat = ijklin + dim1*dim2*dim3*(ibd-1)
                  ug%dble_gridless(iloc,1,1,ibd) = &
                      (ug%dga_wght(iloc)) * data_array(ijkdat)
               enddo
            enddo
          endif
!
          if ( ndim == 2 ) then
!   Specific to Two-dimensional grid
            do iloc = 1, imaxg
               indi  = ug%iga_igl(iloc,1) - fp%var_shape(1,1) + 1
               indj  = ug%iga_igl(iloc,2) - fp%var_shape(1,2) + 1
               ilink = ug%iga_igl(iloc,4)
               ijklin = indi + dim1*(indj-1)
#ifdef DEBUG
               print *, 'indi-j  ilink = ',indi, indj, ilink
               print *, 'ijkdat  et data = ',ijklin, data_array(ijklin)
#endif
               do ibd = 1, il_nbdl
                  ijkdat = ijklin + dim1*dim2*(ibd-1)
                  ug%dble_gridless(iloc,1,1,ibd) = &
                     (ug%dga_wght(iloc)) * data_array(ijkdat)
               enddo
            enddo
          endif
!!!        endif        ! lg_nolink

      elseif (il_side == 1 ) then
!   Target side : ug%dble_gridless is input, and data_array is output
!!!        if ( .NOT.ug%lg_nolink ) then
          if ( ndim == 3 ) then
            do iloc = 1, imaxg
               indi  = ug%iga_igl(iloc,1) - fp%var_shape(1,1) + 1
               indj  = ug%iga_igl(iloc,2) - fp%var_shape(1,2) + 1
               indk  = ug%iga_igl(iloc,3) - fp%var_shape(1,3) + 1
               ilink = ug%iga_igl(iloc,4)
               ijklin = indi + dim1*(indj-1) + dim1*dim2*(indk-1)
#ifdef DEBUG
               print *, 'indi-j-k  ilink = ',indi, indj, indk, ilink
#endif
               do ibd = 1, il_nbdl
                  ijkdat = ijklin + dim1*dim2*dim3*(ibd-1)
                  data_array(ijkdat) = data_array(ijkdat) +  &
                                       ug%dble_gridless(iloc,1,1,ibd)
#ifdef DEBUG
                  print *, 'ijkdat  et data = ',ijkdat, data_array(ijkdat)
#endif
               enddo
            enddo
          endif

          if ( ndim == 2 ) then
            do iloc = 1, imaxg
               indi  = ug%iga_igl(iloc,1) - fp%var_shape(1,1) + 1
               indj  = ug%iga_igl(iloc,2) - fp%var_shape(1,2) + 1
               ilink = ug%iga_igl(iloc,4)
#ifdef DEBUG
               print *, 'indi-j-k  ilink = ',indi, indj, ilink
#endif
               ijklin = indi + dim1*(indj-1)
               do ibd = 1, il_nbdl
                  ijkdat = ijklin + dim1*dim2*(ibd-1)
                  data_array(ijkdat) = data_array(ijkdat) +  &
                                       ug%dble_gridless(iloc,1,1,ibd)
#ifdef DEBUG
                  print *, 'ijkdat  et data = ',ijkdat, data_array(ijkdat)
#endif
               enddo
            enddo
          endif
!!!        endif   ! lg_nolink

      endif     ! il_side

#ifdef DEBUG
      PRINT *, ' local values of function on gridless grid (j = k = 1) '
      PRINT *, (ug%dble_gridless(i,1,1,1),i=1,imaxg)
#endif

      return
      end subroutine psmile_gridless_func_dble