#include "w3macros.h" !/ ------------------------------------------------------------------- / MODULE W3GSRUMD !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 15-Jun-2012 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ 12-Nov-2010 : Change T_NNS, W3NN*, W3SORT, W3ISRT to public. !/ Add W3GFIJ (public). Implement r4 & r8 interfaces. !/ Change to number of search buckets based on !/ dimensions of input grid. ( version 3.14 ) !/ 01-Dec-2010 : Assign cells to buckets based on overlap. The !/ nearest-neighbor bucket search is removed (no longer !/ needed). Add support for tripole grids (JCLO). !/ Add W3GFCD (public). ( version 3.14 ) !/ 06-Dec-2010 : Remove restriction on longitude range. Change ICLO !/ to integer and remove JCLO. Implement support for !/ r4 and r8 source grids. ( version 3.14 ) !/ 15-Jun-2012 : Fixed various format statem,ents that gave compile !/ warnings with Intel compiler on NCEP R&D machine !/ zeus (H. L. Tolman) ( version 4.07 ) !/ ! 1. Purpose : ! ! Search and regrid utilities (data structures and associated ! methods) for logically rectangular grids. ! ! The grid-search-utility (GSU) object can be used for rapid searching ! of the associated grid to identify a grid cell that encloses a target ! point and to compute interpolation weights. The GSU object maintains ! internal pointers to the associated grid coordinate arrays. Rapid ! searching is done using a bucket search algorithm. The search buckets ! are based on the bounding box for the associated grid and an optional ! user defined approximate number of grid cells per search bucket. ! ! Grid cells are identified by the cell's lower-left corner grid point. ! The vertices (grid points) associated with a grid cell are assigned a ! sequential index in a counterclockwise order beginning with the cell's ! lower-left corner grid point. That is, when moving from vertex 1 to ! vertex 2 to vertex 3, etc., the grid cell interior is always to the left. ! Note that though cell will be counterclockwise w.r.t. indices, this does ! not necessarily mean that the cell will be counterclockwise geographically, ! specifically in situation of curvilinear grid. ! ! (x4,y4) (x3,y3) ! _____________________ ! / / ! / / ! / / ! / / ! /____________________/ ! (x1,y1) (x2,y2) ! ! There are two types of index space closure supported for lat/lon grids. ! ! 1) Simple closure: Grid is periodic in the i-index and wraps ! at i=NX+1. In other words, (NX+1,j) => (1,j). ! ! 2) Tripole grid closure: Grid is periodic in the i-index and ! and wraps at i=NX+1 and has closure at j=NY+1. In other words, ! (NX+1,j<=NY) => (1,j) and (i,NY+1) => (MOD(NX-i+1,NX)+1,NY). ! The tripole grid closure requires that NX be even. ! ! A simple interpolation example: ! ! ----------------------------------------------------------- ! ! Define data ! TYPE(T_GSU) :: GSU ! INTEGER :: NX, NY !source grid dimensions ! REAL, POINTER :: XS(:,:), YS(:,:) !source grid coordinates ! REAL :: FS(NX,NY) !source field ! INTEGER :: NT !number of target points ! REAL :: XT(NT), YT(NT), FT(NT) !target coordinates and field ! INTEGER :: IS(4), JS(4) !interpolation points ! REAL :: RW(4) !interpolation weights ! ! ! Setup source grid and field and target points ! < ... > ! ! ! Create grid-search-utility object for source grid ! GSU = W3GSUC( .TRUE., .FALSE., .FALSE., NX, NY, XS, YS ) ! ! ! Interpolate source field to target points ! DO K=1,NT ! FT(K) = 0 ! IF ( W3GRMP( GSU, XT(K), YT(K), IS, JS, RW ) ) THEN ! DO L=1,4 ! FT(K) = FT(K) + RW(L)*FS(IS(L),JS(L)) ! END DO ! END IF ! END DO ! ! ! Destroy grid-search-utility object ! CALL W3GSUD( GSU ) ! ----------------------------------------------------------- ! ! 2. Variables and types : ! ! All module variables and types are scoped private by default. ! The private module variables and types are not listed in this section. ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! MSKC_NONE I.P. Public Named constant identifying a non-masked ! enclosing grid cell ! MSKC_PART I.P. Public Named constant identifying a partially ! masked enclosing grid cell ! MSKC_FULL I.P. Public Named constant identifying a fully ! masked enclosing grid cell ! ICLO_NONE I.P. Public Named constant identifying a grid with ! no closure in index space ! ICLO_SMPL I.P. Public Named constant identifying a grid with ! simple closure: (NX+1,j) => (1,j) ! ICLO_TRPL I.P. Public Named constant identifying a grid with ! tripole closure: (NX+1,j<=NY) => (1,j) ! and (i,NY+1) => (MOD(NX-i+1,NX)+1,NY) ! T_GSU TYPE Public Grid-search-utility type (opaque) ! T_NNS TYPE Public Nearest-neighbor grid-point search type ! ---------------------------------------------------------------- ! ! 3. Subroutines and functions : ! ! All module subroutines and functions are scoped private by default. ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! W3GSUC Func. Public Create grid-search-utility object. ! W3GSUD Subr. Public Destroy grid-search-utility object. ! W3GSUP Subr. Public Print grid-search-utility object to stdout. ! W3GFCL Func. Public Find grid cell that encloses target point (bucket search). ! W3GFCD Func. Public Find grid cell that encloses target point (direct search). ! W3GFPT Func. Public Find grid point that is closest to target point. ! W3GRMP Func. Public Compute interpolation coeff. from grid. ! W3DIST Func. Public Compute distance between two points. ! W3INAN Func. Public Check if input is infinite or NaN. ! W3NNSC Func. Public Create nearest-neighbor-search object. ! W3NNSD Subr. Public Destroy nearest-neighbor-search object. ! W3NNSP Subr. Public Print nearest-neighbor-search object to stdout. ! W3SORT Subr. Public Sort input arrays in increasing order. ! W3ISRT Subr. Public Insert data into array. ! W3CKCL Func. Public Check if point lies within grid cell. ! ---------------------------------------------------------------- ! W3RMBL Subr. Private Compute bilinear interpolation coeff. from cell. ! ---------------------------------------------------------------- ! ! 4. Subroutines and functions used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! STRACE Subr. W3SERVMD Subroutine tracing. ! EXTCDE Subr. W3SERVMD Abort program with exit code. ! ---------------------------------------------------------------- ! ! 5. Remarks : ! ! - The GSU object is an "opaque" object. This means that the ! internals of the object are not accessible outside this module. ! - The burden is upon the user to invoke the destroy method when ! finished with a GSU object. If created GSU objects are ! not properly destroyed, then memory leaks may be introduced. ! - Currently, this module does not work correctly for tripole ! grids. ! ! 6. Switches : ! ! !/S Enable subroutine tracing. ! ! 7. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ Use associated modules !/ USE W3SERVMD, ONLY: EXTCDE !/ !/ Specify default data typing !/ IMPLICIT NONE !/ !/ Specify default accessibility !/ PRIVATE !/ !/ Public module methods !/ PUBLIC W3GSUC PUBLIC W3GSUD PUBLIC W3GSUP PUBLIC W3GFCL PUBLIC W3GFCD PUBLIC W3GFPT PUBLIC W3GFIJ PUBLIC W3GRMP PUBLIC W3DIST PUBLIC W3INAN PUBLIC W3NNSC PUBLIC W3NNSD PUBLIC W3NNSP PUBLIC W3SORT PUBLIC W3ISRT PUBLIC W3CKCL !/ !/ Public return codes !/ INTEGER, PARAMETER, PUBLIC :: MSKC_NONE = 0 INTEGER, PARAMETER, PUBLIC :: MSKC_PART = 1 INTEGER, PARAMETER, PUBLIC :: MSKC_FULL = 2 !/ !/ Public index closure types (for lat/lon grids only) !/ ICLO_NONE : no closure in index space !/ ICLO_SMPL : closure in i-index at i=NX+1: (NX+1,j) => (1,j) !/ ICLO_TRPL : tripole grid closure: (NX+1,j<=NY) => (1,j) and !/ (i,NY+1) => (MOD(NX-i+1,NX)+1,NY) !/ INTEGER, PARAMETER, PUBLIC :: ICLO_NONE = 0 INTEGER, PARAMETER, PUBLIC :: ICLO_SMPL = 1 INTEGER, PARAMETER, PUBLIC :: ICLO_TRPL = 2 !/ !/ Public grid-search-utility type !/ This is an opaque type -- that is, it's internals are private and only !/ accessible to subroutines in this module where the type is declared. !/ TYPE, PUBLIC :: T_GSU PRIVATE TYPE(CLASS_GSU), POINTER :: PTR => NULL() END TYPE T_GSU !/ !/ Private grid-search-utility class !/ TYPE :: CLASS_GSU LOGICAL :: IJG ! grid array ordering flag: T = (NX,NY), F = (NY,NX) LOGICAL :: LLG ! spherical coordinate flag of associated grid INTEGER :: ICLO ! parameter indicating type of index space closure ! this flag must be set by the user LOGICAL :: LCLO ! flag indicating longitudinal periodicity ! this flag is calculated internally ! ICLO != ICLO_NONE => LCLO = T LOGICAL :: L360 ! flag indicating longitude range: ! T = [0:360], F = [-180:180] INTEGER :: GKIND ! kind (precision: 4 or 8) of associated grid INTEGER :: NX, NY ! dimensions of associated grid REAL(4), POINTER :: XG4(:,:), YG4(:,:) ! coordinates of associated grid (r4) REAL(8), POINTER :: XG8(:,:), YG8(:,:) ! coordinates of associated grid (r8) TYPE(T_NNS), POINTER :: NNP ! nearest-neighbor point search indices object INTEGER :: NBX, NBY ! number of buckets in each spatial direction REAL(8) :: DXB, DYB ! spatial extent of each search bucket REAL(8) :: XMIN, YMIN, XMAX, YMAX ! bounding box of search domain TYPE(T_BKT), POINTER :: B(:,:) ! array of search buckets END TYPE CLASS_GSU !/ !/ Private search bucket type !/ TYPE :: T_BKT INTEGER :: N ! number of cells in bucket INTEGER, POINTER :: I(:) ! i-index of cell c INTEGER, POINTER :: J(:) ! j-index of cell c END TYPE T_BKT !/ !/ Public nearest-neighbor grid-point search type !/ TYPE, PUBLIC :: T_NNS INTEGER :: NLVL ! number of nnbr levels INTEGER :: NNBR ! total number of nnbr's INTEGER, POINTER :: N1(:) ! starting nearest-nbr loop index for level l INTEGER, POINTER :: N2(:) ! ending nearest-nbr loop index for level l INTEGER, POINTER :: DI(:) ! i-index delta for nearest-nbr n INTEGER, POINTER :: DJ(:) ! j-index delta for nearest-nbr n END TYPE T_NNS !/ !/ Private module parameters !/ REAL(8), PARAMETER :: PI = 3.14159265358979323846D0 REAL(8), PARAMETER :: PI2 = 2D0*PI REAL(8), PARAMETER :: PI3H = 3D0*PI/2D0 REAL(8), PARAMETER :: D2R = PI/180D0 REAL(8), PARAMETER :: R2D = 1D0/D2R REAL(8), PARAMETER :: D360 = 360D0 REAL(8), PARAMETER :: D270 = 270D0 REAL(8), PARAMETER :: D180 = 180D0 REAL(8), PARAMETER :: D90 = 90D0 REAL(8), PARAMETER :: ZERO = 0.0D0 REAL(8), PARAMETER :: ONE = 1.0D0 REAL(8), PARAMETER :: HALF = 0.5D0 !/ !/ Module Interfaces !/ INTERFACE W3GSUC MODULE PROCEDURE W3GSUC_R4 MODULE PROCEDURE W3GSUC_R8 END INTERFACE W3GSUC INTERFACE W3GFCL MODULE PROCEDURE W3GFCL_R4 MODULE PROCEDURE W3GFCL_R8 END INTERFACE W3GFCL INTERFACE W3GFCD MODULE PROCEDURE W3GFCD_R4 MODULE PROCEDURE W3GFCD_R8 END INTERFACE W3GFCD INTERFACE W3GFPT MODULE PROCEDURE W3GFPT_R4 MODULE PROCEDURE W3GFPT_R8 END INTERFACE W3GFPT INTERFACE W3GFIJ MODULE PROCEDURE W3GFIJ_R4 MODULE PROCEDURE W3GFIJ_R8 END INTERFACE W3GFIJ INTERFACE W3GRMP MODULE PROCEDURE W3GRMP_R4 MODULE PROCEDURE W3GRMP_R8 END INTERFACE W3GRMP INTERFACE W3RMBL MODULE PROCEDURE W3RMBL_R4 MODULE PROCEDURE W3RMBL_R8 END INTERFACE W3RMBL INTERFACE W3DIST MODULE PROCEDURE W3DIST_R4 MODULE PROCEDURE W3DIST_R8 END INTERFACE W3DIST INTERFACE W3CKCL MODULE PROCEDURE W3CKCL_R4 MODULE PROCEDURE W3CKCL_R8 END INTERFACE W3CKCL INTERFACE W3SORT MODULE PROCEDURE W3SORT_R4 MODULE PROCEDURE W3SORT_R8 END INTERFACE W3SORT INTERFACE W3ISRT MODULE PROCEDURE W3ISRT_R4 MODULE PROCEDURE W3ISRT_R8 END INTERFACE W3ISRT INTERFACE W3INAN MODULE PROCEDURE W3INAN_R4 MODULE PROCEDURE W3INAN_R8 END INTERFACE W3INAN !/ CONTAINS !/ ------------------------------------------------------------------- / FUNCTION W3GSUC_R4(IJG, LLG, ICLO, NX, NY, XG, YG, NCB, NNP, DEBUG) & RESULT(GSU) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 06-Dec-2010 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ 12-Nov-2010 : Change to number of search buckets based on !/ dimensions of input grid. ( version 3.14 ) !/ 01-Dec-2010 : Restore NCB optional input. Assign cells to buckets !/ based on overlap. The nearest-neighbor bucket search !/ is removed (no longer needed). Add support for !/ tripole grids (JCLO). ( version 3.14 ) !/ 06-Dec-2010 : Remove restriction on longitude range. Change !/ ICLO to integer and remove JCLO. Implement r4 and r8 !/ input grid versions. ( version 3.14 ) !/ ! 1. Purpose : ! ! Create grid-search-utility (GSU) object for a logically rectangular ! grid defined by the input coordinates. ! Single precision input grid. ! ! 2. Method : ! ! 3. Parameters : ! ! Return parameter ! ---------------------------------------------------------------- ! GSU Type O Created grid-search-utility object. ! ---------------------------------------------------------------- ! ! Parameter list ! ---------------------------------------------------------------- ! IJG Log. I Logical flag indicating ordering of input ! coord. arrays: T = (NX,NY) and F = (NY,NX). ! LLG Log. I Logical flag indicating the coordinate system: ! T = spherical lat/lon (degrees) and F = Cartesian. ! ICLO Int. I Parameter indicating type of index space closure ! NX, NY Int. I Dimensions of input grid. ! XG R.A. I Pointer to array of x-coordinates of input grid. ! YG R.A. I Pointer to array of y-coordinates of input grid. ! NCB Int. I Optional (approximate) number of cells (in each ! direction) per search bucket. (default is NCB_DEFAULT) ! NCB >= 1 is required. NCB = 1 gives most efficient ! searching, but uses more memory. Increasing NCB leads ! to fewer buckets (less memory) but slower searching. ! NNP Int. I Optional maximum number of nearest-neighbor grid ! point search levels. (default is NNP_DEFAULT) ! DEBUG Log. I Optional logical flag to turn on debug mode. ! Default is FALSE. ! ! Internal parameters ! ---------------------------------------------------------------- ! NCB_DEFAULT Int. Default number of grid cells (in each direction) ! per search bucket. ! NNP_DEFAULT Int. Default maximum number of nearest-neighbor grid ! point search levels. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! - Check on correct coordinate system with global grid. ! - Check on association of input grid coordinate array pointers. ! ! 7. Remarks : ! ! - LCLO is calculated internally. ! - ICLO != ICLO_NONE => LCLO = T. ! - Periodic Cartesian grids are not allowed. ! ! 8. Structure : ! ! ----------------------------------------------------------------- ! 1. Test input ! 2. Allocate object and set grid related data and pointers ! 3. Create nearest-neighbor point search object ! 4. Construct bucket search "object" ! 5. Set return parameter ! ----------------------------------------------------------------- ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! !/T8 Enables debugging flag. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Return parameter !/ TYPE(T_GSU) :: GSU !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ LOGICAL, INTENT(IN) :: IJG LOGICAL, INTENT(IN) :: LLG INTEGER, INTENT(IN) :: ICLO INTEGER, INTENT(IN) :: NX, NY REAL(4), POINTER :: XG(:,:), YG(:,:) INTEGER, INTENT(IN), OPTIONAL :: NCB INTEGER, INTENT(IN), OPTIONAL :: NNP LOGICAL, INTENT(IN), OPTIONAL :: DEBUG !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ TYPE(CLASS_GSU), POINTER :: PTR INTEGER, PARAMETER :: NCB_DEFAULT = 1 INTEGER, PARAMETER :: NNP_DEFAULT = 2 LOGICAL :: LDBG, LBC, LPL, LNPL, LSPL INTEGER :: I, J, K, L, N, IC(4), JC(4), IB, JB, NXC, NYC INTEGER :: NS, IB1(2), IB2(2), JB1(2), JB2(2), IBC(4), JBC(4) INTEGER :: ISTEP, ISTAT REAL(8) :: X, Y, XC(4), YC(4) !/ ! -------------------------------------------------------------------- / ! 1. Test input ! SELECT CASE ( ICLO ) CASE ( ICLO_NONE, ICLO_SMPL, ICLO_TRPL ) CONTINUE CASE DEFAULT WRITE(*,'(/1A,1A,1I2/)') 'W3GSUC_R4 ERROR -- ', & 'unsupported ICLO: ',ICLO CALL EXTCDE (1) END SELECT IF ( ICLO.NE.ICLO_NONE .AND. .NOT.LLG ) THEN WRITE(*,'(/1A,1A/)') 'W3GSUC_R4 ERROR -- ', & 'index space closure with cartesian grids is not supported' CALL EXTCDE (1) END IF IF ( ICLO.EQ.ICLO_TRPL .AND. MOD(NX,2).NE.0 ) THEN WRITE(*,'(/1A,1A/)') 'W3GSUC_R4 ERROR -- ', & 'tripole grid closure requires NX even' CALL EXTCDE (1) END IF IF ( .NOT.ASSOCIATED(XG) .OR. .NOT.ASSOCIATED(YG) ) THEN WRITE(*,'(/1A,1A/)') 'W3GSUC_R4 ERROR -- ', & 'input grid coordinate array pointers are not associated' CALL EXTCDE (1) END IF IF ( PRESENT(NCB) ) THEN IF ( NCB .LE. 0 ) THEN WRITE(*,'(/1A,1A/)') 'W3GSUC_R4 ERROR -- ', & 'NCB must be greater than zero' CALL EXTCDE (1) END IF END IF ! IF ( PRESENT(DEBUG) ) THEN LDBG = DEBUG ELSE LDBG = .FALSE. END IF ! ! -------------------------------------------------------------------- / ! 2. Allocate object and set grid related data and pointers ! ALLOCATE(PTR, STAT=ISTAT) IF ( ISTAT .NE. 0 ) THEN WRITE(*,'(/1A,1A/)') 'W3GSUC_R4 ERROR -- ', & 'gsu object allocation failed' CALL EXTCDE (ISTAT) END IF PTR%IJG = IJG PTR%LLG = LLG PTR%ICLO = ICLO PTR%NX = NX PTR%NY = NY PTR%XG4 => XG PTR%YG4 => YG PTR%GKIND = 4 ! ! -------------------------------------------------------------------- / ! 3. Create nearest-neighbor point search object ! IF ( PRESENT(NNP) ) THEN PTR%NNP => W3NNSC(NNP) ELSE PTR%NNP => W3NNSC(NNP_DEFAULT) END IF ! ! -------------------------------------------------------------------- / ! 4. Construct bucket search "object" ! !-----number of cells SELECT CASE ( ICLO ) CASE ( ICLO_NONE ) NXC = NX-1; NYC = NY-1; CASE ( ICLO_SMPL ) NXC = NX; NYC = NY-1; CASE ( ICLO_TRPL ) NXC = NX; NYC = NY; END SELECT ! !-----initialize longitudinal periodicity flag (LCLO) IF ( LLG .AND. ICLO.NE.ICLO_NONE ) THEN PTR%LCLO = .TRUE. ELSE PTR%LCLO = .FALSE. END IF ! !-----check existence of longitudinal branch cut !-----check if source grid includes poles IF ( LDBG ) THEN WRITE(*,'(/A)') 'W3GSUC_R4 - check source grid' END IF LNPL = .FALSE. LSPL = .FALSE. DO I=1,NXC DO J=1,NYC !-------------create list of cell vertices IC(1) = I ; JC(1) = J ; IC(2) = I+1; JC(2) = J ; IC(3) = I+1; JC(3) = J+1; IC(4) = I ; JC(4) = J+1; DO L=1,4 !-----------------i-closure IF ( ICLO.NE.ICLO_NONE ) THEN IF ( IC(L) .LT. 1 ) IC(L) = IC(L) + NX IF ( IC(L) .GT. NX ) IC(L) = IC(L) - NX END IF !-----------------j-closure IF ( ICLO.EQ.ICLO_TRPL ) THEN IF ( JC(L) .GT. NY ) THEN JC(L) = NY IC(L) = MOD(NX-IC(L),NX) + 1 END IF END IF !-----------------copy cell vertex coordinates into local variables IF ( IJG ) THEN XC(L) = XG(IC(L),JC(L)); YC(L) = YG(IC(L),JC(L)); ELSE XC(L) = XG(JC(L),IC(L)); YC(L) = YG(JC(L),IC(L)); END IF END DO !L !-------------check if cell includes a pole or branch cut LPL = .FALSE. LBC = .FALSE. IF ( LLG ) THEN !-----------------count longitudinal branch cut crossings N = 0 DO L=1,4 K = MOD(L,4)+1 IF ( ABS(XC(K)-XC(L)) .GT. D180 ) N = N + 1 END DO !-----------------multiple longitudinal branch cut crossing => cell includes branch cut LBC = N.GT.1 IF ( LBC .AND. LDBG ) & WRITE(*,'(A,8I6)') & 'W3GSUC_R4 -- cell includes branch cut:',IC(:),JC(:) !-----------------single longitudinal branch cut crossing ! or single vertex at 90 degrees => cell includes pole LPL = N.EQ.1 .OR. COUNT(ABS(YC).EQ.D90).EQ.1 IF ( LPL.AND.MINVAL(YC).GT.ZERO ) THEN IF ( LDBG ) & WRITE(*,'(A,8I6)') & 'W3GSUC_R4 -- cell includes N-pole:',IC(:),JC(:) LNPL = .TRUE. END IF IF ( LPL.AND.MAXVAL(YC).LT.ZERO ) THEN IF ( LDBG ) & WRITE(*,'(A,8I6)') & 'W3GSUC_R4 -- cell includes S-pole:',IC(:),JC(:) LSPL = .TRUE. END IF !-----------------longitudinal branch cut crossing => longitudinal closure IF ( N.GT.0 ) PTR%LCLO = .TRUE. END IF !LLG END DO !J END DO !I ! !-----compute domain for search buckets ! if longitudinal periodicity, then force domain in x to [0:360] ! if grid includes north pole, then set ymax = 90 degrees ! if grid includes south pole, then set ymin = -90 degrees PTR%XMIN = MINVAL(XG); PTR%XMAX = MAXVAL(XG); PTR%YMIN = MINVAL(YG); PTR%YMAX = MAXVAL(YG); IF ( PTR%LCLO ) THEN PTR%XMIN = ZERO; PTR%XMAX = D360; END IF IF ( LSPL ) PTR%YMIN = -D90 IF ( LNPL ) PTR%YMAX = D90 PTR%L360 = PTR%XMIN.GE.ZERO ! !-----compute number of search buckets and bucket size IF ( PRESENT(NCB) ) THEN PTR%NBX = MAX(1,NX/NCB) PTR%NBY = MAX(1,NY/NCB) ELSE PTR%NBX = MAX(1,NX/NCB_DEFAULT) PTR%NBY = MAX(1,NY/NCB_DEFAULT) END IF PTR%DXB = (PTR%XMAX-PTR%XMIN)/REAL(PTR%NBX) PTR%DYB = (PTR%YMAX-PTR%YMIN)/REAL(PTR%NBY) ! !-----print debug info IF ( LDBG ) THEN WRITE(*,'(/A,1I2,1L2,1I2)') 'W3GSUC_R4 - ICLO,LCLO,GKIND: ', & PTR%ICLO,PTR%LCLO,PTR%GKIND WRITE(*,'(A,4E14.6)') 'W3GSUC_R4 - grid search domain:', & PTR%XMIN,PTR%YMIN,PTR%XMAX,PTR%YMAX WRITE(*,'(A,2I6)') 'W3GSUC_R4 - number of search buckets:', & PTR%NBX,PTR%NBY WRITE(*,'(A,2E14.6)') 'W3GSUC_R4 - search bucket size:', & PTR%DXB,PTR%DYB END IF ! !-----allocate array of search buckets ALLOCATE(PTR%B(PTR%NBY,PTR%NBX),STAT=ISTAT) IF ( ISTAT .NE. 0 ) THEN WRITE(*,'(/1A,1A/)') 'W3GSUC_R4 ERROR -- ', & 'search bucket array allocation failed' CALL EXTCDE (ISTAT) END IF ! !-----BEGIN ISTEP_LOOP ! first step: compute number of cells in each bucket ! second step: allocate buckets and assign cells to buckets ISTEP_LOOP: DO ISTEP=1,2 ! !-----allocate search bucket cell lists IF ( ISTEP .EQ. 2 ) THEN DO IB=1,PTR%NBX DO JB=1,PTR%NBY NULLIFY(PTR%B(JB,IB)%I) NULLIFY(PTR%B(JB,IB)%J) IF ( PTR%B(JB,IB)%N .GT. 0 ) THEN ALLOCATE(PTR%B(JB,IB)%I(PTR%B(JB,IB)%N),STAT=ISTAT) IF ( ISTAT .NE. 0 ) THEN WRITE(*,'(/1A,2A,3I6/)') 'W3GSUC_R4 ERROR -- ', & 'search bucket cell-i list allocation failed -- ', & 'bucket: ',IB,JB,N CALL EXTCDE (ISTAT) END IF ALLOCATE(PTR%B(JB,IB)%J(PTR%B(JB,IB)%N),STAT=ISTAT) IF ( ISTAT .NE. 0 ) THEN WRITE(*,'(/1A,2A,3I6/)') 'W3GSUC_R4 ERROR -- ', & 'search bucket cell-j list allocation failed -- ', & 'bucket: ',IB,JB,N CALL EXTCDE (ISTAT) END IF END IF END DO END DO END IF !ISTEP.EQ.2 ! !-----build search bucket cell lists PTR%B(:,:)%N = 0 DO I=1,NXC DO J=1,NYC IF ( ICLO.EQ.ICLO_TRPL ) THEN IF ( J.EQ.NYC .AND. I.GT.NX/2+1 ) CYCLE ENDIF !-------------create list of cell vertices IC(1) = I ; JC(1) = J ; IC(2) = I+1; JC(2) = J ; IC(3) = I+1; JC(3) = J+1; IC(4) = I ; JC(4) = J+1; DO L=1,4 !-----------------i-closure IF ( ICLO.NE.ICLO_NONE ) THEN IF ( IC(L) .LT. 1 ) IC(L) = IC(L) + NX IF ( IC(L) .GT. NX ) IC(L) = IC(L) - NX END IF !-----------------j-closure IF ( ICLO.EQ.ICLO_TRPL ) THEN IF ( JC(L) .GT. NY ) THEN JC(L) = NY IC(L) = MOD(NX-IC(L),NX) + 1 END IF END IF !-----------------copy cell vertex coordinates into local variables IF ( IJG ) THEN XC(L) = XG(IC(L),JC(L)); YC(L) = YG(IC(L),JC(L)); ELSE XC(L) = XG(JC(L),IC(L)); YC(L) = YG(JC(L),IC(L)); END IF END DO !L !-------------check if cell includes a pole or branch cut LPL = .FALSE. LBC = .FALSE. IF ( LLG ) THEN !-----------------shift longitudes to appropriate range XC = MOD(XC,D360) IF ( PTR%LCLO .OR. PTR%L360 ) THEN WHERE ( XC.LT.ZERO ) XC = XC + D360 ELSE WHERE ( XC.GT.D180 ) XC = XC - D360 END IF !-----------------count longitudinal branch cut crossings N = 0 DO L=1,4 K = MOD(L,4)+1 IF ( ABS(XC(K)-XC(L)) .GT. D180 ) N = N + 1 END DO !-----------------multiple longitudinal branch cut crossing => cell includes branch cut LBC = N.GT.1 !-----------------single longitudinal branch cut crossing ! or single vertex at 90 degrees => cell includes pole LPL = N.EQ.1 .OR. COUNT(ABS(YC).EQ.D90).EQ.1 END IF !LLG !-------------set bucket id for each cell vertex DO L=1,4 IBC(L) = INT((XC(L)-PTR%XMIN)/PTR%DXB)+1 IF ( .NOT.PTR%LCLO ) IBC(L) = MIN(IBC(L),PTR%NBX) JBC(L) = MIN(INT((YC(L)-PTR%YMIN)/PTR%DYB)+1,PTR%NBY) END DO !L !-------------set bucket overlap bounds IF ( LPL ) THEN !---------------cell includes pole: overlap includes full longitudinal range NS = 1 IB1(1) = 1 IB2(1) = PTR%NBX IF ( MINVAL(YC).GT.ZERO ) THEN JB1(1) = MAX(1,MINVAL(JBC)) JB2(1) = PTR%NBY END IF IF ( MAXVAL(YC).LT.ZERO ) THEN JB1(1) = 1 JB2(1) = MIN(PTR%NBY,MAXVAL(JBC)) END IF IB1(2) = 0 IB2(2) = 0 JB1(2) = 0 JB2(2) = 0 ELSE IF ( LBC ) THEN !---------------cell includes branch cut: split overlap into two sets NS = 2 IB1(1) = PTR%NBX IB2(1) = PTR%NBX IB1(2) = 1 IB2(2) = 1 DO L=1,4 IF ( IBC(L) .GT. PTR%NBX/2 ) THEN IB1(1) = MIN(IB1(1),IBC(L)) ELSE IB2(2) = MAX(IB2(2),IBC(L)) END IF END DO !L JB1(:) = MAX(1,MINVAL(JBC)) JB2(:) = MIN(PTR%NBY,MAXVAL(JBC)) ELSE !---------------default: overlap computed from min/max NS = 1 IB1(1) = MAX(1,MINVAL(IBC)) IB2(1) = MIN(PTR%NBX,MAXVAL(IBC)) JB1(1) = MAX(1,MINVAL(JBC)) JB2(1) = MIN(PTR%NBY,MAXVAL(JBC)) IB1(2) = 0 IB2(2) = 0 JB1(2) = 0 JB2(2) = 0 END IF !-------------debug output IF ( LDBG .AND. ISTEP.EQ.1 ) THEN WRITE(*,'(/A,2I6)') 'W3GSUC_R4 -- BUCKET SORT:',I,J WRITE(*,'(A,2L6,1I6)') 'W3GSUC_R4 -- LBC,LPL:',LBC,LPL WRITE(*,'(A,4I6)') 'W3GSUC_R4 -- IC:',IC(:) WRITE(*,'(A,4I6)') 'W3GSUC_R4 -- JC:',JC(:) WRITE(*,'(A,4E14.6)') 'W3GSUC_R4 -- XC:',XC(:) WRITE(*,'(A,4E14.6)') 'W3GSUC_R4 -- YC:',YC(:) WRITE(*,'(A,4I6)') 'W3GSUC_R4 -- IBC:',IBC(:) WRITE(*,'(A,4I6)') 'W3GSUC_R4 -- JBC:',JBC(:) WRITE(*,'(A,1I6)') 'W3GSUC_R4 -- NS:',NS WRITE(*,'(A,4I6)') 'W3GSUC_R4 -- IB1:',IB1(:) WRITE(*,'(A,4I6)') 'W3GSUC_R4 -- JB1:',JB1(:) WRITE(*,'(A,4I6)') 'W3GSUC_R4 -- IB2:',IB2(:) WRITE(*,'(A,4I6)') 'W3GSUC_R4 -- JB2:',JB2(:) END IF !-------------assign cell to buckets based on overlap DO K=1,NS DO IB=IB1(K),IB2(K) DO JB=JB1(K),JB2(K) PTR%B(JB,IB)%N = PTR%B(JB,IB)%N + 1 IF ( ISTEP .EQ. 2 ) THEN PTR%B(JB,IB)%I(PTR%B(JB,IB)%N) = IC(1) PTR%B(JB,IB)%J(PTR%B(JB,IB)%N) = JC(1) END IF END DO !JB END DO !IB END DO !K END DO !J END DO !I ! !-----END ISTEP_LOOP END DO ISTEP_LOOP ! !-----print debug info IF ( LDBG ) THEN WRITE(*,'(/A,3I6,4E14.6)') 'W3GSUC_R4 - search bucket list:' WRITE(*,'(3A6,4A14)') 'I','J','N','X1','Y1','X2','Y2' DO IB=1,PTR%NBX DO JB=1,PTR%NBY WRITE(*,'(3I6,4E14.6)') IB,JB,PTR%B(JB,IB)%N, & PTR%XMIN+(IB-1)*PTR%DXB,PTR%YMIN+(JB-1)*PTR%DYB, & PTR%XMIN+(IB-0)*PTR%DXB,PTR%YMIN+(JB-0)*PTR%DYB END DO END DO END IF ! ! -------------------------------------------------------------------- / ! 5. Set return parameter ! GSU%PTR => PTR !/ !/ End of W3GSUC_R4 -------------------------------------------------- / !/ END FUNCTION W3GSUC_R4 !/ ------------------------------------------------------------------- / FUNCTION W3GSUC_R8(IJG, LLG, ICLO, NX, NY, XG, YG, NCB, NNP, DEBUG) & RESULT(GSU) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 06-Dec-2010 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ 12-Nov-2010 : Change to number of search buckets based on !/ dimensions of input grid. ( version 3.14 ) !/ 01-Dec-2010 : Restore NCB optional input. Assign cells to buckets !/ based on overlap. The nearest-neighbor bucket search !/ is removed (no longer needed). Add support for !/ tripole grids (JCLO). ( version 3.14 ) !/ 06-Dec-2010 : Remove restriction on longitude range. Change !/ ICLO to integer and remove JCLO. Implement r4 and r8 !/ input grid versions. ( version 3.14 ) !/ ! 1. Purpose : ! ! Create grid-search-utility (GSU) object for a logically rectangular ! grid defined by the input coordinates. ! Double precision input grid. ! ! 2. Method : ! ! 3. Parameters : ! ! Return parameter ! ---------------------------------------------------------------- ! GSU Type O Created grid-search-utility object. ! ---------------------------------------------------------------- ! ! Parameter list ! ---------------------------------------------------------------- ! IJG Log. I Logical flag indicating ordering of input ! coord. arrays: T = (NX,NY) and F = (NY,NX). ! LLG Log. I Logical flag indicating the coordinate system: ! T = spherical lat/lon (degrees) and F = Cartesian. ! ICLO Int. I Parameter indicating type of index space closure ! NX, NY Int. I Dimensions of input grid. ! XG R.A. I Pointer to array of x-coordinates of input grid. ! YG R.A. I Pointer to array of y-coordinates of input grid. ! NCB Int. I Optional (approximate) number of cells (in each ! direction) per search bucket. (default is NCB_DEFAULT) ! NCB >= 1 is required. NCB = 1 gives most efficient ! searching, but uses more memory. Increasing NCB leads ! to fewer buckets (less memory) but slower searching. ! NNP Int. I Optional maximum number of nearest-neighbor grid ! point search levels. (default is NNP_DEFAULT) ! DEBUG Log. I Optional logical flag to turn on debug mode. ! Default is FALSE. ! ! Internal parameters ! ---------------------------------------------------------------- ! NCB_DEFAULT Int. Default number of grid cells (in each direction) ! per search bucket. ! NNP_DEFAULT Int. Default maximum number of nearest-neighbor grid ! point search levels. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! - Check on correct coordinate system with global grid. ! - Check on association of input grid coordinate array pointers. ! ! 7. Remarks : ! ! - LCLO is calculated internally. ! - ICLO != ICLO_NONE => LCLO = T. ! - Periodic Cartesian grids are not allowed. ! ! 8. Structure : ! ! ----------------------------------------------------------------- ! 1. Test input ! 2. Allocate object and set grid related data and pointers ! 3. Create nearest-neighbor point search object ! 4. Construct bucket search "object" ! 5. Set return parameter ! ----------------------------------------------------------------- ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! !/T8 Enables debugging flag. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Return parameter !/ TYPE(T_GSU) :: GSU !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ LOGICAL, INTENT(IN) :: IJG LOGICAL, INTENT(IN) :: LLG INTEGER, INTENT(IN) :: ICLO INTEGER, INTENT(IN) :: NX, NY REAL(8), POINTER :: XG(:,:), YG(:,:) INTEGER, INTENT(IN), OPTIONAL :: NCB INTEGER, INTENT(IN), OPTIONAL :: NNP LOGICAL, INTENT(IN), OPTIONAL :: DEBUG !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ TYPE(CLASS_GSU), POINTER :: PTR INTEGER, PARAMETER :: NCB_DEFAULT = 1 INTEGER, PARAMETER :: NNP_DEFAULT = 2 LOGICAL :: LDBG, LBC, LPL, LNPL, LSPL INTEGER :: I, J, K, L, N, IC(4), JC(4), IB, JB, NXC, NYC INTEGER :: NS, IB1(2), IB2(2), JB1(2), JB2(2), IBC(4), JBC(4) INTEGER :: ISTEP, ISTAT REAL(8) :: X, Y, XC(4), YC(4) !/ ! -------------------------------------------------------------------- / ! 1. Test input ! SELECT CASE ( ICLO ) CASE ( ICLO_NONE, ICLO_SMPL, ICLO_TRPL ) CONTINUE CASE DEFAULT WRITE(*,'(/1A,1A,1I2/)') 'W3GSUC_R8 ERROR -- ', & 'unsupported ICLO: ',ICLO CALL EXTCDE (1) END SELECT IF ( ICLO.NE.ICLO_NONE .AND. .NOT.LLG ) THEN WRITE(*,'(/1A,1A/)') 'W3GSUC_R8 ERROR -- ', & 'index space closure with cartesian grids is not supported' CALL EXTCDE (1) END IF IF ( ICLO.EQ.ICLO_TRPL .AND. MOD(NX,2).NE.0 ) THEN WRITE(*,'(/1A,1A/)') 'W3GSUC_R8 ERROR -- ', & 'tripole grid closure requires NX even' CALL EXTCDE (1) END IF IF ( .NOT.ASSOCIATED(XG) .OR. .NOT.ASSOCIATED(YG) ) THEN WRITE(*,'(/1A,1A/)') 'W3GSUC_R8 ERROR -- ', & 'input grid coordinate array pointers are not associated' CALL EXTCDE (1) END IF IF ( PRESENT(NCB) ) THEN IF ( NCB .LE. 0 ) THEN WRITE(*,'(/1A,1A/)') 'W3GSUC_R8 ERROR -- ', & 'NCB must be greater than zero' CALL EXTCDE (1) END IF END IF ! IF ( PRESENT(DEBUG) ) THEN LDBG = DEBUG ELSE LDBG = .FALSE. END IF ! ! -------------------------------------------------------------------- / ! 2. Allocate object and set grid related data and pointers ! ALLOCATE(PTR, STAT=ISTAT) IF ( ISTAT .NE. 0 ) THEN WRITE(*,'(/1A,1A/)') 'W3GSUC_R8 ERROR -- ', & 'gsu object allocation failed' CALL EXTCDE (ISTAT) END IF PTR%IJG = IJG PTR%LLG = LLG PTR%ICLO = ICLO PTR%NX = NX PTR%NY = NY PTR%XG8 => XG PTR%YG8 => YG PTR%GKIND = 8 ! ! -------------------------------------------------------------------- / ! 3. Create nearest-neighbor point search object ! IF ( PRESENT(NNP) ) THEN PTR%NNP => W3NNSC(NNP) ELSE PTR%NNP => W3NNSC(NNP_DEFAULT) END IF ! ! -------------------------------------------------------------------- / ! 4. Construct bucket search "object" ! !-----number of cells SELECT CASE ( ICLO ) CASE ( ICLO_NONE ) NXC = NX-1; NYC = NY-1; CASE ( ICLO_SMPL ) NXC = NX; NYC = NY-1; CASE ( ICLO_TRPL ) NXC = NX; NYC = NY; END SELECT ! !-----initialize longitudinal periodicity flag (LCLO) IF ( LLG .AND. ICLO.NE.ICLO_NONE ) THEN PTR%LCLO = .TRUE. ELSE PTR%LCLO = .FALSE. END IF ! !-----check existence of longitudinal branch cut !-----check if source grid includes poles IF ( LDBG ) THEN WRITE(*,'(/A)') 'W3GSUC_R8 - check source grid' END IF LNPL = .FALSE. LSPL = .FALSE. DO I=1,NXC DO J=1,NYC !-------------create list of cell vertices IC(1) = I ; JC(1) = J ; IC(2) = I+1; JC(2) = J ; IC(3) = I+1; JC(3) = J+1; IC(4) = I ; JC(4) = J+1; DO L=1,4 !-----------------i-closure IF ( ICLO.NE.ICLO_NONE ) THEN IF ( IC(L) .LT. 1 ) IC(L) = IC(L) + NX IF ( IC(L) .GT. NX ) IC(L) = IC(L) - NX END IF !-----------------j-closure IF ( ICLO.EQ.ICLO_TRPL ) THEN IF ( JC(L) .GT. NY ) THEN JC(L) = NY IC(L) = MOD(NX-IC(L),NX) + 1 END IF END IF !-----------------copy cell vertex coordinates into local variables IF ( IJG ) THEN XC(L) = XG(IC(L),JC(L)); YC(L) = YG(IC(L),JC(L)); ELSE XC(L) = XG(JC(L),IC(L)); YC(L) = YG(JC(L),IC(L)); END IF END DO !L !-------------check if cell includes a pole or branch cut LPL = .FALSE. LBC = .FALSE. IF ( LLG ) THEN !-----------------count longitudinal branch cut crossings N = 0 DO L=1,4 K = MOD(L,4)+1 IF ( ABS(XC(K)-XC(L)) .GT. D180 ) N = N + 1 END DO !-----------------multiple longitudinal branch cut crossing => cell includes branch cut LBC = N.GT.1 IF ( LBC .AND. LDBG ) & WRITE(*,'(A,8I6)') & 'W3GSUC_R8 -- cell includes branch cut:',IC(:),JC(:) !-----------------single longitudinal branch cut crossing ! or single vertex at 90 degrees => cell includes pole LPL = N.EQ.1 .OR. COUNT(ABS(YC).EQ.D90).EQ.1 IF ( LPL.AND.MINVAL(YC).GT.ZERO ) THEN IF ( LDBG ) & WRITE(*,'(A,8I6)') & 'W3GSUC_R8 -- cell includes N-pole:',IC(:),JC(:) LNPL = .TRUE. END IF IF ( LPL.AND.MAXVAL(YC).LT.ZERO ) THEN IF ( LDBG ) & WRITE(*,'(A,8I6)') & 'W3GSUC_R8 -- cell includes S-pole:',IC(:),JC(:) LSPL = .TRUE. END IF !-----------------longitudinal branch cut crossing => longitudinal closure IF ( N.GT.0 ) PTR%LCLO = .TRUE. END IF !LLG END DO !J END DO !I ! !-----compute domain for search buckets ! if longitudinal periodicity, then force domain in x to [0:360] ! if grid includes north pole, then set ymax = 90 degrees ! if grid includes south pole, then set ymin = -90 degrees PTR%XMIN = MINVAL(XG); PTR%XMAX = MAXVAL(XG); PTR%YMIN = MINVAL(YG); PTR%YMAX = MAXVAL(YG); IF ( PTR%LCLO ) THEN PTR%XMIN = ZERO; PTR%XMAX = D360; END IF IF ( LSPL ) PTR%YMIN = -D90 IF ( LNPL ) PTR%YMAX = D90 PTR%L360 = PTR%XMIN.GE.ZERO ! !-----compute number of search buckets and bucket size IF ( PRESENT(NCB) ) THEN PTR%NBX = MAX(1,NX/NCB) PTR%NBY = MAX(1,NY/NCB) ELSE PTR%NBX = MAX(1,NX/NCB_DEFAULT) PTR%NBY = MAX(1,NY/NCB_DEFAULT) END IF PTR%DXB = (PTR%XMAX-PTR%XMIN)/REAL(PTR%NBX) PTR%DYB = (PTR%YMAX-PTR%YMIN)/REAL(PTR%NBY) ! !-----print debug info IF ( LDBG ) THEN WRITE(*,'(/A,1I2,1L2,1I2)') 'W3GSUC_R8 - ICLO,LCLO,GKIND: ', & PTR%ICLO,PTR%LCLO,PTR%GKIND WRITE(*,'(A,4E14.6)') 'W3GSUC_R8 - grid search domain:', & PTR%XMIN,PTR%YMIN,PTR%XMAX,PTR%YMAX WRITE(*,'(A,2I6)') 'W3GSUC_R8 - number of search buckets:', & PTR%NBX,PTR%NBY WRITE(*,'(A,2E14.6)') 'W3GSUC_R8 - search bucket size:', & PTR%DXB,PTR%DYB END IF ! !-----allocate array of search buckets ALLOCATE(PTR%B(PTR%NBY,PTR%NBX),STAT=ISTAT) IF ( ISTAT .NE. 0 ) THEN WRITE(*,'(/1A,1A/)') 'W3GSUC_R8 ERROR -- ', & 'search bucket array allocation failed' CALL EXTCDE (ISTAT) END IF ! !-----BEGIN ISTEP_LOOP ! first step: compute number of cells in each bucket ! second step: allocate buckets and assign cells to buckets ISTEP_LOOP: DO ISTEP=1,2 ! !-----allocate search bucket cell lists IF ( ISTEP .EQ. 2 ) THEN DO IB=1,PTR%NBX DO JB=1,PTR%NBY NULLIFY(PTR%B(JB,IB)%I) NULLIFY(PTR%B(JB,IB)%J) IF ( PTR%B(JB,IB)%N .GT. 0 ) THEN ALLOCATE(PTR%B(JB,IB)%I(PTR%B(JB,IB)%N),STAT=ISTAT) IF ( ISTAT .NE. 0 ) THEN WRITE(*,'(/1A,2A,3I6/)') 'W3GSUC_R8 ERROR -- ', & 'search bucket cell-i list allocation failed -- ', & 'bucket: ',IB,JB,N CALL EXTCDE (ISTAT) END IF ALLOCATE(PTR%B(JB,IB)%J(PTR%B(JB,IB)%N),STAT=ISTAT) IF ( ISTAT .NE. 0 ) THEN WRITE(*,'(/1A,2A,3I6/)') 'W3GSUC_R8 ERROR -- ', & 'search bucket cell-j list allocation failed -- ', & 'bucket: ',IB,JB,N CALL EXTCDE (ISTAT) END IF END IF END DO END DO END IF !ISTEP.EQ.2 ! !-----build search bucket cell lists PTR%B(:,:)%N = 0 DO I=1,NXC DO J=1,NYC IF ( ICLO.EQ.ICLO_TRPL ) THEN IF ( J.EQ.NYC .AND. I.GT.NX/2+1 ) CYCLE ENDIF !-------------create list of cell vertices IC(1) = I ; JC(1) = J ; IC(2) = I+1; JC(2) = J ; IC(3) = I+1; JC(3) = J+1; IC(4) = I ; JC(4) = J+1; DO L=1,4 !-----------------i-closure IF ( ICLO.NE.ICLO_NONE ) THEN IF ( IC(L) .LT. 1 ) IC(L) = IC(L) + NX IF ( IC(L) .GT. NX ) IC(L) = IC(L) - NX END IF !-----------------j-closure IF ( ICLO.EQ.ICLO_TRPL ) THEN IF ( JC(L) .GT. NY ) THEN JC(L) = NY IC(L) = MOD(NX-IC(L),NX) + 1 END IF END IF !-----------------copy cell vertex coordinates into local variables IF ( IJG ) THEN XC(L) = XG(IC(L),JC(L)); YC(L) = YG(IC(L),JC(L)); ELSE XC(L) = XG(JC(L),IC(L)); YC(L) = YG(JC(L),IC(L)); END IF END DO !L !-------------check if cell includes a pole or branch cut LPL = .FALSE. LBC = .FALSE. IF ( LLG ) THEN !-----------------shift longitudes to appropriate range XC = MOD(XC,D360) IF ( PTR%LCLO .OR. PTR%L360 ) THEN WHERE ( XC.LT.ZERO ) XC = XC + D360 ELSE WHERE ( XC.GT.D180 ) XC = XC - D360 END IF !-----------------count longitudinal branch cut crossings N = 0 DO L=1,4 K = MOD(L,4)+1 IF ( ABS(XC(K)-XC(L)) .GT. D180 ) N = N + 1 END DO !-----------------multiple longitudinal branch cut crossing => cell includes branch cut LBC = N.GT.1 !-----------------single longitudinal branch cut crossing ! or single vertex at 90 degrees => cell includes pole LPL = N.EQ.1 .OR. COUNT(ABS(YC).EQ.D90).EQ.1 END IF !LLG !-------------set bucket id for each cell vertex DO L=1,4 IBC(L) = INT((XC(L)-PTR%XMIN)/PTR%DXB)+1 IF ( .NOT.PTR%LCLO ) IBC(L) = MIN(IBC(L),PTR%NBX) JBC(L) = MIN(INT((YC(L)-PTR%YMIN)/PTR%DYB)+1,PTR%NBY) END DO !L !-------------set bucket overlap bounds IF ( LPL ) THEN !---------------cell includes pole: overlap includes full longitudinal range NS = 1 IB1(1) = 1 IB2(1) = PTR%NBX IF ( MINVAL(YC).GT.ZERO ) THEN JB1(1) = MAX(1,MINVAL(JBC)) JB2(1) = PTR%NBY END IF IF ( MAXVAL(YC).LT.ZERO ) THEN JB1(1) = 1 JB2(1) = MIN(PTR%NBY,MAXVAL(JBC)) END IF IB1(2) = 0 IB2(2) = 0 JB1(2) = 0 JB2(2) = 0 ELSE IF ( LBC ) THEN !---------------cell includes branch cut: split overlap into two sets NS = 2 IB1(1) = PTR%NBX IB2(1) = PTR%NBX IB1(2) = 1 IB2(2) = 1 DO L=1,4 IF ( IBC(L) .GT. PTR%NBX/2 ) THEN IB1(1) = MIN(IB1(1),IBC(L)) ELSE IB2(2) = MAX(IB2(2),IBC(L)) END IF END DO !L JB1(:) = MAX(1,MINVAL(JBC)) JB2(:) = MIN(PTR%NBY,MAXVAL(JBC)) ELSE !---------------default: overlap computed from min/max NS = 1 IB1(1) = MAX(1,MINVAL(IBC)) IB2(1) = MIN(PTR%NBX,MAXVAL(IBC)) JB1(1) = MAX(1,MINVAL(JBC)) JB2(1) = MIN(PTR%NBY,MAXVAL(JBC)) IB1(2) = 0 IB2(2) = 0 JB1(2) = 0 JB2(2) = 0 END IF !-------------debug output IF ( LDBG .AND. ISTEP.EQ.1 ) THEN WRITE(*,'(/A,2I6)') 'W3GSUC_R8 -- BUCKET SORT:',I,J WRITE(*,'(A,2L6,1I6)') 'W3GSUC_R8 -- LBC,LPL:',LBC,LPL WRITE(*,'(A,4I6)') 'W3GSUC_R8 -- IC:',IC(:) WRITE(*,'(A,4I6)') 'W3GSUC_R8 -- JC:',JC(:) WRITE(*,'(A,4E14.6)') 'W3GSUC_R8 -- XC:',XC(:) WRITE(*,'(A,4E14.6)') 'W3GSUC_R8 -- YC:',YC(:) WRITE(*,'(A,4I6)') 'W3GSUC_R8 -- IBC:',IBC(:) WRITE(*,'(A,4I6)') 'W3GSUC_R8 -- JBC:',JBC(:) WRITE(*,'(A,1I6)') 'W3GSUC_R8 -- NS:',NS WRITE(*,'(A,4I6)') 'W3GSUC_R8 -- IB1:',IB1(:) WRITE(*,'(A,4I6)') 'W3GSUC_R8 -- JB1:',JB1(:) WRITE(*,'(A,4I6)') 'W3GSUC_R8 -- IB2:',IB2(:) WRITE(*,'(A,4I6)') 'W3GSUC_R8 -- JB2:',JB2(:) END IF !-------------assign cell to buckets based on overlap DO K=1,NS DO IB=IB1(K),IB2(K) DO JB=JB1(K),JB2(K) PTR%B(JB,IB)%N = PTR%B(JB,IB)%N + 1 IF ( ISTEP .EQ. 2 ) THEN PTR%B(JB,IB)%I(PTR%B(JB,IB)%N) = IC(1) PTR%B(JB,IB)%J(PTR%B(JB,IB)%N) = JC(1) END IF END DO !JB END DO !IB END DO !K END DO !J END DO !I ! !-----END ISTEP_LOOP END DO ISTEP_LOOP ! !-----print debug info IF ( LDBG ) THEN WRITE(*,'(/A,3I6,4E14.6)') 'W3GSUC_R8 - search bucket list:' WRITE(*,'(3A6,4A14)') 'I','J','N','X1','Y1','X2','Y2' DO IB=1,PTR%NBX DO JB=1,PTR%NBY WRITE(*,'(3I6,4E14.6)') IB,JB,PTR%B(JB,IB)%N, & PTR%XMIN+(IB-1)*PTR%DXB,PTR%YMIN+(JB-1)*PTR%DYB, & PTR%XMIN+(IB-0)*PTR%DXB,PTR%YMIN+(JB-0)*PTR%DYB END DO END DO END IF ! ! -------------------------------------------------------------------- / ! 5. Set return parameter ! GSU%PTR => PTR !/ !/ End of W3GSUC_R8 -------------------------------------------------- / !/ END FUNCTION W3GSUC_R8 !/ ------------------------------------------------------------------- / SUBROUTINE W3GSUD(GSU) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 30-Oct-2009 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ ! 1. Purpose : ! ! Destroy grid search utility (GSU) object. ! ! 2. Method : ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! GSU Type I Grid-search-utility object. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! - Check on previous creation of grid-search-utility object. ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ TYPE(T_GSU), INTENT(INOUT) :: GSU !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: IB, JB !/ ! ! -------------------------------------------------------------------- / ! IF ( ASSOCIATED(GSU%PTR) ) THEN ! CALL W3NNSD(GSU%PTR%NNP) ! DO IB=1,GSU%PTR%NBX DO JB=1,GSU%PTR%NBY IF ( GSU%PTR%B(JB,IB)%N .GT. 0 ) THEN DEALLOCATE(GSU%PTR%B(JB,IB)%I) NULLIFY(GSU%PTR%B(JB,IB)%I) DEALLOCATE(GSU%PTR%B(JB,IB)%J) NULLIFY(GSU%PTR%B(JB,IB)%J) END IF END DO END DO DEALLOCATE(GSU%PTR%B) NULLIFY(GSU%PTR%B) ! DEALLOCATE(GSU%PTR) NULLIFY(GSU%PTR) ! END IF !/ !/ End of W3GSUD ----------------------------------------------------- / !/ END SUBROUTINE W3GSUD !/ ------------------------------------------------------------------- / SUBROUTINE W3GSUP(GSU, IUNIT, LFULL) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 06-Dec-2010 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ 01-Dec-2010 : Add output of approx memory usage. ( version 3.14 ) !/ 06-Dec-2010 : Change ICLO to int. Remove JCLO. !/ Add GKIND. ( version 3.14 ) !/ ! 1. Purpose : ! ! Print grid-search-utility (GSU) object to IUNIT. ! ! 2. Method : ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! GSU Type I Grid-search-utility object. ! IUNIT Int. I Optional unit for output. Default is stdout. ! LFULL Log. I Optional logical flag to turn on full-output ! mode. Default is FALSE. When full-output ! is enabled the search bucket cell lists and ! nearest-neighbor point search indices are output. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! - Check on previous creation of grid-search-utility object. ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ TYPE(T_GSU), INTENT(IN) :: GSU INTEGER, OPTIONAL, INTENT(IN) :: IUNIT LOGICAL, OPTIONAL, INTENT(IN) :: LFULL !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER, PARAMETER :: NBYTE_PTR=4 INTEGER, PARAMETER :: NBYTE_INT=4 TYPE(CLASS_GSU), POINTER :: PTR INTEGER :: NDST, I, J, K, L, N, IB, JB, NBYTE !/ ! ! -------------------------------------------------------------------- / ! 1. Test input ! IF ( .NOT.ASSOCIATED(GSU%PTR) ) THEN WRITE(*,'(/1A,1A/)') 'W3GSUP ERROR -- ', & 'grid search utility object not created' CALL EXTCDE (1) END IF IF ( PRESENT(IUNIT) ) THEN NDST = IUNIT ELSE NDST = 6 END IF PTR => GSU%PTR ! ! -------------------------------------------------------------------- / ! 2. Compute approximate memory usage ! NBYTE = (NBYTE_INT+NBYTE_PTR*2)*SIZE(PTR%B) DO IB=1,PTR%NBX DO JB=1,PTR%NBY NBYTE = NBYTE + NBYTE_INT*2*PTR%B(JB,IB)%N END DO END DO ! ! -------------------------------------------------------------------- / ! 3. Output ! WRITE(NDST,'(//80A)') ('-',K=1,80) WRITE(NDST,'(A)') 'Report on grid search utility object' WRITE(NDST,'( 80A)') ('-',K=1,80) WRITE(NDST,'(A,1L2)') 'Grid ijg:',PTR%IJG WRITE(NDST,'(A,1L2)') 'Grid llg:',PTR%LLG WRITE(NDST,'(A,1I2)') 'Grid iclo:',PTR%ICLO WRITE(NDST,'(A,1L2)') 'Grid lclo:',PTR%LCLO WRITE(NDST,'(A,1I2)') 'Grid precision:',PTR%GKIND WRITE(NDST,'(A,2I6)') 'Grid nx,ny:',PTR%NX,PTR%NY IF ( PRESENT(LFULL) ) THEN IF ( LFULL ) THEN WRITE(NDST,'( 80A)') ('-',K=1,80) WRITE(NDST,'(A)') 'Nearest-neighbor point search indices' WRITE(NDST,'( 80A)') ('-',K=1,80) CALL W3NNSP(PTR%NNP,NDST) END IF END IF WRITE(NDST,'( 80A)') ('-',K=1,80) WRITE(NDST,'(A)') 'Bucket-search object' WRITE(NDST,'( 80A)') ('-',K=1,80) WRITE(NDST,'(A,4E14.6)') 'Spatial grid search domain: ', & PTR%XMIN,PTR%YMIN,PTR%XMAX,PTR%YMAX WRITE(NDST,'(A,2I6)') 'nbx,nby:',PTR%NBX,PTR%NBY WRITE(NDST,'(A,2E14.6)') 'dxb,dyb:',PTR%DXB,PTR%DYB WRITE(NDST,'(A,1F10.1)') 'Approximate memory usage (MB):', & REAL(NBYTE)/2**20 IF ( PRESENT(LFULL) ) THEN IF ( LFULL ) THEN WRITE(NDST,'( 80A)') ('-',K=1,80) WRITE(NDST,'(A)') 'Search bucket bounds:' WRITE(NDST,'( 80A)') ('-',K=1,80) WRITE(NDST,'(2A4,4A14)') 'IB','JB','X1','Y1','X2','Y2' DO IB=1,PTR%NBX DO JB=1,PTR%NBY WRITE(*,'(2I4,4E14.6)') IB,JB, & PTR%XMIN+(IB-1)*PTR%DXB,PTR%YMIN+(JB-1)*PTR%DYB, & PTR%XMIN+(IB )*PTR%DXB,PTR%YMIN+(JB )*PTR%DYB END DO END DO WRITE(NDST,'( 80A)') ('-',K=1,80) WRITE(NDST,'(A)') 'Number of cells in each search bucket:' WRITE(NDST,'( 80A)') ('-',K=1,80) DO JB=PTR%NBY,1,-1 WRITE(NDST,'(500I4)') (PTR%B(JB,IB)%N,IB=1,PTR%NBX) END DO WRITE(NDST,'( 80A)') ('-',K=1,80) WRITE(NDST,'(A)') 'Search bucket cell lists:' WRITE(NDST,'( 80A)') ('-',K=1,80) WRITE(NDST,'(3A4,A)') 'IB','JB','NC',': ( IC, JC), ...' DO JB=1,PTR%NBY DO IB=1,PTR%NBX WRITE(NDST,'(3I4,A,500(A,I3,A,I3,A))') IB,JB, & PTR%B(JB,IB)%N, ': ', & ( '(',PTR%B(JB,IB)%I(K),',',PTR%B(JB,IB)%J(K),') ', & K=1,PTR%B(JB,IB)%N ) END DO END DO END IF !LFULL END IF !PRESENT(LFULL) WRITE(NDST,'( 80A)') ('-',K=1,80) WRITE(NDST,'( 80A)') ('-',K=1,80) !/ !/ End of W3GSUP ----------------------------------------------------- / !/ END SUBROUTINE W3GSUP !/ ------------------------------------------------------------------- / FUNCTION W3GFCL_R4(GSU, XT, YT, IS, JS, XS, YS, POLE, DEBUG) & RESULT(INGRID) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 06-Dec-2010 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ 12-Nov-2010 : Implement r4 & r8 interfaces. ( version 3.14 ) !/ 01-Dec-2010 : Remove search using nnbr buckets. ( version 3.14 ) !/ 06-Dec-2010 : Remove restriction on longitude range. Change ICLO !/ to integer and remove JCLO. Implement support for !/ r4 and r8 source grids. ( version 3.14 ) !/ ! 1. Purpose : ! ! Find cell in grid, associated with the input grid-search-utility ! object (GSU), that encloses the target point (xt,yt). ! Single precision interface. ! ! 2. Method : ! ! 3. Parameters : ! ! Return parameter ! ---------------------------------------------------------------- ! INGRID Log. O Logical flag indicating if target point lies ! within the source grid domain. ! ---------------------------------------------------------------- ! ! Parameter list ! ---------------------------------------------------------------- ! GSU Type I Grid-search-utility object. ! XT Real I X-coordinate of target point. ! YT Real I Y-coordinate of target point. ! IS,JS I.A. O (I,J) indices of vertices of enclosing grid cell. ! XS,YS R.A. O (X,Y) coord. of vertices of enclosing grid cell. ! POLE Log. O Optional logical flag to indicate whether or not ! the enclosing grid cell includes a pole. ! DEBUG Log. I Optional logical flag to turn on debug mode. ! Default is FALSE. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! - Check on previous creation of grid-search-utility object. ! ! 7. Remarks : ! ! - The target point coordinates may be modified by this routine. ! - The target point longitude will be shifted to the source grid ! longitudinal range. ! - If enclosing cell includes a branch cut, then the coordinates of ! of the cell vertices AND the target point will be adjusted so ! that the branch cut is shifted 180 degrees. ! ! 8. Structure : ! ! ----------------------------------------------------------------- ! 1. Test input ! 2. Initialize search ! 3. Search for enclosing cell in central and nearest nbr buckets ! ----------------------------------------------------------------- ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Return parameter !/ LOGICAL :: INGRID !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ TYPE(T_GSU), INTENT(IN) :: GSU REAL(4), INTENT(INOUT) :: XT REAL(4), INTENT(INOUT) :: YT INTEGER, INTENT(INOUT) :: IS(4), JS(4) REAL(4), INTENT(INOUT) :: XS(4), YS(4) LOGICAL, INTENT(OUT),OPTIONAL :: POLE LOGICAL, INTENT(IN), OPTIONAL :: DEBUG !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ LOGICAL :: LDBG, LPLC INTEGER :: I, J, K, L, N, IB, JB LOGICAL :: IJG, LLG, LCLO, L360 INTEGER :: ICLO, GKIND INTEGER :: NX, NY REAL(4), POINTER :: XG4(:,:), YG4(:,:) REAL(8), POINTER :: XG8(:,:), YG8(:,:) INTEGER :: NBX, NBY REAL(8) :: DXB, DYB, XMIN, XMAX, YMIN, YMAX TYPE(T_BKT), POINTER :: B(:,:) !/ ! ! -------------------------------------------------------------------- / ! 1. Test input ! IF ( .NOT.ASSOCIATED(GSU%PTR) ) THEN WRITE(*,'(/2A/)') 'W3GFCL_R4 ERROR -- ', & 'grid search utility object not created' CALL EXTCDE (1) END IF ! ! -------------------------------------------------------------------- / ! 2. Initialize search ! IF ( PRESENT(DEBUG) ) THEN LDBG = DEBUG ELSE LDBG = .FALSE. END IF ! ! Local pointers to grid search utility object data IJG = GSU%PTR%IJG LLG = GSU%PTR%LLG ICLO = GSU%PTR%ICLO LCLO = GSU%PTR%LCLO L360 = GSU%PTR%L360 GKIND = GSU%PTR%GKIND NX = GSU%PTR%NX; NY = GSU%PTR%NY; IF ( GKIND.EQ.4 ) THEN XG4 => GSU%PTR%XG4; YG4 => GSU%PTR%YG4; ELSE XG8 => GSU%PTR%XG8; YG8 => GSU%PTR%YG8; END IF NBX = GSU%PTR%NBX; NBY = GSU%PTR%NBY; DXB = GSU%PTR%DXB; DYB = GSU%PTR%DYB; XMIN = GSU%PTR%XMIN; YMIN = GSU%PTR%YMIN; XMAX = GSU%PTR%XMAX; YMAX = GSU%PTR%YMAX; B => GSU%PTR%B ! INGRID = .FALSE. ! ! Shift target to appropriate longitude range IF ( LLG ) THEN XT = MOD(XT,REAL(D360,4)) IF ( LCLO .OR. L360 ) THEN IF ( XT.LT.ZERO ) XT = XT + D360 ELSE IF ( XT.GT.D180 ) XT = XT - D360 END IF END IF IF ( LDBG ) WRITE(*,'(/A,2E14.6)') 'W3GFCL_R4 - TARGET POINT:',XT,YT ! ! Target point must lie within search domain IF ( XT.LT.XMIN .OR. XT.GT.XMAX .OR. & YT.LT.YMIN .OR. YT.GT.YMAX ) THEN IF ( LDBG ) WRITE(*,'(A)') & 'W3GFCL_R4 - TARGET POINT OUTSIDE SEARCH DOMAIN' RETURN END IF ! ! Search bucket that contains the target point. IB = INT((XT-XMIN)/DXB)+1; IF ( .NOT.LCLO ) IB = MIN(IB,NBX); JB = INT((YT-YMIN)/DYB)+1; JB = MIN(JB,NBY); ! ! -------------------------------------------------------------------- / ! 3. Search for enclosing cell in bucket ! IF ( LDBG ) & WRITE(*,'(A,3I6,4E14.6)') & 'W3GFCL_R4 - BUCKET SEARCH:',IB,JB,B(JB,IB)%N, & XMIN+(IB-1)*DXB,YMIN+(JB-1)*DYB,XMIN+IB*DXB,YMIN+JB*DYB CELL_LOOP: DO K=1,B(JB,IB)%N !---------setup cell corner indices IS(1) = B(JB,IB)%I(K) ; JS(1) = B(JB,IB)%J(K) ; IS(2) = B(JB,IB)%I(K)+1; JS(2) = B(JB,IB)%J(K) ; IS(3) = B(JB,IB)%I(K)+1; JS(3) = B(JB,IB)%J(K)+1; IS(4) = B(JB,IB)%I(K) ; JS(4) = B(JB,IB)%J(K)+1; !---------setup cell corner coordinates and adjust for periodicity DO L=1,4 IF ( ICLO.NE.ICLO_NONE ) THEN IF ( IS(L) .LT. 1 ) IS(L) = IS(L) + NX IF ( IS(L) .GT. NX ) IS(L) = IS(L) - NX END IF IF ( ICLO.EQ.ICLO_TRPL ) THEN IF ( JS(L) .GT. NY ) THEN JS(L) = NY IS(L) = MOD(NX-IS(L),NX) + 1 END IF END IF IF ( IJG ) THEN IF ( GKIND.EQ.4 ) THEN XS(L) = XG4(IS(L),JS(L)); YS(L) = YG4(IS(L),JS(L)); ELSE XS(L) = XG8(IS(L),JS(L)); YS(L) = YG8(IS(L),JS(L)); END IF ELSE IF ( GKIND.EQ.4 ) THEN XS(L) = XG4(JS(L),IS(L)); YS(L) = YG4(JS(L),IS(L)); ELSE XS(L) = XG8(JS(L),IS(L)); YS(L) = YG8(JS(L),IS(L)); END IF END IF IF ( LLG ) THEN XS(L) = MOD(XS(L),REAL(D360,4)) IF ( LCLO .OR. L360 ) THEN IF ( XS(L).LT.ZERO ) XS(L) = XS(L) + D360 ELSE IF ( XS(L).GT.D180 ) XS(L) = XS(L) - D360 END IF END IF END DO !L IF ( LDBG ) & WRITE(*,'(A,3I6,4(/A,1I1,A,2I6,2E14.6))') & 'W3GFCL_R4 - CHECK CELL:',IB,JB,K, & (' CORNER(',L,'):',IS(L),JS(L),XS(L),YS(L),L=1,4) !---------check if point is enclosed in cell defined by xs(1:4) & ys(1:4) INGRID = W3CKCL(LLG,XT,YT,4,XS,YS,LPLC,LDBG) IF ( LDBG ) WRITE(*,'(A,1L2)')'W3GFCL_R4 - INGRID:',INGRID IF ( INGRID ) THEN !-------------exit search IF ( LDBG ) & WRITE(*,'(A,3I6,4(2I6))') & 'W3GFCL_R4 - ENCLOSING CELL:',IB,JB,K,(IS(L),JS(L),L=1,4) IF ( PRESENT(POLE) ) POLE = LPLC EXIT CELL_LOOP END IF !point in cell END DO CELL_LOOP !/ !/ End of W3GFCL_R4 -------------------------------------------------- / !/ END FUNCTION W3GFCL_R4 !/ ------------------------------------------------------------------- / FUNCTION W3GFCL_R8(GSU, XT, YT, IS, JS, XS, YS, POLE, DEBUG) & RESULT(INGRID) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 06-Dec-2010 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ 12-Nov-2010 : Implement r4 & r8 interfaces. ( version 3.14 ) !/ 01-Dec-2010 : Remove search using nnbr buckets. ( version 3.14 ) !/ 06-Dec-2010 : Remove restriction on longitude range. Change ICLO !/ to integer and remove JCLO. Implement support for !/ r4 and r8 source grids. ( version 3.14 ) !/ ! 1. Purpose : ! ! Find cell in grid, associated with the input grid-search-utility ! object (GSU), that encloses the target point (xt,yt). ! Double precision interface. ! ! 2. Method : ! ! 3. Parameters : ! ! Return parameter ! ---------------------------------------------------------------- ! INGRID Log. O Logical flag indicating if target point lies ! within the source grid domain. ! ---------------------------------------------------------------- ! ! Parameter list ! ---------------------------------------------------------------- ! GSU Type I Grid-search-utility object. ! XT Real I X-coordinate of target point. ! YT Real I Y-coordinate of target point. ! IS,JS I.A. O (I,J) indices of vertices of enclosing grid cell. ! XS,YS R.A. O (X,Y) coord. of vertices of enclosing grid cell. ! POLE Log. O Optional logical flag to indicate whether or not ! the enclosing grid cell includes a pole. ! DEBUG Log. I Optional logical flag to turn on debug mode. ! Default is FALSE. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! - Check on previous creation of grid-search-utility object. ! ! 7. Remarks : ! ! - The target point coordinates may be modified by this routine. ! - The target point longitude will be shifted to the source grid ! longitudinal range. ! - If enclosing cell includes a branch cut, then the coordinates of ! of the cell vertices AND the target point will be adjusted so ! that the branch cut is shifted 180 degrees. ! ! 8. Structure : ! ! ----------------------------------------------------------------- ! 1. Test input ! 2. Initialize search ! 3. Search for enclosing cell in central and nearest nbr buckets ! ----------------------------------------------------------------- ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Return parameter !/ LOGICAL :: INGRID !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ TYPE(T_GSU), INTENT(IN) :: GSU REAL(8), INTENT(INOUT) :: XT REAL(8), INTENT(INOUT) :: YT INTEGER, INTENT(INOUT) :: IS(4), JS(4) REAL(8), INTENT(INOUT) :: XS(4), YS(4) LOGICAL, INTENT(OUT),OPTIONAL :: POLE LOGICAL, INTENT(IN), OPTIONAL :: DEBUG !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ LOGICAL :: LDBG, LPLC INTEGER :: I, J, K, L, N, IB, JB LOGICAL :: IJG, LLG, LCLO, L360 INTEGER :: ICLO, GKIND INTEGER :: NX, NY REAL(4), POINTER :: XG4(:,:), YG4(:,:) REAL(8), POINTER :: XG8(:,:), YG8(:,:) INTEGER :: NBX, NBY REAL(8) :: DXB, DYB, XMIN, XMAX, YMIN, YMAX TYPE(T_BKT), POINTER :: B(:,:) !/ ! ! -------------------------------------------------------------------- / ! 1. Test input ! IF ( .NOT.ASSOCIATED(GSU%PTR) ) THEN WRITE(*,'(/2A/)') 'W3GFCL_R8 ERROR -- ', & 'grid search utility object not created' CALL EXTCDE (1) END IF ! ! -------------------------------------------------------------------- / ! 2. Initialize search ! IF ( PRESENT(DEBUG) ) THEN LDBG = DEBUG ELSE LDBG = .FALSE. END IF ! ! Local pointers to grid search utility object data IJG = GSU%PTR%IJG LLG = GSU%PTR%LLG ICLO = GSU%PTR%ICLO LCLO = GSU%PTR%LCLO L360 = GSU%PTR%L360 GKIND = GSU%PTR%GKIND NX = GSU%PTR%NX; NY = GSU%PTR%NY; IF ( GKIND.EQ.4 ) THEN XG4 => GSU%PTR%XG4; YG4 => GSU%PTR%YG4; ELSE XG8 => GSU%PTR%XG8; YG8 => GSU%PTR%YG8; END IF NBX = GSU%PTR%NBX; NBY = GSU%PTR%NBY; DXB = GSU%PTR%DXB; DYB = GSU%PTR%DYB; XMIN = GSU%PTR%XMIN; YMIN = GSU%PTR%YMIN; XMAX = GSU%PTR%XMAX; YMAX = GSU%PTR%YMAX; B => GSU%PTR%B ! INGRID = .FALSE. ! ! Shift target to appropriate longitude range IF ( LLG ) THEN XT = MOD(XT,REAL(D360,8)) IF ( LCLO .OR. L360 ) THEN IF ( XT.LT.ZERO ) XT = XT + D360 ELSE IF ( XT.GT.D180 ) XT = XT - D360 END IF END IF IF ( LDBG ) WRITE(*,'(/A,2E14.6)') 'W3GFCL_R8 - TARGET POINT:',XT,YT ! ! Target point must lie within search domain IF ( XT.LT.XMIN .OR. XT.GT.XMAX .OR. & YT.LT.YMIN .OR. YT.GT.YMAX ) THEN IF ( LDBG ) WRITE(*,'(A)') & 'W3GFCL_R8 - TARGET POINT OUTSIDE SEARCH DOMAIN' RETURN END IF ! ! Search bucket that contains the target point. IB = INT((XT-XMIN)/DXB)+1; IF ( .NOT.LCLO ) IB = MIN(IB,NBX); JB = INT((YT-YMIN)/DYB)+1; JB = MIN(JB,NBY); ! ! -------------------------------------------------------------------- / ! 3. Search for enclosing cell in bucket ! IF ( LDBG ) & WRITE(*,'(A,3I6,4E14.6)') & 'W3GFCL_R8 - BUCKET SEARCH:',IB,JB,B(JB,IB)%N, & XMIN+(IB-1)*DXB,YMIN+(JB-1)*DYB,XMIN+IB*DXB,YMIN+JB*DYB CELL_LOOP: DO K=1,B(JB,IB)%N !---------setup cell corner indices IS(1) = B(JB,IB)%I(K) ; JS(1) = B(JB,IB)%J(K) ; IS(2) = B(JB,IB)%I(K)+1; JS(2) = B(JB,IB)%J(K) ; IS(3) = B(JB,IB)%I(K)+1; JS(3) = B(JB,IB)%J(K)+1; IS(4) = B(JB,IB)%I(K) ; JS(4) = B(JB,IB)%J(K)+1; !---------setup cell corner coordinates and adjust for periodicity DO L=1,4 IF ( ICLO.NE.ICLO_NONE ) THEN IF ( IS(L) .LT. 1 ) IS(L) = IS(L) + NX IF ( IS(L) .GT. NX ) IS(L) = IS(L) - NX END IF IF ( ICLO.EQ.ICLO_TRPL ) THEN IF ( JS(L) .GT. NY ) THEN JS(L) = NY IS(L) = MOD(NX-IS(L),NX) + 1 END IF END IF IF ( IJG ) THEN IF ( GKIND.EQ.4 ) THEN XS(L) = XG4(IS(L),JS(L)); YS(L) = YG4(IS(L),JS(L)); ELSE XS(L) = XG8(IS(L),JS(L)); YS(L) = YG8(IS(L),JS(L)); END IF ELSE IF ( GKIND.EQ.4 ) THEN XS(L) = XG4(JS(L),IS(L)); YS(L) = YG4(JS(L),IS(L)); ELSE XS(L) = XG8(JS(L),IS(L)); YS(L) = YG8(JS(L),IS(L)); END IF END IF IF ( LLG ) THEN XS(L) = MOD(XS(L),REAL(D360,8)) IF ( LCLO .OR. L360 ) THEN IF ( XS(L).LT.ZERO ) XS(L) = XS(L) + D360 ELSE IF ( XS(L).GT.D180 ) XS(L) = XS(L) - D360 END IF END IF END DO !L IF ( LDBG ) & WRITE(*,'(A,3I6,4(/A,1I1,A,2I6,2E14.6))') & 'W3GFCL_R8 - CHECK CELL:',IB,JB,K, & (' CORNER(',L,'):',IS(L),JS(L),XS(L),YS(L),L=1,4) !---------check if point is enclosed in cell defined by xs(1:4) & ys(1:4) INGRID = W3CKCL(LLG,XT,YT,4,XS,YS,LPLC,LDBG) IF ( LDBG ) WRITE(*,'(A,1L2)')'W3GFCL_R8 - INGRID:',INGRID IF ( INGRID ) THEN !-------------exit search IF ( LDBG ) & WRITE(*,'(A,3I6,4(2I6))') & 'W3GFCL_R8 - ENCLOSING CELL:',IB,JB,K,(IS(L),JS(L),L=1,4) IF ( PRESENT(POLE) ) POLE = LPLC EXIT CELL_LOOP END IF !point in cell END DO CELL_LOOP !/ !/ End of W3GFCL_R8 -------------------------------------------------- / !/ END FUNCTION W3GFCL_R8 !/ ------------------------------------------------------------------- / FUNCTION W3GFCD_R4(GSU, XT, YT, IS, JS, XS, YS, POLE, DEBUG) & RESULT(INGRID) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 06-Dec-2010 | !/ +-----------------------------------+ !/ !/ 01-Dec-2010 : Origination. ( version 3.14 ) !/ 06-Dec-2010 : Remove restriction on longitude range. Change ICLO !/ to integer and remove JCLO. Implement support for !/ r4 and r8 source grids. ( version 3.14 ) !/ ! 1. Purpose : ! ! Find cell in grid, associated with the input grid-search-utility ! object (GSU), that encloses the target point (xt,yt), using direct ! grid search (i.e., no bucket search). ! Single precision interface. ! ! 2. Method : ! ! 3. Parameters : ! ! Return parameter ! ---------------------------------------------------------------- ! INGRID Log. O Logical flag indicating if target point lies ! within the source grid domain. ! ---------------------------------------------------------------- ! ! Parameter list ! ---------------------------------------------------------------- ! GSU Type I Grid-search-utility object. ! XT Real I X-coordinate of target point. ! YT Real I Y-coordinate of target point. ! IS,JS I.A. O (I,J) indices of vertices of enclosing grid cell. ! XS,YS R.A. O (X,Y) coord. of vertices of enclosing grid cell. ! POLE Log. O Optional logical flag to indicate whether or not ! the enclosing grid cell includes a pole. ! DEBUG Log. I Optional logical flag to turn on debug mode. ! Default is FALSE. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! - Check on previous creation of grid-search-utility object. ! ! 7. Remarks : ! ! - The target point coordinates may be modified by this routine. ! - The target point longitude will be shifted to the source grid ! longitudinal range. ! - If enclosing cell includes a branch cut, then the coordinates of ! of the cell vertices AND the target point will be adjusted so ! that the branch cut is shifted 180 degrees. ! ! 8. Structure : ! ! ----------------------------------------------------------------- ! 1. Test input ! 2. Initialize search ! 3. Search for enclosing cell ! ----------------------------------------------------------------- ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Return parameter !/ LOGICAL :: INGRID !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ TYPE(T_GSU), INTENT(IN) :: GSU REAL(4), INTENT(INOUT) :: XT REAL(4), INTENT(INOUT) :: YT INTEGER, INTENT(INOUT) :: IS(4), JS(4) REAL(4), INTENT(INOUT) :: XS(4), YS(4) LOGICAL, INTENT(OUT),OPTIONAL :: POLE LOGICAL, INTENT(IN), OPTIONAL :: DEBUG !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ LOGICAL :: LDBG, LPLC INTEGER :: I, J, K, L, N, IB, JB LOGICAL :: IJG, LLG, LCLO, L360 INTEGER :: ICLO, GKIND INTEGER :: NX, NY, NXC, NYC REAL(4), POINTER :: XG4(:,:), YG4(:,:) REAL(8), POINTER :: XG8(:,:), YG8(:,:) !/ ! ! -------------------------------------------------------------------- / ! 1. Test input ! IF ( .NOT.ASSOCIATED(GSU%PTR) ) THEN WRITE(*,'(/2A/)') 'W3GFCD_R4 ERROR -- ', & 'grid search utility object not created' CALL EXTCDE (1) END IF ! ! -------------------------------------------------------------------- / ! 2. Initialize search ! IF ( PRESENT(DEBUG) ) THEN LDBG = DEBUG ELSE LDBG = .FALSE. END IF ! ! Local pointers to grid search utility object data IJG = GSU%PTR%IJG LLG = GSU%PTR%LLG ICLO = GSU%PTR%ICLO LCLO = GSU%PTR%LCLO L360 = GSU%PTR%L360 GKIND = GSU%PTR%GKIND NX = GSU%PTR%NX; NY = GSU%PTR%NY; IF ( GKIND.EQ.4 ) THEN XG4 => GSU%PTR%XG4; YG4 => GSU%PTR%YG4; ELSE XG8 => GSU%PTR%XG8; YG8 => GSU%PTR%YG8; END IF ! INGRID = .FALSE. ! ! Shift target to appropriate longitude range IF ( LLG ) THEN XT = MOD(XT,REAL(D360,4)) IF ( LCLO .OR. L360 ) THEN IF ( XT.LT.ZERO ) XT = XT + D360 ELSE IF ( XT.GT.D180 ) XT = XT - D360 END IF END IF IF ( LDBG ) WRITE(*,'(/A,2E14.6)') 'W3GFCD_R4 - TARGET POINT:',XT,YT !-----number of cells SELECT CASE ( ICLO ) CASE ( ICLO_NONE ) NXC = NX-1; NYC = NY-1; CASE ( ICLO_SMPL ) NXC = NX; NYC = NY-1; CASE ( ICLO_TRPL ) NXC = NX; NYC = NY; END SELECT ! ! -------------------------------------------------------------------- / ! 3. Search for enclosing cell ! CELL_LOOP: DO I=1,NXC DO J=1,NYC !-------------create list of cell vertices IS(1) = I ; JS(1) = J ; IS(2) = I+1; JS(2) = J ; IS(3) = I+1; JS(3) = J+1; IS(4) = I ; JS(4) = J+1; !-------------setup cell corner coordinates and adjust for periodicity DO L=1,4 IF ( ICLO.NE.ICLO_NONE ) THEN IF ( IS(L) .LT. 1 ) IS(L) = IS(L) + NX IF ( IS(L) .GT. NX ) IS(L) = IS(L) - NX END IF IF ( ICLO.EQ.ICLO_TRPL ) THEN IF ( JS(L) .GT. NY ) THEN JS(L) = NY IS(L) = MOD(NX-IS(L),NX) + 1 END IF END IF IF ( IJG ) THEN IF ( GKIND.EQ.4 ) THEN XS(L) = XG4(IS(L),JS(L)); YS(L) = YG4(IS(L),JS(L)); ELSE XS(L) = XG8(IS(L),JS(L)); YS(L) = YG8(IS(L),JS(L)); END IF ELSE IF ( GKIND.EQ.4 ) THEN XS(L) = XG4(JS(L),IS(L)); YS(L) = YG4(JS(L),IS(L)); ELSE XS(L) = XG8(JS(L),IS(L)); YS(L) = YG8(JS(L),IS(L)); END IF END IF IF ( LLG ) THEN XS(L) = MOD(XS(L),REAL(D360,4)) IF ( LCLO .OR. L360 ) THEN IF ( XS(L).LT.ZERO ) XS(L) = XS(L) + D360 ELSE IF ( XS(L).GT.D180 ) XS(L) = XS(L) - D360 END IF END IF END DO !L IF ( LDBG ) & WRITE(*,'(A,4(/A,1I1,A,2I6,2E14.6))') & 'W3GFCD_R4 - CHECK CELL:', & (' CORNER(',L,'):',IS(L),JS(L),XS(L),YS(L),L=1,4) !-------------check if point is enclosed in cell defined by xs(1:4) & ys(1:4) INGRID = W3CKCL(LLG,XT,YT,4,XS,YS,LPLC,LDBG) IF ( LDBG ) WRITE(*,'(A,1L2)')'W3GFCD_R4 - INGRID:',INGRID IF ( INGRID ) THEN !-----------------exit search IF ( LDBG ) & WRITE(*,'(A,4(2I6))') & 'W3GFCD_R4 - ENCLOSING CELL:',(IS(L),JS(L),L=1,4) IF ( PRESENT(POLE) ) POLE = LPLC EXIT CELL_LOOP END IF !point in cell END DO !J END DO CELL_LOOP !/ !/ End of W3GFCD_R4--------------------------------------------------- / !/ END FUNCTION W3GFCD_R4 !/ ------------------------------------------------------------------- / FUNCTION W3GFCD_R8(GSU, XT, YT, IS, JS, XS, YS, POLE, DEBUG) & RESULT(INGRID) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 06-Dec-2010 | !/ +-----------------------------------+ !/ !/ 01-Dec-2010 : Origination. ( version 3.14 ) !/ 06-Dec-2010 : Remove restriction on longitude range. Change ICLO !/ to integer and remove JCLO. Implement support for !/ r4 and r8 source grids. ( version 3.14 ) !/ ! 1. Purpose : ! ! Find cell in grid, associated with the input grid-search-utility ! object (GSU), that encloses the target point (xt,yt), using direct ! grid search (i.e., no bucket search). ! Double precision interface. ! ! 2. Method : ! ! 3. Parameters : ! ! Return parameter ! ---------------------------------------------------------------- ! INGRID Log. O Logical flag indicating if target point lies ! within the source grid domain. ! ---------------------------------------------------------------- ! ! Parameter list ! ---------------------------------------------------------------- ! GSU Type I Grid-search-utility object. ! XT Real I X-coordinate of target point. ! YT Real I Y-coordinate of target point. ! IS,JS I.A. O (I,J) indices of vertices of enclosing grid cell. ! XS,YS R.A. O (X,Y) coord. of vertices of enclosing grid cell. ! POLE Log. O Optional logical flag to indicate whether or not ! the enclosing grid cell includes a pole. ! DEBUG Log. I Optional logical flag to turn on debug mode. ! Default is FALSE. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! - Check on previous creation of grid-search-utility object. ! ! 7. Remarks : ! ! - The target point coordinates may be modified by this routine. ! - The target point longitude will be shifted to the source grid ! longitudinal range. ! - If enclosing cell includes a branch cut, then the coordinates of ! of the cell vertices AND the target point will be adjusted so ! that the branch cut is shifted 180 degrees. ! ! 8. Structure : ! ! ----------------------------------------------------------------- ! 1. Test input ! 2. Initialize search ! 3. Search for enclosing cell ! ----------------------------------------------------------------- ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Return parameter !/ LOGICAL :: INGRID !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ TYPE(T_GSU), INTENT(IN) :: GSU REAL(8), INTENT(INOUT) :: XT REAL(8), INTENT(INOUT) :: YT INTEGER, INTENT(INOUT) :: IS(4), JS(4) REAL(8), INTENT(INOUT) :: XS(4), YS(4) LOGICAL, INTENT(OUT),OPTIONAL :: POLE LOGICAL, INTENT(IN), OPTIONAL :: DEBUG !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ LOGICAL :: LDBG, LPLC INTEGER :: I, J, K, L, N, IB, JB LOGICAL :: IJG, LLG, LCLO, L360 INTEGER :: ICLO, GKIND INTEGER :: NX, NY, NXC, NYC REAL(4), POINTER :: XG4(:,:), YG4(:,:) REAL(8), POINTER :: XG8(:,:), YG8(:,:) !/ ! ! -------------------------------------------------------------------- / ! 1. Test input ! IF ( .NOT.ASSOCIATED(GSU%PTR) ) THEN WRITE(*,'(/2A/)') 'W3GFCD_R8 ERROR -- ', & 'grid search utility object not created' CALL EXTCDE (1) END IF ! ! -------------------------------------------------------------------- / ! 2. Initialize search ! IF ( PRESENT(DEBUG) ) THEN LDBG = DEBUG ELSE LDBG = .FALSE. END IF ! ! Local pointers to grid search utility object data IJG = GSU%PTR%IJG LLG = GSU%PTR%LLG ICLO = GSU%PTR%ICLO LCLO = GSU%PTR%LCLO L360 = GSU%PTR%L360 GKIND = GSU%PTR%GKIND NX = GSU%PTR%NX; NY = GSU%PTR%NY; IF ( GKIND.EQ.4 ) THEN XG4 => GSU%PTR%XG4; YG4 => GSU%PTR%YG4; ELSE XG8 => GSU%PTR%XG8; YG8 => GSU%PTR%YG8; END IF ! INGRID = .FALSE. ! ! Shift target to appropriate longitude range IF ( LLG ) THEN XT = MOD(XT,REAL(D360,8)) IF ( LCLO .OR. L360 ) THEN IF ( XT.LT.ZERO ) XT = XT + D360 ELSE IF ( XT.GT.D180 ) XT = XT - D360 END IF END IF IF ( LDBG ) WRITE(*,'(/A,2E14.6)') 'W3GFCD_R8 - TARGET POINT:',XT,YT !-----number of cells SELECT CASE ( ICLO ) CASE ( ICLO_NONE ) NXC = NX-1; NYC = NY-1; CASE ( ICLO_SMPL ) NXC = NX; NYC = NY-1; CASE ( ICLO_TRPL ) NXC = NX; NYC = NY; END SELECT ! ! -------------------------------------------------------------------- / ! 3. Search for enclosing cell ! CELL_LOOP: DO I=1,NXC DO J=1,NYC !-------------create list of cell vertices IS(1) = I ; JS(1) = J ; IS(2) = I+1; JS(2) = J ; IS(3) = I+1; JS(3) = J+1; IS(4) = I ; JS(4) = J+1; !-------------setup cell corner coordinates and adjust for periodicity DO L=1,4 IF ( ICLO.NE.ICLO_NONE ) THEN IF ( IS(L) .LT. 1 ) IS(L) = IS(L) + NX IF ( IS(L) .GT. NX ) IS(L) = IS(L) - NX END IF IF ( ICLO.EQ.ICLO_TRPL ) THEN IF ( JS(L) .GT. NY ) THEN JS(L) = NY IS(L) = MOD(NX-IS(L),NX) + 1 END IF END IF IF ( IJG ) THEN IF ( GKIND.EQ.4 ) THEN XS(L) = XG4(IS(L),JS(L)); YS(L) = YG4(IS(L),JS(L)); ELSE XS(L) = XG8(IS(L),JS(L)); YS(L) = YG8(IS(L),JS(L)); END IF ELSE IF ( GKIND.EQ.4 ) THEN XS(L) = XG4(JS(L),IS(L)); YS(L) = YG4(JS(L),IS(L)); ELSE XS(L) = XG8(JS(L),IS(L)); YS(L) = YG8(JS(L),IS(L)); END IF END IF IF ( LLG ) THEN XS(L) = MOD(XS(L),REAL(D360,8)) IF ( LCLO .OR. L360 ) THEN IF ( XS(L).LT.ZERO ) XS(L) = XS(L) + D360 ELSE IF ( XS(L).GT.D180 ) XS(L) = XS(L) - D360 END IF END IF END DO !L IF ( LDBG ) & WRITE(*,'(A,4(/A,1I1,A,2I6,2E14.6))') & 'W3GFCD_R8 - CHECK CELL:', & (' CORNER(',L,'):',IS(L),JS(L),XS(L),YS(L),L=1,4) !-------------check if point is enclosed in cell defined by xs(1:4) & ys(1:4) INGRID = W3CKCL(LLG,XT,YT,4,XS,YS,LPLC,LDBG) IF ( LDBG ) WRITE(*,'(A,1L2)')'W3GFCD_R8 - INGRID:',INGRID IF ( INGRID ) THEN !-----------------exit search IF ( LDBG ) & WRITE(*,'(A,4(2I6))') & 'W3GFCD_R8 - ENCLOSING CELL:',(IS(L),JS(L),L=1,4) IF ( PRESENT(POLE) ) POLE = LPLC EXIT CELL_LOOP END IF !point in cell END DO !J END DO CELL_LOOP !/ !/ End of W3GFCD_R8--------------------------------------------------- / !/ END FUNCTION W3GFCD_R8 !/ ------------------------------------------------------------------- / FUNCTION W3GFPT_R4(GSU, XTIN, YTIN, IX, IY, DEBUG) & RESULT(INGRID) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 01-Dec-2010 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ 12-Nov-2010 : Implement r4 & r8 interfaces. ( version 3.14 ) !/ 01-Dec-2010 : Some cleanup. ( version 3.14 ) !/ ! 1. Purpose : ! ! Find point in grid, associated with the input grid-search-utility ! object (GSU), that is closest to the target point (xtin,ytin). ! Single precision interface. ! ! 2. Method : ! ! 3. Parameters : ! ! Return parameter ! ---------------------------------------------------------------- ! INGRID Log. O Logical flag indicating if target point lies ! within the source grid domain. ! ---------------------------------------------------------------- ! ! Parameter list ! ---------------------------------------------------------------- ! GSU Type I Grid-search-utility object. ! XTIN Real I X-coordinate of target point. ! YTIN Real I Y-coordinate of target point. ! IX,JX I.A. O (I,J) indices of nearest grid point. ! DEBUG Log. I Optional logical flag to turn on debug mode. ! Default is FALSE. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! - Check on previous initialization of grid search utility object. ! ! 7. Remarks : ! ! 8. Structure : ! ! ----------------------------------------------------------------- ! 1. Test input ! 2. Initialize search ! 3. Find enclosing cell and compute closest point ! ----------------------------------------------------------------- ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Return parameter !/ LOGICAL :: INGRID !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ TYPE(T_GSU), INTENT(IN) :: GSU REAL(4), INTENT(IN) :: XTIN REAL(4), INTENT(IN) :: YTIN INTEGER, INTENT(OUT) :: IX, IY LOGICAL, INTENT(IN), OPTIONAL :: DEBUG !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ REAL(8), PARAMETER :: BIG = 1D16 LOGICAL :: LDBG INTEGER :: I, J, K, L REAL(4) :: XT, YT INTEGER :: IS(4), JS(4) REAL(4) :: XS(4), YS(4) REAL(4) :: DD, DMIN LOGICAL :: IJG, LLG !/ ! ! -------------------------------------------------------------------- / ! 1. Test input ! IF ( .NOT.ASSOCIATED(GSU%PTR) ) THEN WRITE(*,'(/2A/)') 'W3GFPT_R4 ERROR -- ', & 'grid search utility object not created' CALL EXTCDE (1) END IF ! ! -------------------------------------------------------------------- / ! 2. Initialize search ! IF ( PRESENT(DEBUG) ) THEN LDBG = DEBUG ELSE LDBG = .FALSE. END IF ! ! Local pointers to grid search utility object data IJG = GSU%PTR%IJG LLG = GSU%PTR%LLG ! INGRID = .FALSE. ! XT = XTIN; YT = YTIN; IF ( LDBG ) WRITE(*,'(/A,2E14.6)') 'W3GFPT_R4 - TARGET POINT:',XT,YT ! ! -------------------------------------------------------------------- / ! 3. Find enclosing cell and compute closest point ! INGRID = W3GFCL(GSU,XT,YT,IS,JS,XS,YS,DEBUG=LDBG) IF ( INGRID ) THEN DMIN = BIG DO L=1,4 DD = W3DIST(LLG,XT,YT,XS(L),YS(L)) IF ( DD .LT. DMIN ) THEN DMIN = DD; IX = IS(L); IY = JS(L); END IF END DO !L ELSE IX = 0; IY = 0; END IF !/ !/ End of W3GFPT_R4 -------------------------------------------------- / !/ END FUNCTION W3GFPT_R4 !/ ------------------------------------------------------------------- / FUNCTION W3GFPT_R8(GSU, XTIN, YTIN, IX, IY, DEBUG) & RESULT(INGRID) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 01-Dec-2010 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ 12-Nov-2010 : Implement r4 & r8 interfaces. ( version 3.14 ) !/ 01-Dec-2010 : Some cleanup. ( version 3.14 ) !/ ! 1. Purpose : ! ! Find point in grid, associated with the input grid-search-utility ! object (GSU), that is closest to the target point (xtin,ytin). ! Double precision interface. ! ! 2. Method : ! ! 3. Parameters : ! ! Return parameter ! ---------------------------------------------------------------- ! INGRID Log. O Logical flag indicating if target point lies ! within the source grid domain. ! ---------------------------------------------------------------- ! ! Parameter list ! ---------------------------------------------------------------- ! GSU Type I Grid-search-utility object. ! XTIN Real I X-coordinate of target point. ! YTIN Real I Y-coordinate of target point. ! IX,JX I.A. O (I,J) indices of nearest grid point. ! DEBUG Log. I Optional logical flag to turn on debug mode. ! Default is FALSE. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! - Check on previous initialization of grid search utility object. ! ! 7. Remarks : ! ! 8. Structure : ! ! ----------------------------------------------------------------- ! 1. Test input ! 2. Initialize search ! 3. Find enclosing cell and compute closest point ! ----------------------------------------------------------------- ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Return parameter !/ LOGICAL :: INGRID !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ TYPE(T_GSU), INTENT(IN) :: GSU REAL(8), INTENT(IN) :: XTIN REAL(8), INTENT(IN) :: YTIN INTEGER, INTENT(OUT) :: IX, IY LOGICAL, INTENT(IN), OPTIONAL :: DEBUG !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ REAL(8), PARAMETER :: BIG = 1D16 LOGICAL :: LDBG INTEGER :: I, J, K, L REAL(8) :: XT, YT INTEGER :: IS(4), JS(4) REAL(8) :: XS(4), YS(4) REAL(8) :: DD, DMIN LOGICAL :: IJG, LLG !/ ! ! -------------------------------------------------------------------- / ! 1. Test input ! IF ( .NOT.ASSOCIATED(GSU%PTR) ) THEN WRITE(*,'(/2A/)') 'W3GFPT_R8 ERROR -- ', & 'grid search utility object not created' CALL EXTCDE (1) END IF ! ! -------------------------------------------------------------------- / ! 2. Initialize search ! IF ( PRESENT(DEBUG) ) THEN LDBG = DEBUG ELSE LDBG = .FALSE. END IF ! ! Local pointers to grid search utility object data IJG = GSU%PTR%IJG LLG = GSU%PTR%LLG ! INGRID = .FALSE. ! XT = XTIN; YT = YTIN; IF ( LDBG ) WRITE(*,'(/A,2E14.6)') 'W3GFPT_R8 - TARGET POINT:',XT,YT ! ! -------------------------------------------------------------------- / ! 3. Find enclosing cell and compute closest point ! INGRID = W3GFCL(GSU,XT,YT,IS,JS,XS,YS,DEBUG=LDBG) IF ( INGRID ) THEN DMIN = BIG DO L=1,4 DD = W3DIST(LLG,XT,YT,XS(L),YS(L)) IF ( DD .LT. DMIN ) THEN DMIN = DD; IX = IS(L); IY = JS(L); END IF END DO !L ELSE IX = 0; IY = 0; END IF !/ !/ End of W3GFPT_R8 -------------------------------------------------- / !/ END FUNCTION W3GFPT_R8 !/ ------------------------------------------------------------------- / FUNCTION W3GFIJ_R4(GSU, XTIN, YTIN, IX, JX, DEBUG) RESULT(INGRID) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 01-Dec-2010 | !/ +-----------------------------------+ !/ !/ 12-Nov-2010 : Origination. ( version 3.14 ) !/ 01-Dec-2010 : Some cleanup. ( version 3.14 ) !/ ! 1. Purpose : ! ! Compute coordinates ( ix, jx ) of target point ( xtin, ytin ) in ! source grid index space from source grid associated with the input ! grid search utility object (GSU). ! Single precision interface. ! ! 2. Method : ! ! 3. Parameters : ! ! Return parameter ! ---------------------------------------------------------------- ! INGRID Log. O Logical flag indicating if target point lies ! within the source grid domain. ! ---------------------------------------------------------------- ! ! Parameter list ! ---------------------------------------------------------------- ! GSU Type I Grid-search-utility object. ! XTIN Real I X-coordinate of target point. ! YTIN Real I Y-coordinate of target point. ! IX Real O X-coordinate of target point in source grid ! index space. ! JX Real O Y-coordinate of target point in source grid ! index space. ! DEBUG Log. I Optional logical flag to turn on debug mode. ! Default is FALSE. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! - Check on previous initialization of grid search utility object. ! - Check on appropriate input of optional arguments. ! ! 7. Remarks : ! ! 8. Structure : ! ! ----------------------------------------------------------------- ! 1. Test input ! 2. Initialize search ! 3. Find enclosing cell and compute relative coordinates ! ----------------------------------------------------------------- ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Return parameter !/ LOGICAL :: INGRID !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ TYPE(T_GSU), INTENT(IN) :: GSU REAL(4), INTENT(IN) :: XTIN REAL(4), INTENT(IN) :: YTIN REAL(4), INTENT(OUT) :: IX REAL(4), INTENT(OUT) :: JX LOGICAL, INTENT(IN) , OPTIONAL :: DEBUG !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ LOGICAL :: LDBG, POLE INTEGER :: IS(4), JS(4) REAL(4) :: XT, YT, XS(4), YS(4) !/ ! ! -------------------------------------------------------------------- / ! 1. Test input ! IF ( .NOT.ASSOCIATED(GSU%PTR) ) THEN WRITE(*,'(/2A/)') 'W3GFIJ_R4 ERROR -- ', & 'grid search utility object not created' CALL EXTCDE (1) END IF ! ! -------------------------------------------------------------------- / ! 2. Initialize search ! IF ( PRESENT(DEBUG) ) THEN LDBG = DEBUG ELSE LDBG = .FALSE. END IF ! XT = XTIN; YT = YTIN; IF ( LDBG ) WRITE(*,'(/A,2E14.6)') 'W3GFIJ_R4 - TARGET POINT:',XT,YT ! ! -------------------------------------------------------------------- / ! 3. Find enclosing cell and compute point location ! INGRID = W3GFCL(GSU,XT,YT,IS,JS,XS,YS,POLE=POLE,DEBUG=LDBG) IF ( .NOT. INGRID ) RETURN ! IF ( .NOT.POLE ) THEN !---------non-pole cell: compute relative location CALL W3RMBL(XT,YT,XS,YS,IX=IX,JX=JX,DEBUG=LDBG) IX = REAL(IS(1),4)+IX; JX = REAL(JS(1),4)+JX; ELSE !---------pole cell: set to center of cell IX = REAL(IS(1),4)+HALF; JX = REAL(JS(1),4)+HALF; ENDIF !/ !/ End of W3GFIJ_R4 -------------------------------------------------- / !/ END FUNCTION W3GFIJ_R4 !/ ------------------------------------------------------------------- / FUNCTION W3GFIJ_R8(GSU, XTIN, YTIN, IX, JX, DEBUG) RESULT(INGRID) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 01-Dec-2010 | !/ +-----------------------------------+ !/ !/ 12-Nov-2010 : Origination. ( version 3.14 ) !/ 01-Dec-2010 : Some cleanup. ( version 3.14 ) !/ ! 1. Purpose : ! ! Compute coordinates ( ix, jx ) of target point ( xtin, ytin ) in ! source grid index space from source grid associated with the input ! grid search utility object (GSU). ! Double precision interface. ! ! 2. Method : ! ! 3. Parameters : ! ! Return parameter ! ---------------------------------------------------------------- ! INGRID Log. O Logical flag indicating if target point lies ! within the source grid domain. ! ---------------------------------------------------------------- ! ! Parameter list ! ---------------------------------------------------------------- ! GSU Type I Grid-search-utility object. ! XTIN Real I X-coordinate of target point. ! YTIN Real I Y-coordinate of target point. ! IX Real O X-coordinate of target point in source grid ! index space. ! JX Real O Y-coordinate of target point in source grid ! index space. ! DEBUG Log. I Optional logical flag to turn on debug mode. ! Default is FALSE. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! - Check on previous initialization of grid search utility object. ! - Check on appropriate input of optional arguments. ! ! 7. Remarks : ! ! 8. Structure : ! ! ----------------------------------------------------------------- ! 1. Test input ! 2. Initialize search ! 3. Find enclosing cell and compute relative coordinates ! ----------------------------------------------------------------- ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Return parameter !/ LOGICAL :: INGRID !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ TYPE(T_GSU), INTENT(IN) :: GSU REAL(8), INTENT(IN) :: XTIN REAL(8), INTENT(IN) :: YTIN REAL(8), INTENT(OUT) :: IX REAL(8), INTENT(OUT) :: JX LOGICAL, INTENT(IN) , OPTIONAL :: DEBUG !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ LOGICAL :: LDBG, POLE INTEGER :: IS(4), JS(4) REAL(8) :: XT, YT, XS(4), YS(4) !/ ! ! -------------------------------------------------------------------- / ! 1. Test input ! IF ( .NOT.ASSOCIATED(GSU%PTR) ) THEN WRITE(*,'(/2A/)') 'W3GFIJ_R8 ERROR -- ', & 'grid search utility object not created' CALL EXTCDE (1) END IF ! ! -------------------------------------------------------------------- / ! 2. Initialize search ! IF ( PRESENT(DEBUG) ) THEN LDBG = DEBUG ELSE LDBG = .FALSE. END IF ! XT = XTIN; YT = YTIN; IF ( LDBG ) WRITE(*,'(/A,2E14.6)') 'W3GFIJ_R8 - TARGET POINT:',XT,YT ! ! -------------------------------------------------------------------- / ! 3. Find enclosing cell and compute point location ! INGRID = W3GFCL(GSU,XT,YT,IS,JS,XS,YS,POLE=POLE,DEBUG=LDBG) IF ( .NOT. INGRID ) RETURN ! IF ( .NOT.POLE ) THEN !---------non-pole cell: compute relative location CALL W3RMBL(XT,YT,XS,YS,IX=IX,JX=JX,DEBUG=LDBG) IX = REAL(IS(1),8)+IX; JX = REAL(JS(1),8)+JX; ELSE !---------pole cell: set to center of cell IX = REAL(IS(1),8)+HALF; JX = REAL(JS(1),8)+HALF; ENDIF !/ !/ End of W3GFIJ_R8 -------------------------------------------------- / !/ END FUNCTION W3GFIJ_R8 !/ ------------------------------------------------------------------- / FUNCTION W3GRMP_R4(GSU, XTIN, YTIN, IS, JS, RW, & MASK, MSKC, NNBR, DEBUG) RESULT(INGRID) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 15-Jun-2012 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ 12-Nov-2010 : Implement r4 & r8 interfaces. ( version 3.14 ) !/ 01-Dec-2010 : Some cleanup. ( version 3.14 ) !/ 06-Dec-2010 : Remove restriction on longitude range. Change ICLO !/ to integer and remove JCLO. Implement support for !/ r4 and r8 source grids. ( version 3.14 ) !/ 15-Jun-2012 : Fixing format statement that gave warning with !/ Intell compiler (H. L. Tolman). ( version 4.07 ) !/ ! 1. Purpose : ! ! Compute remapping for target point ( xtin, ytin ) from source grid ! associated with the input grid search utility object (GSU). ! The indices of the source points used for remapping are returned in ! is(1:4) and js(1:4). The remapping weights are returned in rw(1:4). ! Single precision interface. ! ! 2. Method : ! ! 3. Parameters : ! ! Return parameter ! ---------------------------------------------------------------- ! INGRID Log. O Logical flag indicating if target point lies ! within the source grid domain. ! ---------------------------------------------------------------- ! ! Parameter list ! ---------------------------------------------------------------- ! GSU Type I Grid-search-utility object. ! XTIN Real I X-coordinate of target point. ! YTIN Real I Y-coordinate of target point. ! IS,JS I.A. O (I,J) indices of vertices of enclosing grid cell. ! RW R.A. O Array of interpolation weights. ! MASK L.A. I Optional logical mask for source grid. ! MSKC Int. O Optional output integer parameter indicating how ! the enclosing cell is masked. Possible values ! are MSKC_NONE, MSKC_PART and MSKC_FULL. ! MSKC is required when MASK is specified. ! NNBR Int. I/O Optional integer parameter indicating the number ! of nearest-neighbor non-masked points used for ! distance-weighted averaging. ! Input: Requested number of nearest-neighbor ! non-masked points (0 < NNBR <= 4). ! Output: Actual number of nearest-neighbor ! non-masked points used. ! DEBUG Log. I Optional logical flag to turn on debug mode. ! Default is FALSE. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! - Check on previous initialization of grid search utility object. ! - Check on appropriate input of optional arguments. ! ! 7. Remarks : ! ! 8. Structure : ! ! ----------------------------------------------------------------- ! 1. Test input ! 2. Initialize search ! 3. Find enclosing cell and compute remapping weights ! - if enclosing cell does not includes a pole, then ! compute bilinear remapping ! - if enclosing cell includes a pole, then ! compute distance weighted remapping ! 4. Handle case of target point located within a partially masked cell. ! 5. Handle case of target point located within a fully masked cell. ! ----------------------------------------------------------------- ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Return parameter !/ LOGICAL :: INGRID !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ TYPE(T_GSU), INTENT(IN) :: GSU REAL(4), INTENT(IN) :: XTIN REAL(4), INTENT(IN) :: YTIN INTEGER, INTENT(OUT) :: IS(4) INTEGER, INTENT(OUT) :: JS(4) REAL(4), INTENT(OUT) :: RW(4) LOGICAL, INTENT(IN) , OPTIONAL :: MASK(:,:) INTEGER, INTENT(OUT) , OPTIONAL :: MSKC INTEGER, INTENT(INOUT), OPTIONAL :: NNBR LOGICAL, INTENT(IN) , OPTIONAL :: DEBUG !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ REAL(8), PARAMETER :: BIG = 1D16 REAL(8), PARAMETER :: SMALL = 1D-6 LOGICAL :: LDBG, POLE INTEGER :: I, J, K, L, IB, JB, IBC, JBC LOGICAL :: M, MSK(4) INTEGER :: LVL, N, NS, ICC, JCC REAL(4) :: XT, YT, XS(4), YS(4) REAL(4) :: X, Y, D(4), DD, DMIN, DSUM LOGICAL :: IJG, LLG, LCLO INTEGER :: ICLO, GKIND INTEGER :: NX, NY REAL(4), POINTER :: XG4(:,:), YG4(:,:) REAL(8), POINTER :: XG8(:,:), YG8(:,:) TYPE(T_NNS), POINTER :: NNP !/ ! ! -------------------------------------------------------------------- / ! 1. Test input ! IF ( .NOT.ASSOCIATED(GSU%PTR) ) THEN WRITE(*,'(/2A/)') 'W3GRMP_R4 ERROR -- ', & 'grid search utility object not created' CALL EXTCDE (1) END IF ! IF ( PRESENT(MASK) ) THEN IF ( .NOT.PRESENT(MSKC) ) THEN WRITE(*,'(/2A/)') 'W3GRMP_R4 ERROR -- ', & 'MSKC must be specified with MASK' CALL EXTCDE (1) END IF IF ( PRESENT(NNBR) ) THEN IF ( .NOT.ASSOCIATED(GSU%PTR%NNP) ) THEN WRITE(*,'(/3A/)') 'W3GRMP_R4 ERROR -- ', & 'MASK and NNBR input specified, ', & 'but grid point-search object not created' CALL EXTCDE (1) END IF IF ( NNBR .LE. 0 .OR. NNBR .GT. 4 ) THEN WRITE(*,'(/2A/)') 'W3GRMP_R4 ERROR -- ', & 'NNBR must be >= 1 AND <= 4' CALL EXTCDE (1) END IF END IF END IF ! ! -------------------------------------------------------------------- / ! 2. Initialize search ! IF ( PRESENT(DEBUG) ) THEN LDBG = DEBUG ELSE LDBG = .FALSE. END IF ! ! Local pointers to grid search utility object data IJG = GSU%PTR%IJG LLG = GSU%PTR%LLG ICLO = GSU%PTR%ICLO LCLO = GSU%PTR%LCLO GKIND = GSU%PTR%GKIND NX = GSU%PTR%NX; NY = GSU%PTR%NY; IF ( GKIND.EQ.4 ) THEN XG4 => GSU%PTR%XG4; YG4 => GSU%PTR%YG4; ELSE XG8 => GSU%PTR%XG8; YG8 => GSU%PTR%YG8; END IF NNP => GSU%PTR%NNP ! RW = ZERO; ! XT = XTIN; YT = YTIN; IF ( LDBG ) WRITE(*,'(/A,2E14.6)') 'W3GRMP_R4 - TARGET POINT:',XT,YT ! ! -------------------------------------------------------------------- / ! 3. Find enclosing cell and compute remapping ! INGRID = W3GFCL(GSU,XT,YT,IS,JS,XS,YS,POLE=POLE,DEBUG=LDBG) IF ( .NOT. INGRID ) RETURN ! IF ( .NOT.POLE ) THEN !---------non-pole cell: compute bilinear remapping CALL W3RMBL(XT,YT,XS,YS,RW=RW,DEBUG=LDBG) IF ( LDBG ) THEN WRITE(*,'(A,2E14.6)') 'W3GRMP_R4 - BILINEAR (TGT):',XT,YT DO L=1,4 WRITE(*,'(A,3I6,E14.6)') 'W3GRMP_R4 - BILINEAR (SRC):', & L,IS(L),JS(L),RW(L) END DO END IF !LDBG ELSE !---------pole cell: compute distance-weighted remapping DSUM = ZERO DO L=1,4 D(L) = W3DIST(LLG,XT,YT,XS(L),YS(L)) DSUM = DSUM + ONE/(D(L)+SMALL) END DO RW(1:4) = ONE/(D(1:4)+SMALL)/DSUM IF ( LDBG ) THEN WRITE(*,'(A,2E14.6)') 'W3GRMP_R4 - DISTWGHT (TGT):',XT,YT DO L=1,4 WRITE(*,'(A,3I6,E14.6)') 'W3GRMP_R4 - DISTWGHT (SRC):', & L,IS(L),JS(L),RW(L) END DO END IF !LDBG ENDIF ! IF ( .NOT.PRESENT(MASK) ) RETURN ! ! -------------------------------------------------------------------- / ! 4. Handle case of target point located within a partially masked cell. ! !-----copy cell mask values according to array ordering IF ( IJG ) THEN DO L=1,4 MSK(L) = MASK(IS(L),JS(L)) END DO ELSE DO L=1,4 MSK(L) = MASK(JS(L),IS(L)) END DO END IF ! !-----adjust weights for a partially masked cell DSUM = ZERO NS = 4 DO L=1,4 IF ( MSK(L) ) THEN NS = NS - 1 RW(L) = ZERO END IF DSUM = DSUM + RW(L) END DO IF ( NS .EQ. 4 ) THEN MSKC = MSKC_NONE RETURN END IF IF ( NS .GT. 0 .AND. DSUM .GT. SMALL ) THEN RW = RW / DSUM IF ( LDBG ) & WRITE(*,'(A,2E14.6,4(2I6,E14.6))') & 'W3GRMP_R4 - PARTIAL MASKED CELL:', & XT,YT,(IS(L),JS(L),RW(L),L=1,4) MSKC = MSKC_PART RETURN ELSE MSKC = MSKC_FULL IF ( .NOT.PRESENT(NNBR) ) RETURN END IF ! ! -------------------------------------------------------------------- / ! 5. Handle case of target point located within a fully masked cell. ! ! Choose closest point in enclosing land cell to be the central point DMIN = BIG DO L=1,4 DD = W3DIST(LLG,XT,YT,XS(L),YS(L)) IF ( DD .LT. DMIN ) THEN DMIN = DD; ICC = IS(L); JCC = JS(L); END IF END DO !L ! ! Search nearest-neighbor source points for closest nnbr un-masked ! points and compute distance-weighted average remapping. IF ( LDBG ) & WRITE(*,'(A,2I6)') & 'W3GRMP_R4 - BEGIN POINT NNBR SEARCH:',ICC,JCC NS = 0; D(:) = BIG; LEVEL_LOOP: DO LVL=0,NNP%NLVL NNBR_LOOP: DO N=NNP%N1(LVL),NNP%N2(LVL) I = ICC + NNP%DI(N); J = JCC + NNP%DJ(N); IF ( ICLO.EQ.ICLO_NONE ) THEN IF ( I.LT.1 .OR. I.GT.NX ) CYCLE NNBR_LOOP IF ( J.LT.1 .OR. J.GT.NY ) CYCLE NNBR_LOOP END IF IF ( ICLO.NE.ICLO_NONE ) THEN IF ( I .LT. 1 ) I = I + NX IF ( I .GT. NX ) I = I - NX END IF IF ( ICLO.EQ.ICLO_TRPL ) THEN IF ( J .GT. NY ) THEN J = NY I = MOD(NX-I,NX) + 1 END IF END IF IF ( IJG ) THEN M = MASK(I,J) ELSE M = MASK(J,I) END IF IF ( LDBG ) & WRITE(*,'(A,4I6,1L6)') & 'W3GRMP_R4 - POINT NNBR SEARCH:',LVL,N,I,J,M !-------------if masked point, then skip IF ( M ) CYCLE NNBR_LOOP !-------------compute distance IF ( IJG ) THEN IF ( GKIND.EQ.4 ) THEN X = XG4(I,J); Y = YG4(I,J); ELSE X = XG8(I,J); Y = YG8(I,J); END IF ELSE IF ( GKIND.EQ.4 ) THEN X = XG4(J,I); Y = YG4(J,I); ELSE X = XG8(J,I); Y = YG8(J,I); END IF END IF DD = W3DIST(LLG,XT,YT,X,Y) !-------------still need nnbr points IF ( NS .LT. NNBR ) THEN !-----------------add to list NS = NS + 1 IS(NS) = I; JS(NS) = J; D(NS) = DD; !-----------------once list is full sort according to increasing distance IF ( NS .EQ. NNBR ) CALL W3SORT(NS,IS,JS,D) !---------------we have found nnbr points ELSE !list is full !-----------------insert into list if the newest point is closer CALL W3ISRT(I,J,DD,NS,IS,JS,D) END IF !list is full IF ( LDBG ) & WRITE(*,'(A,I2,I3,I6,4(2I6,E14.6))') & 'W3GRMP_R4 - POINT NNBR LIST:', & LVL,N,NS,(IS(L),JS(L),D(L),L=1,NS) END DO NNBR_LOOP !---------if we have found nnbr_rqd points, then exit the search IF ( NS .EQ. NNBR ) EXIT LEVEL_LOOP END DO LEVEL_LOOP NNBR = NS ! ! If zero unmasked points found, then return nnbr=0 as error indicator IF ( NNBR .EQ. 0 ) RETURN ! ! Compute distance-weighted remapping for nnbr points DSUM = ZERO DO L=1,NNBR DSUM = DSUM + ONE/(D(L)+SMALL) END DO RW(1:NNBR) = ONE/(D(1:NNBR)+SMALL)/DSUM IF ( LDBG ) THEN WRITE(*,'(A,2E14.6,I6)') & 'W3GRMP_R4 - FULLY MASKED CELL (TGT):',XT,YT,NNBR DO L=1,NNBR WRITE(*,'(A,3I6,E14.6)') & 'W3GRMP_R4 - FULLY MASKED CELL (SRC):', & L,IS(L),JS(L),RW(L) END DO END IF !LDBG !/ !/ End of W3GRMP_R4 -------------------------------------------------- / !/ END FUNCTION W3GRMP_R4 !/ ------------------------------------------------------------------- / FUNCTION W3GRMP_R8(GSU, XTIN, YTIN, IS, JS, RW, & MASK, MSKC, NNBR, DEBUG) RESULT(INGRID) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 15-Jun-2012 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ 12-Nov-2010 : Implement r4 & r8 interfaces. ( version 3.14 ) !/ 01-Dec-2010 : Some cleanup. ( version 3.14 ) !/ 06-Dec-2010 : Remove restriction on longitude range. Change ICLO !/ to integer and remove JCLO. Implement support for !/ r4 and r8 source grids. ( version 3.14 ) !/ 15-Jun-2012 : Fixing format statement that gave warning with !/ Intell compiler (H. L. Tolman). ( version 4.07 ) !/ ! 1. Purpose : ! ! Compute remapping for target point ( xtin, ytin ) from source grid ! associated with the input grid search utility object (GSU). ! The indices of the source points used for remapping are returned in ! is(1:4) and js(1:4). The remapping weights are returned in rw(1:4). ! Double precision interface. ! ! 2. Method : ! ! 3. Parameters : ! ! Return parameter ! ---------------------------------------------------------------- ! INGRID Log. O Logical flag indicating if target point lies ! within the source grid domain. ! ---------------------------------------------------------------- ! ! Parameter list ! ---------------------------------------------------------------- ! GSU Type I Grid-search-utility object. ! XTIN Real I X-coordinate of target point. ! YTIN Real I Y-coordinate of target point. ! IS,JS I.A. O (I,J) indices of vertices of enclosing grid cell. ! RW R.A. O Array of interpolation weights. ! MASK L.A. I Optional logical mask for source grid. ! MSKC Int. O Optional output integer parameter indicating how ! the enclosing cell is masked. Possible values ! are MSKC_NONE, MSKC_PART and MSKC_FULL. ! MSKC is required when MASK is specified. ! NNBR Int. I/O Optional integer parameter indicating the number ! of nearest-neighbor non-masked points used for ! distance-weighted averaging. ! Input: Requested number of nearest-neighbor ! non-masked points (0 < NNBR <= 4). ! Output: Actual number of nearest-neighbor ! non-masked points used. ! DEBUG Log. I Optional logical flag to turn on debug mode. ! Default is FALSE. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! - Check on previous initialization of grid search utility object. ! - Check on appropriate input of optional arguments. ! ! 7. Remarks : ! ! 8. Structure : ! ! ----------------------------------------------------------------- ! 1. Test input ! 2. Initialize search ! 3. Find enclosing cell and compute remapping weights ! - if enclosing cell does not includes a pole, then ! compute bilinear remapping ! - if enclosing cell includes a pole, then ! compute distance weighted remapping ! 4. Handle case of target point located within a partially masked cell. ! 5. Handle case of target point located within a fully masked cell. ! ----------------------------------------------------------------- ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Return parameter !/ LOGICAL :: INGRID !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ TYPE(T_GSU), INTENT(IN) :: GSU REAL(8), INTENT(IN) :: XTIN REAL(8), INTENT(IN) :: YTIN INTEGER, INTENT(OUT) :: IS(4) INTEGER, INTENT(OUT) :: JS(4) REAL(8), INTENT(OUT) :: RW(4) LOGICAL, INTENT(IN) , OPTIONAL :: MASK(:,:) INTEGER, INTENT(OUT) , OPTIONAL :: MSKC INTEGER, INTENT(INOUT), OPTIONAL :: NNBR LOGICAL, INTENT(IN) , OPTIONAL :: DEBUG !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ REAL(8), PARAMETER :: BIG = 1D16 REAL(8), PARAMETER :: SMALL = 1D-6 LOGICAL :: LDBG, POLE INTEGER :: I, J, K, L, IB, JB, IBC, JBC LOGICAL :: M, MSK(4) INTEGER :: LVL, N, NS, ICC, JCC REAL(8) :: XT, YT, XS(4), YS(4) REAL(8) :: X, Y, D(4), DD, DMIN, DSUM LOGICAL :: IJG, LLG, LCLO INTEGER :: ICLO, GKIND INTEGER :: NX, NY REAL(4), POINTER :: XG4(:,:), YG4(:,:) REAL(8), POINTER :: XG8(:,:), YG8(:,:) TYPE(T_NNS), POINTER :: NNP !/ ! ! -------------------------------------------------------------------- / ! 1. Test input ! IF ( .NOT.ASSOCIATED(GSU%PTR) ) THEN WRITE(*,'(/2A/)') 'W3GRMP_R8 ERROR -- ', & 'grid search utility object not created' CALL EXTCDE (1) END IF ! IF ( PRESENT(MASK) ) THEN IF ( .NOT.PRESENT(MSKC) ) THEN WRITE(*,'(/2A/)') 'W3GRMP_R8 ERROR -- ', & 'MSKC must be specified with MASK' CALL EXTCDE (1) END IF IF ( PRESENT(NNBR) ) THEN IF ( .NOT.ASSOCIATED(GSU%PTR%NNP) ) THEN WRITE(*,'(/3A/)') 'W3GRMP_R8 ERROR -- ', & 'MASK and NNBR input specified, ', & 'but grid point-search object not created' CALL EXTCDE (1) END IF IF ( NNBR .LE. 0 .OR. NNBR .GT. 4 ) THEN WRITE(*,'(/2A/)') 'W3GRMP_R8 ERROR -- ', & 'NNBR must be >= 1 AND <= 4' CALL EXTCDE (1) END IF END IF END IF ! ! -------------------------------------------------------------------- / ! 2. Initialize search ! IF ( PRESENT(DEBUG) ) THEN LDBG = DEBUG ELSE LDBG = .FALSE. END IF ! ! Local pointers to grid search utility object data IJG = GSU%PTR%IJG LLG = GSU%PTR%LLG ICLO = GSU%PTR%ICLO LCLO = GSU%PTR%LCLO GKIND = GSU%PTR%GKIND NX = GSU%PTR%NX; NY = GSU%PTR%NY; IF ( GKIND.EQ.4 ) THEN XG4 => GSU%PTR%XG4; YG4 => GSU%PTR%YG4; ELSE XG8 => GSU%PTR%XG8; YG8 => GSU%PTR%YG8; END IF NNP => GSU%PTR%NNP ! RW = ZERO; ! XT = XTIN; YT = YTIN; IF ( LDBG ) WRITE(*,'(/A,2E14.6)') 'W3GRMP_R8 - TARGET POINT:',XT,YT ! ! -------------------------------------------------------------------- / ! 3. Find enclosing cell and compute remapping ! INGRID = W3GFCL(GSU,XT,YT,IS,JS,XS,YS,POLE=POLE,DEBUG=LDBG) IF ( .NOT. INGRID ) RETURN ! IF ( .NOT.POLE ) THEN !---------non-pole cell: compute bilinear remapping CALL W3RMBL(XT,YT,XS,YS,RW=RW,DEBUG=LDBG) IF ( LDBG ) THEN WRITE(*,'(A,2E14.6)') 'W3GRMP_R8 - BILINEAR (TGT):',XT,YT DO L=1,4 WRITE(*,'(A,3I6,E14.6)') 'W3GRMP_R8 - BILINEAR (SRC):', & L,IS(L),JS(L),RW(L) END DO END IF !LDBG ELSE !---------pole cell: compute distance-weighted remapping DSUM = ZERO DO L=1,4 D(L) = W3DIST(LLG,XT,YT,XS(L),YS(L)) DSUM = DSUM + ONE/(D(L)+SMALL) END DO RW(1:4) = ONE/(D(1:4)+SMALL)/DSUM IF ( LDBG ) THEN WRITE(*,'(A,2E14.6)') 'W3GRMP_R8 - DISTWGHT (TGT):',XT,YT DO L=1,4 WRITE(*,'(A,3I6,E14.6)') 'W3GRMP_R8 - DISTWGHT (SRC):', & L,IS(L),JS(L),RW(L) END DO END IF !LDBG ENDIF ! IF ( .NOT.PRESENT(MASK) ) RETURN ! ! -------------------------------------------------------------------- / ! 4. Handle case of target point located within a partially masked cell. ! !-----copy cell mask values according to array ordering IF ( IJG ) THEN DO L=1,4 MSK(L) = MASK(IS(L),JS(L)) END DO ELSE DO L=1,4 MSK(L) = MASK(JS(L),IS(L)) END DO END IF ! !-----adjust weights for a partially masked cell DSUM = ZERO NS = 4 DO L=1,4 IF ( MSK(L) ) THEN NS = NS - 1 RW(L) = ZERO END IF DSUM = DSUM + RW(L) END DO IF ( NS .EQ. 4 ) THEN MSKC = MSKC_NONE RETURN END IF IF ( NS .GT. 0 .AND. DSUM .GT. SMALL ) THEN RW = RW / DSUM IF ( LDBG ) & WRITE(*,'(A,2E14.6,4(2I6,E14.6))') & 'W3GRMP_R8 - PARTIAL MASKED CELL:', & XT,YT,(IS(L),JS(L),RW(L),L=1,4) MSKC = MSKC_PART RETURN ELSE MSKC = MSKC_FULL IF ( .NOT.PRESENT(NNBR) ) RETURN END IF ! ! -------------------------------------------------------------------- / ! 5. Handle case of target point located within a fully masked cell. ! ! Choose closest point in enclosing land cell to be the central point DMIN = BIG DO L=1,4 DD = W3DIST(LLG,XT,YT,XS(L),YS(L)) IF ( DD .LT. DMIN ) THEN DMIN = DD; ICC = IS(L); JCC = JS(L); END IF END DO !L ! ! Search nearest-neighbor source points for closest nnbr un-masked ! points and compute distance-weighted average remapping. IF ( LDBG ) & WRITE(*,'(A,2I6)') & 'W3GRMP_R8 - BEGIN POINT NNBR SEARCH:',ICC,JCC NS = 0; D(:) = BIG; LEVEL_LOOP: DO LVL=0,NNP%NLVL NNBR_LOOP: DO N=NNP%N1(LVL),NNP%N2(LVL) I = ICC + NNP%DI(N); J = JCC + NNP%DJ(N); IF ( ICLO.EQ.ICLO_NONE ) THEN IF ( I.LT.1 .OR. I.GT.NX ) CYCLE NNBR_LOOP IF ( J.LT.1 .OR. J.GT.NY ) CYCLE NNBR_LOOP END IF IF ( ICLO.NE.ICLO_NONE ) THEN IF ( I .LT. 1 ) I = I + NX IF ( I .GT. NX ) I = I - NX END IF IF ( ICLO.EQ.ICLO_TRPL ) THEN IF ( J .GT. NY ) THEN J = NY I = MOD(NX-I,NX) + 1 END IF END IF IF ( IJG ) THEN M = MASK(I,J) ELSE M = MASK(J,I) END IF IF ( LDBG ) & WRITE(*,'(A,4I6,1L6)') & 'W3GRMP_R8 - POINT NNBR SEARCH:',LVL,N,I,J,M !-------------if masked point, then skip IF ( M ) CYCLE NNBR_LOOP !-------------compute distance IF ( IJG ) THEN IF ( GKIND.EQ.4 ) THEN X = XG4(I,J); Y = YG4(I,J); ELSE X = XG8(I,J); Y = YG8(I,J); END IF ELSE IF ( GKIND.EQ.4 ) THEN X = XG4(J,I); Y = YG4(J,I); ELSE X = XG8(J,I); Y = YG8(J,I); END IF END IF DD = W3DIST(LLG,XT,YT,X,Y) !-------------still need nnbr points IF ( NS .LT. NNBR ) THEN !-----------------add to list NS = NS + 1 IS(NS) = I; JS(NS) = J; D(NS) = DD; !-----------------once list is full sort according to increasing distance IF ( NS .EQ. NNBR ) CALL W3SORT(NS,IS,JS,D) !---------------we have found nnbr points ELSE !list is full !-----------------insert into list if the newest point is closer CALL W3ISRT(I,J,DD,NS,IS,JS,D) END IF !list is full IF ( LDBG ) & WRITE(*,'(A,I2,I3,I6,4(2I6,E14.6))') & 'W3GRMP_R8 - POINT NNBR LIST:', & LVL,N,NS,(IS(L),JS(L),D(L),L=1,NS) END DO NNBR_LOOP !---------if we have found nnbr_rqd points, then exit the search IF ( NS .EQ. NNBR ) EXIT LEVEL_LOOP END DO LEVEL_LOOP NNBR = NS ! ! If zero unmasked points found, then return nnbr=0 as error indicator IF ( NNBR .EQ. 0 ) RETURN ! ! Compute distance-weighted remapping for nnbr points DSUM = ZERO DO L=1,NNBR DSUM = DSUM + ONE/(D(L)+SMALL) END DO RW(1:NNBR) = ONE/(D(1:NNBR)+SMALL)/DSUM IF ( LDBG ) THEN WRITE(*,'(A,2E14.6,I6)') & 'W3GRMP_R8 - FULLY MASKED CELL (TGT):',XT,YT,NNBR DO L=1,NNBR WRITE(*,'(A,3I6,E14.6)') & 'W3GRMP_R8 - FULLY MASKED CELL (SRC):', & L,IS(L),JS(L),RW(L) END DO END IF !LDBG !/ !/ End of W3GRMP_R8 -------------------------------------------------- / !/ END FUNCTION W3GRMP_R8 !/ ------------------------------------------------------------------- / FUNCTION W3NNSC(NLVL) RESULT(NNS) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 30-Oct-2009 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ ! 1. Purpose : ! ! Create nearest-neighbor (NNBR) search object. ! ! 2. Method : ! ! Notation ! ( L, N): L = NNBR level; N = NNBR sequential index ! {DI, DJ}: DI = I-index delta; DJ = J-index delta ! ! --------------------------------------------------- ! | ( 2,21) | ( 2,20) | ( 2,19) | ( 2,18) | ( 2,17) | ! | {-2,+2} | {-1,+2} | { 0,+2} | {+1,+2} | {+2,+2} | ! --------------------------------------------------- ! | ( 2,22) | ( 1, 7) | ( 1, 6) | ( 1, 5) | ( 2,16) | ! | {-2,+1} | {-1,+1} | { 0,+1} | {+1,+1} | {+2,+1} | ! --------------------------------------------------- ! | ( 2,23) | ( 1, 8) | ( 0, 0) | ( 1, 4) | ( 2,15) | ! | {-2, 0} | {-1, 0} | { 0, 0} | {+1, 0} | {+2, 0} | ! --------------------------------------------------- ! | ( 2,24) | ( 1, 1) | ( 1, 2) | ( 1, 3) | ( 2,14) | ! | {-2,-1} | {-1,-1} | { 0,-1} | {+1,-1} | {+2,-1} | ! --------------------------------------------------- ! | ( 2, 9) | ( 2,10) | ( 2,11) | ( 2,12) | ( 2,13) | ! | {-2,-2} | {-1,-2} | { 0,-2} | {+1,-2} | {+2,-2} | ! --------------------------------------------------- ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Return parameter !/ TYPE(T_NNS), POINTER :: NNS !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ INTEGER, INTENT(IN) :: NLVL !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: I, J, K, L, N !/ ! ! -------------------------------------------------------------------- / ! !-----allocate object ALLOCATE(NNS) !-----initialize sizes NNS%NLVL = NLVL NNS%NNBR = (2*NLVL+1)**2 !-----allocate arrays ALLOCATE(NNS%N1(0:NNS%NLVL)) ALLOCATE(NNS%N2(0:NNS%NLVL)) ALLOCATE(NNS%DI(0:NNS%NNBR-1)) ALLOCATE(NNS%DJ(0:NNS%NNBR-1)) !-----compute index deltas for nearest-neighbor searches N = 0 !-----central point L = 0 NNS%N1(L) = 0; NNS%N2(L) = (2*L+1)**2-1; NNS%DI(N) = 0; NNS%DJ(N) = 0; !-----loop over levels DO L=1,NNS%NLVL !---------nnbr loop bounds NNS%N1(L) = (2*L-1)**2; NNS%N2(L) = (2*L+1)**2-1; !---------bottom-layer J = -L DO I=-L,L-1 N = N + 1 NNS%DI(N) = I; NNS%DJ(N) = J; END DO !---------right-layer I = L DO J=-L,L-1 N = N + 1 NNS%DI(N) = I; NNS%DJ(N) = J; END DO !---------top-layer J = L DO I=L,-L+1,-1 N = N + 1 NNS%DI(N) = I; NNS%DJ(N) = J; END DO !---------left-layer I = -L DO J=L,-L+1,-1 N = N + 1 NNS%DI(N) = I; NNS%DJ(N) = J; END DO END DO !loop over levels !/ !/ End of W3NNSC ----------------------------------------------------- / !/ END FUNCTION W3NNSC !/ ------------------------------------------------------------------- / SUBROUTINE W3NNSD(NNS) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 30-Oct-2009 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ ! 1. Purpose : ! ! Destroy nearest-neighbor (NNBR) search object. ! ! 2. Method : ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ TYPE(T_NNS), POINTER :: NNS !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ !/ ! ! -------------------------------------------------------------------- / ! IF ( ASSOCIATED(NNS) ) THEN NNS%NLVL = 0 NNS%NNBR = 0 IF ( ASSOCIATED(NNS%N1) ) THEN DEALLOCATE(NNS%N1); NULLIFY(NNS%N1); END IF IF ( ASSOCIATED(NNS%N2) ) THEN DEALLOCATE(NNS%N2); NULLIFY(NNS%N2); END IF IF ( ASSOCIATED(NNS%DI) ) THEN DEALLOCATE(NNS%DI); NULLIFY(NNS%DI); END IF IF ( ASSOCIATED(NNS%DJ) ) THEN DEALLOCATE(NNS%DJ); NULLIFY(NNS%DJ); END IF DEALLOCATE(NNS) NULLIFY(NNS) END IF !/ !/ End of W3NNSD ----------------------------------------------------- / !/ END SUBROUTINE W3NNSD !/ ------------------------------------------------------------------- / SUBROUTINE W3NNSP(NNS, IUNIT) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 30-Oct-2009 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ ! 1. Purpose : ! ! Print nearest-neighbor (NNBR) search object to IUNIT. ! ! 2. Method : ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! NNBR Type I Nearest-neighbor search object. ! IUNIT Int. I Optional unit for output. Default is stdout. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ TYPE(T_NNS), INTENT(IN) :: NNS INTEGER, OPTIONAL, INTENT(IN) :: IUNIT !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: NDST, K, L, N !/ ! ! -------------------------------------------------------------------- / ! IF ( PRESENT(IUNIT) ) THEN NDST = IUNIT ELSE NDST = 6 END IF ! WRITE(NDST,'(A,2I6)') 'nlvl,nnbr:',NNS%NLVL,NNS%NNBR DO L=0,NNS%NLVL DO N=NNS%N1(L),NNS%N2(L) WRITE(NDST,'(A,4I6)') 'l,n,di,dj:',L,N,NNS%DI(N),NNS%DJ(N) END DO END DO !/ !/ End of W3NNSP ----------------------------------------------------- / !/ END SUBROUTINE W3NNSP !/ ------------------------------------------------------------------- / SUBROUTINE W3RMBL_R4(XT, YT, XS, YS, RW, IX, JX, DEBUG) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 01-Dec-2010 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ 12-Nov-2010 : Add IX,IY optional arguments. !/ Implement r4 & r8 interfaces. ( version 3.14 ) !/ 01-Dec-2010 : Add check for target point coincident with a cell !/ vertex. Change to error exit when unable to !/ determine local (i,j). ( version 3.14 ) !/ ! 1. Purpose : ! ! Bilinear remapping for target point ( xt, yt ) in a cell defined ! by the source points ( xs(1:4), ys(1:4) ). Remapping weights are ! returned in rw(1:4). It is the caller's responsibility to ensure ! that the target point is located within the input cell and that the ! cell corner points are properly defined. ! ! (xs4,ys4) (xs3,ys3) ! _____________________ ! / / ! / x / ! / (xtin,ytin) / ! / / ! /____________________/ ! (xs1,ys1) (xs2,ys2) ! ! In spherical coordinates it is assumed that the longitude range of ! the target point is the same as that of the cell vertices. It is ! also assumed that the cell does not includes a branch cut. ! Single precision interface. ! ! 2. Method : ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! 7. Remarks : ! ! - Implementation is based on SCRIP. ! - In the case of spherical coordinates, the method results in ! bogus weights if enclosing cell contains a pole. ! ! 8. Structure : ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ REAL(4), INTENT(IN) :: XT REAL(4), INTENT(IN) :: YT REAL(4), INTENT(IN) :: XS(4) REAL(4), INTENT(IN) :: YS(4) REAL(4), INTENT(OUT), OPTIONAL :: RW(4) REAL(4), INTENT(OUT), OPTIONAL :: IX REAL(4), INTENT(OUT), OPTIONAL :: JX LOGICAL, INTENT(IN) , OPTIONAL :: DEBUG !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER, PARAMETER :: MAX_ITER = 10 REAL(8), PARAMETER :: CONVERGE = 1D-6 LOGICAL :: LDBG INTEGER :: K, ITER REAL(8) :: DXT, DX1, DX2, DX3, DXP, DYT, DY1, DY2, DY3, DYP REAL(8) :: IGUESS, JGUESS, MAT1, MAT2, MAT3, MAT4, DELI, DELJ, DET !/ ! ! -------------------------------------------------------------------- / ! IF ( PRESENT(DEBUG) ) THEN LDBG = DEBUG ELSE LDBG = .FALSE. END IF ! !-----handle point coincident with a cell vertex DO K=1,4 IF ( XT.EQ.XS(K) .AND. YT.EQ.YS(K) ) THEN SELECT CASE ( K ) CASE ( 1 ) IGUESS = ZERO; JGUESS = ZERO; CASE ( 2 ) IGUESS = ONE; JGUESS = ZERO; CASE ( 3 ) IGUESS = ONE; JGUESS = ONE; CASE ( 4 ) IGUESS = ZERO; JGUESS = ONE; END SELECT IF ( LDBG ) & WRITE(*,'(A,I3,2E14.6)') 'W3RMBL_R4 - COINCIDENT:', & K,IGUESS,JGUESS IF ( PRESENT(RW) ) THEN RW(1) = (ONE-IGUESS)*(ONE-JGUESS) RW(2) = IGUESS*(ONE-JGUESS) RW(3) = IGUESS*JGUESS RW(4) = (ONE-IGUESS)*JGUESS END IF IF ( PRESENT(IX) .AND. PRESENT(JX) ) THEN IX = IGUESS JX = JGUESS END IF RETURN END IF END DO ! !-----set iteration parameters and initial guess IF ( PRESENT(RW) ) RW = ZERO IGUESS = HALF JGUESS = HALF DYT = YT - YS(1) DY1 = YS(2) - YS(1) DY2 = YS(4) - YS(1) DY3 = YS(3) - YS(2) - DY2 DXT = XT - XS(1) DX1 = XS(2) - XS(1) DX2 = XS(4) - XS(1) DX3 = XS(3) - XS(2) - DX2 !-----iterate to find (i,j) for bilinear approximation ITER_LOOP: DO ITER=1,MAX_ITER DYP = DYT - DY1*IGUESS - DY2*JGUESS - DY3*IGUESS*JGUESS DXP = DXT - DX1*IGUESS - DX2*JGUESS - DX3*IGUESS*JGUESS MAT1 = DY1 + DY3*JGUESS MAT2 = DY2 + DY3*IGUESS MAT3 = DX1 + DX3*JGUESS MAT4 = DX2 + DX3*IGUESS DET = MAT1*MAT4 - MAT2*MAT3 DELI = (DYP*MAT4 - MAT2*DXP)/DET DELJ = (MAT1*DXP - DYP*MAT3)/DET IF ( LDBG ) & WRITE(*,'(A,I3,4E14.6)') 'W3RMBL_R4 - ITER:', & ITER,IGUESS,JGUESS,DELI,DELJ IF ( ABS(DELI) < CONVERGE .AND. & ABS(DELJ) < CONVERGE ) EXIT ITER_LOOP IGUESS = IGUESS + DELI JGUESS = JGUESS + DELJ END DO ITER_LOOP !-----if successful in finding (i,j), then compute weights IF ( ITER .LE. MAX_ITER ) THEN IF ( PRESENT(RW) ) THEN RW(1) = (ONE-IGUESS)*(ONE-JGUESS) RW(2) = IGUESS*(ONE-JGUESS) RW(3) = IGUESS*JGUESS RW(4) = (ONE-IGUESS)*JGUESS END IF IF ( PRESENT(IX) .AND. PRESENT(JX) ) THEN IX = IGUESS JX = JGUESS END IF ELSE WRITE(*,'(/A)') & 'W3RMBL_R4 -- ERROR: exceeded max iteration count' WRITE(*,'(A,2E14.6)') 'W3RMBL_R4 - DEST POINT COORDS: ',XT,YT DO K=1,4 WRITE(*,'(A,I1,A,2E14.6)') & 'W3RMBL_R4 - SRC POINT ',K,': ',XS(K),YS(K) END DO WRITE(*,'(A,2E14.6)') 'W3RMBL_R4 - CURRENT I,J: ',IGUESS,JGUESS CALL EXTCDE (1) END IF !(ITER.LE.MAX_ITER) !/ !/ End of W3RMBL_R4 -------------------------------------------------- / !/ END SUBROUTINE W3RMBL_R4 !/ ------------------------------------------------------------------- / SUBROUTINE W3RMBL_R8(XT, YT, XS, YS, RW, IX, JX, DEBUG) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 01-Dec-2010 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ 12-Nov-2010 : Add IX,IY optional arguments. !/ Implement r4 & r8 interfaces. ( version 3.14 ) !/ 01-Dec-2010 : Add check for target point coincident with a cell !/ vertex. Change to error exit when unable to !/ determine local (i,j). ( version 3.14 ) !/ ! 1. Purpose : ! ! Bilinear remapping for target point ( xt, yt ) in a cell defined ! by the source points ( xs(1:4), ys(1:4) ). Remapping weights are ! returned in rw(1:4). It is the caller's responsibility to ensure ! that the target point is located within the input cell and that the ! cell corner points are properly defined. ! ! (xs4,ys4) (xs3,ys3) ! _____________________ ! / / ! / x / ! / (xtin,ytin) / ! / / ! /____________________/ ! (xs1,ys1) (xs2,ys2) ! ! In spherical coordinates it is assumed that the longitude range of ! the target point is the same as that of the cell vertices. It is ! also assumed that the cell does not includes a branch cut. ! Double precision interface. ! ! 2. Method : ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! 7. Remarks : ! ! - Implementation is based on SCRIP. ! - In the case of spherical coordinates, the method results in ! bogus weights if enclosing cell contains a pole. ! ! 8. Structure : ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ REAL(8), INTENT(IN) :: XT REAL(8), INTENT(IN) :: YT REAL(8), INTENT(IN) :: XS(4) REAL(8), INTENT(IN) :: YS(4) REAL(8), INTENT(OUT), OPTIONAL :: RW(4) REAL(8), INTENT(OUT), OPTIONAL :: IX REAL(8), INTENT(OUT), OPTIONAL :: JX LOGICAL, INTENT(IN) , OPTIONAL :: DEBUG !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER, PARAMETER :: MAX_ITER = 10 REAL(8), PARAMETER :: CONVERGE = 1D-6 LOGICAL :: LDBG INTEGER :: K, ITER REAL(8) :: DXT, DX1, DX2, DX3, DXP, DYT, DY1, DY2, DY3, DYP REAL(8) :: IGUESS, JGUESS, MAT1, MAT2, MAT3, MAT4, DELI, DELJ, DET !/ ! ! -------------------------------------------------------------------- / ! IF ( PRESENT(DEBUG) ) THEN LDBG = DEBUG ELSE LDBG = .FALSE. END IF ! !-----handle point coincident with a cell vertex DO K=1,4 IF ( XT.EQ.XS(K) .AND. YT.EQ.YS(K) ) THEN SELECT CASE ( K ) CASE ( 1 ) IGUESS = ZERO; JGUESS = ZERO; CASE ( 2 ) IGUESS = ONE; JGUESS = ZERO; CASE ( 3 ) IGUESS = ONE; JGUESS = ONE; CASE ( 4 ) IGUESS = ZERO; JGUESS = ONE; END SELECT IF ( LDBG ) & WRITE(*,'(A,I3,2E14.6)') 'W3RMBL_R8 - COINCIDENT:', & K,IGUESS,JGUESS IF ( PRESENT(RW) ) THEN RW(1) = (ONE-IGUESS)*(ONE-JGUESS) RW(2) = IGUESS*(ONE-JGUESS) RW(3) = IGUESS*JGUESS RW(4) = (ONE-IGUESS)*JGUESS END IF IF ( PRESENT(IX) .AND. PRESENT(JX) ) THEN IX = IGUESS JX = JGUESS END IF RETURN END IF END DO ! !-----set iteration parameters and initial guess IF ( PRESENT(RW) ) RW = ZERO IGUESS = HALF JGUESS = HALF DYT = YT - YS(1) DY1 = YS(2) - YS(1) DY2 = YS(4) - YS(1) DY3 = YS(3) - YS(2) - DY2 DXT = XT - XS(1) DX1 = XS(2) - XS(1) DX2 = XS(4) - XS(1) DX3 = XS(3) - XS(2) - DX2 !-----iterate to find (i,j) for bilinear approximation ITER_LOOP: DO ITER=1,MAX_ITER DYP = DYT - DY1*IGUESS - DY2*JGUESS - DY3*IGUESS*JGUESS DXP = DXT - DX1*IGUESS - DX2*JGUESS - DX3*IGUESS*JGUESS MAT1 = DY1 + DY3*JGUESS MAT2 = DY2 + DY3*IGUESS MAT3 = DX1 + DX3*JGUESS MAT4 = DX2 + DX3*IGUESS DET = MAT1*MAT4 - MAT2*MAT3 DELI = (DYP*MAT4 - MAT2*DXP)/DET DELJ = (MAT1*DXP - DYP*MAT3)/DET IF ( LDBG ) & WRITE(*,'(A,I3,4E14.6)') 'W3RMBL_R8 - ITER:', & ITER,IGUESS,JGUESS,DELI,DELJ IF ( ABS(DELI) < CONVERGE .AND. & ABS(DELJ) < CONVERGE ) EXIT ITER_LOOP IGUESS = IGUESS + DELI JGUESS = JGUESS + DELJ END DO ITER_LOOP !-----if successful in finding (i,j), then compute weights IF ( ITER .LE. MAX_ITER ) THEN IF ( PRESENT(RW) ) THEN RW(1) = (ONE-IGUESS)*(ONE-JGUESS) RW(2) = IGUESS*(ONE-JGUESS) RW(3) = IGUESS*JGUESS RW(4) = (ONE-IGUESS)*JGUESS END IF IF ( PRESENT(IX) .AND. PRESENT(JX) ) THEN IX = IGUESS JX = JGUESS END IF ELSE WRITE(*,'(/A)') & 'W3RMBL_R8 -- ERROR: exceeded max iteration count' WRITE(*,'(A,2E14.6)') 'W3RMBL_R8 - DEST POINT COORDS: ',XT,YT DO K=1,4 WRITE(*,'(A,I1,A,2E14.6)') & 'W3RMBL_R8 - SRC POINT ',K,': ',XS(K),YS(K) END DO WRITE(*,'(A,2E14.6)') 'W3RMBL_R8 - CURRENT I,J: ',IGUESS,JGUESS CALL EXTCDE (1) END IF !(ITER.LE.MAX_ITER) !/ !/ End of W3RMBL_R8 -------------------------------------------------- / !/ END SUBROUTINE W3RMBL_R8 !/ ------------------------------------------------------------------- / FUNCTION W3DIST_R4(LLG, XT, YT, XS, YS) RESULT(DIST) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 12-Nov-2010 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ 14-Jun-2010 : Fix for ACOS argument > 1. ( version 3.14 ) !/ 12-Nov-2010 : Implement r4 & r8 interfaces. ( version 3.14 ) !/ ! 1. Purpose : ! ! Compute distance between two points. If spherical grid, then ! distance is the angle (in degrees) between the two points. ! Single precision interface. ! ! 2. Method : ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! !/T8 Enables NaN check. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Return parameter !/ REAL(4) :: DIST !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ LOGICAL, INTENT(IN) :: LLG REAL(4), INTENT(IN) :: XT, YT REAL(4), INTENT(IN) :: XS, YS !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ REAL(8) :: DX, DY, ARGD !/ ! ! -------------------------------------------------------------------- / ! !-----compute displacements DX = XT - XS DY = YT - YS IF ( LLG ) THEN !spherical coordinates !---------check for longitudinal branch cut crossing IF ( ABS(DX) .GT. D270 ) THEN DX = DX - SIGN(D360,DX) END IF !---------compute angular distance (min required for rare ! situation of acos(1+small) generating NaN) ARGD = MIN( ONE, COS(YT*D2R)*COS(YS*D2R)*COS(DX*D2R) & + SIN(YT*D2R)*SIN(YS*D2R) ) DIST = R2D*ACOS( ARGD ) ELSE !cartesian coordinates !---------compute cartesian distance DIST = SQRT( DX**2 + DY**2 ) END IF !cartesian coordinates !/ !/ End of W3DIST_R4 -------------------------------------------------- / !/ END FUNCTION W3DIST_R4 !/ ------------------------------------------------------------------- / FUNCTION W3DIST_R8(LLG, XT, YT, XS, YS) RESULT(DIST) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 12-Nov-2010 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ 14-Jun-2010 : Fix for ACOS argument > 1. ( version 3.14 ) !/ 12-Nov-2010 : Implement r4 & r8 interfaces. ( version 3.14 ) !/ ! 1. Purpose : ! ! Compute distance between two points. If spherical grid, then ! distance is the angle (in degrees) between the two points. ! Double precision interface. ! ! 2. Method : ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! !/T8 Enables NaN check. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Return parameter !/ REAL(8) :: DIST !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ LOGICAL, INTENT(IN) :: LLG REAL(8), INTENT(IN) :: XT, YT REAL(8), INTENT(IN) :: XS, YS !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ REAL(8) :: DX, DY, ARGD !/ ! ! -------------------------------------------------------------------- / ! !-----compute displacements DX = XT - XS DY = YT - YS IF ( LLG ) THEN !spherical coordinates !---------check for longitudinal branch cut crossing IF ( ABS(DX) .GT. D270 ) THEN DX = DX - SIGN(D360,DX) END IF !---------compute angular distance (min required for rare ! situation of acos(1+small) generating NaN) ARGD = MIN( ONE, COS(YT*D2R)*COS(YS*D2R)*COS(DX*D2R) & + SIN(YT*D2R)*SIN(YS*D2R) ) DIST = R2D*ACOS( ARGD ) ELSE !cartesian coordinates !---------compute cartesian distance DIST = SQRT( DX**2 + DY**2 ) END IF !cartesian coordinates !/ !/ End of W3DIST_R8 -------------------------------------------------- / !/ END FUNCTION W3DIST_R8 !/ ------------------------------------------------------------------- / FUNCTION W3CKCL_R4(LLG, XT, YT, NS, XS, YS, POLE, DEBUG) & RESULT(INCELL) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 06-Dec-2010 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ 12-Nov-2010 : Add subcell check for grid cell that includes a pole. !/ Implement r4 & r8 interfaces. ( version 3.14 ) !/ 06-Dec-2010 : Remove restriction on longitude range. Change ICLO !/ to integer and remove JCLO. ( version 3.14 ) !/ ! 1. Purpose : ! ! Check if point lies within grid cell. ! ! 2. Method : ! ! Calculates cross products for vertex to vertex (i.e. cell side) ! vs vertex to target. If all cross products have the same sign, ! the point is considered to be within the cell. Since they can ! be "all positive" *or* "all negative", there are no pre-conditions ! that the order of specification of the vertices be clockwise vs. ! counter-clockwise geographically. The logical variable POLE is ! set to true if the grid cell includes a pole. ! Single precision interface. ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! 7. Remarks : ! ! - For LL grids, this method assumes that the longitudes of point ! and grid cell vertices lie in the same range (i.e., both in [0:360] ! or [-180:180]). If the longitudes are not in the same range, then ! this method may result in a false positive. The burden is upon the ! caller to ensure that the longitude range of the point is the same ! as that of the grid cell vertices. ! - If enclosing cell includes a branch cut, then the coordinates of ! of the cell vertices AND the target point will be adjusted so ! that the branch cut is shifted 180 degrees. ! - If the enclosing cell includes a pole, then the cross-product check ! is performed for each quadrilateral subcell (with two cell ! vertices at the pole). ! ! 8. Structure : ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Return parameter !/ LOGICAL :: INCELL !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ LOGICAL, INTENT(IN) :: LLG REAL(4), INTENT(INOUT) :: XT, YT INTEGER, INTENT(IN) :: NS REAL(4), INTENT(INOUT) :: XS(NS), YS(NS) LOGICAL, INTENT(OUT):: POLE LOGICAL, INTENT(IN), OPTIONAL :: DEBUG !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ REAL(8), PARAMETER :: SMALL = 1D-6 LOGICAL :: LDBG, LSBC, BCUT INTEGER :: I, J, K, N REAL(8) :: XXT, YYT, XXS(NS), YYS(NS) REAL(8) :: V1X, V1Y, V2X, V2Y, S90 REAL(8) :: CROSS REAL(8) :: SIGN1 !/ ! ! -------------------------------------------------------------------- / ! INCELL = .TRUE. ! !-----must have >= 3 points to be a cell IF ( NS .LT. 3 ) THEN INCELL = .FALSE. RETURN END IF ! IF ( PRESENT(DEBUG) ) THEN LDBG = DEBUG ELSE LDBG = .FALSE. END IF ! !-----copy into locals XXT = XT; XXS = XS; YYT = YT; YYS = YS; ! !-----check if cell includes a pole or branch cut IF ( LLG ) THEN N = 0 !---------count longitudinal branch cut crossings DO I=1,NS J = MOD(I,NS) + 1 IF ( ABS(XXS(J)-XXS(I)) .GT. D180 ) N = N + 1 END DO !---------multiple longitudinal branch cut crossing => cell includes branch cut BCUT = N.GT.1 IF ( BCUT .AND. LDBG ) & WRITE(*,'(A)') 'W3CKCL_R4 - CELL INCLUDES A BRANCH CUT' !---------single longitudinal branch cut crossing ! or single vertex at 90 degrees => cell includes pole POLE = N.EQ.1 .OR. COUNT(ABS(YYS).EQ.D90).EQ.1 IF ( POLE .AND. LDBG ) & WRITE(*,'(A)') 'W3CKCL_R4 - CELL INCLUDES A POLE' ELSE POLE = .FALSE. BCUT = .FALSE. END IF ! !-----handle cell that includes a pole IF ( POLE ) THEN S90 = D90; IF ( MAXVAL(YS).LT.ZERO ) S90 = -D90; !---------perform cross-product check for each subcell SUBCELL_LOOP: DO I=1,NS LSBC = .TRUE. J = MOD(I,NS) + 1 DO K=1,4 SELECT CASE (K) CASE (1) !---------------------vector from (xi,yi) to (xj,yj) V1X = XXS(J) - XXS(I) V1Y = YYS(J) - YYS(I) !---------------------vector from (xi,yi) to (xt,yt) V2X = XXT - XXS(I) V2Y = YYT - YYS(I) CASE (2) !---------------------vector from (xj,yj) to (xj,90) V1X = XXS(J) - XXS(J) V1Y = S90 - YYS(J) !---------------------vector from (xj,yj) to (xt,yt) V2X = XXT - XXS(J) V2Y = YYT - YYS(J) CASE (3) !---------------------vector from (xj,90) to (xi,90) V1X = XXS(I) - XXS(J) V1Y = S90 - S90 !---------------------vector from (xj,90) to (xt,yt) V2X = XXT - XXS(J) V2Y = YYT - S90 CASE (4) !---------------------vector from (xi,90) to (xi,yi) V1X = XXS(I) - XXS(I) V1Y = YYS(I) - S90 !---------------------vector from (xi,90) to (xt,yt) V2X = XXT - XXS(I) V2Y = YYT - S90 END SELECT !-----------------check for longitudinal branch cut crossing IF ( ABS(V1X) .GT. D180 ) THEN V1X = V1X - SIGN(D360,V1X) END IF IF ( ABS(V2X) .GT. D180 ) THEN V2X = V2X - SIGN(D360,V2X) END IF !-----------------cross product CROSS = V1X*V2Y - V1Y*V2X IF ( LDBG ) & WRITE(*,'(A,3(I1,A),5E14.6)') 'W3CKCL_R4 - CROSS(', & I,',',J,',',K,'):',V1X,V1Y,V2X,V2Y,CROSS !-----------------if sign of cross product is not "unanimous" among the ! subcell sides, then target is outside the subcell IF ( K .EQ. 1 ) THEN SIGN1 = SIGN(ONE,CROSS) ELSE IF ( SIGN(ONE,CROSS) .NE. SIGN1 ) THEN LSBC = .FALSE. CYCLE SUBCELL_LOOP END IF END IF END DO !K IF ( LSBC ) RETURN END DO SUBCELL_LOOP INCELL = .FALSE. RETURN END IF !POLE ! !-----shift branch cut if necessary IF ( BCUT ) THEN IF ( MINVAL(XXS) .GE. ZERO ) THEN WHERE ( XXS .GT. D180 ) XXS = XXS - D360 IF ( XXT .GT. D180 ) XXT = XXT - D360 ELSE WHERE ( XXS .LT. ZERO ) XXS = XXS + D360 IF ( XXT .LT. ZERO ) XXT = XXT + D360 END IF IF ( LDBG ) THEN WRITE(*,'(A,2E14.6,4(/A,1I1,A,2E14.6))') & 'W3CKCL_R4 - SHIFT BRANCH CUT:',XXT,YYT, & (' CORNER(',K,'):',XXS(K),YYS(K),K=1,4) END IF END IF ! !-----check if target point lies outside cell bounding box IF ( XXT.LT.MINVAL(XXS) .OR. XXT.GT.MAXVAL(XXS) .OR. & YYT.LT.MINVAL(YYS) .OR. YYT.GT.MAXVAL(YYS) ) THEN IF ( LDBG ) THEN WRITE(*,'(A)') & 'W3CKCL_R4 - TARGET POINT LIES OUTSIDE CELL BOUNDING BOX' WRITE(*,'(A,2E14.6)') 'W3CKCL_R4 - TARGET: ',XXT,YYT WRITE(*,'(A,4E14.6)') 'W3CKCL_R4 - SOURCE: ', & MINVAL(XXS),MAXVAL(XXS),MINVAL(YYS),MAXVAL(YYS) END IF INCELL = .FALSE. RETURN END IF ! !-----perform cross-product cell check CORNER_LOOP: DO I=1,NS J = MOD(I,NS) + 1 !---------if target point is coincident a cell vertex, then ! exit cross-product check (flag as in cell) IF ( ABS(XXT-XXS(I)).LT.SMALL .AND. & & ABS(YYT-YYS(I)).LT.SMALL ) EXIT CORNER_LOOP !---------vector from (xi,yi) to (xj,yj) V1X = XXS(J) - XXS(I) V1Y = YYS(J) - YYS(I) !---------vector from (xi,yi) to (xt,yt) V2X = XXT - XXS(I) V2Y = YYT - YYS(I) !---------cross product CROSS = V1X*V2Y - V1Y*V2X IF ( LDBG ) & WRITE(*,'(A,2(I1,A),5E14.6)') 'W3CKCL_R4 - CROSS(', & I,',',J,'):',V1X,V1Y,V2X,V2Y,CROSS !---------if sign of cross product is not "unanimous" among the cell sides, ! then target is outside the cell IF ( I .EQ. 1 ) THEN SIGN1 = SIGN(ONE,CROSS) ELSE IF ( SIGN(ONE,CROSS) .NE. SIGN1 ) THEN INCELL = .FALSE. RETURN END IF END IF END DO CORNER_LOOP ! !-----return branch cut shifted coordinates IF ( BCUT ) THEN XT = XXT; XS = XXS; END IF !/ !/ End of W3CKCL_R4 -------------------------------------------------- / !/ END FUNCTION W3CKCL_R4 !/ ------------------------------------------------------------------- / FUNCTION W3CKCL_R8(LLG, XT, YT, NS, XS, YS, POLE, DEBUG) & RESULT(INCELL) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 06-Dec-2010 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ 12-Nov-2010 : Add subcell check for grid cell that includes a pole. !/ Implement r4 & r8 interfaces. ( version 3.14 ) !/ 06-Dec-2010 : Remove restriction on longitude range. Change ICLO !/ to integer and remove JCLO. ( version 3.14 ) !/ ! 1. Purpose : ! ! Check if point lies within grid cell. ! ! 2. Method : ! ! Calculates cross products for vertex to vertex (i.e. cell side) ! vs vertex to target. If all cross products have the same sign, ! the point is considered to be within the cell. Since they can ! be "all positive" *or* "all negative", there are no pre-conditions ! that the order of specification of the vertices be clockwise vs. ! counter-clockwise geographically. The logical variable POLE is ! set to true if the grid cell includes a pole. ! Double precision interface. ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! 7. Remarks : ! ! - For LL grids, this method assumes that the longitudes of point ! and grid cell vertices lie in the same range (i.e., both in [0:360] ! or [-180:180]). If the longitudes are not in the same range, then ! this method may result in a false positive. The burden is upon the ! caller to ensure that the longitude range of the point is the same ! as that of the grid cell vertices. ! - If enclosing cell includes a branch cut, then the coordinates of ! of the cell vertices AND the target point will be adjusted so ! that the branch cut is shifted 180 degrees. ! - If the enclosing cell includes a pole, then the cross-product check ! is performed for each quadrilateral subcell (with two cell ! vertices at the pole). ! ! 8. Structure : ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Return parameter !/ LOGICAL :: INCELL !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ LOGICAL, INTENT(IN) :: LLG REAL(8), INTENT(INOUT) :: XT, YT INTEGER, INTENT(IN) :: NS REAL(8), INTENT(INOUT) :: XS(NS), YS(NS) LOGICAL, INTENT(OUT):: POLE LOGICAL, INTENT(IN), OPTIONAL :: DEBUG !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ REAL(8), PARAMETER :: SMALL = 1D-6 LOGICAL :: LDBG, LSBC, BCUT INTEGER :: I, J, K, N REAL(8) :: XXT, YYT, XXS(NS), YYS(NS) REAL(8) :: V1X, V1Y, V2X, V2Y, S90 REAL(8) :: CROSS REAL(8) :: SIGN1 !/ ! ! -------------------------------------------------------------------- / ! INCELL = .TRUE. ! !-----must have >= 3 points to be a cell IF ( NS .LT. 3 ) THEN INCELL = .FALSE. RETURN END IF ! IF ( PRESENT(DEBUG) ) THEN LDBG = DEBUG ELSE LDBG = .FALSE. END IF ! !-----copy into locals XXT = XT; XXS = XS; YYT = YT; YYS = YS; ! !-----check if cell includes a pole or branch cut IF ( LLG ) THEN N = 0 !---------count longitudinal branch cut crossings DO I=1,NS J = MOD(I,NS) + 1 IF ( ABS(XXS(J)-XXS(I)) .GT. D180 ) N = N + 1 END DO !---------multiple longitudinal branch cut crossing => cell includes branch cut BCUT = N.GT.1 IF ( BCUT .AND. LDBG ) & WRITE(*,'(A)') 'W3CKCL_R8 - CELL INCLUDES A BRANCH CUT' !---------single longitudinal branch cut crossing ! or single vertex at 90 degrees => cell includes pole POLE = N.EQ.1 .OR. COUNT(ABS(YYS).EQ.D90).EQ.1 IF ( POLE .AND. LDBG ) & WRITE(*,'(A)') 'W3CKCL_R8 - CELL INCLUDES A POLE' ELSE POLE = .FALSE. BCUT = .FALSE. END IF ! !-----handle cell that includes a pole IF ( POLE ) THEN S90 = D90; IF ( MAXVAL(YS).LT.ZERO ) S90 = -D90; !---------perform cross-product check for each subcell SUBCELL_LOOP: DO I=1,NS LSBC = .TRUE. J = MOD(I,NS) + 1 DO K=1,4 SELECT CASE (K) CASE (1) !---------------------vector from (xi,yi) to (xj,yj) V1X = XXS(J) - XXS(I) V1Y = YYS(J) - YYS(I) !---------------------vector from (xi,yi) to (xt,yt) V2X = XXT - XXS(I) V2Y = YYT - YYS(I) CASE (2) !---------------------vector from (xj,yj) to (xj,90) V1X = XXS(J) - XXS(J) V1Y = S90 - YYS(J) !---------------------vector from (xj,yj) to (xt,yt) V2X = XXT - XXS(J) V2Y = YYT - YYS(J) CASE (3) !---------------------vector from (xj,90) to (xi,90) V1X = XXS(I) - XXS(J) V1Y = S90 - S90 !---------------------vector from (xj,90) to (xt,yt) V2X = XXT - XXS(J) V2Y = YYT - S90 CASE (4) !---------------------vector from (xi,90) to (xi,yi) V1X = XXS(I) - XXS(I) V1Y = YYS(I) - S90 !---------------------vector from (xi,90) to (xt,yt) V2X = XXT - XXS(I) V2Y = YYT - S90 END SELECT !-----------------check for longitudinal branch cut crossing IF ( ABS(V1X) .GT. D180 ) THEN V1X = V1X - SIGN(D360,V1X) END IF IF ( ABS(V2X) .GT. D180 ) THEN V2X = V2X - SIGN(D360,V2X) END IF !-----------------cross product CROSS = V1X*V2Y - V1Y*V2X IF ( LDBG ) & WRITE(*,'(A,3(I1,A),5E14.6)') 'W3CKCL_R4 - CROSS(', & I,',',J,',',K,'):',V1X,V1Y,V2X,V2Y,CROSS !-----------------if sign of cross product is not "unanimous" among the ! subcell sides, then target is outside the subcell IF ( K .EQ. 1 ) THEN SIGN1 = SIGN(ONE,CROSS) ELSE IF ( SIGN(ONE,CROSS) .NE. SIGN1 ) THEN LSBC = .FALSE. CYCLE SUBCELL_LOOP END IF END IF END DO !K IF ( LSBC ) RETURN END DO SUBCELL_LOOP INCELL = .FALSE. RETURN END IF !POLE ! !-----shift branch cut if necessary IF ( BCUT ) THEN IF ( MINVAL(XXS) .GE. ZERO ) THEN WHERE ( XXS .GT. D180 ) XXS = XXS - D360 IF ( XXT .GT. D180 ) XXT = XXT - D360 ELSE WHERE ( XXS .LT. ZERO ) XXS = XXS + D360 IF ( XXT .LT. ZERO ) XXT = XXT + D360 END IF IF ( LDBG ) THEN WRITE(*,'(A,2E14.6,4(/A,1I1,A,2E14.6))') & 'W3CKCL_R8 - SHIFT BRANCH CUT:',XXT,YYT, & (' CORNER(',K,'):',XXS(K),YYS(K),K=1,4) END IF END IF ! !-----check if target point lies outside cell bounding box IF ( XXT.LT.MINVAL(XXS) .OR. XXT.GT.MAXVAL(XXS) .OR. & YYT.LT.MINVAL(YYS) .OR. YYT.GT.MAXVAL(YYS) ) THEN IF ( LDBG ) THEN WRITE(*,'(A)') & 'W3CKCL_R8 - TARGET POINT LIES OUTSIDE CELL BOUNDING BOX' WRITE(*,'(A,2E14.6)') 'W3CKCL_R8 - TARGET: ',XXT,YYT WRITE(*,'(A,4E14.6)') 'W3CKCL_R8 - SOURCE: ', & MINVAL(XXS),MAXVAL(XXS),MINVAL(YYS),MAXVAL(YYS) END IF INCELL = .FALSE. RETURN END IF ! !-----perform cross-product cell check CORNER_LOOP: DO I=1,NS J = MOD(I,NS) + 1 !---------if target point is coincident a cell vertex, then ! exit cross-product check (flag as in cell) IF ( ABS(XXT-XXS(I)).LT.SMALL .AND. & & ABS(YYT-YYS(I)).LT.SMALL ) EXIT CORNER_LOOP !---------vector from (xi,yi) to (xj,yj) V1X = XXS(J) - XXS(I) V1Y = YYS(J) - YYS(I) !---------vector from (xi,yi) to (xt,yt) V2X = XXT - XXS(I) V2Y = YYT - YYS(I) !---------cross product CROSS = V1X*V2Y - V1Y*V2X IF ( LDBG ) & WRITE(*,'(A,2(I1,A),5E14.6)') 'W3CKCL_R8 - CROSS(', & I,',',J,'):',V1X,V1Y,V2X,V2Y,CROSS !---------if sign of cross product is not "unanimous" among the cell sides, ! then target is outside the cell IF ( I .EQ. 1 ) THEN SIGN1 = SIGN(ONE,CROSS) ELSE IF ( SIGN(ONE,CROSS) .NE. SIGN1 ) THEN INCELL = .FALSE. RETURN END IF END IF END DO CORNER_LOOP ! !-----return branch cut shifted coordinates IF ( BCUT ) THEN XT = XXT; XS = XXS; END IF !/ !/ End of W3CKCL_R8 -------------------------------------------------- / !/ END FUNCTION W3CKCL_R8 !/ ------------------------------------------------------------------- / SUBROUTINE W3SORT_R4(N, I, J, D) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 12-Nov-2010 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ 12-Nov-2010 : Implement r4 & r8 interfaces. ( version 3.14 ) !/ ! 1. Purpose : ! ! Sort input arrays in increasing order according to input array D. ! Single precision interface. ! ! 2. Method : ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ INTEGER, INTENT(IN) :: N INTEGER, INTENT(INOUT) :: I(N) INTEGER, INTENT(INOUT) :: J(N) REAL(4), INTENT(INOUT) :: D(N) !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: K, L, IM, JM REAL(4) :: DM !/ ! ! -------------------------------------------------------------------- / ! DO K=1, N-1 DO L=K+1, N IF ( D(L) .LT. D(K) ) THEN IM = I(K); JM = J(K); DM = D(K); I(K) = I(L); J(K) = J(L); D(K) = D(L); I(L) = IM; J(L) = JM; D(L) = DM; END IF END DO !L END DO !K !/ !/ End of W3SORT_R4 -------------------------------------------------- / !/ END SUBROUTINE W3SORT_R4 !/ ------------------------------------------------------------------- / SUBROUTINE W3SORT_R8(N, I, J, D) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 12-Nov-2010 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ 12-Nov-2010 : Implement r4 & r8 interfaces. ( version 3.14 ) !/ ! 1. Purpose : ! ! Sort input arrays in increasing order according to input array D. ! Double precision interface. ! ! 2. Method : ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ INTEGER, INTENT(IN) :: N INTEGER, INTENT(INOUT) :: I(N) INTEGER, INTENT(INOUT) :: J(N) REAL(8), INTENT(INOUT) :: D(N) !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: K, L, IM, JM REAL(8) :: DM !/ ! ! -------------------------------------------------------------------- / ! DO K=1, N-1 DO L=K+1, N IF ( D(L) .LT. D(K) ) THEN IM = I(K); JM = J(K); DM = D(K); I(K) = I(L); J(K) = J(L); D(K) = D(L); I(L) = IM; J(L) = JM; D(L) = DM; END IF END DO !L END DO !K !/ !/ End of W3SORT_R8 -------------------------------------------------- / !/ END SUBROUTINE W3SORT_R8 !/ ------------------------------------------------------------------- / SUBROUTINE W3ISRT_R4(II, JJ, DD, N, I, J, D) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 12-Nov-2010 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ 12-Nov-2010 : Implement r4 & r8 interfaces. ( version 3.14 ) !/ ! 1. Purpose : ! ! Insert DD data into D at location where DD < D(K). ! Single precision interface. ! ! 2. Method : ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ INTEGER, INTENT(IN) :: II INTEGER, INTENT(IN) :: JJ REAL(4), INTENT(IN) :: DD INTEGER, INTENT(IN) :: N INTEGER, INTENT(INOUT) :: I(N) INTEGER, INTENT(INOUT) :: J(N) REAL(4), INTENT(INOUT) :: D(N) !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: K, L !/ ! ! -------------------------------------------------------------------- / ! K_LOOP: DO K=1,N IF ( DD .LT. D(K) ) THEN !---------right-shift list (>= k) DO L=N,K+1,-1 I(L) = I(L-1); J(L) = J(L-1); D(L) = D(L-1); END DO !L !---------insert point into list at k I(K) = II; J(K) = JJ; D(K) = DD; EXIT K_LOOP END IF !dd.lt.d(k) END DO K_LOOP !/ !/ End of W3ISRT_R4 -------------------------------------------------- / !/ END SUBROUTINE W3ISRT_R4 !/ ------------------------------------------------------------------- / SUBROUTINE W3ISRT_R8(II, JJ, DD, N, I, J, D) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 12-Nov-2010 | !/ +-----------------------------------+ !/ !/ 30-Oct-2009 : Origination. ( version 3.14 ) !/ 12-Nov-2010 : Implement r4 & r8 interfaces. ( version 3.14 ) !/ ! 1. Purpose : ! ! Insert DD data into D at location where DD < D(K). ! Double precision interface. ! ! 2. Method : ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ INTEGER, INTENT(IN) :: II INTEGER, INTENT(IN) :: JJ REAL(8), INTENT(IN) :: DD INTEGER, INTENT(IN) :: N INTEGER, INTENT(INOUT) :: I(N) INTEGER, INTENT(INOUT) :: J(N) REAL(8), INTENT(INOUT) :: D(N) !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: K, L !/ ! ! -------------------------------------------------------------------- / ! K_LOOP: DO K=1,N IF ( DD .LT. D(K) ) THEN !---------right-shift list (>= k) DO L=N,K+1,-1 I(L) = I(L-1); J(L) = J(L-1); D(L) = D(L-1); END DO !L !---------insert point into list at k I(K) = II; J(K) = JJ; D(K) = DD; EXIT K_LOOP END IF !dd.lt.d(k) END DO K_LOOP !/ !/ End of W3ISRT_R8 -------------------------------------------------- / !/ END SUBROUTINE W3ISRT_R8 !/ ------------------------------------------------------------------- / FUNCTION W3INAN_R4(X) RESULT(INAN) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 14-Jun-2010 | !/ +-----------------------------------+ !/ !/ 14-Jun-2010 : Origination. ( version 3.14 ) !/ ! 1. Purpose : ! ! Return TRUE if input is infinite or NaN (not a number). ! Single precision interface. ! ! 2. Method : ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Return parameter !/ LOGICAL :: INAN !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ REAL(4), INTENT(IN) :: X !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ ! ! -------------------------------------------------------------------- / ! !-----return true if X is NaN or +Inf or -Inf INAN = .NOT. ( X .GE. -HUGE(X) .AND. X .LE. HUGE(X) ) !/ !/ End of W3INAN_R4 -------------------------------------------------- / !/ END FUNCTION W3INAN_R4 !/ ------------------------------------------------------------------- / FUNCTION W3INAN_R8(X) RESULT(INAN) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | T. J. Campbell, NRL | !/ | FORTRAN 90 | !/ | Last update : 14-Jun-2010 | !/ +-----------------------------------+ !/ !/ 14-Jun-2010 : Origination. ( version 3.14 ) !/ ! 1. Purpose : ! ! Return TRUE if input is infinite or NaN (not a number). ! Double precision interface. ! ! 2. Method : ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ Return parameter !/ LOGICAL :: INAN !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ REAL(8), INTENT(IN) :: X !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ ! ! -------------------------------------------------------------------- / ! !-----return true if X is NaN or +Inf or -Inf INAN = .NOT. ( X .GE. -HUGE(X) .AND. X .LE. HUGE(X) ) !/ !/ End of W3INAN_R8 -------------------------------------------------- / !/ END FUNCTION W3INAN_R8 !/ ------------------------------------------------------------------- / !/ !/ End of module W3GSRUMD -------------------------------------------- / !/ END MODULE W3GSRUMD