#include "swancpp.h"
MODULE interp_swan_mod
#if defined WRF_COUPLING || defined NESTING
!
!=======================================================================
! This module contains some original functions as well as !
! some routines modified from a ROMS subroutine interpolate.F and !
! contains several all purpuse generic routines: !
! !
! Routines: !
! !
! compute_angle computes swan grid angle.
! swan_ref_init Initialize swan refinement. !
! shindices Finds model grid cell for any datum. !
! stry_range Binary search of model grid cell for any datum. !
! sinside Closed polygon datum search. !
! !
!=======================================================================
!
USE mod_coupler_kinds
implicit none
CONTAINS
# ifdef WRF_COUPLING
SUBROUTINE compute_angler (ng)
!
!=======================================================================
! !
! This routine computes the anlge of the swan grid. !
! !
!=======================================================================
!
USE M_GENARR
USE M_PARALL
USE SWCOMM2
USE SWCOMM4
!
! Imported.
!
integer, intent(in) :: ng
!
! Locals.
!
integer :: i, j, count
real, allocatable :: xlon(:,:), ylat(:,:)
real, allocatable :: xlonb(:,:), ylatb(:,:)
real :: cff1, cff2, cff3, angler, DEG2RAD
!
! Compute SWAN grid angles for rotation of winds from WRF.
!
ALLOCATE(M_GENARR_MOD(ng)%CosAngler_G(MXCGL*MYCGL))
CosAngler=>M_GENARR_MOD(ng)%CosAngler_G
ALLOCATE(M_GENARR_MOD(ng)%SinAngler_G(MXCGL*MYCGL))
SinAngler=>M_GENARR_MOD(ng)%SinAngler_G
ALLOCATE(xlon(MXCGL,MYCGL))
ALLOCATE(ylat(MXCGL,MYCGL))
ALLOCATE(xlonb(MXCGL+1,MYCGL))
ALLOCATE(ylatb(MXCGL+1,MYCGL))
!
! xlon and xlat are the locations at the center points.
! xlonb and ylatb are offset by 1 and used to compute the angle
! at the center points.
!
DO j=1,MYCGL
DO i=1,MXCGL
cff1=XOFFS_G(ng)
cff2=PARALL_MOD(ng)%XGRDGL_G(i,j)
xlon(i,j)=cff1+cff2
cff1=YOFFS_G(ng)
cff2=PARALL_MOD(ng)%YGRDGL_G(i,j)
ylat(i,j)=cff1+cff2
END DO
END DO
!
DO i=2,MXCGL
DO j=1,MYCGL
xlonb(i,j)=0.5*(xlon(i-1,j)+xlon(i,j))
ylatb(i,j)=0.5*(ylat(i-1,j)+ylat(i,j))
END DO
END DO
DO j=1,MYCGL
xlonb(1,j)=xlonb(2,j)
xlonb(MXCGL+1,j)=xlonb(MXCGL,j)
ylatb(1,j)=ylatb(2,j)
ylatb(MXCGL+1,j)=ylatb(MXCGL,j)
END DO
!
! Now compute the angles.
!
IF (KSPHER.eq.0) THEN ! Cartesian
DEG2RAD=1.
ELSE
DEG2RAD=3.14159/180.
END IF
count=0
DO j=1,MYCGL
DO i=1,MXCGL
count=count+1
cff1=(xlonb(i+1,j)-xlonb(i,j))*DEG2RAD
cff2=(ylatb(i+1,j)-ylatb(i,j))*DEG2RAD
cff3=COS(ylat(i,j)*DEG2RAD)
IF (KSPHER.eq.0) THEN ! Cartesian
angler=ATAN2(cff2,cff1)
ELSE
angler=ATAN2(cff2,cff1*cff3)
END IF
CosAngler(count)=COS(angler)
SinAngler(count)=SIN(angler)
END DO
END DO
!
DEALLOCATE(xlon, ylat, xlonb, ylatb)
RETURN
END SUBROUTINE compute_angler
# endif
# ifdef NESTING
SUBROUTINE swan_ref_init (ng, ngp)
!
!=======================================================================
! !
! This routine determines the bounding cells for each child grid, !
! determines the indices of the child cell in the parent grid, !
! and associates a parent tile for each bounding point. !
! !
!=======================================================================
!
USE M_PARALL
USE M_SREFINED
USE SWCOMM2
USE SWCOMM3
USE M_GENARR
integer, intent(in) :: ng, ngp
! Locals.
integer :: i, j, count, IX, IY, MyError
integer :: Xmin, Xmax, Ymin, Ymax
INTEGER IARRL(5), IARRC(5,1:NPROC)
real, allocatable :: Ipos(:,:), Jpos(:,:), angler(:,:)
real :: IJspv
logical :: rectangular
!
! Create arrays of boundary points for current ng grid because
! it is a child grid.
!
ALLOCATE (REFINED_MOD(ng)%BOUNDHINDX_G(Numspec(ng)))
BOUNDHINDX =>REFINED_MOD(ng)%BOUNDHINDX_G
BOUNDHINDX=0.
ALLOCATE (REFINED_MOD(ng)%BOUNDHINDY_G(Numspec(ng)))
BOUNDHINDY =>REFINED_MOD(ng)%BOUNDHINDY_G
BOUNDHINDY=0.
ALLOCATE (REFINED_MOD(ng)%BOUNDHINDPF_G(Numspec(ng)))
BOUNDHINDPF =>REFINED_MOD(ng)%BOUNDHINDPF_G
ALLOCATE (REFINED_MOD(ng)%BOUNDHINDPR_G(Numspec(ng)))
BOUNDHINDPR =>REFINED_MOD(ng)%BOUNDHINDPR_G
IF (NPROC.gt.1) THEN
BOUNDHINDPF=0
BOUNDHINDPR=0
ELSE
BOUNDHINDPF=1
BOUNDHINDPR=1
END IF
ALLOCATE (REFINED_MOD(ng)%BOUNDHINDIX_G(Numspec(ng)))
BOUNDHINDIX =>REFINED_MOD(ng)%BOUNDHINDIX_G
BOUNDHINDIX=0
ALLOCATE (REFINED_MOD(ng)%BOUNDHINDIY_G(Numspec(ng)))
BOUNDHINDIY =>REFINED_MOD(ng)%BOUNDHINDIY_G
BOUNDHINDIY=0
!
! Compute number of boundary array points to be passed.
!
ac2size(ng)=MSC*MDC*Numspec(ng)
!
! Determine the horizontal indices of each child
! point in the parent grid.
!
IJspv=9999.
rectangular=.false.
allocate(Ipos(MXCGL,MYCGL),Jpos(MXCGL,MYCGL))
allocate(angler(MXCGL_G(ngp),MYCGL_G(ngp)))
angler=0.
!
CALL shindices(ng, 1, MXCGL_G(ngp), 1, MYCGL_G(ngp), &
& 1, MXCGL_G(ngp), 1, MYCGL_G(ngp), &
& angler, &
& PARALL_MOD(ngp)%XGRDGL_G+XOFFS_G(ngp), &
& PARALL_MOD(ngp)%YGRDGL_G+YOFFS_G(ngp), &
& 1, MXCGL, 1, MYCGL, &
& 1, MXCGL, 1, MYCGL, &
& XGRDGL+XOFFS_G(ng), YGRDGL+YOFFS_G(ng), Ipos, Jpos,&
& IJspv, rectangular, Xmin, Xmax, Ymin, Ymax)
!
! Fill the boundary arrays with the hindices.
!
count=0
DO i=1,MXCGL
count=count+1
IX=i
IY=1
BOUNDHINDX(count)=Ipos(IX,IY)
BOUNDHINDY(count)=Jpos(IX,IY)
BOUNDHINDIX(count)=IX
BOUNDHINDIY(count)=IY
END DO
DO i=1,MXCGL
count=count+1
IX=i
IY=MYCGL
BOUNDHINDX(count)=Ipos(IX,IY)
BOUNDHINDY(count)=Jpos(IX,IY)
BOUNDHINDIX(count)=IX
BOUNDHINDIY(count)=IY
END DO
DO i=2,MYCGL-1
count=count+1
IX=1
IY=i
BOUNDHINDX(count)=Ipos(IX,IY)
BOUNDHINDY(count)=Jpos(IX,IY)
BOUNDHINDIX(count)=IX
BOUNDHINDIY(count)=IY
END DO
DO i=2,MYCGL-1
count=count+1
IX=MXCGL
IY=i
BOUNDHINDX(count)=Ipos(IX,IY)
BOUNDHINDY(count)=Jpos(IX,IY)
BOUNDHINDIX(count)=IX
BOUNDHINDIY(count)=IY
END DO
!
! Determine which parent tile will give AC2 for each child point.
!
IARRL(1) = MXF_G(ngp)
IARRL(2) = MXL_G(ngp)
IARRL(3) = MYF_G(ngp)
IARRL(4) = MYL_G(ngp)
IARRL(5) = MCGRD_G(ngp)
CALL SWGATHER (IARRC, 5*NPROC, IARRL, 5, SWINT )
!
CALL MPI_BCAST(IARRC,5*NPROC,SWINT,0, &
& WAV_COMM_WORLD,MyError)
!
! BOUNDHINDPF = which processor is going to fill the AC2 for this point.
! BOUNDHINDPR = which processor is able to read this data to fill its AC2 array.
!
IF (NPROC.gt.1) THEN
DO i=1,Numspec(ng)
IF (MXCGL_G(ngp).gt.MYCGL_G(ngp)) THEN
DO j=1,NPROC
IF (BOUNDHINDPF(i).eq.0) THEN
IF ((BOUNDHINDX(i).ge.IARRC(1,j)).and. &
& (BOUNDHINDX(i).le.IARRC(2,j)-3)) THEN
BOUNDHINDPF(i)=j
ENDIF
ENDIF
END DO
ELSE
DO j=1,NPROC
IF (BOUNDHINDPF(i).eq.0) THEN
IF ((BOUNDHINDY(i).ge.IARRC(3,j)).and. &
& (BOUNDHINDY(i).le.IARRC(4,j)-3)) THEN
BOUNDHINDPF(i)=j
ENDIF
ENDIF
END DO
ENDIF
IF (MXCGL.gt.MYCGL) THEN
IF ((BOUNDHINDIX(i).ge.MXF).and. &
& (BOUNDHINDIX(i).le.MXL)) THEN
BOUNDHINDPR(i)=INODE
ENDIF
ELSE
IF ((BOUNDHINDIY(i).ge.MYF).and. &
& (BOUNDHINDIY(i).le.MYL)) THEN
BOUNDHINDPR(i)=INODE
ENDIF
ENDIF
END DO
END IF
deallocate(Ipos,Jpos)
RETURN
END SUBROUTINE swan_ref_init
SUBROUTINE shindices (ng, LBi, UBi, LBj, UBj, &
& Is, Ie, Js, Je, &
& angler, Xgrd, Ygrd, &
& LBm, UBm, LBn, UBn, &
& Ms, Me, Ns, Ne, &
& Xpos, Ypos, Ipos, Jpos, &
& IJspv, rectangular, Xmin, Xmax, Ymin, Ymax)
!
!=======================================================================
! !
! Given any geographical locations Xpos and Ypos, this routine finds !
! the corresponding array cell indices (Ipos, Jpos) of gridded data !
! Xgrd and Ygrd containing each requested location. This indices are !
! usually used for interpolation. !
! !
! On Input: !
! !
! ng Nested grid number. !
! LBi I-dimension Lower bound of gridded data. !
! UBi I-dimension Upper bound of gridded data. !
! LBj J-dimension Lower bound of gridded data. !
! UBj J-dimension Upper bound of gridded data. !
! Is Starting gridded data I-index to search. !
! Ie Ending gridded data I-index to search. !
! Js Starting gridded data J-index to search. !
! Je Ending gridded data J-index to search. !
! angler Gridded data angle between X-axis and true EAST !
! (radians). !
! Xgrd Gridded data X-locations (usually, longitude). !
! Ygrd Gridded data Y-locations (usually, latitude). !
! LBm I-dimension Lower bound of requested locations. !
! UBm I-dimension Upper bound of requested locations. !
! LBn J-dimension Lower bound of requested locations. !
! UBn J-dimension Upper bound of requested locations. !
! Ms Starting requested locations I-index to search. !
! Me Ending requested locations I-index to search. !
! Ns Starting requested locations J-index to search. !
! Ne Ending requested locations J-index to search. !
! Xpos Requested X-locations to process (usually longitude).!
! Ypos Requested Y-locations to process (usually latitude). !
! IJspv Unbounded special value to assign. !
! rectangular Logical switch indicating that gridded data has a !
! plaid distribution. !
! !
! On Output: !
! !
! Ipos Fractional I-cell index containing locations in data. !
! Jpos Fractional J-cell index containing locations in data. !
! Xmin minimal cell index in x-dir containing data !
! Xmax maximum cell index in x-dir containing data !
! Ymin minimal cell index in y-dir containing data !
! Ymax maximum cell index in y-dir containing data !
! !
! Calls: Try_Range !
! !
!=======================================================================
!
USE mod_coupler_kinds
USE SWCOMM3
USE SWCOMM4
!
! Imported variable declarations.
!
logical, intent(in) :: rectangular
integer, intent(in) :: ng
integer, intent(in) :: LBi, UBi, LBj, UBj, Is, Ie, Js, Je
integer, intent(in) :: LBm, UBm, LBn, UBn, Ms, Me, Ns, Ne
real, intent(in) :: IJspv
!
real, intent(in) :: angler(LBi:UBi,LBj:UBj)
real, intent(in) :: Xgrd(LBi:UBi,LBj:UBj)
real, intent(in) :: Ygrd(LBi:UBi,LBj:UBj)
real, intent(in) :: Xpos(LBm:UBm,LBn:UBn)
real, intent(in) :: Ypos(LBm:UBm,LBn:UBn)
real, intent(out) :: Ipos(LBm:UBm,LBn:UBn)
real, intent(out) :: Jpos(LBm:UBm,LBn:UBn)
integer, intent(out), optional :: Xmin, Xmax, Ymin, Ymax
!
! Local variable declarations.
!
logical :: found, foundi, foundj, foundxy
integer :: Imax, Imin, Jmax, Jmin, i, i0, j, j0, mp, np
real :: aa2, ang, bb2, diag2, dxx, dyy, phi
real :: xfac, xpp, yfac, ypp
!
!-----------------------------------------------------------------------
! Determine grid cell indices containing requested position points.
! Then, interpolate to fractional cell position.
!-----------------------------------------------------------------------
!
Xmin=-9999
Xmax=-9999
Ymin=-9999
Ymax=-9999
foundxy=.TRUE.
DO np=Ns,Ne
DO mp=Ms,Me
Ipos(mp,np)=IJspv
Jpos(mp,np)=IJspv
!
! The gridded data has a plaid distribution so the search is trivial.
!
IF (rectangular) THEN
foundi=.FALSE.
I_LOOP : DO i=LBi,UBi-1
IF ((Xgrd(i ,1).le.Xpos(mp,np)).and. &
& (Xgrd(i+1,1).gt.Xpos(mp,np))) THEN
Imin=i
foundi=.TRUE.
EXIT I_LOOP
END IF
END DO I_LOOP
foundj=.FALSE.
J_LOOP : DO j=LBj,UBj-1
IF ((Ygrd(1,j ).le.Ypos(mp,np)).and. &
& (Ygrd(1,j+1).gt.Ypos(mp,np))) THEN
Jmin=j
foundj=.TRUE.
EXIT J_LOOP
END IF
END DO J_LOOP
found=foundi.and.foundj
!
! Check each position to find if it falls inside the whole domain.
! Once it is stablished that it inside, find the exact cell to which
! it belongs by successively dividing the domain by a half (binary
! search).
!
ELSE
found=stry_range(ng, LBi, UBi, LBj, UBj, &
& Xgrd, Ygrd, &
& Is, Ie, Js, Je, &
& Xpos(mp,np), Ypos(mp,np))
IF (found) THEN
Imin=Is
Imax=Ie
Jmin=Js
Jmax=Je
DO while (((Imax-Imin).gt.1).or.((Jmax-Jmin).gt.1))
IF ((Imax-Imin).gt.1) THEN
i0=(Imin+Imax)/2
found=stry_range(ng, LBi, UBi, LBj, UBj, &
& Xgrd, Ygrd, &
& Imin, i0, Jmin, Jmax, &
& Xpos(mp,np), Ypos(mp,np))
IF (found) THEN
Imax=i0
ELSE
Imin=i0
END IF
END IF
IF ((Jmax-Jmin).gt.1) THEN
j0=(Jmin+Jmax)/2
found=stry_range(ng, LBi, UBi, LBj, UBj, &
& Xgrd, Ygrd, &
& Imin, Imax, Jmin, j0, &
& Xpos(mp,np), Ypos(mp,np))
IF (found) THEN
Jmax=j0
ELSE
Jmin=j0
END IF
END IF
END DO
found=(Is.le.Imin).and.(Imin.le.Ie).and. &
& (Is.le.Imax).and.(Imax.le.Ie).and. &
& (Js.le.Jmin).and.(Jmin.le.Je).and. &
& (Js.le.Jmax).and.(Jmax.le.Je)
END IF
END IF
!
! Knowing the correct cell, calculate the exact indices, accounting
! for a possibly rotated grid. If spherical, convert all positions
! to meters first.
!
IF (found) THEN
IF (KSPHER.gt.0) THEN
yfac=REARTH*DEGRAD
xfac=yfac*COS(Ypos(mp,np)*DEGRAD)
xpp=(Xpos(mp,np)-Xgrd(Imin,Jmin))*xfac
ypp=(Ypos(mp,np)-Ygrd(Imin,Jmin))*yfac
ELSE
xfac=1.0_m8
yfac=1.0_m8
xpp=Xpos(mp,np)-Xgrd(Imin,Jmin)
ypp=Ypos(mp,np)-Ygrd(Imin,Jmin)
END IF
!
! Use Law of Cosines to get cell parallelogram "shear" angle.
!
diag2=((Xgrd(Imin+1,Jmin)-Xgrd(Imin,Jmin+1))*xfac)**2+ &
& ((Ygrd(Imin+1,Jmin)-Ygrd(Imin,Jmin+1))*yfac)**2
aa2=((Xgrd(Imin,Jmin)-Xgrd(Imin+1,Jmin))*xfac)**2+ &
& ((Ygrd(Imin,Jmin)-Ygrd(Imin+1,Jmin))*yfac)**2
bb2=((Xgrd(Imin,Jmin)-Xgrd(Imin,Jmin+1))*xfac)**2+ &
& ((Ygrd(Imin,Jmin)-Ygrd(Imin,Jmin+1))*yfac)**2
phi=ASIN((diag2-aa2-bb2)/(2.0_m8*SQRT(aa2*bb2)))
!
! Transform float position into curvilinear coordinates. Assume the
! cell is rectanglar, for now.
!
ang=angler(Imin,Jmin)
dxx=xpp*COS(ang)+ypp*SIN(ang)
dyy=ypp*COS(ang)-xpp*SIN(ang)
!
! Correct for parallelogram.
!
dxx=dxx+dyy*TAN(phi)
dyy=dyy/COS(phi)
!
! Scale with cell side lengths to translate into cell indices.
!
dxx=MIN(MAX(0.0,dxx/SQRT(aa2)),1.0)
dyy=MIN(MAX(0.0,dyy/SQRT(bb2)),1.0)
Ipos(mp,np)=REAL(Imin)+dxx
Jpos(mp,np)=REAL(Jmin)+dyy
!
! Set min and max indices
!
IF (foundxy) THEN
Xmin=mp
Ymin=np
foundxy=.FALSE.
ENDIF
Xmax=MAX(Xmax, mp)
Ymax=MAX(Ymax, np)
END IF
END DO
END DO
RETURN
END SUBROUTINE shindices
LOGICAL FUNCTION stry_range (ng, LBi, UBi, LBj, UBj, Xgrd, Ygrd, &
& Imin, Imax, Jmin, Jmax, Xo, Yo)
!
!=======================================================================
! !
! Given a grided domain with matrix coordinates Xgrd and Ygrd, this !
! function finds if the point (Xo,Yo) is inside the box defined by !
! the requested corners (Imin,Jmin) and (Imax,Jmax). It will return !
! logical switch try_range=.TRUE. if (Xo,Yo) is inside, otherwise !
! it will return false. !
! !
! Calls: inside !
! !
!=======================================================================
!
!
! Imported variable declarations.
!
integer, intent(in) :: ng, LBi, UBi, LBj, UBj
integer, intent(in) :: Imin, Imax, Jmin, Jmax
# ifdef ASSUMED_SHAPE
real, intent(in) :: Xgrd(LBi:,LBj:)
real, intent(in) :: Ygrd(LBi:,LBj:)
# else
real, intent(in) :: Xgrd(LBi:UBi,LBj:UBj)
real, intent(in) :: Ygrd(LBi:UBi,LBj:UBj)
# endif
real, intent(in) :: Xo, Yo
!
! Local variable declarations.
!
integer :: Nb, i, j, shft, ic
real, dimension(2*(Jmax-Jmin+Imax-Imin)+1) :: Xb, Yb
!
!-----------------------------------------------------------------------
! Define closed polygon.
!-----------------------------------------------------------------------
!
! Note that the last point (Xb(Nb),Yb(Nb)) does not repeat first
! point (Xb(1),Yb(1)). Instead, in function inside, it is implied
! that the closing segment is (Xb(Nb),Yb(Nb))-->(Xb(1),Yb(1)). In
! fact, function inside sets Xb(Nb+1)=Xb(1) and Yb(Nb+1)=Yb(1).
!
Nb=2*(Jmax-Jmin+Imax-Imin)
shft=1-Imin
DO i=Imin,Imax-1
Xb(i+shft)=Xgrd(i,Jmin)
Yb(i+shft)=Ygrd(i,Jmin)
END DO
shft=1-Jmin+Imax-Imin
DO j=Jmin,Jmax-1
Xb(j+shft)=Xgrd(Imax,j)
Yb(j+shft)=Ygrd(Imax,j)
END DO
shft=1+Jmax-Jmin+2*Imax-Imin
DO i=Imax,Imin+1,-1
Xb(shft-i)=Xgrd(i,Jmax)
Yb(shft-i)=Ygrd(i,Jmax)
END DO
shft=1+2*Jmax-Jmin+2*(Imax-Imin)
DO j=Jmax,Jmin+1,-1
Xb(shft-j)=Xgrd(Imin,j)
Yb(shft-j)=Ygrd(Imin,j)
END DO
!
!-----------------------------------------------------------------------
! Check if point (Xo,Yo) is inside of the defined polygon.
!-----------------------------------------------------------------------
!
stry_range=sinside(Nb, Xb, Yb, Xo, Yo)
RETURN
END FUNCTION stry_range
LOGICAL FUNCTION sinside (Nb, Xb, Yb, Xo, Yo)
!
!=======================================================================
! !
! Given the vectors Xb and Yb of size Nb, defining the coordinates !
! of a closed polygon, this function find if the point (Xo,Yo) is !
! inside the polygon. If the point (Xo,Yo) falls exactly on the !
! boundary of the polygon, it still considered inside. !
! !
! This algorithm does not rely on the setting of Xb(Nb)=Xb(1) and !
! Yb(Nb)=Yb(1). Instead, it assumes that the last closing segment !
! is (Xb(Nb),Yb(Nb)) --> (Xb(1),Yb(1)). !
! !
! Reference: !
! !
! Reid, C., 1969: A long way from Euclid. Oceanography EMR, !
! page 174. !
! !
! Algorithm: !
! !
! The decision whether the point is inside or outside the polygon !
! is done by counting the number of crossings from the ray (Xo,Yo) !
! to (Xo,-infinity), hereafter called meridian, by the boundary of !
! the polygon. In this counting procedure, a crossing is counted !
! as +2 if the crossing happens from "left to right" or -2 if from !
! "right to left". If the counting adds up to zero, then the point !
! is outside. Otherwise, it is either inside or on the boundary. !
! !
! This routine is a modified version of the Reid (1969) algorithm, !
! where all crossings were counted as positive and the decision is !
! made based on whether the number of crossings is even or odd. !
! This new algorithm may produce different results in cases where !
! Xo accidentally coinsides with one of the (Xb(k),k=1:Nb) points. !
! In this case, the crossing is counted here as +1 or -1 depending !
! of the sign of (Xb(k+1)-Xb(k)). Crossings are not counted if !
! Xo=Xb(k)=Xb(k+1). Therefore, if Xo=Xb(k0) and Yo>Yb(k0), and if !
! Xb(k0-1) < Xb(k0) < Xb(k0+1), the crossing is counted twice but !
! with weight +1 (for segments with k=k0-1 and k=k0). Similarly if !
! Xb(k0-1) > Xb(k0) > Xb(k0+1), the crossing is counted twice with !
! weight -1 each time. If, on the other hand, the meridian only !
! touches the boundary, that is, for example, Xb(k0-1) < Xb(k0)=Xo !
! and Xb(k0+1) < Xb(k0)=Xo, then the crossing is counted as +1 for !
! segment k=k0-1 and -1 for segment k=k0, resulting in no crossing. !
! !
! Note 1: (Explanation of the logical condition) !
! !
! Suppose that there exist two points (x1,y1)=(Xb(k),Yb(k)) and !
! (x2,y2)=(Xb(k+1),Yb(k+1)), such that, either (x1 < Xo < x2) or !
! (x1 > Xo > x2). Therefore, meridian x=Xo intersects the segment !
! (x1,y1) -> (x2,x2) and the ordinate of the point of intersection !
! is: !
! !
! y1*(x2-Xo) + y2*(Xo-x1) !
! y = ----------------------- !
! x2-x1 !
! !
! The mathematical statement that point (Xo,Yo) either coinsides !
! with the point of intersection or lies to the north (Yo>=y) from !
! it is, therefore, equivalent to the statement: !
! !
! Yo*(x2-x1) >= y1*(x2-Xo) + y2*(Xo-x1), if x2-x1 > 0 !
! or !
! Yo*(x2-x1) <= y1*(x2-Xo) + y2*(Xo-x1), if x2-x1 < 0 !
! !
! which, after noting that Yo*(x2-x1) = Yo*(x2-Xo + Xo-x1) may be !
! rewritten as: !
! !
! (Yo-y1)*(x2-Xo) + (Yo-y2)*(Xo-x1) >= 0, if x2-x1 > 0 !
! or !
! (Yo-y1)*(x2-Xo) + (Yo-y2)*(Xo-x1) <= 0, if x2-x1 < 0 !
! !
! and both versions can be merged into essentially the condition !
! that (Yo-y1)*(x2-Xo)+(Yo-y2)*(Xo-x1) has the same sign as x2-x1. !
! That is, the product of these two must be positive or zero. !
! !
!=======================================================================
!
! Imported variable declarations.
!
integer, intent(in) :: Nb
real, intent(in) :: Xo, Yo
# ifdef ASSUMED_SHAPE
real, intent(inout) :: Xb(:), Yb(:)
# else
real, intent(inout) :: Xb(Nb+1), Yb(Nb+1)
# endif
!
! Local variable declarations.
!
integer, parameter :: Nstep =128
integer :: crossings, i, inc, k, kk, nc
integer, dimension(Nstep) :: Sindex
real :: dx1, dx2, dxy
!
!-----------------------------------------------------------------------
! Find intersections.
!-----------------------------------------------------------------------
!
! Set crossings counter and close the contour of the polygon.
!
crossings=0
Xb(Nb+1)=Xb(1)
Yb(Nb+1)=Yb(1)
!
! The search is optimized. First select the indices of segments
! where Xb(k) is different from Xb(k+1) and Xo falls between them.
! Then, further investigate these segments in a separate loop.
! Doing it in two stages takes less time because the first loop is
! pipelined.
!
DO kk=0,Nb-1,Nstep
nc=0
DO k=kk+1,MIN(kk+Nstep,Nb)
IF (((Xb(k+1)-Xo)*(Xo-Xb(k)).ge.0.0_m8).and. &
& (Xb(k).ne.Xb(k+1))) THEN
nc=nc+1
Sindex(nc)=k
END IF
END DO
DO i=1,nc
k=Sindex(i)
IF (Xb(k).ne.Xb(k+1)) THEN
dx1=Xo-Xb(k)
dx2=Xb(k+1)-Xo
dxy=dx2*(Yo-Yb(k))-dx1*(Yb(k+1)-Yo)
inc=0
IF ((Xb(k).eq.Xo).and.(Yb(k).eq.Yo)) THEN
crossings=1
goto 10
ELSE IF (((dx1.eq.0.0_m8).and.(Yo.ge.Yb(k ))).or. &
& ((dx2.eq.0.0_m8).and.(Yo.ge.Yb(k+1)))) THEN
inc=1
ELSE IF ((dx1*dx2.gt.0.0_m8).and. &
& ((Xb(k+1)-Xb(k))*dxy.ge.0.0_m8)) THEN ! see note 1
inc=2
END IF
IF (Xb(k+1).gt.Xb(k)) THEN
crossings=crossings+inc
ELSE
crossings=crossings-inc
END IF
END IF
END DO
END DO
!
! Determine if point (Xo,Yo) is inside of closed polygon.
!
10 IF (crossings.eq.0) THEN
sinside=.FALSE.
ELSE
sinside=.TRUE.
END IF
RETURN
END FUNCTION sinside
# endif
#endif
END MODULE interp_swan_mod