prismtrs_intersection.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_intersection
!
! !INTERFACE!
      subroutine prismtrs_intersection(location,intrsct_lat,intrsct_lon,lcoinc, &
                             beglat, beglon, endlat, endlon, begseg, &
                             lbegin, lrevers, &
                             num_srch_cells, id_nbcorners, &
                             srch_corner_lat, srch_corner_lon, srch_add)

!-----------------------------------------------------------------------
!
!     intent(in): 
!
!-----------------------------------------------------------------------
     !
        USE PRISMDrv, dummy_INTERFACE => PRISMTrs_intersection
        IMPLICIT NONE

        INTEGER, INTENT(in) :: num_srch_cells, id_nbcorners
        DOUBLE PRECISION, DIMENSION(id_nbcorners, num_srch_cells), INTENT(in) ::
           srch_corner_lat, srch_corner_lon
        INTEGER, DIMENSION(num_srch_cells), INTENT(in) :: srch_add
        !
      logical, intent(in) :: 
          lbegin,  ! flag for first integration along this segment
          lrevers ! flag whether segment integrated in reverse

      DOUBLE PRECISION, intent(in) :: 
          beglat, beglon,  ! beginning lat/lon endpoints for segment
          endlat, endlon   ! ending    lat/lon endpoints for segment

      DOUBLE PRECISION, dimension(2), intent(inout) :: 
          begseg ! begin lat/lon of full segment

!-----------------------------------------------------------------------
!
!     intent(out): 
!
!-----------------------------------------------------------------------

      integer, intent(out) :: 
             location  ! address in destination array containing this
                        ! segment

      logical, intent(out) :: 
             lcoinc    ! flag segments which are entirely coincident
                        ! with a grid line

      DOUBLE PRECISION, intent(out) :: 
          intrsct_lat, intrsct_lon ! lat/lon coords of next intersect.
!
! !DESCRIPTION
!
! SUBROUTINE PRISMTrs_intersection finds the next intersection of a 
!     destination grid 
!     line with the line segment given by beglon, endlon, etc.
!     a coincidence flag is returned if the segment is entirely 
!     coincident with an ocean grid line.  the cells in which to search
!     for an intersection must have already been restricted in the
!     calling routine.
!     This routine comes from the SCRIP 1.4 library available 
!     from Los Alamos National Laboratory.
!
! !REVISED HISTORY
!   Date      Programmer   Description
! ----------  ----------   -----------
! 21/03/2007  S. Valcke    Duplicated from the SCRIP library
!
! EOP
!----------------------------------------------------------------------
! $Id: prismtrs_intersection.F90 2082 2009-10-21 13:31:19Z hanke $
! $Author: hanke $
!----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!     local variables
!
!-----------------------------------------------------------------------
      CHARACTER(LEN=len_cvs_string), SAVE  :: mycvs = 
     '$Id: prismtrs_intersection.F90 2082 2009-10-21 13:31:19Z hanke $'
      integer :: n, next_n, cell, srch_corners
! for test of non-convexe cell
      INTEGER :: next2_n, neg, pos  

      integer, save :: 
          last_loc  ! save location when crossing threshold

      logical :: 
          loutside  ! flags points outside grid

      logical, save :: 
          lthresh = .false.  ! flags segments crossing threshold bndy

      DOUBLE PRECISION :: 
          lon1, lon2,       ! local longitude variables for segment
          lat1, lat2,       ! local latitude  variables for segment
          grdlon1, grdlon2,  ! local longitude variables for grid cell
          grdlat1, grdlat2,  ! local latitude  variables for grid cell
          vec1_lat, vec1_lon,  ! vectors and cross products used
          vec2_lat, vec2_lon,  ! during grid search
          cross_product, 
          eps, offset,        ! small offset away from intersect
          s1, s2, determ,     ! variables used for linear solve to
          mat1, mat2, mat3, mat4, rhs1, rhs2,  ! find intersection
          rl_halfpi, rl_v2lonmpi2, rl_v2lonppi2

      DOUBLE PRECISION, save :: 
          intrsct_lat_off, intrsct_lon_off ! lat/lon coords offset 
                                            ! for next search

