psmile_set_userdef.F90

Go to the documentation of this file.

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

      subroutine psmile_set_userdef (var_id, il_side, chan_id, ierror)

!
! !USES:
!
      Use PRISM_Constants
      Use PSMILe, dummy_interface => psmile_set_userdef
      Use PSMILe_smioc
      Use prism
      use psmile_user_data, only : psmile_store_data_intern_field, &
                                   psmile_store_data_intern_points
!
      Implicit none
!
      INCLUDE 'netcdf.inc'
!
! !INPUT PARAMETERS:
!
      Integer, intent (In)                :: var_id

!     Index to the global Field array (identify the transient variable)

      Integer, intent (In)                :: il_side

!     = 0 on source il_side,    = 1 on target il_side of the transient
!
      Integer, intent (In)                :: chan_id

!     last index of input or output channel using a User-Defined Interpolation
!
! !OUTPUT PARAMETERS:
!
      Integer, Intent (Out)               :: ierror

!     Returns the error code of prism_enddef;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !LOCAL VARIABLES
!
      Character(len=max_name)      :: cl_file_name
      Character(len=16), parameter :: point_name_2 = 'grid-point'

      Integer                  :: i, j, k, il_fileid, il_varid, il_dimid
      Integer                  :: nlinks,  ilink, il_celldim
      Integer                  :: il_grid_id
      Integer                  :: indi_min, indi_max, indj_min, indj_max
      Integer                  :: indk_min, indk_max
      Integer                  :: il_val11, il_val12, il_val13
      Integer                  :: il_nb_ppp, iloc, il_blocks
      Integer                  :: il_nbr_blocks, il_ndim
      Integer                  :: isbd, il_deb
      Integer                  :: iming, imaxg, jming, jmaxg, kming, kmaxg
      Integer                  :: indi, indj, indk
      Integer                  :: gr_id, meth_id, mk_id
      Integer                  :: lenbuf
!
!     Variables for gridless grid, mask and function
!
      Integer                  :: method_id_2
      Integer                  :: grid_id_2
      Integer                  :: mask_id_2
      Integer                  :: grid_type_2
      Integer                  :: id_comp_id
      Integer                  :: ass_var_id_2
      Integer                  :: userdef_id_2
      Integer                  :: valid_shape_pr(2,4)
      Integer                  :: actual_shape_pr(2,4)
      Integer, allocatable     :: il_offspr(:,:)
      Integer, allocatable     :: il_extepr(:,:)
      Integer                  :: il_var_nodims(2)
      Integer                  :: id_var_type
      Integer                  :: nbr_fields

      Character(len=max_name)  :: var_name_gl
      Character(len=max_name)  :: grid_name_gl
      Character(len=max_name)  :: cl_buffer

      Logical                  :: ll_novar
      Logical, parameter       :: new_points = .TRUE.
      Logical, parameter       :: new_mask = .TRUE.
      Logical, allocatable     :: gprmask (:,:,:)
!
!     local array for the links when reading the weight and address file
      Type(PSMILe_Link), pointer :: ila_links(:)
      Integer                    :: ierrp(2)
!
!     Pointers on global variables
!
      Type(Grid),POINTER            :: gp
      Type(Userdef),POINTER         :: up
      Type(GridFunction),POINTER    :: fp



#ifdef VERBOSE
      print *, trim(ch_id), ': psmile_set_userdef : start'
      call PSMILe_Flushstd
#endif /* VERBOSE */
!
!   Initialization
!
      ierror = 0
!
!  Get a handle in the Userdefs array

      call psmile_get_userdef_handle (userdef_id_2, ierror)
      up => Userdefs(userdef_id_2)
!
      up%var_id     = var_id               ! geographical field id
      up%ig_transi_side = il_side             !  source   .OR.   target
!
      fp => Fields(var_id)                 ! geographical field structure
!  
!  Get the name of the file containing the weights and addresses
!
      if ( il_side == 0 ) then
         cl_file_name = &
         trim(fp%Taskout(chan_id)%interp%arg10%cg_file_name)
      elseif ( il_side == 1 ) then
         cl_file_name = &
         trim(fp%Taskin%In_channel(chan_id)%interp%arg10%cg_file_name)
      endif
#ifdef DEBUG
      PRINT *, trim(ch_id), ': psmile_set_userdef : cl_file_name = ',trim(cl_file_name)
#endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!! 1. Open and read  user-defined weight-and-address file
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
      CALL psmile_ud_hdlerr(NF_OPEN(TRIM(cl_file_name), NF_NOWRITE, il_fileid))
