prismtrs_pole_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_pole_intersection
!
! !INTERFACE
      subroutine 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)

!-----------------------------------------------------------------------
!
!     this routine 
!
!-----------------------------------------------------------------------

        !
      USE PRISMDrv, dummy_INTERFACE => PRISMTrs_pole_intersection
      IMPLICIT NONE
!
! !PARAMETERS:
!
      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
        !
      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

      logical, intent(in) :: 
              lrevers   ! flag true if segment integrated in reverse
!
! !RETURN VALUE
!
      integer, intent(inout) :: 
              location  ! address in destination array containing this
                        ! segment -- also may contain last location on
                        ! entry

      logical, intent(out) :: 
             lcoinc    ! flag segment coincident with grid line

      logical, intent(inout) :: 
             lthresh   ! flag segment crossing threshold boundary

      DOUBLE PRECISION, intent(out) :: 
          intrsct_lat, intrsct_lon ! lat/lon coords of next intersect.
!
! !DESCRIPTION
!
! SUBROUTINE "PRISMTrs_pole_intersection" is identical to the intersection 
! routine except that a coordinate transformation (using a Lambert azimuthal
! equivalent projection) is performed to treat polar cells more accurately.
! 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_pole_intersection.F90 2082 2009-10-21 13:31:19Z hanke $
! $Author: hanke $
!----------------------------------------------------------------------
!
!----------------------------------------------------------------------
! Local declarations
!-----------------------------------------------------------------------
  CHARACTER(LEN=len_cvs_string), SAVE  :: mycvs = 
     '$Id: prismtrs_pole_intersection.F90 2082 2009-10-21 13:31:19Z hanke $'
      integer  :: n, next_n, cell, srch_corners

      logical  :: loutside ! flags points outside grid

      DOUBLE PRECISION :: pi4, rns,  ! north/south conversion
          x1, x2,       ! local x variables for segment
          y1, y2,       ! local y variables for segment
          begx, begy,   ! beginning x,y variables for segment
          endx, endy,   ! beginning x,y variables for segment
          begsegx, begsegy,   ! beginning x,y variables for segment
          grdx1, grdx2,  ! local x variables for grid cell
          grdy1, grdy2,  ! local y variables for grid cell
          vec1_y, vec1_x,  ! vectors and cross products used
          vec2_y, vec2_x,  ! during grid search
          cross_product, eps,  ! eps=small offset away from intersect
          s1, s2, determ,     ! variables used for linear solve to
          mat1, mat2, mat3, mat4, rhs1, rhs2  ! find intersection

      DOUBLE PRECISION, dimension(:,:), allocatable :: 
           srch_corner_x,  ! x of each corner of srch cells
           srch_corner_y   ! y of each corner of srch cells

      !***
      !*** save last intersection to avoid roundoff during coord
      !*** transformation
      !***

      logical, save :: luse_last = .false.

      DOUBLE PRECISION, save :: 
           intrsct_x, intrsct_y  ! x,y for intersection

      !***
      !*** variables necessary if segment manages to hit pole
      !***

      integer, save :: 
           avoid_pole_count = 0  ! count attempts to avoid pole

      DOUBLE PRECISION, save :: 
           avoid_pole_offset = tiny  ! endpoint offset to avoid pole

#ifdef VERBOSE
  PRINT *, '| | | | | | | Enter PRISMTrs_pole_intersection'
  call psmile_flushstd
#endif
!-----------------------------------------------------------------------
!
!     initialize defaults, flags, etc.
!
!-----------------------------------------------------------------------

      if (.not. lthresh) location = 0
      lcoinc = .false.
      intrsct_lat = endlat
      intrsct_lon = endlon

      loutside = .false.
      s1 = zero

