prism_enddef.F90

Go to the documentation of this file.

!-----------------------------------------------------------------------
! Copyright 2006-2010, CERFACS, Toulouse, France.
! Copyright 2006-2010, SGI Germany, Munich, Germany.
! Copyright 2006-2010, NEC Europe Ltd., London, UK.
! All rights reserved. Use is subject to OASIS4 license terms.
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: PRISM_Enddef
!
! !INTERFACE:

      subroutine prism_enddef (ierror)

!
! !USES:
!
      use PRISM, dummy_interface => prism_enddef
      use psmile_user_data, only : psmile_apply_user_data, psmile_merge_fields
      use PSMILe
      use PSMILe_SMIOC, only : sga_smioc_comp, transient
      use psmile_timer, only : psmile_timer_stop

      implicit none
!
! !OUTPUT PARAMETERS:
!
      integer, Intent (Out)               :: ierror

!     Returns the error code of prism_enddef;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !LOCAL VARIABLES
!
!     a_comps (1, *) = Global component id of component
!     a_comps (2, *) = Number of grids for the component
!     a_comps (3, *) = Local component id
!
!     Second dimension: Number_of_Comps_allocated is sufficent
!                       since only active components of actual
!                       process are stored in "acomps".
!
      integer, parameter           :: nd_acomps = 3
      integer                      :: a_comps (nd_acomps, noComponents)
      Integer                      :: n_active
      Integer                      :: my_icomp0_coll_comps
      Integer                      :: icomp, lastc
      Integer                      :: use_Grid(Number_of_Grids_allocated)
!
      Type (Transient), pointer    :: sga_smioc_transi(:)
      Logical                      :: found

      Integer                      :: i, j, n
!
      logical, save                :: first = .true.
!
!  ... for userdef interpolations
!
      integer                      :: nbr_in, nbr_out, il_side, il_i, il_o
      integer                      :: il_interp_meth
      integer                      :: il_userdef_id, il_ass_var_id

!  ... for intersections found

      integer                      :: ninter, nmyint, nnull, max_num_intersect
      integer, allocatable         :: num_intersect_per_grid(:), num_dummy_intersect_per_grid(:)

!  ... for communication

      integer                      :: lastag
      integer                      :: grid_id
      logical                      :: get_halo
      integer                      :: color
      integer                      :: key
!
!  ... for error parameters

      Integer, parameter           :: nerrp = 3
      Integer                      :: ierrp (nerrp)
!
#ifdef TIMING
!  ... for timing prism_enddef

      Double Precision            :: tic, toc
#endif
!
! !DESCRIPTION:
!
! Subroutine "prism_enddef" finish the definition phase of the
! PSMILe library.
!
! !SEE ALSO:
!
!     psmile_get_act_comp.F90
!
! !REVISION HISTORY:
!   Date      Programmer   Description
! ----------  ----------   -----------
! 01.12.02    H. Ritzdorf  created
!
!EOP
!----------------------------------------------------------------------
!
! $Id: prism_enddef.F90 2846 2011-01-04 12:02:30Z hanke $
! $Author: hanke $
!
  Character(len=len_cvs_string), save :: mycvs = 
      '$Id: prism_enddef.F90 2846 2011-01-04 12:02:30Z hanke $'
!
!----------------------------------------------------------------------

#ifdef VERBOSE
      print 9990, trim(ch_id)
      call psmile_flushstd
#endif /* VERBOSE */

      if ( .not. first ) then
         ierror = -1
#ifdef VERBOSE
         print *, trim(ch_id), ': prism_enddef: It is not allowed to call prism_enddef twice.'
         print *, trim(ch_id), ': prism_enddef: eof ierror =', ierror
#endif /* VERBOSE */
         return
      endif

#ifdef TIMING
      call MPI_Barrier ( comm_psmile, ierror )
      tic=MPI_Wtime()
#endif
!
!  Initialization
!
      ierror = 0
      lastag = PSMILe_Enddef_Tag
      first  = .false.
      use_Grid = PSMILe_Undef
      get_halo = .false.
!
!  Generate psmile data structures from data provided by the user which is stored in psmile_user_data
!
      call psmile_apply_user_data(ierror)