#ifdef DEBUGX
      INTEGER, PARAMETER :: nbr_tgt_cells = 13
      INTEGER, DIMENSION(nbr_tgt_cells)  :: grid_tgt  ! EPIO target number of the target cell under investigation
      INTEGER :: ii
      DOUBLE PRECISION, PARAMETER :: dbl_rad2deg = 360.0/6.2831853
#endif

#ifdef DEBUGX
  PRINT *, '| | | | | | | Enter PRISMTrs_intersection'
  call psmile_flushstd
#endif
!-----------------------------------------------------------------------
!
!     initialize defaults, flags, etc.
!
!-----------------------------------------------------------------------
#ifdef DEBUGX
  grid_tgt(1)= 7215
  grid_tgt(2)= 7216
  grid_tgt(3)= 7251
  grid_tgt(4)= 7252
  grid_tgt(5)= 7253
  grid_tgt(6)= 7254
  grid_tgt(7)= 7279
  grid_tgt(8)= 7280
  grid_tgt(9)= 7281
  grid_tgt(10)= 7282
  grid_tgt(11)= 7290
  grid_tgt(12)= 7291
  grid_tgt(13)= 7292
#endif

      location = 0
      lcoinc = .false.
      intrsct_lat = endlat
      intrsct_lon = endlon

      if (num_srch_cells == 0) return

      if (beglat > north_thresh .or. beglat < south_thresh) then

        if (lthresh) location = last_loc
        call prismtrs_pole_intersection(location, &
                     intrsct_lat,intrsct_lon,lcoinc,lthresh, &
                     beglat, beglon, endlat, endlon, begseg, lrevers, &
                     num_srch_cells, id_nbcorners, &
                     srch_corner_lat, srch_corner_lon, srch_add)
        if (lthresh) then
          last_loc = location
          intrsct_lat_off = intrsct_lat
          intrsct_lon_off = intrsct_lon
        endif
        return

      endif

      loutside = .false.
      if (lbegin) then
        lat1 = beglat
        lon1 = beglon
      else
        lat1 = intrsct_lat_off
        lon1 = intrsct_lon_off
      endif
      lat2 = endlat
      lon2 = endlon
      if ((lon2-lon1) > three*dble_pih) then
        lon2 = lon2 - dble_pi2
      else if ((lon2-lon1) < -three*dble_pih) then
        lon2 = lon2 + dble_pi2
      endif
      s1 = zero
!-----------------------------------------------------------------------
!
!     search for location of this segment in ocean grid using cross
!     product method to determine whether a point is enclosed by a cell
!
!-----------------------------------------------------------------------