!
#ifdef DEBUG
      PRINT *, trim(ch_id), ': psmile_set_userdef : il_fileid = ',il_fileid
#endif
!  Get dimension of array src_ind1 i.e. the number of links
      CALL psmile_ud_hdlerr(NF_INQ_VARID   (il_fileid, 'src_ind1', il_varid))
      CALL psmile_ud_hdlerr(NF_INQ_VARDIMID(il_fileid, il_varid, il_dimid))
      CALL psmile_ud_hdlerr(NF_INQ_DIMLEN  (il_fileid, il_dimid, nlinks))
      up%ig_nb_links = nlinks
#ifdef DEBUG
      PRINT *, trim(ch_id), ': psmile_set_userdef : nlinks =', nlinks
#endif
!
!  Allocate array for links ( same number on any side )
      ALLOCATE (ila_links(1:nlinks), STAT=ierror )
      IF (ierror /= 0) THEN
         PRINT *, 'Error in allocating weight-and-address arrays = ', ierror
         CALL PSMILe_Flushstd
         CALL PSMILe_Abort
      ENDIF
      do k = 1, nlinks
         ila_links(k)%cell_id(1) = 0
         ila_links(k)%cell_id(2) = 0
         ila_links(k)%cell_id(3) = 0
         ila_links(k)%weight = 0.0D0
      enddo
!
!  Get dimension of cell_id array defined in the file by the user(depends on side)
!
      il_celldim = 3
      if ( il_side == 0 ) then
         CALL psmile_ud_hdlerr2( NF_INQ_VARID(il_fileid, 'src_ind3', il_varid), ll_novar)
         if ( ll_novar ) il_celldim = 2
         CALL psmile_ud_hdlerr2( NF_INQ_VARID(il_fileid, 'src_ind2', il_varid), ll_novar)
         if ( ll_novar ) il_celldim = 1
      elseif (il_side == 1 ) then
         CALL psmile_ud_hdlerr2( NF_INQ_VARID(il_fileid, 'tgt_ind3', il_varid), ll_novar)
         if ( ll_novar ) il_celldim = 2
         CALL psmile_ud_hdlerr2( NF_INQ_VARID(il_fileid, 'tgt_ind2', il_varid), ll_novar)
         if ( ll_novar ) il_celldim = 1
      endif
!     keep geographical cell dimension in Userdef structure
      up%ig_celldim = il_celldim

#ifdef DEBUG
      PRINT *, trim(ch_id), ': psmile_set_userdef : cell_dim =', il_celldim
#endif
!
!  Read links depending on the side (source links .OR. target links)
!
      if ( il_side == 0 ) then
!        Get source links in array ila_links on the source side
         CALL psmile_ud_hdlerr(NF_INQ_VARID   (il_fileid, 'src_ind1', il_varid))
         CALL psmile_ud_hdlerr(NF_GET_VAR_INT (il_fileid, il_varid,     &
                     ila_links(:)%cell_id(1)))
         if ( il_celldim >= 2 ) then
            CALL psmile_ud_hdlerr(NF_INQ_VARID(il_fileid, 'src_ind2', il_varid))
            CALL psmile_ud_hdlerr(NF_GET_VAR_INT (il_fileid, il_varid,     &
                        ila_links(:)%cell_id(2)))
         endif
         if ( il_celldim == 3 ) then         
            CALL psmile_ud_hdlerr(NF_INQ_VARID(il_fileid, 'src_ind3', il_varid))
            CALL psmile_ud_hdlerr(NF_GET_VAR_INT (il_fileid, il_varid,     &
                        ila_links(:)%cell_id(3)))
         endif
         CALL psmile_ud_hdlerr(NF_INQ_VARID(il_fileid, 'weight', il_varid))
         CALL psmile_ud_hdlerr(NF_GET_VAR_DOUBLE (il_fileid, il_varid,  &
                     ila_links(:)%weight))
!
      elseif ( il_side == 1 ) then