!
!  Check existence of Userdef Interpolations and fill Userdefs structures
!
#ifdef DEBUG
            print *, trim(ch_id),': prism_enddef: Fields allocated = ',Number_of_Fields_allocated
#endif
      do i = 1, Number_of_Fields_allocated

         if ( Fields(i)%status  == PSMILe_status_defined .and. &
              Fields(i)%used_for_coupling ) then

!  First check if we are on the source or on the target side (or both)
!     for each Field "used_in_coupling" : find
!     nbr_in  = number of "In channels"
!     nbr_out = number of "Taskout structures"
            nbr_in  = 0
            nbr_out = 0
!
            nbr_in = Fields(i)%Taskin%nbr_inchannels
            if ( Associated(Fields(i)%Taskout) ) then
               nbr_out = size (Fields(i)%Taskout)
            endif
#ifdef DEBUG
            print *, trim(ch_id),': prism_enddef: index Field = ',i
            print *, trim(ch_id),': prism_enddef: nbr_in nbr_out = ',nbr_in,nbr_out
#endif /* DEBUG */

! check if a User_Defined interpolation is done on this transient field
! look for interpolation method (PSMILe_user3D).
!
! Loop on all In-channels
         do il_i = 1, nbr_in
            il_side = 1
!  Initialize "userdef" variables to PSMILe_undef for all Fields
            Fields(i)%Taskin%In_channel(il_i)%assoc_var_id  = PSMILe_undef
            Fields(i)%Taskin%In_channel(il_i)%userdef_id    = PSMILe_undef
            Grids(Methods(Fields(i)%method_id)%grid_id)%assoc_grid_id = PSMILe_undef
!
            il_interp_meth = Fields(i)%Taskin%In_channel(il_i)%interp%interp_meth(1)
            if ( il_interp_meth == PSMILe_user3D ) then
!
!  Read weight and addresses file and define associated gridless grid and function
               call psmile_set_userdef(i, il_side, il_i, ierror)
#ifdef DEBUG
      print *, trim(ch_id), ': prism_enddef: index Field = ',i
      print *, trim(ch_id), ': prism_enddef: index channel in ',il_i
      print *, trim(ch_id), ': prism_enddef: assoc_var_id = ',  &
                               Fields(i)%Taskin%In_channel(il_i)%assoc_var_id
      print *, trim(ch_id), ': prism_enddef: userdef_id = ',    &
                               Fields(i)%Taskin%In_channel(il_i)%userdef_id
      print *, trim(ch_id), ': prism_enddef: assoc_grid_id = ', &
                               Grids(Methods(Fields(i)%method_id)%grid_id)%assoc_grid_id
#endif /* DEBUG */
            endif
         enddo        !  In channels
!
!  Loop on all channels "Out"
         do il_o = 1, nbr_out
            il_side = 0
!  Initialize "userdef" variables to PSMILe_undef for all Fields
            Fields(i)%Taskout(il_o)%assoc_var_id  = PSMILe_undef
            Fields(i)%Taskout(il_o)%userdef_id    = PSMILe_undef
            Grids(Methods(Fields(i)%method_id)%grid_id)%assoc_grid_id = PSMILe_undef
!
            il_interp_meth = Fields(i)%Taskout(il_o)%interp%interp_meth(1)
            if ( il_interp_meth == PSMILe_user3D ) then
!
!  Read weight and addresses file and define associated gridless grid and function
               call psmile_set_userdef(i, il_side, il_o, ierror)
#ifdef DEBUG
      print *, trim(ch_id), ': prism_enddef: index Field = ',i
      print *, trim(ch_id), ': prism_enddef: index channel out ',il_o
      print *, trim(ch_id), ': prism_enddef: assoc_var_id = ',  &
                               Fields(i)%Taskout(il_o)%assoc_var_id
      print *, trim(ch_id), ': prism_enddef: userdef_id = ',    &
                               Fields(i)%Taskout(il_o)%userdef_id
      print *, trim(ch_id), ': prism_enddef: assoc_grid_id = ', &
                               Grids(Methods(Fields(i)%method_id)%grid_id)%assoc_grid_id
