prismtrs_remap_conserv.F90

Go to the documentation of this file.

!------------------------------------------------------------------------
! Copyright 2006-2010, CERFACS, Toulouse, France.
! All rights reserved. Use is subject to OASIS4 license terms.
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: PRISMTrs_remap_conserv
!
! !INTERFACE
subroutine prismtrs_remap_conserv (     &
   grid1_size, grid2_size,              &
   grid1_corners, grid2_corners,        &
   id_src_lonlatz_size,                 &
   id_index_size1,                      &
   grid1_corner_lat, grid1_corner_lon,  &
   grid2_corner_lat, grid2_corner_lon,  &
   ida_nbsrccells_pertgtpt,             &
   ida_index_array,                     &
   ida_srcepio_add,                     &
   id_epio_id,                          &
   id_interp_id,                        &
   id_idim,                             &
   num_wts,                             &
   id_err)

!
! !USES:
  USE PRISMDrv, dummy_INTERFACE => PRISMTrs_remap_conserv
  IMPLICIT NONE
!
! !PARAMETERS:
!
  INTEGER :: grid1_size, grid2_size ! Source and target epio size
  !
  INTEGER :: grid1_corners, grid2_corners
  INTEGER :: id_src_lonlatz_size ! Size of array containing source corner info
  INTEGER :: id_index_size1 ! Size of array containing info of all source cells
  ! contributing to all epio target cells (sum of ida_nbsrccells_pertgtpt(:))
  DOUBLE PRECISION, DIMENSION (id_src_lonlatz_size)       :: grid1_corner_lat
  DOUBLE PRECISION, DIMENSION (id_src_lonlatz_size)       :: grid1_corner_lon
  DOUBLE PRECISION, DIMENSION (grid2_size, grid2_corners) :: grid2_corner_lat
  DOUBLE PRECISION, DIMENSION (grid2_size, grid2_corners) :: grid2_corner_lon
  ! Number of source cells contributing to the rremapping of each target cell
  INTEGER, DIMENSION(grid2_size)                          :: ida_nbsrccells_pertgtpt
  ! Index in id_src_lonlatz_size of each contributing source cell corner
  INTEGER, DIMENSION(id_index_size1,grid1_corners)        :: ida_index_array
  ! Index in grid1_size of each contributing source cell
  INTEGER, DIMENSION(id_index_size1)                      :: ida_srcepio_add
  INTEGER :: id_epio_id   ! EPIO Id
  INTEGER :: id_interp_id ! Interpolation Id
  INTEGER :: id_idim      ! Dimension concerned with the current interpolation
  INTEGER :: num_wts      ! Number of weights for 2D conservative remapping
  INTEGER :: id_err, il_add
!
! !RETURN VALUE
!
! !DESCRIPTION
!
! SUBROUTINE "PRISMTrs_remap_conserv" computes for weights with a
! 2d conservative remapping method, derived from the SCRIP 1.4
! library available from Los Alamos National Laboratory.
!
! !REVISED HISTORY
!   Date      Programmer   Description
! ----------  ----------   -----------
! 14/12/2006  S. Valcke    Created (based on SCRIP library)
!
! EOP
!----------------------------------------------------------------------
! $Id: prismtrs_remap_conserv.F90 2082 2009-10-21 13:31:19Z hanke $
! $Author: hanke $
!----------------------------------------------------------------------
!
! Local declarations
! 
  CHARACTER(LEN=len_cvs_string), SAVE  :: mycvs = 
     '$Id: prismtrs_remap_conserv.F90 2082 2009-10-21 13:31:19Z hanke $'
! ---------------------------------------------------------------------
!
! max number of subsegments per segment to prevent infinite loop
  INTEGER, PARAMETER ::  max_subseg = 10000
!
  INTEGER  ::   
     grid2_add,  ! current linear address for grid2 cell in tgt epio
     grid1_add,  ! current linear address for grid1 cell in src epio
     n, nwgt,    ! generic counters
     npart,      ! generic counters
     corner,     ! corner of cell that segment starts from
     next_corn,  ! corner of cell that segment ends on
     num_subseg, ! number of subsegments
     il_loop, il_count, il_bou, i, il_grid1,  ! loop counter
     num_srch_cells,  ! Number of search cell for each source and target cell
     il_nbrsrc_cells, ! Number of source cells treated in grid1 sweep
     il_cumul,   ! index in ida_index_array 
     il_action,  ! defines action in prismtrs_store_link_cnsrv: 
                  ! 0 if first call, 2 if last call, 1 otherwise
     il_stride, & ! stride in lat/lon corner array for Gaussian Red src grid
     nb_corner_grid1cell, &  ! real number of corners for source grid
     il_grid1_new   ! new index of neighbour cells when cells ovelarp to
                    ! keep only non masked cells
                    ! (for example when orca is source) 

#ifdef DEBUGX
  INTEGER, PARAMETER :: nbr_src_cells = 10
  INTEGER, DIMENSION(nbr_src_cells) :: grid_src ! number of source cells linked to grid_tgt
  INTEGER  :: grid_tgt  ! EPIO target number of the target cell under investigation
  INTEGER :: ii
  DOUBLE PRECISION, PARAMETER :: dbl_rad2deg = 360.0/6.2831853
#endif

  LOGICAL  :: 
     lcoinc,     ! flag for coincident segments
     lrevers,    ! flag for reversing direction of segment
     lbegin       ! flag for first integration of a segment

  LOGICAL ::  ll_gaussred ! flag true if source grid is Gaussian Reduced
  LOGICAL ::  ll_debug    ! for debug outputs

  LOGICAL, DIMENSION(:), ALLOCATABLE :: 
     srch_mask   ! mask for restricting searches

  DOUBLE PRECISION ::                     
     intrsct_lat, intrsct_lon,        ! lat/lon of next intersect
     beglat, endlat, beglon, endlon,  ! endpoints of current seg.
     norm_factor                       ! factor for normalizing wts

  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: 
     grid1cell_corner_lat, grid1cell_corner_lon ! corner of a given source cell
 
  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE ::    
     grid2_centroid_lat, grid2_centroid_lon ,      ! centroid coords on each grid
     grid1_centroid_lat, grid1_centroid_lon

  DOUBLE PRECISION, DIMENSION(grid1_size) ::   
     grid1_center_lat, grid1_center_lon,       ! center coords on each grid
     grid1_frac,                               ! intersected part of the cell
     grid1_larea                                ! cell area as estimated here

  DOUBLE PRECISION, DIMENSION(grid2_size) ::   
     grid2_center_lat, grid2_center_lon,       ! center coords on each grid
     grid2_frac,                               ! intersected part of the cell
     grid2_larea                                ! cell area as estimated here

  DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE ::    
     srch_corner_lat,   ! lat of each corner of srch cells
     srch_corner_lon     ! lon of each corner of srch cells

  DOUBLE PRECISION, DIMENSION(2) :: begseg ! begin lat/lon for full segment

  DOUBLE PRECISION, DIMENSION(6) :: lcl_weights ! local wgt array

  INTEGER, DIMENSION(id_index_size1) :: ila_grid2_assoc
  LOGICAL, DIMENSION(grid1_size)     :: ll_grid1done ! .true. if associated source grid has been treated
  INTEGER, DIMENSION(:), ALLOCATABLE :: srch_add,search_add_temp

  !
  ! ---------------------------------------------------------------------
  !