!-----------------------------------------------------------------------
!
!     convert coordinates
!
!-----------------------------------------------------------------------

      allocate(srch_corner_x(size(srch_corner_lat,DIM=1), &
                             size(srch_corner_lat,DIM=2)), &
               srch_corner_y(size(srch_corner_lat,DIM=1), &
                            size(srch_corner_lat,DIM=2)))

      if (beglat > zero) then
        pi4 = quart*dble_pi
        rns = one
      else
        pi4 = -quart*dble_pi
        rns = -one
      endif

      if (luse_last) then
        x1 = intrsct_x
        y1 = intrsct_y
      else
        x1 = rns*two*sin(pi4 - half*beglat)*cos(beglon)
        y1 =     two*sin(pi4 - half*beglat)*sin(beglon)
        luse_last = .true.
      endif
      x2 = rns*two*sin(pi4 - half*endlat)*cos(endlon)
      y2 =     two*sin(pi4 - half*endlat)*sin(endlon)
      srch_corner_x = rns*two*sin(pi4 - half*srch_corner_lat)* &
                              cos(srch_corner_lon)
      srch_corner_y =     two*sin(pi4 - half*srch_corner_lat)* &
                              sin(srch_corner_lon)

      begx = x1
      begy = y1
      endx = x2
      endy = y2
      begsegx = rns*two*sin(pi4 - half*begseg(1))*cos(begseg(2))
      begsegy =     two*sin(pi4 - half*begseg(1))*sin(begseg(2))
      intrsct_x = endx
      intrsct_y = endy

!-----------------------------------------------------------------------
!
!     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) == location) then
              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_x = srch_corner_x(next_n,cell) -  &
                    srch_corner_x(n     ,cell)
            vec1_y = srch_corner_y(next_n,cell) -  &
                    srch_corner_y(n     ,cell)
            vec2_x = x1 - srch_corner_x(n,cell)
            vec2_y = y1 - srch_corner_y(n,cell)

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

            if (vec2_x == 0 .and. vec2_y == 0) then
              x1 = x1 + 1.d-10*(x2-x1)
              y1 = y1 + 1.d-10*(y2-y1)
              vec2_x = x1 - srch_corner_x(n,cell)
              vec2_y = y1 - srch_corner_y(n,cell)
            endif

            cross_product = vec1_x*vec2_y - vec2_x*vec1_y

            !***
            !*** if the cross product for a side is zero, the point 
            !***   lies exactly on the side or the length of a side
            !***   is zero.  if the length is zero set det > 0.
            !***   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_x /= zero .or. vec1_y /= 0) then
                vec2_x = x2 - x1
                vec2_y = y2 - y1
                cross_product = vec1_x*vec2_y - vec2_x*vec1_y
              else
                cross_product = one
              endif

              if (cross_product == zero) then
                lcoinc = .true.
                cross_product = vec1_x*vec2_x + vec1_y*vec2_y
                if (lrevers) cross_product = -cross_product
              endif
            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
              x2 = begx
              y2 = begy
              location = 0
              eps  = -tiny
            else
              eps  = tiny
            endif

            exit srch_loop
          endif

          !***
          !*** otherwise move on to next cell
          !***

        end do cell_loop

        !***
        !*** if 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
        x1 = begx + s1*(x2 - begx)
        y1 = begy + s1*(y2 - begy)

        !***
        !*** reached the end of the segment and still outside the grid
        !*** return no intersection
        !***

        if (s1 >= one) then
          deallocate(srch_corner_x, srch_corner_y)
          luse_last = .false.
          return
        endif

      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

        grdy1 = srch_corner_y(n     ,cell)
        grdy2 = srch_corner_y(next_n,cell)
        grdx1 = srch_corner_x(n     ,cell)
        grdx2 = srch_corner_x(next_n,cell)

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

        mat1 = x2 - x1
        mat2 = grdx1 - grdx2
        mat3 = y2 - y1
        mat4 = grdy1 - grdy2
        rhs1 = grdx1 - x1
        rhs2 = grdy1 - y1

        determ = mat1*mat4 - mat2*mat3

        !***
        !*** if the determinant is zero, the segments are either 
        !***   parallel or coincident or one segment has zero length.  
        !***   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 using entire segment
            !*** for consistency between sweeps
            !***

            if (.not. loutside) then
              mat1 = x2 - begsegx
              mat3 = y2 - begsegy
              rhs1 = grdx1 - begsegx
              rhs2 = grdy1 - begsegy
            else 
              mat1 = x2 - endx
              mat3 = y2 - endy
              rhs1 = grdx1 - endx
              rhs2 = grdy1 - endy
            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

              if (.not. loutside) then
                intrsct_x = begsegx + s1*mat1
                intrsct_y = begsegy + s1*mat3
              else 
                intrsct_x = endx + s1*mat1
                intrsct_y = endy + s1*mat3
              endif

              !***
              !*** convert back to lat/lon coordinates
              !***

              intrsct_lon = rns*atan2(intrsct_y,intrsct_x)
              if (intrsct_lon < zero)  &
                intrsct_lon = intrsct_lon + dble_pi2

              if (abs(intrsct_x) > 1.d-10) then
                intrsct_lat = (pi4 -  &
                  asin(rns*half*intrsct_x/cos(intrsct_lon)))*two
              else if (abs(intrsct_y) > 1.d-10) then
                intrsct_lat = (pi4 -  &
                  asin(half*intrsct_y/sin(intrsct_lon)))*two
              else
                intrsct_lat = two*pi4
              endif

              !***
              !*** add offset in transformed space for next pass.
              !***

              if (s1 - eps/determ < one) then
                intrsct_x = intrsct_x - mat1*(eps/determ)
                intrsct_y = intrsct_y - mat3*(eps/determ)
              else
                if (.not. loutside) then
                  intrsct_x = endx
                  intrsct_y = endy
                  intrsct_lat = endlat
                  intrsct_lon = endlon
                else 
                  intrsct_x = begsegx
                  intrsct_y = begsegy
                  intrsct_lat = begseg(1)
                  intrsct_lon = begseg(2)
                endif
              endif

              exit intrsct_loop
            endif
          endif
        endif

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

      end do intrsct_loop
      ! call timer_stop(13)

      deallocate(srch_corner_x, srch_corner_y)