!        Get target links in array ila_links on the target side
         CALL psmile_ud_hdlerr(NF_INQ_VARID(il_fileid, 'tgt_ind1', il_varid))
         CALL psmile_ud_hdlerr(NF_GET_VAR_INT (il_fileid, il_varid,     &
                     ila_links(:)%cell_id(1)))
         if ( il_celldim >= 2 ) then
            CALL psmile_ud_hdlerr(NF_INQ_VARID(il_fileid, 'tgt_ind2', il_varid))
            CALL psmile_ud_hdlerr(NF_GET_VAR_INT (il_fileid, il_varid,     &
                        ila_links(:)%cell_id(2)))
         endif
         if ( il_celldim == 3 ) then
            CALL psmile_ud_hdlerr(NF_INQ_VARID(il_fileid, 'tgt_ind3', il_varid))
            CALL psmile_ud_hdlerr(NF_GET_VAR_INT (il_fileid, il_varid,     &
                        ila_links(:)%cell_id(3)))
         endif
      endif
!
!  Close Netcdf file
         CALL psmile_ud_hdlerr(NF_CLOSE(il_fileid))
#ifdef DEBUGX
    PRINT *,'links(:)%cell_id(1)= ',ila_links(:)%cell_id(1)
    PRINT *,'links(:)%cell_id(2)= ',ila_links(:)%cell_id(2)
    PRINT *,'links(:)%cell_id(3)= ',ila_links(:)%cell_id(3)
    if (il_side == 0 ) then
       PRINT *,'links(:)%weight = ',ila_links(:)%weight
    endif
#endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!! 2. Define and declare intermediate gridless exchange grid
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  2.1  Get offsets and extents of geographical grid 
!
      il_grid_id = Methods(fp%method_id)%grid_id
      gp => Grids(il_grid_id)

      il_blocks = size(gp%partition(:,:),DIM=1)
      IF (il_blocks /= 1) THEN
          PRINT*, 'psmile_set_userdef does not support more than one block yet'
          CALL PSMILe_Abort
      ENDIF
!
!  2.2  store limits of geographical grid partition
!
      indi_min = gp%partition(1,1)+1
      indi_max = gp%partition(1,1)+gp%extent(1,1)
      if ( il_celldim >= 2 ) then
         indj_min = gp%partition(1,2)+1
         indj_max = gp%partition(1,2)+gp%extent(1,2)
      endif
      if ( il_celldim == 3 ) then
         indk_min = gp%partition(1,3)+1
         indk_max = gp%partition(1,3)+gp%extent(1,3)
      endif
  !
      il_val11 = fp%var_shape(1,1)
      if ( il_celldim >= 2 ) then
         il_val12 = fp%var_shape(1,2)
      endif
      if ( il_celldim == 3 ) then
         il_val13 = fp%var_shape(1,3)
      endif

#ifdef DEBUG
  PRINT *," i min max, j min max, k min max "
  PRINT *, indi_min,indi_max,indj_min,indj_max,indk_min,indk_max
  PRINT *, il_val11, il_val12, il_val13
#endif
!
! 2.3  Compute the number of points in my gridlessgrid partition :
!      = number of links for which the src index in geographic grid
!        is within the limits of the geographic partition for this PE
!      Tests on indexes depend on the number of dimensions of the geo-grid
!
      il_nb_ppp = 0           ! initial (and may be final !) value
!
!   Specific to One-dimensional grid
      if ( il_celldim == 1 ) then
         do i = 1, nlinks
            indi = ila_links(i)%cell_id(1)
            IF (indi .GE. indi_min .AND. indi .LE. indi_max ) THEN
               il_nb_ppp = il_nb_ppp + 1
            ENDIF
         enddo
!   Specific to Bi-dimensional grid
      elseif ( il_celldim == 2 ) then
         do i = 1, nlinks
            indi = ila_links(i)%cell_id(1)
            indj = ila_links(i)%cell_id(2)
            IF (indi .GE. indi_min .AND. indi .LE. indi_max &
               .AND. indj .GE. indj_min .AND. indj .LE. indj_max) THEN
               il_nb_ppp = il_nb_ppp + 1
            ENDIF
         enddo
!   Specific to Three-dimensional grid
      elseif ( il_celldim == 3 ) then
         do i = 1, nlinks
            indi = ila_links(i)%cell_id(1)
            indj = ila_links(i)%cell_id(2)
            indk = ila_links(i)%cell_id(3)
            IF (indi .GE. indi_min .AND. indi .LE. indi_max &
               .AND. indj .GE. indj_min .AND. indj .LE. indj_max &
               .AND. indk .GE. indk_min .AND. indk .LE. indk_max) THEN
               il_nb_ppp = il_nb_ppp + 1
            ENDIF
         enddo
      endif
      up%ig_nb_ppp = il_nb_ppp
#ifdef DEBUG
      PRINT *, "il_nb_ppp = ", il_nb_ppp