#endif /* DEBUG */
            endif
         enddo                !  Out channnels
!
         endif             ! Field is used for coupling
!
      end do            ! loop on Field index (var_id)
!
!  Control Grids
!
      call psmile_control_grids (ierror)
      if (ierror /= 0) return
!
!  Check existence of matching prism_def_var call for every SMIOC field name
!
      do n = 1, Number_of_Comps_allocated
         if ( Comps(n)%status == PSMILe_status_defined ) then
            sga_smioc_transi => sga_smioc_comp(n)%sga_smioc_transi

            do i = 1, size(sga_smioc_transi)

               found = .false.
               do j = 1, Number_of_Fields_allocated
                  if ( Fields(j)%status == PSMILe_status_defined ) then
                     if ( trim(sga_smioc_transi(i)%cg_local_name) == trim(adjustl(Fields(j)%local_name)) ) &
                        found = .true.
                  endif
               enddo
               if ( .not. found ) then
                  ierror =  PRISM_Warn_NoDefVar
                  ierrp (1) = j
                  print *, trim(ch_id), ': transient with name ', &
                       trim(sga_smioc_transi(i)%cg_local_name), ' failed.'
                   call psmile_warning ( ierror,  trim(sga_smioc_transi(i)%cg_local_name), &
                        ierrp(1), 1, __FILE__, __LINE__ )
               endif

            enddo

         endif
      enddo
!
!  For all components: Get number of grids to be coupled and
!                      get sorted list of components
!
      call psmile_get_act_comps (a_comps, nd_acomps, n_act_comp, ierror)
      if (ierror /= 0) return