!      call timer_start(12)
      srch_corners = size(srch_corner_lat,DIM=1)
      srch_loop: do

        !***
        !*** if last segment crossed threshold, use that location
        !***

        if (lthresh) then
          do cell=1,num_srch_cells
            if (srch_add(cell) == last_loc) then
              location = last_loc
              eps = tiny
              exit srch_loop
            endif
          end do
        endif

        !***
        !*** otherwise normal search algorithm
        !***

        cell_loop: DO cell=1,num_srch_cells
          corner_loop: do n=1,srch_corners
            next_n = MOD(n,srch_corners) + 1

            !***
            !*** here we take the cross product of the vector making 
            !*** up each cell side with the vector formed by the vertex
            !*** and search point.  if all the cross products are 
            !*** positive, the point is contained in the cell.
            !***

            vec1_lat = srch_corner_lat(next_n,cell) -  &
                      srch_corner_lat(n     ,cell)
            vec1_lon = srch_corner_lon(next_n,cell) -  &
                      srch_corner_lon(n     ,cell)
            vec2_lat = lat1 - srch_corner_lat(n,cell)
            vec2_lon = lon1 - srch_corner_lon(n,cell)

            !***
            !*** if endpoint coincident with vertex, offset
            !*** the endpoint
            !***

            if (vec2_lat == 0 .and. vec2_lon == 0) then
              lat1 = lat1 + 1.d-10*(lat2-lat1)
              lon1 = lon1 + 1.d-10*(lon2-lon1)
              vec2_lat = lat1 - srch_corner_lat(n,cell)
              vec2_lon = lon1 - srch_corner_lon(n,cell)
            ENDIF

            !***
            !*** check for 0,2pi crossings
            !***

            if (vec1_lon >  dble_pi) then
              vec1_lon = vec1_lon - dble_pi2
            else if (vec1_lon < -dble_pi) then
              vec1_lon = vec1_lon + dble_pi2
            endif
            if (vec2_lon >  dble_pi) then
              vec2_lon = vec2_lon - dble_pi2
            else if (vec2_lon < -dble_pi) then
              vec2_lon = vec2_lon + dble_pi2
            endif

            cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat

            !***
            !*** if the cross product for a side is zero, the point 
            !***   lies exactly on the side or the side is degenerate
            !***   (zero length).  if degenerate, set the cross 
            !***   product to a positive number.  otherwise perform 
            !***   another cross product between the side and the 
            !***   segment itself. 
            !*** if this cross product is also zero, the line is 
            !***   coincident with the cell boundary - perform the 
            !***   dot product and only choose the cell if the dot 
            !***   product is positive (parallel vs anti-parallel).
            !***

            if (cross_product == zero) then
              if (vec1_lat /= zero .or. vec1_lon /= zero) then
                vec2_lat = lat2 - lat1
                vec2_lon = lon2 - lon1

                if (vec2_lon >  dble_pi) then
                  vec2_lon = vec2_lon - dble_pi2
                else if (vec2_lon < -dble_pi) then
                  vec2_lon = vec2_lon + dble_pi2
                endif

                cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat
              else
                cross_product = one
              endif

              if (cross_product == zero) then
                lcoinc = .true.
                cross_product = vec1_lon*vec2_lon + vec1_lat*vec2_lat
                if (lrevers) cross_product = -cross_product
              endif
            endif

#ifdef DEBUGX
            DO ii = 1,nbr_tgt_cells
              IF (srch_add(cell) == grid_tgt(ii) ) THEN
                  WRITE(88,*) 'lrevers :',lrevers
                  WRITE(88,1110) srch_add(cell),n,next_n
                  WRITE(88,1111) dbl_rad2deg*lon1,dbl_rad2deg*lat1
                  WRITE(88,1112) dbl_rad2deg*lon2,dbl_rad2deg*lat2
                  WRITE(88,1113) cross_product
              ENDIF
            ENDDO
!
 1110 FORMAT ('    grid2 cell, corners = ', 2X, I6, 2X, I2, 2X, I2)
 1111 format ('    lon1, lat1 = ', 2X, F12.4, 2X, F12.4)
 1112 format ('    lon2, lat2 = ', 2X, F12.4, 2X, F12.4)
 1113 format ('    cross_product  = ', 2X, F12.6)