#endif
!
! Case where there is no link for this PE : 
      IF (il_nb_ppp .EQ. 0) THEN
        PRINT *, " Warning : No gridless point for this PE"
        CALL PSMILe_Flushstd
        gp%used_for_coupling = .false.
        fp%used_for_coupling = .false.
!  If no gridless point : no gridless grid for this PE
#ifdef VERBOSE
      print *, trim(ch_id), ': psmile_set_userdef : eof, ierror = ', ierror
      call PSMILe_Flushstd
#endif /* VERBOSE */

        return
      ENDIF
!
!  2.4  Allocate  arrays in global Userdefs structure  for the gridless grid
!
      ALLOCATE (up%iga_igl(1:il_nb_ppp,4), STAT=ierror)
      ALLOCATE (up%dga_wght(1:il_nb_ppp), STAT=ierror)
!
!  2.5  Select links that are relevant for this PE partition
!       indi, indj, indk are the indices of the geographical points
!       ilink : index of the link in the w&a file (or grless grid) 1 to nlinks
!       iloc  : local index in the local partition (1 to il_nb_ppp)
!
!   Case where at least one usr-defined link concerns this PE partition
      iloc = 0
!   Specific to One-dimensional grid
      if ( il_celldim == 1 ) then
         do ilink = 1, nlinks
            indi = ila_links(ilink)%cell_id(1)
            IF (indi .GE. indi_min .AND. indi .LE. indi_max ) THEN
               iloc = iloc + 1
               up%iga_igl(iloc,1) = indi + il_val11 - indi_min
               up%iga_igl(iloc,4) = ilink
               up%dga_wght(iloc)  = ila_links(ilink)%weight
            ENDIF
         enddo
!   Specific to Bi-dimensional grid
      elseif ( il_celldim == 2 ) then
         do ilink = 1, nlinks
            indi = ila_links(ilink)%cell_id(1)
            indj = ila_links(ilink)%cell_id(2)
            IF (indi .GE. indi_min .AND. indi .LE. indi_max &
               .AND. indj .GE. indj_min .AND. indj .LE. indj_max) THEN
               iloc = iloc + 1
               up%iga_igl(iloc,1) = indi + il_val11 - indi_min
               up%iga_igl(iloc,2) = indj + il_val12 - indj_min
               up%iga_igl(iloc,4) = ilink
               up%dga_wght(iloc)  = ila_links(ilink)%weight
            ENDIF
         enddo
!   Specific to Three-dimensional grid
      elseif ( il_celldim == 3 ) then
         do ilink = 1, nlinks
            indi = ila_links(ilink)%cell_id(1)
            indj = ila_links(ilink)%cell_id(2)
            indk = ila_links(ilink)%cell_id(3)
            IF (indi .GE. indi_min .AND. indi .LE. indi_max &
               .AND. indj .GE. indj_min .AND. indj .LE. indj_max &
               .AND. indk .GE. indk_min .AND. indk .LE. indk_max) THEN
               iloc = iloc + 1
               up%iga_igl(iloc,1) = indi + il_val11 - indi_min
               up%iga_igl(iloc,2) = indj + il_val12 - indj_min
               up%iga_igl(iloc,3) = indk + il_val13 - indk_min
               up%iga_igl(iloc,4) = ilink
               up%dga_wght(iloc)  = ila_links(ilink)%weight
            ENDIF
         enddo
      endif
!!!      endif                ! lg_nolink
!
!   Deallocate ila_links. Local info is kept in up%iga_igl and up%dga_wght
!
      DEALLOCATE ( ila_links, STAT=ierror )
      IF (ierror > 0) THEN
          ierrp (1) = ierror
          ierror = 14

          CALL PSMILe_Error ( ierror, 'ila_links', &
             ierrp, 1, __FILE__, __LINE__ )
          RETURN
      ENDIF


#ifdef DEBUG
      PRINT *," iga_igl(:,1:4) = "
      do i = 1, up%ig_nb_ppp
        PRINT *, up%iga_igl(i,:)
      enddo
#endif
!
! 2.6  find the number of subdomains of the gridlessgrid in my PE partition
!
!      1 subdomain = suite of contiguous indices in gridlessgrid (ilink values)
!
      il_nbr_blocks = 1  ! initial number of sub-domains of gridless grid
      il_ndim = 3        ! il_ndim is always 3 for gridless grid