!
!     Component info structure
!
      Allocate (comp_infos(1:n_act_comp), STAT = ierror)
      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = n_act_comp
         call psmile_error ( PRISM_Error_Alloc, 'comp_infos', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif
!
!     comp_infos(1:n_act_comp)%global_comp_id = a_comps (1, 1:n_act_comp)
!     comp_infos(1:n_act_comp)%comp_id = a_comps (3, 1:n_act_comp)
!     comp_infos(1:n_act_comp)%size = Comps(comp_infos(1:n_act_comp)%comp_id)%size
!
!===> Prepare search for reduced Gauss grids
!
      do grid_id = 1, Number_of_Grids_allocated
      if ( Grids(grid_id)%grid_type == PRISM_Gaussreduced_regvrt ) then
         call psmile_gauss_setup ( grid_id, ierror )
         if (ierror /= 0) return
      endif
      enddo
!
!  For all components: Get extensions of grids
!
         do i = 1, n_act_comp
         call psmile_enddef_comp (a_comps(3,i), a_comps (1,i), &
                                  a_comps(2,i), comp_infos(i), ierror)
         if (ierror /= 0) return
         end do
!
!  Collect the data for all applications
!
!   my_icomp0_coll_comps = 
!     Index 0 of component information of actual application
!     in collected component information "all_comp_infos";
!     i.e. the component information of actual application is stored in
!     all_comp_infos(my_icomp0_coll_comps+1:
!                    my_icomp0_coll_comps+n_active)
!     and n_active = Number of active components in actual application
!
      call psmile_enddef_appl (lastag, my_icomp0_coll_comps, &
                               n_active, ierror)
      if (ierror /= 0) return

#ifdef TIMING
      call MPI_Barrier ( comm_psmile, ierror )
      toc=MPI_Wtime()
      print *, trim(ch_id), ': prism_enddef: Time for initial procedures ', toc-tic
#endif

#ifdef __PSMILE_WITH_IO
!
!  Complete I/O for meta data here which will allow us to perform
!  debugging with a play-back tool based on the debug output.
!  JL
#ifdef PRISM_ASSERTION
      IF (PRISM_noCompsPerProc /= n_act_comp) THEN
          CALL psmile_assert ( __FILE__, __LINE__, &
                'Call to psmile_io_init_pelist supposes that PRISM_noCompsPerProc = n_act_comp')
      ENDIF
#endif
      do i = 1, PRISM_noCompsPerProc
         call psmile_io_init_pelist (i, comp_infos(i),  ierror)
         if (ierror /= 0) return         
      enddo

      call psmile_enddef_metadata(ierror)

      if ( ierror .ne. 0 ) then
        ierrp(1)=ierror
        call psmile_error(PRISM_Error_IO_Meta, &
                         'psmile_enddef_metadata', &
                              ierrp, 1, __FILE__, __LINE__ )
      endif

      if ( Appl%stand_alone ) then
#ifdef VERBOSE
         print *, trim(ch_id), ': prism_enddef: Stand alone application. We return'
         print *, trim(ch_id), ': prism_enddef: eof ierror =', ierror
#endif /* VERBOSE */
         return
      endif


#endif

#ifdef TIMING
      call MPI_Barrier ( comm_psmile, ierror )
      tic=MPI_Wtime()
#endif
!
!  Look for intersections of own grids with other components
!  and send partner data on the extension of the intersection
!
      ninter = 0
      nmyint = 0
      nnull  = 0
!
      lastc = 0

#define CONS_REMAP_DEADLOCK_WORKAROUND
#ifdef CONS_REMAP_DEADLOCK_WORKAROUND
!     computes the maximum number of intersections
      max_num_intersect = 1
      do i = 1,  Number_of_coll_comps

         ! the computation gives a size which is way to big, but I currently do
         ! not know how to caluculate the real maximum number of intersections
         ! this part is a workaround, so it might be better to soft the actual problem
         max_num_intersect = max_num_intersect * &
            SUM(all_comp_infos(i)%Number_of_grids_vector(1:all_comp_infos(i)%size))
      enddo

      call psmile_flushstd
      allocate (paction%intersect_ranks(max_num_intersect), stat = ierror)
      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = max_num_intersect
         call psmile_error ( PRISM_Error_Alloc, 'intersect_ranks', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif
      paction%intersect_ranks = -1
#endif

      allocate (num_intersect_per_grid(Number_of_Grids_allocated), &
                num_dummy_intersect_per_grid(Number_of_Grids_allocated), stat = ierror)
      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = Number_of_Grids_allocated
         call psmile_error ( PRISM_Error_Alloc, 'num_intersect_per_grid, num_dummy_intersect_per_grid', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif

      num_intersect_per_grid(:)        = 0
      num_dummy_intersect_per_grid (:) = 0

      do i = 1, n_act_comp
!
!===> ... Get index of comp_infos(i) in all_comp_infos
!
        do icomp = lastc+1, n_active
           if (all_comp_infos(my_icomp0_coll_comps+icomp)%global_comp_id ==   &
                                            comp_infos(i)%global_comp_id) exit
        end do
!
        if (icomp > n_active) then
           ierror = PRISM_Error_Internal
           ierrp (1) = comp_infos(i)%global_comp_id
           ierrp (2) = lastc
           ierrp (3) = n_active

           call psmile_error ( ierror, 'Global Comp Id not found in all_comp_infos', &
                               ierrp, 3, __FILE__, __LINE__ )
        endif

        lastc = icomp
!
!===> ... Look for intersections for component "comp_infos(i)"
!
        call psmile_find_intersect (comp_infos(i), my_icomp0_coll_comps+icomp, &
                                    num_intersect_per_grid, &
                                    num_dummy_intersect_per_grid, &
                                    ninter, nmyint, nnull, lastag, ierror)
        if (ierror /= 0) return
      end do
!
!-----------------------------------------------------------------------
!     Get the halo regions
!-----------------------------------------------------------------------
!
      color = 0
      key   = 0

      do n = 1, Number_of_Comps_allocated

        if ( Comps(n)%status /= PSMILe_status_free ) then

        do i = 1, Number_of_Fields_allocated
           if ( Fields(i)%smioc_loc /= PSMILe_Undef .and. Fields(i)%comp_id == n ) then
              if ( associated(Fields(i)%Taskout) ) then
                 grid_id = Methods(Fields(i)%method_id)%grid_id
                 if ( Associated (Grids(grid_id)%partition) ) then
                    if ( ( Grids(grid_id)%grid_type == PRISM_Irrlonlat_regvrt &
#ifdef TODO
! see psmile_get_halo_indices.F90
                      .or. Grids(grid_id)%grid_type == PRISM_Irrlonlatvrt
#endif
                         )  .and. &
                         use_Grid(grid_id) == PSMILe_Undef ) then
                         use_Grid(grid_id) = Grids(grid_id)%global_grid_id
                         get_halo = .true.
                         color    = 1
                    endif
                 endif
              endif
           endif
        enddo

        call MPI_Comm_split (Comps(n)%comm, color, key, Comps(n)%comm_halo, ierror)

        if ( get_halo ) then

           call psmile_get_halo_indices ( n, use_Grid, ierror )

           call psmile_get_halo_points ( n, ierror )

        endif

        endif

      enddo ! n-loop
!
!-----------------------------------------------------------------------
!     Receive extension of the intersection and
!     send request for coordinates
!-----------------------------------------------------------------------
!
      call psmile_get_intersect (ninter, nmyint, nnull, &
                                 num_intersect_per_grid, &
                                 num_dummy_intersect_per_grid, &
                                 lastag, ierror)
      if (ierror /= 0) return
!
#ifdef TIMING
      call MPI_Barrier ( comm_psmile, ierror )
      toc=MPI_Wtime()
      print *, trim(ch_id), ': prism_enddef: Time for search', toc-tic
#endif
!
!-----------------------------------------------------------------------
!     Free memory allocated
!-----------------------------------------------------------------------
!
      call psmile_mg_clean (ierror)
      if (ierror /= 0) return
!
!     Note: This memory is freed by freeing corresponnding
!           entries of "all_comp_infos"
!
!     do i = 1, n_act_comp
!        Deallocate (comp_info(i)%Number_of_Grids_vector)
!        Deallocate (comp_info(i)%all_extents)
!        Deallocate (comp_info(i)%all_extent_infos)
!        Deallocate (comp_info(i)%psmile_ranks)
!     end do
!
#define CONS_REMAP_DEADLOCK_WORKAROUND
#ifdef CONS_REMAP_DEADLOCK_WORKAROUND
      Deallocate (paction%intersect_ranks)
#endif
      Deallocate (num_intersect_per_grid,num_dummy_intersect_per_grid)
      Deallocate (comp_infos)
!
      if (Number_of_coll_comps > 0) then
         do i = 1, Number_of_coll_comps
            Deallocate (all_comp_infos(i)%Number_of_Grids_vector)
            Deallocate (all_comp_infos(i)%all_extents)
            Deallocate (all_comp_infos(i)%all_extent_infos)
            Deallocate (all_comp_infos(i)%psmile_ranks)
         end do
!
         Deallocate (all_comp_infos)
         Number_of_coll_comps = 0
      endif

!
!-----------------------------------------------------------------------
!     Get data from restart file
!-----------------------------------------------------------------------
!
      call psmile_get_restart ( ierror )

      if ( ierror .ne. 0 ) then
        ierrp(1)=ierror
        call psmile_error(999, &
                         'psmile_get_restart', &
                              ierrp, 1, __FILE__, __LINE__ )
      endif

#ifdef DEBUG
!
! Print out information/status about grids, fields and methods 
!
      do i = 1, Number_of_Grids_allocated
        call psmile_print_grid_info (i)
      enddo

      do i = 1, Number_of_Fields_allocated
        call psmile_print_field_info (i)
      enddo

      do i = 1, Number_of_Methods_allocated
        call psmile_print_method_info (i)
      enddo
#endif
!
!-----------------------------------------------------------------------
!     Merge fields
!-----------------------------------------------------------------------
!
      call psmile_merge_fields(ierror)
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id), ierror
      call psmile_flushstd
#endif /* VERBOSE */
!
!  time profile measurements (optional)
!
#ifdef PROFILE
      call psmile_timer_stop(2)
#endif
!
!  Formats:
!
#ifdef VERBOSE
9990 format (1x, a, ': prism_enddef: start')
9980 format (1x, a, ': prism_enddef: eof ierror =', i3)
#endif /* VERBOSE */
!
      end subroutine prism_enddef