#endif

            !***
            !*** if cross product is less than zero, this cell
            !*** doesn't work
            !***
            if (cross_product < zero) exit corner_loop
          end do corner_loop

          !***
          !*** if cross products all positive, we found the location
          !***

          if (n > srch_corners) then
            location = srch_add(cell)

            !***
            !*** if the beginning of this segment was outside the
            !*** grid, invert the segment so the intersection found
            !*** will be the first intersection with the grid
            !***
            if (loutside) then
               !*** do a test to see if the cell really is outside 
               !*** or if it is a non-convexe cell 
               neg=0
               pos=0
               do n=1,srch_corners
                  next_n = MOD(n,srch_corners) + 1
                  next2_n = MOD(next_n,srch_corners) + 1

                  vec1_lat = srch_corner_lat(next_n,cell) - &
                      srch_corner_lat(n,cell)
                  vec1_lon = srch_corner_lon(next_n,cell) - &
                      srch_corner_lon(n,cell)
                  vec2_lat = srch_corner_lat(next2_n,cell) - &
                      srch_corner_lat(next_n,cell)
                  vec2_lon =  srch_corner_lon(next2_n,cell) - & 
                      srch_corner_lon(next_n,cell)
                  
                  if (vec1_lon > three*dble_pih) then
                     vec1_lon = vec1_lon - dble_pi2
                  else if (vec1_lon < -three*dble_pih) then
                     vec1_lon = vec1_lon + dble_pi2
                  endif
                  
                  if (vec2_lon > three*dble_pih) then
                     vec2_lon = vec2_lon - dble_pi2
                  else if (vec2_lon < -three*dble_pih) then
                     vec2_lon = vec2_lon + dble_pi2
                  endif
                  
                  cross_product= &
                      vec1_lat*vec2_lon - vec2_lat*vec1_lon

                  if (cross_product < zero) then
                     neg=neg+1
                  else if (cross_product > zero) then
                     pos=pos+1
                  endif
               enddo
               !*** the cell is non-convexe if not all cross_products 
               !*** have the same signe
               if (neg/=0 .and. pos/=0) then
                  loutside=.false.
                  print*, 'The mesh ',srch_add(cell),' is non-convex'
                  print*, 'srch_corner_lat=',srch_corner_lat(:,cell)
                  print*, 'srch_corner_lon=',srch_corner_lon(:,cell)
               endif
            endif
            if (loutside) then
              lat2 = beglat
              lon2 = beglon
              location = 0
              eps  = -tiny
            else
              eps  = tiny
            endif
            exit srch_loop
          endif

          !***
          !*** otherwise move on to next cell
          !***
        end do cell_loop

        !***
        !*** if still no cell found, the point lies outside the grid.
        !***   take some baby steps along the segment to see if any
        !***   part of the segment lies inside the grid.  
        !***

        loutside = .true.
        s1 = s1 + baby
        lat1 = beglat + s1*(endlat - beglat)
        lon1 = beglon + s1*(lon2   - beglon)

        !***
        !*** reached the end of the segment and still outside the grid
        !*** return no intersection
        !***
        if (s1 >= one) return

      end do srch_loop
!      call timer_stop(12)

!-----------------------------------------------------------------------
!
!     now that a cell is found, search for the next intersection.
!     loop over sides of the cell to find intersection with side
!     must check all sides for coincidences or intersections
!
!-----------------------------------------------------------------------