!
      do iloc=2, il_nb_ppp
         if ( (up%iga_igl(iloc,4) - up%iga_igl(iloc-1,4)) .GT. 1) then
!      A new subdomain in the gridless grid indices is found
            il_nbr_blocks = il_nbr_blocks + 1
         endif
      enddo
#ifdef DEBUG
      PRINT *, "il_nbr_blocks = ", il_nbr_blocks
#endif
!
! 2.7   For the Gridless grid : allocate arrays for offsets and extents 
!       of the subdomains with exact dimensions
!
      ALLOCATE (il_offspr(il_nbr_blocks,il_ndim), STAT=ierror)
      ALLOCATE (il_extepr(il_nbr_blocks,il_ndim), STAT=ierror)

#ifdef DEBUG
      IF ( ierror /= 0 ) then
        PRINT *, 'Error allocating partition il_offspr or il_extepr = ', ierror
        CALL PSMILe_Flushstd
        CALL PSMILe_Abort
      ENDIF
#endif
!
!  2.8  compute index limits for this PE (offsets and extent of each subdomain)
!
      isbd = 1            ! Subdomain (block) index
      il_deb = 1          ! Keep the index value of the first point in subdomain
      il_offspr(isbd,1) = up%iga_igl(1,4) - 1           ! offset of first block
      il_extepr(isbd,1) = 1                             ! min number for extent
      do iloc=2, il_nb_ppp
         if ( (up%iga_igl(iloc,4) - up%iga_igl(iloc-1,4)) .GT. 1) then
!        A new block in the gridless grid indices is found
           il_extepr(isbd,1) = iloc - 1 - il_deb + 1   ! previous block extent
           il_deb = iloc                               ! a new bloc begin now
           isbd = isbd + 1                             ! increment block index
           il_offspr(isbd,1) = up%iga_igl(iloc,4) - 1  ! offset of new block
         endif
         il_extepr(isbd,1) = iloc - il_deb + 1         ! current block extent
      enddo
!
!  2.9  Gridless grid : index limits for this PE
!
      valid_shape_pr(1,1) = 1
      valid_shape_pr(2,1) = il_nb_ppp
      valid_shape_pr(1:2,2:4) = 1
!
      il_extepr(:,2:3) = 1
      il_offspr(:,2:3) = 0

      actual_shape_pr(1:2,1:4) = valid_shape_pr(1:2,1:4)
#ifdef DEBUG
      PRINT *, 'il_extepr = ',il_extepr
      PRINT *, 'il_offspr = ',il_offspr
      PRINT*, 'actual_shape_pr(:,1) :',actual_shape_pr(1,1),actual_shape_pr(2,1)
      PRINT*, 'actual_shape_pr(:,2) :',actual_shape_pr(1,2),actual_shape_pr(2,2)
      PRINT*, 'actual_shape_pr(:,3) :',actual_shape_pr(1,3),actual_shape_pr(2,3)
#endif
!
!  Build gridless grid
!
      grid_type_2 = PRISM_Gridless
      id_comp_id  = fp%comp_id
!  Uses transient local_name when generating the name of the gridless grid
!  (Matches the name generated by prismdrv_get_all_grids)
!
      call put_udef_suffix ( fp%local_name, grid_name_gl, chan_id, il_side )
#ifdef DEBUG
   PRINT *, ' psmile_set_userdef,  grid_name_gl = ',grid_name_gl
#endif

      CALL psmile_def_grid ( grid_id_2 , grid_name_gl, id_comp_id, valid_shape_pr, &
                            grid_type_2, ierror )
!
!  Keep values in Userdef structure
!
      up%igl_grid_id = grid_id_2
      up%status = PSMILe_Status_defined
!
#ifdef DEBUG
      IF ( ierror /= 0 ) PRINT *,'Error in gridless psmile_def_grid = ', ierror
#endif
!
      CALL psmile_def_partition ( grid_id_2, il_nbr_blocks, il_offspr, il_extepr, &
                                 ierror )
#ifdef DEBUG
      IF ( ierror /= 0 ) PRINT *,'Error in gridless psmile_def_partition = ',ierror
#endif
!
      CALL psmile_set_points_gridless (method_id_2, point_name_2, grid_id_2, new_points, ierror)

      ! make this method know to psmile_user_data, this way it will not be ignored in psmile_merge_fields
      call psmile_store_data_intern_points (method_id_2)
!
#ifdef DEBUG
      IF ( ierror /= 0 ) PRINT *,'Error in psmile_set_points_gridless = ', ierror