#ifdef VERBOSE
  PRINT *, '| | | | | | | Enter PRISMTrs_remap_conserv'
  call psmile_flushstd
#endif
  !
  nb_corner_grid1cell = grid1_corners
  !
#ifdef DEBUGX
  grid_tgt = 4382
  grid_src(1) = 18300
  grid_src(2) = 18381
  grid_src(3) = 18382
  grid_src(4) = 10583
  grid_src(5) = 10579
  grid_src(6) = 18379
  grid_src(7) = 18380
  grid_src(8) = 18383
  grid_src(9) = 18310
  grid_src(10) = 18305
!!$  grid_src(11) = 10610
!!$  grid_src(12) = 10554
!!$  grid_src(13) = 10558
!!$  grid_src(14) = 10548
#endif
  !
  ! Check if the source grid is Gaussian Reduced
  !
  ll_gaussred = .false.
  il_stride = 0

  IF (Drv_Epios(id_epio_id)%src_grid_type == PRISM_Gaussreduced_regvrt .OR. &
     Drv_Epios(id_epio_id)%src_grid_type == PRISM_Gaussreduced_sigmavrt) THEN
      ll_gaussred = .TRUE.
      il_stride = Drv_Epios(id_epio_id)%gaussred_stride
      nb_corner_grid1cell = 4

#ifdef DEBUG
      IF (il_stride == 0) THEN
          PRINT *, '| | | | | | Pb in prismtrs_remap_conserv: Stride in '
          PRINT *, '| | | | | | lat/lon corner array for Gaussian Red = 0'
          call psmile_abort
      ENDIF
#endif
  ENDIF
  !
  ! Allocate array that will contain the corner of a given source cell
  !
  ALLOCATE (grid1cell_corner_lat(nb_corner_grid1cell))
  ALLOCATE (grid1cell_corner_lon(nb_corner_grid1cell))
  !
  !  initialize centroid arrays
  !
      allocate( grid1_centroid_lat(grid1_size), &
                grid1_centroid_lon(grid1_size), &
                grid2_centroid_lat(grid2_size), &
                grid2_centroid_lon(grid2_size))

      grid1_centroid_lat = zero
      grid1_centroid_lon = zero
      grid2_centroid_lat = zero
      grid2_centroid_lon = zero
  !
      ! Initialize cell areas and fraction for local calculation
      grid1_larea = zero
      grid2_larea = zero
      grid1_frac  = zero
      grid2_frac  = zero
      !
      ! For now the center of the grid cells will be calculated
      grid1_center_lat = zero
      grid1_center_lon = zero
      grid2_center_lat = zero
      grid2_center_lon = zero
      !
      lcl_weights = zero
  ! integrate around each cell on grid1
  !
#ifdef VERBOSE
  PRINT *,'grid1 sweep '
  call psmile_flushstd
#endif

!  DO grid1_add = 1,grid1_size
        !***
        !*** create search arrays for each source mesh
        !***
!
      allocate(srch_mask(grid2_size))
!
  srch_mask (:) = .FALSE.
  ila_grid2_assoc(:) = 0
  il_count = 1
  !
  ! Define grid2_add associated to each source cell of ida_index_array
  ! and calculate center lat and lon
  !
  DO grid2_add = 1,grid2_size
    ila_grid2_assoc(il_count : il_count+ida_nbsrccells_pertgtpt(grid2_add)-1) = grid2_add
    il_count = il_count + ida_nbsrccells_pertgtpt(grid2_add)
    grid2_center_lat(grid2_add) = SUM(grid2_corner_lat(grid2_add,:))/grid2_corners
    grid2_center_lon(grid2_add) = SUM(grid2_corner_lon(grid2_add,:))/grid2_corners
  ENDDO
  !
  ll_grid1done(:) = .FALSE.
  il_nbrsrc_cells = 0
  il_action = 0

  PRINT *,'grid1 sweep '

  DO il_loop = 1, id_index_size1
    grid1_add = ida_srcepio_add(il_loop)
    !
    ! Define all target cells associated to each source cell of ida_index_array
    ! Each source cell may appear more than once in ida_index_array
    ! ll_grid1done indicates if this source cell was already treated or not
    !
    IF (.NOT. ll_grid1done(grid1_add)) THEN
        ll_grid1done(grid1_add) = .TRUE.
        ! Number of source cells treated
        il_nbrsrc_cells = il_nbrsrc_cells + 1
        !
#ifdef DEBUGX
        DO ii=1,nbr_src_cells
          IF (grid1_add == grid_src(ii)) THEN
              WRITE(88,*) 'grid1_add=', grid1_add
          ENDIF
        ENDDO