!      call timer_start(13)
      intrsct_loop: do n=1,srch_corners
        next_n = mod(n,srch_corners) + 1

        grdlon1 = srch_corner_lon(n     ,cell)
        grdlon2 = srch_corner_lon(next_n,cell)
        grdlat1 = srch_corner_lat(n     ,cell)
        grdlat2 = srch_corner_lat(next_n,cell)

        !***
        !*** set up linear system to solve for intersection
        !***

        mat1 = lat2 - lat1
        mat2 = grdlat1 - grdlat2
        mat3 = lon2 - lon1
        mat4 = grdlon1 - grdlon2
        rhs1 = grdlat1 - lat1
        rhs2 = grdlon1 - lon1

        if (mat3 >  dble_pi) then
          mat3 = mat3 - dble_pi2
        else if (mat3 < -dble_pi) then
          mat3 = mat3 + dble_pi2
        endif
        if (mat4 >  dble_pi) then
          mat4 = mat4 - dble_pi2
        else if (mat4 < -dble_pi) then
          mat4 = mat4 + dble_pi2
        endif
        if (rhs2 >  dble_pi) then
          rhs2 = rhs2 - dble_pi2
        else if (rhs2 < -dble_pi) then
          rhs2 = rhs2 + dble_pi2
        endif

        determ = mat1*mat4 - mat2*mat3

        !***
        !*** if the determinant is zero, the segments are either 
        !***   parallel or coincident.  coincidences were detected 
        !***   above so do nothing.
        !*** if the determinant is non-zero, solve for the linear 
        !***   parameters s for the intersection point on each line 
        !***   segment.
        !*** if 0<s1,s2<1 then the segment intersects with this side.
        !***   return the point of intersection (adding a small
        !***   number so the intersection is off the grid line).
        !***

        if (abs(determ) > 1.e-30) then

          s1 = (rhs1*mat4 - mat2*rhs2)/determ
          s2 = (mat1*rhs2 - rhs1*mat3)/determ

          if (s2 >= zero .and. s2 <= one .and. &
             s1 >  zero .and. s1 <= one) then

            !***
            !*** recompute intersection based on full segment
            !*** so intersections are consistent for both sweeps
            !***

            if (.not. loutside) then
              mat1 = lat2 - begseg(1)
              mat3 = lon2 - begseg(2)
              rhs1 = grdlat1 - begseg(1)
              rhs2 = grdlon1 - begseg(2)
            else
              mat1 = begseg(1) - endlat
              mat3 = begseg(2) - endlon
              rhs1 = grdlat1 - endlat
              rhs2 = grdlon1 - endlon
            endif

            if (mat3 >  dble_pi) then
              mat3 = mat3 - dble_pi2
            else if (mat3 < -dble_pi) then
              mat3 = mat3 + dble_pi2
            endif
            if (rhs2 >  dble_pi) then
              rhs2 = rhs2 - dble_pi2
            else if (rhs2 < -dble_pi) then
              rhs2 = rhs2 + dble_pi2
            endif

            determ = mat1*mat4 - mat2*mat3

            !***
            !*** sometimes due to roundoff, the previous 
            !*** determinant is non-zero, but the lines
            !*** are actually coincident.  if this is the
            !*** case, skip the rest.
            !***

            if (determ /= zero) then
              s1 = (rhs1*mat4 - mat2*rhs2)/determ
              s2 = (mat1*rhs2 - rhs1*mat3)/determ

              offset = s1 + eps/determ
              if (offset > one) offset = one

              if (.not. loutside) then
                intrsct_lat = begseg(1) + mat1*s1
                intrsct_lon = begseg(2) + mat3*s1
                intrsct_lat_off = begseg(1) + mat1*offset
                intrsct_lon_off = begseg(2) + mat3*offset
              else
                intrsct_lat = endlat + mat1*s1
                intrsct_lon = endlon + mat3*s1
                intrsct_lat_off = endlat + mat1*offset
                intrsct_lon_off = endlon + mat3*offset
              endif
              exit intrsct_loop
            endif

          endif
        endif

        !***
        !*** no intersection this side, move on to next side
        !***

      end do intrsct_loop
!      call timer_stop(13)

!-----------------------------------------------------------------------
!
!     if the segment crosses a pole threshold, reset the intersection
!     to be the threshold latitude.  only check if this was not a
!     threshold segment since sometimes coordinate transform can end
!     up on other side of threshold again.
!
!-----------------------------------------------------------------------

      if (lthresh) then
        if (intrsct_lat < north_thresh .or. intrsct_lat > south_thresh) &
           lthresh = .false.
      else if (lat1 > zero .and. intrsct_lat > north_thresh) then
        intrsct_lat = north_thresh + tiny
        intrsct_lat_off = north_thresh + eps*mat1
        s1 = (intrsct_lat - begseg(1))/mat1
        intrsct_lon     = begseg(2) + s1*mat3
        intrsct_lon_off = begseg(2) + (s1+eps)*mat3
        last_loc = location
        lthresh = .true.
      else if (lat1 < zero .and. intrsct_lat < south_thresh) then
        intrsct_lat = south_thresh - tiny
        intrsct_lat_off = south_thresh + eps*mat1
        s1 = (intrsct_lat - begseg(1))/mat1
        intrsct_lon     = begseg(2) + s1*mat3
        intrsct_lon_off = begseg(2) + (s1+eps)*mat3
        last_loc = location
        lthresh = .true.
      endif

#ifdef DEBUGX
  PRINT *, '| | | | | | | Quit PRISMTrs_intersection'
  call psmile_flushstd
#endif

!-----------------------------------------------------------------------

      end subroutine prismtrs_intersection

!***********************************************************************