#endif
!
      iming = valid_shape_pr(1,1)
      imaxg = valid_shape_pr(2,1)
      jming = valid_shape_pr(1,2)
      jmaxg = valid_shape_pr(2,2)
      kming = valid_shape_pr(1,3)
      kmaxg = valid_shape_pr(2,3)
!
! ... allocate mask
!
      ALLOCATE ( gprmask(iming:imaxg,jming:jmaxg,kming:kmaxg), STAT=ierror )
!
#ifdef DEBUG
      IF (ierror /= 0) THEN
          PRINT *, 'Error in allocating gprmask arrays = ', ierror
          CALL PSMILe_Flushstd
          CALL PSMILe_Abort
      ENDIF
#endif
!
! Definition du masque
!
      DO k = kming, kmaxg
        DO j = jming, jmaxg
          DO i = iming, imaxg
            gprmask(i,j,k)=.TRUE.
          ENDDO
        ENDDO
      ENDDO
!
      CALL  psmile_set_mask ( mask_id_2, grid_id_2, valid_shape_pr, &
                             gprmask, new_mask, ierror )
#ifdef DEBUG
      IF ( ierror /= 0 ) PRINT *, 'Error in psmile_set_mask = ', ierror
#endif
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!! 3. declares the associates variable : gridless grid function
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!     Gridless functions is always 3dim + bundles
      il_var_nodims(1) = 3
!     Detect bundle number from Geographical Function
      if ( fp%transi_type == PSMILe_bundle ) then
          nbr_fields = fp%var_shape(2,gp%n_dim+1)
          il_var_nodims(1) = 4
          il_var_nodims(2) = nbr_fields
      else
          nbr_fields = 1
          il_var_nodims(2) = 0
      endif
      up%ig_nbr_fields = nbr_fields

!!! Set by def_var in fp%var_shape for the "gridfunction"
      actual_shape_pr(1,4) = 1
      actual_shape_pr(2,4) = nbr_fields
!
      id_var_type = fp%dataType
      call put_udef_suffix ( fp%local_name, var_name_gl, chan_id, il_side )
!
#ifdef DEBUG
   print *,' psmile_set_userdef: var_name_gl nbr_bundles = ',var_name_gl, il_var_nodims(2)
#endif
!
      CALL psmile_def_var (ass_var_id_2, var_name_gl, grid_id_2, method_id_2, &
           mask_id_2, il_var_nodims, actual_shape_pr, id_var_type, ierror )


      ! make this field know to psmile_user_data, this way it will not be ignored in psmile_merge_fields
      call psmile_store_data_intern_field (ass_var_id_2)
!
!    4. Updates geographical grid and field structures with associated userdef values
!
        gp%assoc_grid_id = grid_id_2
        gp%used_for_coupling = .false.
!
        if  ( il_side == 0 ) then
!       Geographic transient
           fp%Taskout(chan_id)%assoc_var_id  = ass_var_id_2
           fp%Taskout(chan_id)%userdef_id    = userdef_id_2
        elseif ( il_side == 1 ) then
!       Geographic transient
           fp%Taskin%In_channel(chan_id)%assoc_var_id  = ass_var_id_2
           fp%Taskin%In_channel(chan_id)%userdef_id    = userdef_id_2
        endif
!
#ifdef VERBOSE
      print *, trim(ch_id), ': psmile_set_userdef : eof, ierror = ', ierror
      call PSMILe_Flushstd
#endif /* VERBOSE */
!
      return
!
      END SUBROUTINE psmile_set_userdef
!
!
      SUBROUTINE psmile_ud_hdlerr2(istatus, lnotfound)
      integer     :: istatus
      logical     :: lnotfound
      INCLUDE 'netcdf.inc'
      lnotfound = .false.
      IF (istatus .NE. NF_NOERR) THEN
         if ( istatus .EQ. NF_ENOTVAR ) then
            lnotfound = .true.
            return
         endif 
         PRINT *, NF_STRERROR(istatus)
         STOP 'netcdf stopped'
      ENDIF
      return
      END SUBROUTINE psmile_ud_hdlerr2

      SUBROUTINE psmile_ud_hdlerr(istatus)
      INTEGER        :: istatus
      INCLUDE 'netcdf.inc'
      IF (istatus .NE. NF_NOERR) THEN
         PRINT *, NF_STRERROR(istatus)
         STOP 'netcdf stopped'
      ENDIF
      RETURN
      END SUBROUTINE psmile_ud_hdlerr