#endif

        ! Define the corners of the source cell
        !
        IF (ll_gaussred) THEN 
            grid1cell_corner_lon(1) = &
               grid1_corner_lon(ida_index_array(il_loop,1))
            grid1cell_corner_lon(2) = &
               grid1_corner_lon(ida_index_array(il_loop,1)+il_stride)
            grid1cell_corner_lon(3) = &
               grid1_corner_lon(ida_index_array(il_loop,1)+il_stride)
            grid1cell_corner_lon(4) = &
               grid1_corner_lon(ida_index_array(il_loop,1))
            grid1cell_corner_lat(1) = &
               grid1_corner_lat(ida_index_array(il_loop,1))
            grid1cell_corner_lat(2) = &
               grid1_corner_lat(ida_index_array(il_loop,1))
            grid1cell_corner_lat(3) = &
               grid1_corner_lat(ida_index_array(il_loop,1)+il_stride)
            grid1cell_corner_lat(4) = &
               grid1_corner_lat(ida_index_array(il_loop,1)+il_stride)
        ELSE
            grid1cell_corner_lat(:) = grid1_corner_lat(ida_index_array(il_loop,:))
            grid1cell_corner_lon(:) = grid1_corner_lon(ida_index_array(il_loop,:))
        ENDIF
        !
        ! Calculate lat and lon for centers of source grid cells 
        ! Hypothesis: it is an average of the corner lats and lons
        !
        DO il_bou = 1, nb_corner_grid1cell
          grid1_center_lat(grid1_add) = grid1_center_lat(grid1_add) + &
             grid1cell_corner_lat(il_bou)
          grid1_center_lon(grid1_add) = grid1_center_lon(grid1_add) + &
             grid1cell_corner_lon(il_bou)
        ENDDO
        grid1_center_lat(grid1_add) = grid1_center_lat(grid1_add)/nb_corner_grid1cell
        grid1_center_lon(grid1_add) = grid1_center_lon(grid1_add)/nb_corner_grid1cell
        !
        num_srch_cells = 0
        srch_mask(:) = .FALSE.
        DO il_bou = 1, id_index_size1
          !
          ! Identify all target cells for grid1_add source cell 
          IF (ida_srcepio_add(il_bou) == grid1_add) THEN
              num_srch_cells = num_srch_cells + 1
              srch_mask(ila_grid2_assoc(il_bou)) = .TRUE.
          ENDIF
        ENDDO
        !
        IF (num_srch_cells .GT. grid2_size) THEN
            PRINT *, '| | | | | | Number of search cell on target grid'
            PRINT *, '| | | | | | greater than epio target size ' 
            PRINT *, '| | | | | | Calling psmile_abort in prismtrs_remap_conserv'
            call psmile_abort 
        ENDIF
        !
        ALLOCATE(srch_add(num_srch_cells), &
           srch_corner_lat(grid2_corners,num_srch_cells), &
           srch_corner_lon(grid2_corners,num_srch_cells))
        srch_add(:) = 0
        srch_corner_lat(:,:) = zero
        srch_corner_lon(:,:) = zero
        !
        n = 0
        !
        ! Store lon and lat of target cells associated to grid1_add source cell
        !
        gather1: DO grid2_add = 1, grid2_size 
          if (srch_mask(grid2_add)) then
            n = n+1
            srch_add(n) = grid2_add
            srch_corner_lat(:,n) = grid2_corner_lat(grid2_add,:)
            srch_corner_lon(:,n) = grid2_corner_lon(grid2_add,:)
          endif
        end do gather1
        !
#ifdef DEBUGX
        DO ii=1,nbr_src_cells
          IF (grid1_add == grid_src(ii)) THEN
              WRITE(88,*)'    '
              WRITE(88,*)'  ** Grid1 cell and associated search cells **'
              WRITE(88,*) 'grid1_add=', grid1_add
              DO corner = 1, nb_corner_grid1cell
                WRITE(88,1111) corner, &
                   dbl_rad2deg*grid1cell_corner_lon(corner), &
                   dbl_rad2deg*grid1cell_corner_lat(corner)
              ENDDO
              WRITE(88,*) '  num_srch_cells=', num_srch_cells
              WRITE(88,*) '    '
              IF (num_srch_cells .NE. 0) & 
                 WRITE(88,*) '  srch_add(:)=', srch_add(:)
              DO n=1, num_srch_cells
                DO corner = 1,grid2_corners
                  WRITE(88,1112) n, dbl_rad2deg*srch_corner_lon(corner,n), &
                               dbl_rad2deg*srch_corner_lat(corner,n)
                END DO
              END DO
              WRITE(88,*)'    ***************************************'
              WRITE(88,*)'    '
          ENDIF
        ENDDO

 1111 format ('   grid1 corner, lon, lat = ', I2, 2X, F12.4, 2X, F12.4)
 1112 format ('     srch cell, lon, lat = ', I2, 2X, F12.4, 2X, F12.4) 
#endif
        !
        DO corner = 1,nb_corner_grid1cell
          next_corn = mod(corner,nb_corner_grid1cell) + 1

          !***
          !*** define endpoints of the current segment
          !***
          beglat = grid1cell_corner_lat(corner)
          beglon = grid1cell_corner_lon(corner)
          endlat = grid1cell_corner_lat(next_corn)
          endlon = grid1cell_corner_lon(next_corn)

          lrevers = .false.

          !***
          !*** to ensure exact path taken during both
          !*** sweeps, always integrate segments in the same 
          !*** direction (SW to NE).
          !***

          if ((endlat < beglat) .or.  &
            (endlat == beglat .and. endlon < beglon)) then 
              beglat = grid1cell_corner_lat(next_corn)
              beglon = grid1cell_corner_lon(next_corn)
              endlat = grid1cell_corner_lat(corner)
              endlon = grid1cell_corner_lon(corner)
              lrevers = .true.
#ifdef DEBUGX
        DO ii=1,nbr_src_cells
          IF (grid1_add == grid_src(ii)) THEN
              WRITE(88, *) ' sweep1 LREVERS TRUE'
              WRITE(88,*) 'grid1_add=', grid1_add
          ENDIF
        ENDDO
#endif
          endif
          begseg(1) = beglat
          begseg(2) = beglon
          lbegin = .true.
          num_subseg = 0

          !***
          !*** if this is a constant-longitude segment, skip the rest 
          !*** since the line integral contribution will be zero.
          !***
#ifdef DEBUGX
        DO ii=1,nbr_src_cells
          IF (grid1_add == grid_src(ii)) THEN
              IF (endlon .EQ. beglon) THEN
                  WRITE(88,1110) grid1_add,corner,next_corn
                  WRITE(88,1113) dbl_rad2deg*beglon, dbl_rad2deg*beglat 
                  WRITE(88,1114) dbl_rad2deg*endlon, dbl_rad2deg*endlat
                  WRITE(88, *)'  sweep1 endlon == beglon; skip segment'
                  WRITE(88,*) '             '
              ENDIF
          ENDIF
        ENDDO
 1110 FORMAT ('    grid1_add, corners = ', 2X, I6, 2X, I2, 2X, I2)
 1113 format ('   endlon == beglon;  beglon, beglat = ', &
              2X, F12.4, 2X, F12.4)
 1114 format ('   endlon == beglon;  endlon, endlat = ', &
              2X, F12.4, 2X, F12.4)
#endif

          if (endlon /= beglon) then
          !***
          !*** integrate along this segment, detecting intersections 
          !*** and computing the line integral for each sub-segment
          !***

              DO WHILE (beglat /= endlat .OR. beglon /= endlon)
              !***
              !*** prevent infinite loops if integration gets stuck
              !*** near cell or threshold boundary
              !***
                num_subseg = num_subseg + 1
                if (num_subseg > max_subseg) THEN
                    print*, 'ERROR=>integration stalled:' 
                    print*, 'num_subseg exceeded limit'
                    print*, '=>Verify corners in grids.nc, especially'
                    print*, 'if calculated by OASIS routine corners'           
                    PRINT *, 'integration stalled: num_subseg exceeded limit'
                    call psmile_abort
                endif

                !***
                !*** find next intersection of this segment with a grid
                !*** line on grid 2.
                !***