!-----------------------------------------------------------------------
!
!     if segment manages to cross over pole, shift the beginning 
!     endpoint in order to avoid hitting pole directly
!     (it is ok for endpoint to be pole point)
!
!-----------------------------------------------------------------------

      if (abs(intrsct_x) < 1.e-10 .and. abs(intrsct_y) < 1.e-10 .and. &
         (endx /= zero .and. endy /=0)) then
        if (avoid_pole_count > 2) then
           avoid_pole_count = 0
           avoid_pole_offset = 10.*avoid_pole_offset
        endif

        cross_product = begsegx*(endy-begsegy) - begsegy*(endx-begsegx)
        intrsct_lat = begseg(1)
        if (cross_product*intrsct_lat > zero) then
          intrsct_lon = beglon    + avoid_pole_offset
          begseg(2)   = begseg(2) + avoid_pole_offset
        else
          intrsct_lon = beglon    - avoid_pole_offset
          begseg(2)   = begseg(2) - avoid_pole_offset
        endif

        avoid_pole_count = avoid_pole_count + 1
        luse_last = .false.
      else
        avoid_pole_count = 0
        avoid_pole_offset = tiny
      endif

!-----------------------------------------------------------------------
!
!     if the segment crosses a pole threshold, reset the intersection
!     to be the threshold latitude and do not reuse x,y intersect
!     on next entry.  only check if did not cross threshold last
!     time - sometimes the coordinate transformation can place a
!     segment on the other side of the threshold again
!
!-----------------------------------------------------------------------

      if (lthresh) then
        if (intrsct_lat > north_thresh .or. intrsct_lat < south_thresh) &
         lthresh = .false.
      else if (beglat > zero .and. intrsct_lat < north_thresh) then
        mat4 = endlat - begseg(1)
        mat3 = endlon - begseg(2)
        if (mat3 >  dble_pi) mat3 = mat3 - dble_pi2
        if (mat3 < -dble_pi) mat3 = mat3 + dble_pi2
        intrsct_lat = north_thresh - tiny
        s1 = (north_thresh - begseg(1))/mat4
        intrsct_lon = begseg(2) + s1*mat3
        luse_last = .false.
        lthresh = .true.
      else if (beglat < zero .and. intrsct_lat > south_thresh) then
        mat4 = endlat - begseg(1)
        mat3 = endlon - begseg(2)
        if (mat3 >  dble_pi) mat3 = mat3 - dble_pi2
        if (mat3 < -dble_pi) mat3 = mat3 + dble_pi2
        intrsct_lat = south_thresh + tiny
        s1 = (south_thresh - begseg(1))/mat4
        intrsct_lon = begseg(2) + s1*mat3
        luse_last = .false.
        lthresh = .true.
      endif

      !***
      !*** if reached end of segment, do not use x,y intersect 
      !*** on next entry
      !***

      if (intrsct_lat == endlat .and. intrsct_lon == endlon) then
        luse_last = .false.
      endif
#ifdef VERBOSE
  PRINT *, '| | | | | | | Quit PRISMTrs_pole_intersection'
  call psmile_flushstd
#endif
!-----------------------------------------------------------------------

      end subroutine PRISMTrs_pole_intersection