!                call timer_start(2)

#ifdef DEBUGX
        DO ii=1,nbr_src_cells
          IF (grid1_add == grid_src(ii)) THEN
              WRITE(88,1110) grid1_add,corner,next_corn
              WRITE(88,*) '             '
              WRITE(88,1115) dbl_rad2deg*beglon, dbl_rad2deg*beglat 
              WRITE(88,1116) dbl_rad2deg*endlon, dbl_rad2deg*endlat
              WRITE(88,*) '             '
          ENDIF
        ENDDO

 1115 format ('   avant intersection;  beglon, beglat = ', &
              2X, F12.4, 2X, F12.4)
 1116 format ('   avant intersection;  endlon, endlat = ', &
              2X, F12.4, 2X, F12.4)
#endif

                call prismtrs_intersection(grid2_add,intrsct_lat,intrsct_lon, &
                    lcoinc, beglat, beglon, endlat, endlon, begseg, &
                    lbegin, lrevers, &
                    num_srch_cells, grid2_corners, &
                    srch_corner_lat, srch_corner_lon, srch_add)

#ifdef DEBUGX
        DO ii=1,nbr_src_cells
!          IF (grid1_add == grid_src(ii) .AND. grid2_add==grid_tgt) THEN
          IF (grid1_add == grid_src(ii) ) THEN
              WRITE(88,*) 'grid1_add=', grid1_add
              WRITE(88,*) ' After call intersection, grid2_add', &
                                grid2_add
              WRITE(88,1117) dbl_rad2deg*beglon, dbl_rad2deg*beglat 
              WRITE(88,1118) dbl_rad2deg*endlon, dbl_rad2deg*endlat
              WRITE(88,1119) dbl_rad2deg*intrsct_lon, dbl_rad2deg*intrsct_lat
              WRITE(88,*) '   '
          ENDIF
        ENDDO

 1117 format('   après intersection;  beglon, beglat =           ', &
              2X, F12.4, 2X, F12.4)
 1118 format ('   après intersection;  endlon, endlat =          ', &
              2X, F12.4, 2X, F12.4)
 1119 format ('   après intersection; intrsct_lon, intrsct_lat = ', &
                    2X, F12.4, 2X, F12.4)
#endif

!                call timer_stop(2)
                lbegin = .false.

                !***
                !*** compute line integral for this subsegment.
                !***

                ! call timer_start(3)
                if (grid2_add /= 0) THEN
                    !
                    call prismtrs_line_integral(lcl_weights, &
                               beglon, intrsct_lon, beglat, intrsct_lat, &
                               grid1_center_lon(grid1_add), &
                               grid2_center_lon(grid2_add), &
                               num_wts)

#ifdef DEBUGX
        DO ii=1,nbr_src_cells
!          IF (grid1_add == grid_src(ii) .AND. grid2_add==grid_tgt) THEN
          IF (grid1_add == grid_src(ii)) THEN
              WRITE(88,*) 'grid1_add=', grid1_add
              WRITE(88,*) '  A1) WEIGHTS for this subsegment =', &
                                   lcl_weights(1)
              WRITE(88,*) '     '
          ENDIF
        ENDDO
#endif

                else
                    call prismtrs_line_integral(lcl_weights, &
                               beglon, intrsct_lon, beglat, intrsct_lat, &
                               grid1_center_lon(grid1_add), &
                               grid1_center_lon(grid1_add), &
                               num_wts)

#ifdef DEBUGX
        DO ii=1,nbr_src_cells
          IF (grid1_add == grid_src(ii)) THEN
              WRITE(88,*) 'grid1_add=', grid1_add
              WRITE(88,*) '  B1) WEIGHTS for this subsegment =', &
                                   lcl_weights(1)
              WRITE(88,*) '     ' 
          ENDIF
        ENDDO
#endif

                endif
                ! call timer_stop(3)
                !***
                !*** if integrating in reverse order, change
                !*** sign of lcl_weights
                !***

                if (lrevers) then
                    lcl_weights = -lcl_weights
#ifdef DEBUGX
        DO ii=1,nbr_src_cells
!          IF (grid1_add == grid_src(ii) .AND. grid2_add==grid_tgt) THEN
          IF (grid1_add == grid_src(ii)) THEN
              WRITE(88,*) 'grid1_add=', grid1_add
              WRITE(88,*) '  LREVERS; WEIGHTS for this subsegment =', &
                              lcl_weights(1)
              WRITE(88,*) '     '
          ENDIF
        ENDDO
#endif
                endif
                !***
                !*** store the appropriate addresses and lcl_weights. 
                !*** also add contributions to cell areas and centroids.
                !***


                if (grid2_add /= 0) then
                    if (Drv_Epios(id_epio_id)%src_mask_pointer(grid1_add) == 1) then
                        ! call timer_start(4)
                        call prismtrs_store_link_cnsrv(grid1_add, grid2_add, & 
                            lcl_weights, il_action, id_epio_id, num_wts)

#ifdef DEBUGX
        DO ii=1,nbr_src_cells
!          IF (grid1_add == grid_src(ii) .AND. grid2_add==grid_tgt) THEN
          IF (grid1_add == grid_src(ii)) THEN
              WRITE(88,*) '      after store_link_cnsrv'
              WRITE(88,1120) grid1_add, grid2_add,dbl_rad2deg*beglon, dbl_rad2deg*beglat, &
                             dbl_rad2deg*intrsct_lon,dbl_rad2deg*intrsct_lat,lcl_weights(1)
          ENDIF
        ENDDO

 1120 format ('      STORE add1,add2,blon,blat,ilon,ilat,WEIGHTS=', &
               1X,I6,1X,I6,1X,F12.8,1X,F12.8,1X,F12.8,1X,F12.8,1X,E15.8)
#endif

                        il_action = 1
                        ! call timer_stop(4)
                        grid1_frac(grid1_add) = &
                           grid1_frac(grid1_add)+lcl_weights(1)
                        grid2_frac(grid2_add) = &
                           grid2_frac(grid2_add)+lcl_weights(num_wts+1)
                    endif
                endif

                grid1_larea(grid1_add) = &
                   grid1_larea(grid1_add) + lcl_weights(1)
                 
                grid1_centroid_lat(grid1_add) = &
                   grid1_centroid_lat(grid1_add) + lcl_weights(2)
                grid1_centroid_lon(grid1_add) = &
                   grid1_centroid_lon(grid1_add) + lcl_weights(3)

                !***
                !*** reset beglat and beglon for next subsegment.
                !***

                beglat = intrsct_lat
                beglon = intrsct_lon

              END DO

          endif

          !***
          !*** end of segment
          !***

        END DO

        !***
        !*** finished with this cell: deallocate search array and
        !*** start on next cell
        deallocate(srch_add, srch_corner_lat, srch_corner_lon)
    ENDIF                        
  END DO
#ifdef VERBOSE
  PRINT *, 'Number of source cells in epio = ', grid1_size
  PRINT *, 'Number of source cells in grid1 sweep = ', il_nbrsrc_cells
#endif
  IF (il_nbrsrc_cells /= grid1_size) THEN 
      PRINT *, '| | |WARNING: Number of source cells in grid1 sweep '
      PRINT *, '| | |not equal to number of source cells in epio '
  ENDIF

      deallocate(srch_mask)
 
      print *,'grid1 end sweep '
!-----------------------------------------------------------------------
!
!     integrate around each cell on grid2
!
!-----------------------------------------------------------------------

      print *,'grid2 sweep '
      il_cumul = 0
      DO grid2_add = 1,grid2_size

        il_grid1_new = 0

#ifdef DEBUGX
        IF (grid2_add == grid_tgt) THEN 
            WRITE(88,*) 'grid2_add=', grid2_add
        ENDIF
#endif

        num_srch_cells = ida_nbsrccells_pertgtpt (grid2_add)

        ALLOCATE(srch_add(num_srch_cells),search_add_temp(num_srch_cells),&
                 srch_corner_lat(nb_corner_grid1cell,num_srch_cells), &
                 srch_corner_lon(nb_corner_grid1cell,num_srch_cells))

        gather2: DO il_grid1 = 1, num_srch_cells
            il_cumul = il_cumul + 1
            srch_add(il_grid1) = ida_srcepio_add (il_cumul)
#ifdef DEBUGX
            IF (grid2_add == grid_tgt) THEN
                WRITE(88,*), 'srch_add(il_grid1) :', grid2_add,il_grid1,srch_add(il_grid1)
            ENDIF
#endif
            IF (Drv_Epios(id_epio_id)%src_mask_pointer(srch_add(il_grid1)) == 1) THEN
            il_grid1_new = il_grid1_new + 1
            search_add_temp(il_grid1_new) = srch_add(il_grid1)
            IF (ll_gaussred) THEN
                srch_corner_lon(1,il_grid1_new) = &
                   grid1_corner_lon(ida_index_array(il_cumul,1))
                srch_corner_lon(2,il_grid1_new) = &
                   grid1_corner_lon(ida_index_array(il_cumul,1)+il_stride)
                srch_corner_lon(3,il_grid1_new) = &
                   grid1_corner_lon(ida_index_array(il_cumul,1)+il_stride)
                srch_corner_lon(4,il_grid1_new) = &
                   grid1_corner_lon(ida_index_array(il_cumul,1))
                srch_corner_lat(1,il_grid1_new) = &
                   grid1_corner_lat(ida_index_array(il_cumul,1))
                srch_corner_lat(2,il_grid1_new) = &
                   grid1_corner_lat(ida_index_array(il_cumul,1))
                srch_corner_lat(3,il_grid1_new) = &
                   grid1_corner_lat(ida_index_array(il_cumul,1)+il_stride)
                srch_corner_lat(4,il_grid1_new) = &
                   grid1_corner_lat(ida_index_array(il_cumul,1)+il_stride)  
            ELSE
                srch_corner_lat(:,il_grid1_new) = grid1_corner_lat(ida_index_array(il_cumul,:))
                srch_corner_lon(:,il_grid1_new) = grid1_corner_lon(ida_index_array(il_cumul,:))
            ENDIF
        ENDIF
        end do gather2

        num_srch_cells = il_grid1_new
        srch_add(1:num_srch_cells)=search_add_temp(1:il_grid1_new)

#ifdef DEBUGX
        IF (grid2_add == grid_tgt) THEN 
            WRITE(88,*)'    '
            WRITE(88,*)'  ** Grid2 cell and associated search cells **'
            DO corner = 1,grid2_corners 
              WRITE(88,1131) corner, &
                 dbl_rad2deg*grid2_corner_lon(grid2_add,corner), &
                 dbl_rad2deg*grid2_corner_lat(grid2_add,corner)
            ENDDO
            WRITE(88,*) '  num_srch_cells=', num_srch_cells
            WRITE(88,*) '    '
            WRITE(88,*) '  srch_add(:)=', srch_add(:)
            DO n=1, num_srch_cells
              DO corner = 1,nb_corner_grid1cell
                WRITE(88,1132) n, dbl_rad2deg*srch_corner_lon(corner,n), &
                              dbl_rad2deg*srch_corner_lat(corner,n)
              END DO
            END DO
            WRITE(88,*)'    ***************************************'
            WRITE(88,*)'    '
        ENDIF
 1131 format ('   grid2 corner, lon, lat = ', I2, 2X, F12.4, 2X, F12.4)
 1132 format ('     srch cell, lon, lat = ', I2, 2X, F12.4, 2X, F12.4) 
#endif

        DO corner = 1,grid2_corners
          next_corn = mod(corner,grid2_corners) + 1

          beglat = grid2_corner_lat(grid2_add,corner)
          beglon = grid2_corner_lon(grid2_add,corner)
          endlat = grid2_corner_lat(grid2_add,next_corn)
          endlon = grid2_corner_lon(grid2_add,next_corn)
          lrevers = .false.

          !***
          !*** to ensure exact path taken during both
          !*** sweeps, always integrate in the same direction
          !***

          if ((endlat < beglat) .or. &
              (endlat == beglat .and. endlon < beglon)) then 
            beglat = grid2_corner_lat(grid2_add,next_corn)
            beglon = grid2_corner_lon(grid2_add,next_corn)
            endlat = grid2_corner_lat(grid2_add,corner)
            endlon = grid2_corner_lon(grid2_add,corner)
            lrevers = .true.
#ifdef DEBUGX
            IF (grid2_add == grid_tgt) WRITE(88, *) ' sweep2 LREVERS TRUE'
#endif
          endif
          begseg(1) = beglat
          begseg(2) = beglon
          lbegin = .true.

          !***
          !*** if this is a constant-longitude segment, skip the rest 
          !*** since the line integral contribution will be zero.
          !***

#ifdef DEBUGX
          IF (grid2_add == grid_tgt) THEN
              IF (endlon .EQ. beglon) THEN
                  WRITE(88,1113) dbl_rad2deg*beglon, dbl_rad2deg*beglat 
                  WRITE(88,1114) dbl_rad2deg*endlon, dbl_rad2deg*endlat
                  WRITE(88, *)'  sweep2 endlon == beglon; skip segment'
                  WRITE(88,*) '             '
              ENDIF
          ENDIF
#endif

          if (endlon /= beglon) then
          num_subseg = 0

          !***
          !*** integrate along this segment, detecting intersections 
          !*** and computing the line integral for each sub-segment
          !***

          do while (beglat /= endlat .or. beglon /= endlon)

            !***
            !*** prevent infinite loops if integration gets stuck
            !*** near cell or threshold boundary
            !***

            num_subseg = num_subseg + 1
            if (num_subseg > max_subseg) THEN
                print*, 'ERROR=>integration stalled:' 
                print*, 'num_subseg exceeded limit'
                print*, '=>Verify corners in grids.nc, especially'
                print*, 'if calculated by OASIS routine corners' 
                PRINT *, 'integration stalled: num_subseg exceeded limit'
                call psmile_abort
            endif

            !***
            !*** find next intersection of this segment with a line 
            !*** on grid 2.
            !***

            !call timer_start(6)

#ifdef DEBUGX
            IF (grid2_add == grid_tgt) THEN
                WRITE(88,*) '             '
                WRITE(88,1115) dbl_rad2deg*beglon, dbl_rad2deg*beglat 
                WRITE(88,1116) dbl_rad2deg*endlon, dbl_rad2deg*endlat
                WRITE(88,*) '             '
            ENDIF
#endif
            call prismtrs_intersection(grid1_add,intrsct_lat,intrsct_lon,lcoinc, &
                              beglat, beglon, endlat, endlon, begseg, &
                              lbegin, lrevers, &
                              num_srch_cells, nb_corner_grid1cell, &
                              srch_corner_lat, srch_corner_lon, srch_add)

#ifdef DEBUGX
            IF (grid2_add == grid_tgt) THEN
                WRITE(88,*) ' After call intersection, grid1_add', &
                   grid1_add
                WRITE(88,1117) dbl_rad2deg*beglon, dbl_rad2deg*beglat 
                WRITE(88,1118) dbl_rad2deg*endlon, dbl_rad2deg*endlat
                WRITE(88,1119) dbl_rad2deg*intrsct_lon, dbl_rad2deg*intrsct_lat
                WRITE(88,*) '   '
            ENDIF
#endif

            !call timer_stop(6)
            lbegin = .false.

            !***
            !*** compute line integral for this subsegment.
            !***

            ! call timer_start(7)
            if (grid1_add /= 0) then

                ! Check if grid1 center was calculated 
                IF (.NOT. ll_grid1done(grid1_add)) THEN
                    PRINT *, &
                       '| | | | Center of grid1 mesh not known - cannot call prismtrs_line_integral'
                    call psmile_abort
                ENDIF

              call prismtrs_line_integral(lcl_weights, &
                              beglon, intrsct_lon, beglat, intrsct_lat, &
                              grid1_center_lon(grid1_add), &
                              grid2_center_lon(grid2_add), &
                              num_wts)

#ifdef DEBUGX
              IF (grid2_add == grid_tgt) THEN
                  WRITE(88,*) '  A2) WEIGHTS for this subsegment =', &
                                  lcl_weights(1)
                  WRITE(88,*) '     ' 
              ENDIF
#endif

            else
              call prismtrs_line_integral(lcl_weights, &
                             beglon, intrsct_lon, beglat, intrsct_lat, &
                             grid2_center_lon(grid2_add), &
                             grid2_center_lon(grid2_add), &
                             num_wts)

#ifdef DEBUGX
             IF (grid2_add == grid_tgt) THEN
                 WRITE(88,*) '  B2) WEIGHTS for this subsegment =', &
                     lcl_weights(1)
                  WRITE(88,*) '     ' 
              ENDIF
#endif

            endif
            ! call timer_stop(7)

            if (lrevers) then
              lcl_weights = -lcl_weights
#ifdef DEBUGX
              IF (grid2_add == grid_tgt) THEN
                WRITE(88,*) '  LREVERS; WEIGHTS for this subsegment =', &
                   lcl_weights(1)
                WRITE(88,*) '     '
              ENDIF
#endif
            endif
            !***
            !*** store the appropriate addresses and weights. 
            !*** also add contributions to cell areas and centroids.
            !*** if there is a coincidence, do not store weights
            !*** because they have been captured in the previous loop.
            !*** the grid1 mask is the master mask
            !***

#ifdef DEBUGX
            IF (lcoinc .AND. grid2_add == grid_tgt) &
               WRITE(88,*) '  LCOINC is TRUE; weight not stored'
#endif

            if (.not. lcoinc .and. grid1_add /= 0) then
              if (Drv_Epios(id_epio_id)%src_mask_pointer(grid1_add) == 1) then
                  ! call timer_start(8)
                  call prismtrs_store_link_cnsrv(grid1_add, grid2_add, & 
                     lcl_weights, il_action, id_epio_id, num_wts)

#ifdef DEBUGX
                  IF (grid2_add == grid_tgt) THEN
                      WRITE(88,*) '      after store_link_cnsrv'
                      WRITE(88,1120) grid1_add, grid2_add,dbl_rad2deg*beglon, dbl_rad2deg*beglat, &
                                     dbl_rad2deg*intrsct_lon,dbl_rad2deg*intrsct_lat,lcl_weights(1)
                  ENDIF
#endif

                  il_action = 1
                  ! call timer_stop(8)
                  grid1_frac(grid1_add) = &
                     grid1_frac(grid1_add) + lcl_weights(1)
                  grid2_frac(grid2_add) = &
                     grid2_frac(grid2_add) + lcl_weights(num_wts+1)
              endif

            endif

            grid2_larea(grid2_add) = &
               grid2_larea(grid2_add) + lcl_weights(num_wts+1)
            grid2_centroid_lat(grid2_add) = &
               grid2_centroid_lat(grid2_add) + lcl_weights(num_wts+2)
            grid2_centroid_lon(grid2_add) = &
               grid2_centroid_lon(grid2_add) + lcl_weights(num_wts+3)

            !***
            !*** reset beglat and beglon for next subsegment.
            !***

            beglat = intrsct_lat
            beglon = intrsct_lon

          end DO

          END if

          !***
          !*** end of segment
          !***

        END DO

        !***
        !*** finished with this cell: deallocate search array and
        !*** start on next cell

        DEALLOCATE(srch_add, srch_corner_lat, srch_corner_lon, search_add_temp)

#ifdef DEBUGX
        IF (grid2_add == grid_tgt) THEN
        DO n=1,Drv_Epios(id_epio_id)%num_links_map1
          IF (Drv_Epios(id_epio_id)%grid2_add_map1(n) == grid_tgt) THEN
              WRITE(88,*) 'grid1, grid2, weight= ', &
                 Drv_Epios(id_epio_id)%grid1_add_map1(n), Drv_Epios(id_epio_id)%grid2_add_map1(n), &
                 Drv_Epios(id_epio_id)%wts_map1(1,n)
          ENDIF
        END DO
    ENDIF
#endif

    END DO

      print *,'grid2 end sweep '

!-----------------------------------------------------------------------
!
!     correct for situations where N/S pole not explicitly included in
!     grid (i.e. as a grid corner point). if pole is missing from only
!     one grid, need to correct only the area and centroid of that 
!     grid.  if missing from both, do complete weight calculation.
!
!-----------------------------------------------------------------------
      !*** North Pole
      lcl_weights(1) =  dble_pi2
      lcl_weights(2) =  dble_pi*dble_pi
      lcl_weights(3) =  zero
      lcl_weights(4) =  dble_pi2
      lcl_weights(5) =  dble_pi*dble_pi
      lcl_weights(6) =  zero

      grid1_add = 0
      pole_loop1: do n=1,grid1_size
        if (grid1_larea(n) < -three*dble_pih .and. &
            grid1_center_lat(n) > zero) then
          grid1_add = n
          exit pole_loop1
        endif
      end do pole_loop1

      grid2_add = 0
      pole_loop2: do n=1,grid2_size
        if (grid2_larea(n) < -three*dble_pih .and. &
            grid2_center_lat(n) > zero) then
          grid2_add = n
          exit pole_loop2
        endif
      end do pole_loop2

      if (grid1_add /=0) then
        grid1_larea(grid1_add) = grid1_larea(grid1_add) + lcl_weights(1)
        grid1_centroid_lat(grid1_add) = &
        grid1_centroid_lat(grid1_add) + lcl_weights(2)
        grid1_centroid_lon(grid1_add) = &
        grid1_centroid_lon(grid1_add) + lcl_weights(3)
      endif

      if (grid2_add /=0) then
        grid2_larea(grid2_add) = grid2_larea(grid2_add) + &
                                        lcl_weights(num_wts+1)
        grid2_centroid_lat(grid2_add) = &
        grid2_centroid_lat(grid2_add) + lcl_weights(num_wts+2)
        grid2_centroid_lon(grid2_add) = &
        grid2_centroid_lon(grid2_add) + lcl_weights(num_wts+3)
      endif

      if (grid1_add /= 0 .and. grid2_add /=0) then
          call prismtrs_store_link_cnsrv(grid1_add, grid2_add, & 
             lcl_weights, il_action, id_epio_id, num_wts)
          il_action = 1
        grid1_frac(grid1_add) = &
           grid1_frac(grid1_add) + lcl_weights(1)
        grid2_frac(grid2_add) = &
           grid2_frac(grid2_add) + lcl_weights(num_wts+1)
      endif

      !*** South Pole
      lcl_weights(1) =  dble_pi2
      lcl_weights(2) = -dble_pi*dble_pi
      lcl_weights(3) =  zero
      lcl_weights(4) =  dble_pi2
      lcl_weights(5) = -dble_pi*dble_pi
      lcl_weights(6) =  zero

      grid1_add = 0
      pole_loop3: do n=1,grid1_size
        if (grid1_larea(n) < -three*dble_pih .and. &
            grid1_center_lat(n) < zero) then
          grid1_add = n
          exit pole_loop3
        endif
      end do pole_loop3

      grid2_add = 0
      pole_loop4: do n=1,grid2_size
        if (grid2_larea(n) < -three*dble_pih .and. &
            grid2_center_lat(n) < zero) then
          grid2_add = n
          exit pole_loop4
        endif
      end do pole_loop4

      if (grid1_add /=0) then
        grid1_larea(grid1_add) = grid1_larea(grid1_add) + lcl_weights(1)
        grid1_centroid_lat(grid1_add) = &
        grid1_centroid_lat(grid1_add) + lcl_weights(2)
        grid1_centroid_lon(grid1_add) = &
        grid1_centroid_lon(grid1_add) + lcl_weights(3)
      endif

      if (grid2_add /=0) then
        grid2_larea(grid2_add) = grid2_larea(grid2_add) + &
                                        lcl_weights(num_wts+1)
        grid2_centroid_lat(grid2_add) = &
        grid2_centroid_lat(grid2_add) + lcl_weights(num_wts+2)
        grid2_centroid_lon(grid2_add) = &
        grid2_centroid_lon(grid2_add) + lcl_weights(num_wts+3)
      endif

      if (grid1_add /= 0 .and. grid2_add /=0) then
        call prismtrs_store_link_cnsrv(grid1_add, grid2_add, lcl_weights, 1, &
           id_epio_id, num_wts)
        grid1_frac(grid1_add) = grid1_frac(grid1_add) + & 
                                lcl_weights(1)
        grid2_frac(grid2_add) = grid2_frac(grid2_add) + &
                                lcl_weights(num_wts+1)
      endif

!-----------------------------------------------------------------------
!
!     finish centroid computation
!
!-----------------------------------------------------------------------

      where (grid1_larea /= zero)
        grid1_centroid_lat = grid1_centroid_lat/grid1_larea
        grid1_centroid_lon = grid1_centroid_lon/grid1_larea
      end where

      where (grid2_larea /= zero)
        grid2_centroid_lat = grid2_centroid_lat/grid2_larea
        grid2_centroid_lon = grid2_centroid_lon/grid2_larea
      end where
!
!-----------------------------------------------------------------------
!
!     include centroids in weights and normalize using destination
!     area if requested
!
!-----------------------------------------------------------------------

      do n=1,Drv_Epios(id_epio_id)%num_links_map1

        grid1_add = Drv_Epios(id_epio_id)%grid1_add_map1(n)
        grid2_add = Drv_Epios(id_epio_id)%grid2_add_map1(n)
        do nwgt=1,num_wts
          lcl_weights(nwgt) = Drv_Epios(id_epio_id)%wts_map1(nwgt,n)
        end do
#ifdef DEBUGX
        IF (grid2_add == grid_tgt) THEN
            WRITE(88,*), 'Source and weight before normalisation for target point ', grid_tgt
            WRITE(88,*), grid1_add, lcl_weights(1)
        ENDIF
#endif

        select case(Drv_Interps(id_interp_id)%arg4(id_idim))
        case (PSMILe_destarea)
          if (grid2_larea(grid2_add) /= zero) then
            !if (luse_grid2_area) then
            !  norm_factor = one/grid2_area_in(grid2_add)
            !else
              norm_factor = one/grid2_larea(grid2_add)
            !endif
          else
            norm_factor = zero
          endif
        case (PSMILe_fracarea)
          if (grid2_frac(grid2_add) /= zero) then
            !if (luse_grid2_area) then
            !  norm_factor = grid2_larea(grid2_add)/
     !&                     (grid2_frac(grid2_add)*
     !&                      grid2_larea_in(grid2_add))
            !else
              norm_factor = one/grid2_frac(grid2_add)
            !endif
          else
            norm_factor = zero
          endif
        case (PSMILe_none)
          norm_factor = one
        end select

#ifdef DEBUGX
        IF (grid2_add == grid_tgt) THEN
            WRITE(88,*), 'Normalisation coef', grid_tgt,norm_factor
        ENDIF
#endif

        Drv_Epios(id_epio_id)%wts_map1(1,n)= lcl_weights(1)*norm_factor
        Drv_Epios(id_epio_id)%wts_map1(2,n)=(lcl_weights(2) - lcl_weights(1)* &
           grid1_centroid_lat(grid1_add))* norm_factor
        Drv_Epios(id_epio_id)%wts_map1(3,n)=(lcl_weights(3) - lcl_weights(1)* &
           grid1_centroid_lon(grid1_add))* norm_factor

#ifdef DEBUGX
        IF (grid2_add == grid_tgt) THEN
            WRITE(88,*), 'Source and weight after normalisation for target point ', grid_tgt
            WRITE(88,*), grid1_add, Drv_Epios(id_epio_id)%wts_map1(1,n)
            WRITE(88,*), ''
        ENDIF
#endif
      end do
      print *, 'Total number of links = ',Drv_Epios(id_epio_id)%num_links_map1

      where (grid1_larea /= zero) &
         grid1_frac = &
         grid1_frac/grid1_larea
      where (grid2_larea /= zero) &
         grid2_frac = &
         grid2_frac/grid2_larea

!
!-----------------------------------------------------------------------
!
!     perform some error checking on final weights
!
!-----------------------------------------------------------------------

      !grid2_centroid_lat = zero
      !grid2_centroid_lon = zero
#ifdef DEBUG
      do n=1,grid1_size
        if (grid1_larea(n) < -.01) then
          print *,'Grid 1 area error: n, area, mask =' &
                ,n,grid1_larea(n), &
                Drv_Epios(id_epio_id)%src_mask_pointer(n)
        endif
        if (grid1_centroid_lat(n) < -dble_pih-.01 .or. &
            grid1_centroid_lat(n) >  dble_pih+.01) then
            print *,'Grid 1 centroid lat error: n, centroid_lat, mask=' &
                ,n,grid1_centroid_lat(n), &
                Drv_Epios(id_epio_id)%src_mask_pointer(n)
        endif
        grid1_centroid_lat(n) = zero
        grid1_centroid_lon(n) = zero
      end do

      do n=1,grid2_size
        if (grid2_larea(n) < -.01) then
            PRINT *,'Grid 2 area error:  n, area, mask =' &
            ,n,grid2_larea(n), &
            Drv_Epios(id_epio_id)%tgt_mask_pointer(n)
        endif
        if (grid2_centroid_lat(n) < -dble_pih-.01 .or. &
            grid2_centroid_lat(n) >  dble_pih+.01) then
          PRINT *,'Grid 2 centroid lat error: n, centroid_lat, mask =' &
                ,n,grid2_centroid_lat(n), &
                Drv_Epios(id_epio_id)%tgt_mask_pointer(n)
        endif
        grid2_centroid_lat(n) = zero
        grid2_centroid_lon(n) = zero
      end do

! Rene: print statements prohibit vectorisation, better put this in a VERBOSE statement
!   to spped up producton runs. Stripmining of loop to avoid problems on the SX.

      do npart = 1, Drv_Epios(id_epio_id)%num_links_map1, 5000
      do n = npart, min(Drv_Epios(id_epio_id)%num_links_map1,npart+5000-1)

        grid1_add = Drv_Epios(id_epio_id)%grid1_add_map1(n)
        grid2_add = Drv_Epios(id_epio_id)%grid2_add_map1(n)
        
        if (Drv_Epios(id_epio_id)%wts_map1(1,n) < -.01) then
          print *,'Map 1 weight < 0 '
          PRINT *, &
              'grid1_add, grid2_add, wts_map1, grid1_mask, grid2_mask', &
              grid1_add, grid2_add, Drv_Epios(id_epio_id)%wts_map1(1,n), & 
              Drv_Epios(id_epio_id)%src_mask_pointer(grid1_add), &
              Drv_Epios(id_epio_id)%tgt_mask_pointer(grid2_add) 
        endif
        if (Drv_Interps(id_interp_id)%arg4(id_idim) /= PSMILe_none .and. Drv_Epios(id_epio_id)%wts_map1(1,n) > 1.01) then
          print *,'Map 1 weight > 1 '
          PRINT *, &
              'grid1_add, grid2_add, wts_map1, grid1_mask, grid2_mask', &
              grid1_add, grid2_add, Drv_Epios(id_epio_id)%wts_map1(1,n), & 
              Drv_Epios(id_epio_id)%src_mask_pointer(grid1_add), &
              Drv_Epios(id_epio_id)%tgt_mask_pointer(grid2_add) 
        endif
        grid2_centroid_lat(grid2_add) = &
        grid2_centroid_lat(grid2_add) + Drv_Epios(id_epio_id)%wts_map1(1,n)

      end do
      end do

      do n=1,grid2_size
        select case(Drv_Interps(id_interp_id)%arg4(id_idim))
        case (PSMILe_destarea)
          norm_factor = grid2_frac(n)
        case (PSMILe_fracarea)
          norm_factor = one
        case (PSMILe_none)
          ! if (luse_grid2_area) then
          !  norm_factor = grid2_area_in(grid2_add)
          ! else
            norm_factor = grid2_larea(n)
          ! endif
        end select
        if (abs(grid2_centroid_lat(n)-norm_factor) > .01) then
          print *, &
      'Error sum wts map1:grid2_add,grid2_centroid_lat,norm_factor, mask=' &
      ,n,grid2_centroid_lat(n), &
      norm_factor,Drv_Epios(id_epio_id)%tgt_mask_pointer(n)  
        endif
      end do
#endif
      ! Call prismtrs_store_link_cnsrv juste to deallocate local save arrays
      !
      il_action = il_action + 1
      IF (il_action == 2) THEN
          call prismtrs_store_link_cnsrv(grid1_add, grid2_add, lcl_weights, &
             il_action, id_epio_id, num_wts)
      ENDIF
      !
!-----------------------------------------------------------------------
!
!     deallocate allocated arrays
!
!----------------------------------------------------------------------- 

      deallocate (grid1_centroid_lat, grid1_centroid_lon, &
                  grid2_centroid_lat, grid2_centroid_lon)
      DEALLOCATE (grid1cell_corner_lat, grid1cell_corner_lon)

!-----------------------------------------------------------------------
!
#ifdef VERBOSE
  PRINT *, '| | | | | | | Quit PRISMTrs_remap_conserv'
  call psmile_flushstd
#endif
END SUBROUTINE PRISMTrs_remap_conserv