! MODULES containting data structures and subroutines for ! spatial and temporal tracking (part of wave partitioning). ! ! Contents of this file ! ! PARTDAT data structures used for spatial/temporal tracking ! ! MATFUNC simple data maniputation (based on Matlab intr. functions) ! contains: ! UNIQUE removes duplicate reals from an vector ! SORT sorts the vector in ascending or descending order ! SETDIFF returns elements in vector1 that are not in vector2 ! INTERSECT returns elements that are mutual in vector1 and vector2 ! UNION returns the union of vector1 and vector2 ! func. LENGTH finds no. of indices in vector not filled with blank entries. ! func. FINDFIRST returns index of first instance of a search value in vector ! func. STD computes standard deviation ! ! PARTFUNC auxiliary subroutines and functions for tracking ! contains: ! findWay find direction and no. steps in spatial search spiral ! findNext find next point on spatial search spiral ! findSys find all neighboring wave systems for given grid point ! combineWaveSystems combine wave systems, then remove small and low-energy systems ! printFinalSys output the final output systems for this time step ! combineSys combine wave systems ! combinePartitionsV2 combine two partitions that have been assigned to the same system ! checkPoint find whether a lon/lat point has already been processed ! func. mean_angleV2 compute the mean direction from array of directions ! findIJV3 Find indices of system "a" that lie over or next to system "b" ! ! PARTMAIN main subroutine for spatial/temporal tracking ! contains: ! waveTracking_NWS_V2 main subroutine of spatial and temporal tracking algorithm ! spiralTrackV3 performs the spatial spiral tracking for a given time step ! timeTrackingV2 performs the time tracking of all wave systems ! ! MODULE PARTDAT ! ! 0. Authors ! ! 40.81: Andre van der Westhuysen and Jeff Hanson ! ! 1. Updates ! ! 40.81, Jul. 11: Origination, based on Matlab code by Jeff Hanson ! and Eve-Marie Devaliere ! ! 2. Purpose ! ! Data structures for spatial and temporal tracking of wave systems (partitions) ! IMPLICIT NONE ! ! NOTE: In Fortran 90/95 derived types cannot contain allocatable arrays. ! The same functionality is achieved here using pointers (pointing ! to unnamed allocatable arrays). Can be replaced by allocatable arrays ! when transitioning to the Fortran 2003 standard. ! ! param Der. type structure of basic spectrally partitioned results at a geo point ! hs Real arr array of sign. wave height partitions ! tp Real arr array of peak period partitions ! dir Real arr array of mena direction partitions ! ipart Int arr array of partition indices ! sys Int arr array of system indices to which a given partition has been assigned ! ngbrSys Int arr array of system indices of neighboring grid points ! checked Int 0 = geo point not checked yet (in SUBROUTINE findSys) ! 1 = geo point has been checked ! -1 = geo land point (i.e. no partitioning data found for this point) ! -2 = geo land point, second passing. ! TYPE param REAL :: hs(10) 40.PAR !Make allocatable REAL :: tp(10) REAL :: dir(10) INTEGER :: ipart(10) INTEGER :: sys(10) INTEGER :: ngbrSys(50) INTEGER :: checked END TYPE param ! ! wind Der. type structure containing wind-related parameters ! wdir Real wind direction at grid point (Nautical or Cartesian, invariant) ! wspd Real wind speed at grid point ! TYPE wind REAL :: wdir REAL :: wspd END TYPE wind ! ! dat1d Der. type 1d data structure for reading "partRes" input file with raw ! partition data (Redundant if this data is stored in memory) ! lat1d Real arr 1d array of latitudes of input partitioned data ! lon1d Real arr 1d array of longitudes of input partitioned data ! par1d type(param) arr 1d array of partitioned parameter structures ! wnd1d type(wind) arr 1d array of wind parameter structures ! TYPE dat1d REAL, POINTER :: lat1d(:) REAL, POINTER :: lon1d(:) TYPE(param), POINTER :: par1d(:) TYPE(wind), POINTER :: wnd1d(:) END TYPE dat1d ! ! dat2d Der. type 2d data structure for storing raw partitioned data ! lat Real arr 2d array of latitudes of input partitioned data ! lon Real arr 2d array of longitudes of input partitioned data ! par type(param) arr 2d array of partitioned parameter structures ! wnd type(wind) arr 2d array of wind parameter structures ! maxi Int size of 2d array of raw partitioned data in i dimension ! maxj Int size of 2d array of raw partitioned data in j dimension ! TYPE dat2d REAL, POINTER :: lat(:,:) REAL, POINTER :: lon(:,:) TYPE(param), POINTER :: par(:,:) TYPE(wind), POINTER :: wnd(:,:) INTEGER :: maxi INTEGER :: maxj END TYPE dat2d ! ! neighbr Der. type structure for storing data of neighboring grid point ! par type(param) partitioned parameter structure at neighboring grid point ! i Int i index of neighboring grid point ! j Int j index of neighboring grid point ! TYPE neighbr TYPE(param) :: par INTEGER :: i INTEGER :: j END TYPE neighbr ! ! mtchsys Der. type structure for storing data of matched systems ! sysVal Int arr array of indices of matched systems ! tpVal Real arr array of peak period values of matched systems ! TYPE mtchsys INTEGER :: sysVal(50) REAL :: tpVal(50) END TYPE mtchsys ! ! system Der. type structure for storing spatially tracked systems (one time level) ! hs Real arr sign wave height field assoc with wave system (in 1d array) ! tp Real arr peak period field assoc with wave system (in 1d array) ! dir Real arr mean direction field assoc with wave system (in 1d array) ! i Int arr i index of geo grid point in wave system (in 1d array) ! j Int arr j index of geo grid point in wave system (in 1d array) ! lat Real arr latitudes of grid point in wave system (in 1d array) ! lon Real arr longitudes of grid point in wave system (in 1d array) ! sysInd Int index of current wave system ! hsMean Real spatial mean sign wave height of current wave system ! tpMean Real spatial mean peak period of current wave system ! dirMean Real spatial mean wave direction of current wave system ! nPoints Int total number of grid points in current wave system ! ngbr Int arr indices of neighboring wave systems ! grp Int time-tracked group that system is assigned to ! TYPE system REAL, POINTER :: hs(:) REAL, POINTER :: tp(:) REAL, POINTER :: dir(:) INTEGER, POINTER :: i(:) INTEGER, POINTER :: j(:) REAL, POINTER :: lat(:) REAL, POINTER :: lon(:) INTEGER :: sysInd REAL :: hsMean REAL :: tpMean REAL :: dirMean INTEGER :: nPoints INTEGER :: ngbr(1000) INTEGER :: grp END TYPE system ! ! timsys Der. type structure for storing time-tracked systems (all time levels) ! sys type(system) arr array of all spatially+temporally tracked systems at given ! time level ! TYPE timsys TYPE(system), POINTER :: sys(:) END TYPE timsys END MODULE PARTDAT MODULE MATFUNC CONTAINS !*********************************************************************** ! SUBROUTINE UNIQUE (INARRAY,INSIZE,OUTARRAY,OUTSIZE) ! * !*********************************************************************** ! IMPLICIT NONE ! ! 0. Authors ! ! 40.81: Andre van der Westhuysen ! ! 2. Purpose ! ! Replicates the Matlab intrinsic function UNIQUE: ! (i) Removes duplicate reals from an array. ! (ii) Sort the resulting array in ascending order. Based on code by John Mahaffy. ! ! 3. Method ! ! Based on http://rosettacode.org/wiki/Remove_duplicate_elements#Fortran ! ! 4. Argument variables ! ! INARRAY REAL ARR input Input array ! INSIZE INTEGER input Size of input array ! OUTARRAY REAL ARR output Output array ! OUTSIZE INTEGER output Size of output array (number of unique elements) ! INTEGER :: INSIZE, OUTSIZE, I, J REAL :: INARRAY(INSIZE), TEMP(INSIZE) REAL, ALLOCATABLE :: OUTARRAY(:) LOGICAL :: FOUND INTENT (IN) INARRAY, INSIZE INTENT (OUT) OUTARRAY, OUTSIZE ! ! 5. Local variables ! ! INARRAY - array of values within which to find unique values ! IY - array to be carried with X (all swaps of X elements are ??? EDIT! ! matched in IY . After the sort IY(J) contains the original ! postition of the value X(J) in the unsorted X array. ! N - number of values in array X to be sorted ! INTEGER N INTEGER IY(INSIZE) REAL TEMP2 INTEGER JMAX, ITEMP ! ! 6. Subroutines used ! ! - ! ! 7. Subroutines calling ! ! waveTracking_NWS_V2 ! findSys ! printFinalSys ! combineSys ! UNION ! ! 8. Error messages ! ! - ! ! 9. Remarks ! ! 10. Structure ! ! 11. Source text ! ! NOTE: This routine can be improved further by removing "IF ( FOUND ) CYCLE" ! Initialize TEMP array with dummy values TEMP(:) = 9999. OUTSIZE = 1 TEMP(1) = INARRAY(1) DO I=2,SIZE(INARRAY) FOUND = .FALSE. DO J=1,OUTSIZE IF (TEMP(J) == INARRAY(I)) THEN ! Found a match so start looking again FOUND = .TRUE. END IF END DO ! No match found so add it to the output IF ( FOUND ) CYCLE OUTSIZE = OUTSIZE + 1 TEMP(OUTSIZE) = INARRAY(I) END DO ! Allocate OUTARRAY to the size of the unique array ALLOCATE(OUTARRAY(OUTSIZE)) OUTARRAY = TEMP(1:OUTSIZE) ! Sort OUTARRAY in ascending order N = SIZE(OUTARRAY) JMAX=N-1 DO 200 I=1,N-1 TEMP2=1.E38 DO 100 J=1,JMAX IF(OUTARRAY(J).LT.OUTARRAY(J+1)) GO TO 100 TEMP2=OUTARRAY(J) OUTARRAY(J)=OUTARRAY(J+1) OUTARRAY(J+1)=TEMP2 ITEMP=IY(J) IY(J)=IY(J+1) IY(J+1)=ITEMP 100 CONTINUE JMAX=JMAX-1 200 CONTINUE RETURN ! * end of subroutine UNIQUE * END SUBROUTINE UNIQUE !*********************************************************************** ! SUBROUTINE SORT (INARRAY,INSIZE,OUTARRAY,IY,DIRECTION) ! * !*********************************************************************** ! IMPLICIT NONE ! ! 0. Authors ! ! 40.81: Andre van der Westhuysen ! ! 2. Purpose ! ! Replicates the Matlab intrinsic function SORT: ! Sorts the array INARRAY in ascending (Direction = 'A') or ! descending (Direciton = 'D') order. The sorted array is ! stored in OUTARRAY, and the sorted array of the original ! indices is stored in IY. Based on code by John Mahaffy. ! ! 3. Method ! ! Based on http://rosettacode.org/wiki/Remove_duplicate_elements#Fortran ! ! 4. Argument variables ! ! INARRAY REAL ARR input Input array ! INSIZE INTEGER input Size of input array ! OUTARRAY REAL ARR output Sorted output array ! IY INTEGER ARR output Sorted array of the original indices CHARACTER :: DIRECTION *1 INTEGER :: INSIZE, I, J INTEGER :: IY(INSIZE) REAL :: INARRAY(INSIZE), OUTARRAY(INSIZE) REAL :: TEMP(INSIZE) LOGICAL :: FOUND INTENT (IN) INARRAY, INSIZE, DIRECTION INTENT (OUT) OUTARRAY, IY ! ! 5. Local variables ! ! INARRAY - array of values to be sorted ! IY - array to be carried with X (all swaps of X elements are ??? EDIT! ! matched in IY . After the sort IY(J) contains the original ! postition of the value X(J) in the unsorted X array. ! N - number of values in array X to be sorted INTEGER N REAL TEMP2 INTEGER JMAX, ITEMP ! ! 6. Subroutines used ! ! - ! ! 7. Subroutines calling ! ! printFinalSys ! combineSys ! timeTrackingV2 ! ! 8. Error messages ! ! - ! ! 9. Remarks ! ! 10. Structure ! ! 11. Source text ! ! Sort OUTARRAY in as/decending order OUTARRAY = INARRAY N = INSIZE JMAX=N-1 DO I = 1, N IY(I) = I END DO IF (DIRECTION .EQ. 'A') THEN DO 200 I=1,N-1 ! TEMP2=1.E38 DO 100 J=1,JMAX IF(OUTARRAY(J).LE.OUTARRAY(J+1)) GO TO 100 TEMP2=OUTARRAY(J) OUTARRAY(J)=OUTARRAY(J+1) OUTARRAY(J+1)=TEMP2 ITEMP=IY(J) IY(J)=IY(J+1) IY(J+1)=ITEMP 100 CONTINUE JMAX=JMAX-1 200 CONTINUE ELSE IF (DIRECTION .EQ. 'D') THEN DO 201 I=1,N-1 ! TEMP2=1.E38 DO 101 J=1,JMAX IF(OUTARRAY(J).GE.OUTARRAY(J+1)) GO TO 101 TEMP2=OUTARRAY(J) OUTARRAY(J)=OUTARRAY(J+1) OUTARRAY(J+1)=TEMP2 ITEMP=IY(J) IY(J)=IY(J+1) IY(J+1)=ITEMP 101 CONTINUE JMAX=JMAX-1 201 CONTINUE END IF RETURN ! * end of subroutine SORT * END SUBROUTINE SORT !*********************************************************************** ! SUBROUTINE SETDIFF (INARRAY1, INSIZE1, INARRAY2, INSIZE2, & OUTARRAY, OUTSIZE) ! * !*********************************************************************** ! IMPLICIT NONE ! ! 0. Authors ! ! 40.81: Andre van der Westhuysen ! ! 2. Purpose ! ! Replicates the Matlab intrinsic function SETDIFF: ! (i) Returns the elements in INARRAY1 that are not in INARRAY2. ! (ii) Sort the resulting array in ascending order. ! ! 3. Method ! ! 4. Argument variables ! ! INARRAY1 REAL ARR input Input array ! INSIZE1 INTEGER input Size of input array ! INARRAY2 REAL ARR input Input array ! INSIZE2 INTEGER input Size of input array ! OUTARRAY REAL ARR output Output array ! OUTSIZE INTEGER output Size of output array (number of unique elements) INTEGER :: INSIZE1, INSIZE2, OUTSIZE, I, J REAL :: INARRAY1(INSIZE1), INARRAY2(INSIZE2) REAL :: TEMP(INSIZE1) REAL, ALLOCATABLE :: OUTARRAY(:) LOGICAL :: FOUND INTENT (IN) INARRAY1, INSIZE1, INARRAY2, INSIZE2 INTENT (OUT) OUTARRAY, OUTSIZE ! ! 5. Local variables ! ! INARRAY1, INARRAY2 - arrays to find the difference between ! IY - array to be carried with X (all swaps of X elements are ??? EDIT! ! matched in IY . After the sort IY(J) contains the original ! postition of the value X(J) in the unsorted X array. ! N - number of values in array X to be sorted INTEGER N INTEGER IY(INSIZE1) REAL TEMP2 INTEGER JMAX, ITEMP ! ! 6. Subroutines used ! ! - ! ! 7. Subroutines calling ! ! printFinalSys ! combineSys ! timeTrackingV2 ! ! 8. Error messages ! ! - ! ! 9. Remarks ! ! 10. Structure ! ! 11. Source text ! ! NOTE: This routine can be improved further by removing "IF ( FOUND ) CYCLE" ! Initialize TEMP array with dummy values TEMP(:) = 9999. OUTSIZE = 0 DO I=1,SIZE(INARRAY1) ! FOUND = .FALSE. DO J=1,SIZE(INARRAY2) IF (INARRAY2(J).EQ.INARRAY1(I)) THEN ! Found a match ! FOUND = .TRUE. GOTO 90 END IF END DO IF (.NOT.ANY(TEMP(:).EQ.INARRAY1(I))) THEN OUTSIZE = OUTSIZE + 1 TEMP(OUTSIZE) = INARRAY1(I) END IF 90 CONTINUE END DO ! Allocate OUTARRAY to the size of the unique array ALLOCATE(OUTARRAY(OUTSIZE)) OUTARRAY = TEMP(1:OUTSIZE) ! Sort OUTARRAY in ascending order N = SIZE(OUTARRAY) JMAX=N-1 DO 200 I=1,N-1 TEMP2=1.E38 DO 100 J=1,JMAX IF(OUTARRAY(J).LT.OUTARRAY(J+1)) GO TO 100 TEMP2=OUTARRAY(J) OUTARRAY(J)=OUTARRAY(J+1) OUTARRAY(J+1)=TEMP2 ITEMP=IY(J) IY(J)=IY(J+1) IY(J+1)=ITEMP 100 CONTINUE JMAX=JMAX-1 200 CONTINUE RETURN ! * end of subroutine SETDIFF * END SUBROUTINE SETDIFF !*********************************************************************** ! SUBROUTINE INTERSECT (INARRAY1 ,INSIZE1 ,INARRAY2 ,INSIZE2 , & OUTARRAY ,OUTSIZE ,IND1 ,IND2 ) ! * !*********************************************************************** ! IMPLICIT NONE ! ! 0. Authors ! ! 40.81: Andre van der Westhuysen ! ! 2. Purpose ! ! Replicates the Matlab intrinsic function INTERSECT: ! (i) Returns the elements that are mutual in INARRAY1 and INARRAY2. ! (ii) Sort the resulting array in ascending order. ! ! 3. Method ! ! 4. Argument variables ! ! INARRAY1 REAL ARR input Input array ! INSIZE1 INTEGER input Size of input array ! INARRAY2 REAL ARR input Input array ! INSIZE2 INTEGER input Size of input array ! OUTARRAY REAL ARR output Output array ! OUTSIZE INTEGER output Size of output array (number of unique elements) INTEGER :: INSIZE1, INSIZE2, OUTSIZE REAL :: INARRAY1(INSIZE1), INARRAY2(INSIZE2) REAL, ALLOCATABLE :: OUTARRAY(:) REAL, ALLOCATABLE :: IND1(:), IND2(:) INTENT (IN) INARRAY1, INSIZE1, INARRAY2, INSIZE2 INTENT (OUT) OUTARRAY, OUTSIZE, IND1, IND2 ! ! 5. Local variables ! ! INARRAY1, INARRAY2 - arrays for which to find the intersection ! IY - array to be carried with X (all swaps of X elements are ??? EDIT! ! matched in IY . After the sort IY(J) contains the original ! postition of the value X(J) in the unsorted X array. ! N - number of values in array X to be sorted REAL :: TEMP(MIN(INSIZE1,INSIZE2)) REAL :: TEMP2 INTEGER :: N, NN, I, J, JMAX ! ! 6. Subroutines used ! ! - ! ! 7. Subroutines calling ! ! findIJV3 ! ! 8. Error messages ! ! - ! ! 9. Remarks ! ! 10. Structure ! ! 11. Source text ! ! NOTE: This routine can be improved further by removing "IF ( FOUND ) CYCLE" ! Initialize TEMP array with dummy values TEMP(:) = 9999. OUTSIZE = 0 DO I=1,SIZE(INARRAY1) DO J=1,SIZE(INARRAY2) IF (INARRAY2(J).EQ.INARRAY1(I)) THEN ! Found a match IF (.NOT.ANY(TEMP(:).EQ.INARRAY1(I))) THEN OUTSIZE = OUTSIZE + 1 TEMP(OUTSIZE) = INARRAY1(I) GOTO 90 END IF END IF END DO 90 CONTINUE END DO ! Allocate OUTARRAY to the size of the unique array ALLOCATE(OUTARRAY(OUTSIZE)) ALLOCATE(IND1(OUTSIZE)) ALLOCATE(IND2(OUTSIZE)) OUTARRAY = TEMP(1:OUTSIZE) ! Sort OUTARRAY in ascending order N = SIZE(OUTARRAY) JMAX=N-1 DO 200 I=1,N-1 TEMP2=1.E38 DO 100 J=1,JMAX IF(OUTARRAY(J).LT.OUTARRAY(J+1)) GO TO 100 TEMP2=OUTARRAY(J) OUTARRAY(J)=OUTARRAY(J+1) OUTARRAY(J+1)=TEMP2 100 CONTINUE JMAX=JMAX-1 200 CONTINUE ! Find orginal indices of the values in OUTARRAY ! in INARRAY1 and INARRAY2 DO N = 1, OUTSIZE DO NN = 1, INSIZE1 IF (INARRAY1(NN).EQ.OUTARRAY(N)) IND1(N)=NN END DO DO NN = 1, INSIZE2 IF (INARRAY2(NN).EQ.OUTARRAY(N)) IND2(N)=NN END DO END DO RETURN ! * end of subroutine INTERSECT * END SUBROUTINE INTERSECT !*********************************************************************** ! SUBROUTINE UNION (INARRAY1, INSIZE1, INARRAY2, INSIZE2, & OUTARRAY, OUTSIZE) ! * !*********************************************************************** ! IMPLICIT NONE ! ! 0. Authors ! ! 40.81: Andre van der Westhuysen ! ! 2. Purpose ! ! Replicates the Matlab intrinsic function UNION: ! (i) Returns the union of INARRAY1 and INARRAY2. ! (ii) Sort the resulting array in ascending order. ! ! 3. Method ! ! ! 4. Argument variables ! ! INARRAY1 REAL ARR input Input array ! INSIZE1 INTEGER input Size of input array ! INARRAY2 REAL ARR input Input array ! INSIZE2 INTEGER input Size of input array ! OUTARRAY REAL ARR output Output array ! OUTSIZE INTEGER output Size of output array INTEGER :: INSIZE1, INSIZE2, OUTSIZE REAL :: INARRAY1(INSIZE1), INARRAY2(INSIZE2) REAL, ALLOCATABLE :: OUTARRAY(:) INTENT (IN) INARRAY1, INSIZE1, INARRAY2, INSIZE2 INTENT (OUT) OUTARRAY, OUTSIZE ! ! 5. Local variables ! REAL :: COMBINE(INSIZE1+INSIZE2) ! ! 6. Subroutines used ! ! UNIQUE ! ! 7. Subroutines calling ! ! combineSys ! ! 8. Error messages ! ! - ! ! 9. Remarks ! ! 10. Structure ! ! ! 11. Source text ! IF ((INSIZE1+INSIZE2).GT.0) THEN COMBINE = (/INARRAY1,INARRAY2/) CALL UNIQUE(REAL(COMBINE),SIZE(COMBINE),OUTARRAY,OUTSIZE) ELSE ALLOCATE(OUTARRAY(0)) OUTSIZE = 0 END IF RETURN ! * end of subroutine UNION * END SUBROUTINE UNION !*********************************************************************** ! INTEGER FUNCTION LENGTH(ARRAY,ARRAYSIZE,VAL) ! * !*********************************************************************** ! IMPLICIT NONE ! ! 0. Authors ! ! 40.81: Andre van der Westhuysen ! ! 2. Purpose ! ! Replicates the Matlab intrinsic function LENGTH. Find largest ! index in ARRAY with a value not equal to the filler value VAL. ! E.g. If VAL = 9999. and ARRAY = [X X X X 9999. 9999. 9999.], ! the function returns 4. ! ! Argument variables REAL :: ARRAY(ARRAYSIZE) REAL :: VAL INTEGER :: ARRAYSIZE ! Local variables REAL :: FIELD INTEGER :: I I = 1 FIELD = ARRAY(I) DO WHILE (FIELD.NE.VAL) I = I+1 IF (I.GT.SIZE(ARRAY)) EXIT FIELD = ARRAY(I) END DO LENGTH = I-1 RETURN ! * end of function LENGTH * END FUNCTION LENGTH !*********************************************************************** ! INTEGER FUNCTION FINDFIRST(ARRAY,ARRAYSIZE,VAL) ! * !*********************************************************************** ! IMPLICIT NONE ! ! 0. Authors ! ! 40.81: Andre van der Westhuysen ! ! 2. Purpose ! ! Based on the Matlab intrinsic function FIND. Fast ! algorithm to find the *first* index IND in ARRAY ! for which ARRAY(IND) = VAL. Use only when there are ! no duplicates in ARRAY! ! ! Argument variables REAL :: ARRAY(ARRAYSIZE) REAL :: VAL INTEGER :: ARRAYSIZE ! Local variables INTEGER :: IND IND = 1 DO WHILE (IND.LE.ARRAYSIZE) IF ( ARRAY(IND).EQ.VAL ) EXIT IND = IND + 1 END DO IF (IND.GT.ARRAYSIZE) THEN ! WRITE(103,*) '*** In FINDFIRST: IND not found!' FINDFIRST = 0 ELSE FINDFIRST = IND ENDIF RETURN ! * end of function FINDFIRST * END FUNCTION FINDFIRST !*********************************************************************** ! REAL FUNCTION STD(ARRAY,N) ! * !*********************************************************************** ! IMPLICIT NONE ! ! 0. Authors ! ! 40.81: Andre van der Westhuysen ! ! 2. Purpose ! ! Replicates the Matlab intrinsic function STD. ! Computes standard deviation. ! ! Argument variables ! ARRAY REAL Input array for which to compute the std dev. ! N INT Size of ARRAY ! REAL :: ARRAY(N) INTEGER :: N ! ! Local variables ! REAL :: MN IF (N.GT.1) THEN MN = SUM(ARRAY)/N STD = SQRT( 1/(REAL(N)-1)*SUM( (ARRAY(:)-MN)**2 ) ) ELSE STD = 0. END IF RETURN ! * end of function STD * END FUNCTION STD END MODULE MATFUNC MODULE PARTFUNC CONTAINS !*********************************************************************** ! SUBROUTINE findWay (way ,horizStepCount,vertStepCount , & vertBorder ,horizBorder ,stepCount ) ! * !*********************************************************************** ! IMPLICIT NONE ! ! 0. Authors ! ! 40.81: Andre van der Westhuysen and Jeff Hanson ! ! 1. Updates ! ! 40.81, Jul. 11: Origination, based on Matlab code by Jeff Hanson ! and Eve-Marie Devaliere ! ! 2. Purpose ! ! From the direction (way) we were going before, find which direction we ! are going now and how many 'steps' we need to take ! ! 3. Method ! ! - ! ! 4. Argument variables ! ! way Char in/out Direction of spiral search ! vertBorder Int input ! horizBorder Int input ! stepCount Int output Number of steps to go in the selected direction (way) ! CHARACTER :: way *1 40.PAR INTEGER :: horizStepCount, vertStepCount, 40.PAR & vertBorder, horizBorder, stepCount 40.PAR INTENT (IN) vertBorder, horizBorder INTENT (OUT) stepCount INTENT (IN OUT) way ! ! 5. Local variables ! ! ! 6. Subroutines used ! ! - ! ! 7. Subroutines calling ! ! spiralTrackV3 ! ! 8. Error messages ! ! - ! ! 9. Remarks ! ! 10. Structure ! ! - ! ! 11. Source text ! SELECT CASE (way) CASE ('R') way='D' vertStepCount=vertStepCount+1 IF (horizBorder.EQ.1) THEN horizStepCount=horizStepCount-1 END IF stepCount=vertStepCount CASE ('D') way='L' horizStepCount=horizStepCount+1 IF (vertBorder.EQ.1) THEN vertStepCount=vertStepCount-1 END IF stepCount=horizStepCount CASE ('L') way='U' vertStepCount=vertStepCount+1 IF (horizBorder.EQ.1) THEN horizStepCount=horizStepCount-1 END IF stepCount=vertStepCount CASE ('U') way='R' horizStepCount=horizStepCount+1 IF (vertBorder.EQ.1) THEN vertStepCount=vertStepCount-1 END IF stepCount=horizStepCount CASE DEFAULT WRITE(103,*) 'In spaTack:findWay should NOT go here!' END SELECT ! WRITE(103,*) 'In Subroutine findway, way= ', ! & way,' step count:',stepCount RETURN END SUBROUTINE findWay !*********************************************************************** ! SUBROUTINE findNext (i ,j ,maxI ,maxJ , & way ,vertBorder ,horizBorder ) ! * !*********************************************************************** ! IMPLICIT NONE ! ! 0. Authors ! ! 40.81: Andre van der Westhuysen and Jeff Hanson ! ! 1. Updates ! ! 40.81, Jul. 11: Origination, based on Matlab code by Jeff Hanson ! and Eve-Marie Devaliere ! ! 2. Purpose ! ! Find next point on spatial search spiral ! ! 3. Method ! ! - ! ! 4. Argument variables ! ! i,j Int in/out Current grid indices ! maxI, maxJ Int input Maximum indices of wave field ! way Char input Direction of spiral search ! vertBorder Int output Flag indicating that vert domain edge has been hit ! horizBorder Int output Flag indicating that hor domain edge has been hit ! CHARACTER :: way 40.PAR INTEGER :: i, j, maxI, maxJ, vertBorder, horizBorder 40.PAR INTENT (IN) maxI, maxJ, way INTENT (IN OUT) i, j INTENT (OUT) vertBorder, horizBorder ! ! 5. Local variables ! ! - ! ! 6. Subroutines used ! ! - ! ! 7. Subroutines calling ! ! spiralTrackV3 ! ! 8. Error messages ! ! - ! ! 9. Remarks ! ! 10. Structure ! ! - ! ! 11. Source text ! vertBorder=0 horizBorder=0 SELECT CASE (way) CASE ('R') IF (i.LT.maxI) THEN i=i+1 ELSE ! Need to tell findWay that if we hit the border we don't increment stepCount... horizBorder=1 END IF CASE ('D') IF (j.GT.1) THEN j=j-1 ELSE vertBorder=1 END IF CASE ('L') IF (i.GT.1) THEN i=i-1 ELSE horizBorder=1 END IF CASE ('U') IF (j.LT.maxJ) THEN j=j+1 ELSE vertBorder=1 END IF END SELECT ! WRITE(103,*) 'In find next, vertBorder:',vertBorder, ! & ' and horiz border=',horizBorder RETURN END SUBROUTINE findNext !*********************************************************************** ! SUBROUTINE findSys (i ,j ,wsdat ,maxSys , & ngbrExt ,maxI ,maxJ ,freqKnob , & dirKnob ) ! * !*********************************************************************** ! USE PARTDAT USE MATFUNC ! IMPLICIT NONE ! ! 0. Authors ! ! 40.81: Andre van der Westhuysen and Jeff Hanson ! ! 1. Updates ! ! 40.81, Jul. 11: Origination, based on Matlab code by Jeff Hanson ! and Eve-Marie Devaliere ! ! 2. Purpose ! ! Find all wave systems that neighbour the grid point (i,j), and ! match these with the systems at (i,j). ! ! 3. Method ! ! For the given point (i,j), find all wave systems at neighbouring grid ! points within the reach specified by ngbrExt. ! ! 4. Argument variables ! ! i,j Int input Current grid indices ! maxI, maxJ Int input Maximum indices of wave field ! wsdat Type(dat2d) in/out Input data structure to be spiral tracked ! maxSys Int in/out Maximum number of systems identified ! TYPE(dat2d) :: wsdat 40.PAR INTEGER :: i, j, maxI, maxJ, ngbrExt, maxSys 40.PAR REAL :: freqKnob ,dirKnob 40.PAR INTENT (IN) i, j, maxI, maxJ, ngbrExt, freqKnob ,dirKnob INTENT (IN OUT) wsdat, maxSys ! ! 5. Local variables ! ! tmpsys TYPE(system) Temporary instance of the wave system variable ! nngbr Int Number of neighbours found ! TYPE(system), ALLOCATABLE :: tmpsys(:) 40.PAR TYPE(neighbr) :: ngbr(50) TYPE(mtchsys) :: match LOGICAL :: found INTEGER :: counter, ii, jj, nngbr, startCount, endCount, l, & nout, maxS, s, p, n, countAll, ind, minInd, & npart, pp INTEGER :: allFullSys(50) REAL, ALLOCATABLE :: realarr(:) INTEGER, ALLOCATABLE :: allSys(:) REAL :: hsAll(50),tpAll(50),dirAll(50) REAL :: absDir, T, deltaPer, deltaDir, temp ! ! 6. Subroutines used ! ! UNIQUE ! combinePartitionsV2 ! ! 7. Subroutines calling ! ! spiralTrackV3 ! ! 8. Error messages ! ! - ! ! 9. Remarks ! ! 10. Structure ! ! ! 11. Source text ! ! WRITE(103,*) 'findSys: i,j,maxSys =',i,j,maxSys ! First find the checked neighbor counter=1 DO ii=(i-ngbrExt), (i+ngbrExt) DO jj=(j-ngbrExt), (j+ngbrExt) IF ( (ii.GT.0).AND.(jj.GT.0).AND. & (jj.LE.maxJ).AND.(ii.LE.maxI) ) THEN IF ( wsdat%par(ii,jj)%checked.EQ.1 ) THEN ngbr(counter)%par = wsdat%par(ii,jj) 40.81 !Added the par field to maintain the data structure ngbr(counter)%i = ii ngbr(counter)%j = jj counter=counter+1 END IF END IF END DO END DO ! New variable nngbr nngbr=counter-1 IF (nngbr.GT.0) THEN allFullSys(:) = 0 startCount=1 l=1 DO WHILE (l.LE.nngbr) endCount = startCount + LENGTH(REAL(ngbr(l)%par%sys), & SIZE(ngbr(l)%par%sys),REAL(9999))-1 allFullSys(startCount:endCount) = ngbr(l)%par%sys ! WRITE(103,*) 'In findSys:startCount,endCount,allFullSys(:)', ! & startCount,endCount,allFullSys(:) startCount=endCount+1 l=l+1 END DO ! allSys=unique(allFullSys(allFullSys>0)); ??? Test for >0 necessary? sys never gets this value! CALL UNIQUE (REAL(allFullSys),endCount,realarr,nout) 40.81 Can I do this ??? ALLOCATE(allSys(nout)) allSys = INT(realarr) 40.81 Can I do this ??? maxS = MAXVAL(allSys) IF (maxSys.LT.maxS) THEN WRITE(103,*) 'In findSys: updated maxSys from ', & maxSys,' to ',maxS maxSys=maxS END IF ! Initiate sys num ALLOCATE( tmpsys(SIZE(allSys)) ) ! Clear the wsdat%par(i,j)%sys field, new values assigned below. ! System info temporarily stored in allSys wsdat%par(i,j)%sys(1:10) = 9999 DO s=1, SIZE(allSys) hsAll(:) = 0. tpAll(:) = 0. dirAll(:) = 0. n=1 countAll=0 ! DO WHILE ~isempty(ngbr{n}) DO WHILE (n.LE.nngbr) ! Calculate mean of common neighbor wave system ! for every neigbor wave system ! ind=ngbr{n}.sys==allSys(s); found = .FALSE. DO ind = 1, SIZE(ngbr(n)%par%sys) 40.81 !Optimize this??? IF ( ngbr(n)%par%sys(ind).EQ.allSys(s) ) THEN 40.81 !Put sys under par to maintain structure found = .TRUE. EXIT END IF END DO IF (found) THEN countAll=countAll+1 hsAll(countAll)=ngbr(n)%par%hs(ind) tpAll(countAll)=ngbr(n)%par%tp(ind) dirAll(countAll)=ngbr(n)%par%dir(ind) ELSE n=n+1 CYCLE END IF n=n+1 END DO tmpsys(s)%HsMean = SUM(hsAll(1:countAll))/countAll tmpsys(s)%tpMean = SUM(tpAll(1:countAll))/countAll tmpsys(s)%dirMean = & mean_angleV2(dirAll(1:countAll),countAll) END DO ! Find the partition at current (i,j) point that matches previously ! identified wave systems if any... wsdat%par(i,j)%ngbrSys(1:SIZE(allSys)) = allSys npart = LENGTH(REAL(wsdat%par(i,j)%ipart), & SIZE(wsdat%par(i,j)%ipart),REAL(0)) DO p = 1, npart IF ( (wsdat%par(i,j)%hs(p).EQ.0.).AND. & (wsdat%par(i,j)%tp(p).EQ.0.) ) THEN wsdat%par(i,j)%sys(p)=-1 CYCLE END IF ind=0 40.81 !Replaced 'index' by 'ind' match%sysVal(:) = 9999 match%tpVal(:) = 9999. ! Cycle through the neighbouring systems identified above DO s=1,SIZE(allSys) absDir = ABS(wsdat%par(i,j)%dir(p)-tmpsys(s)%dirMean) IF (absDir.GT.180) THEN absDir = 360 - absDir END IF ! Calculate delta dir and freq as a function of the partition ! dir and freq T = wsdat%par(i,j)%tp(p) deltaPer = -0.06*T + 2 + freqKnob deltaDir = -T + (25 + 10*dirKnob) IF ( (ABS(T-tmpsys(s)%tpMean).LT.deltaPer).AND. & (absDir.LT.deltaDir) ) THEN ind=ind+1 match%sysVal(ind) = allSys(s) match%tpVal(ind) = ABS(wsdat%par(i,j)%tp(p) & - tmpsys(s)%tpMean) END IF END DO ! if isfield(match,'sysVal') IF (ind.GT.0) THEN IF (ind.EQ.1) THEN wsdat%par(i,j)%sys(p) = match%sysVal(1) ELSE ! Take the closest match ! minInd = MINLOC(match%tpVal(:)) 40.81 Optimize ??? temp = 9999. DO pp = 1,ind IF (temp.GT.match%tpVal(pp)) THEN temp = match%tpVal(pp) minInd = pp END IF END DO wsdat%par(i,j)%sys(p) = match%sysVal(minInd) 40.81 !The index of the system is swapped - the remaing info stays the same! END IF END IF ! match = 0. 40.81 !Probably not necessary END DO END IF ! Now check if 2 partitions have been associated to the same wave system, if ! so combine them npart = LENGTH(REAL(wsdat%par(i,j)%ipart), & SIZE(wsdat%par(i,j)%ipart),REAL(0)) ! WRITE(103,*) 'npart, wsdat%par(i,j)%ipart',npart,wsdat%par(i,j)%ipart ! WRITE(103,*) 'wsdat%par(i,j)%sys(:) =',wsdat%par(i,j)%sys(:) DO p = 1, (npart-1) 40.81 !Could probably be optimized! DO pp = (p+1), npart IF (wsdat%par(i,j)%sys(p).EQ.wsdat%par(i,j)%sys(pp)) THEN ! There is at least one duplicate, so combine systems WRITE(103,*) 'Call combinePartitionsV2...' CALL combinePartitionsV2(wsdat%par(i,j)) END IF END DO END DO ! Now that we have associated any possible partition to an existing ! wave system, we check if any wave system is free. If so give it a ! new wave system number npart = LENGTH(REAL(wsdat%par(i,j)%ipart), & SIZE(wsdat%par(i,j)%ipart),REAL(0)) DO p = 1, npart IF (wsdat%par(i,j)%sys(p).EQ.9999) THEN WRITE(103,*) 'In findSys: NEW wave system for ',i,j maxSys = maxSys + 1 wsdat%par(i,j)%sys(p) = maxSys END IF END DO wsdat%par(i,j)%checked=1 IF (ALLOCATED(allSys)) DEALLOCATE(allSys) IF (ALLOCATED(realarr)) DEALLOCATE(realarr) IF (ALLOCATED(tmpsys)) DEALLOCATE(tmpsys) RETURN END SUBROUTINE findSys !*********************************************************************** ! SUBROUTINE combineWaveSystems (wsdat ,maxSys ,maxPts , & maxI ,maxJ ,freqKnob , & dirKnob ,hsKnob ,combine , & sys ) ! * !*********************************************************************** ! USE PARTDAT USE MATFUNC ! IMPLICIT NONE ! ! 0. Authors ! ! 40.81: Andre van der Westhuysen and Jeff Hanson ! ! 1. Updates ! ! 40.81, Jul. 11: Origination, based on Matlab code by Jeff Hanson ! and Eve-Marie Devaliere ! ! 2. Purpose ! ! Combine wave systems. Then remove small and low-energy systems from set, ! based on the parameters maxPts and maxHgt. ! ! 3. Method ! ! ! 4. Argument variables ! ! wsdat Type(dat2d) output Combined wave system data structure ! sys Type(system) output Final set of tracked systems, for one time level ! maxI, maxJ Int input Maximum indices of wave field ! maxSys Int input Maximum number of systems identified ! maxPts Int input Number of points req for valid system ! hsKnob Real input Parameter for identifying valid system ! combine Int input Toggle: 1=combine systems; 0=do not combine TYPE(dat2d) :: wsdat 40.PAR TYPE(system), POINTER :: sys(:) 40.PAR INTEGER :: maxSys, maxPts, maxI, maxJ, combine 40.PAR REAL :: freqKnob ,dirKnob, hsKnob INTENT (IN) maxPts, maxI, maxJ, hsKnob, combine INTENT (IN OUT) wsdat, maxSys 40.81 !In the Matlab code maxSys is only input ??? INTENT (OUT) sys ! ! 5. Local variables ! ! nSys Int Number of wave systems (for checking iterative combining loop) ! LOGICAL :: found INTEGER, ALLOCATABLE :: sysOut(:) INTEGER :: iter, ok, nSys, mS, s, so, ss, ind, leng, & iw, jw REAL :: dev, hsCmp, maxHgt, temp(5) ! ! 6. Subroutines used ! ! printFinalSys ! combineSys ! ! 7. Subroutines calling ! ! spiralTrackV3 ! ! 8. Error messages ! ! - ! ! 9. Remarks ! ! 10. Structure ! ! ! 11. Source text ! !global combine hsKnob display maxI maxJ; WRITE(103,*) 'maxSys,maxPts,maxI,maxJ,hsKnob,combine =', & maxSys,maxPts,maxI,maxJ,hsKnob,combine IF (combine.EQ.1) THEN ! Combine wave systems WRITE(103,*) 'Calling printFinalSys...' ! sys=printFinalSys(wsdat,maxSys,1,sys); CALL printFinalSys (wsdat,maxSys,maxI,maxJ,1,sys) iter=0 ok=0 ! Keep on combining wave systems until all possible combining ! has been carried out (based on the combining criteria) DO WHILE (ok.EQ.0) iter = iter+1 ! No of systems before combining nSys = maxSys WRITE(103,'(A,A,I3,A,I3,A)') 'Calling combineSys for ', & 'iteration',iter,' (maxSys =',nSys,').' ! [wsdat,sys]=combineSys(wsdat,sys); CALL combineSys (wsdat,sys,maxSys,maxI,maxJ, & freqKnob,dirKnob) ! No of systems after combining IF (maxSys.EQ.nSys) ok = 1 END DO ! % sys=printFinalSys(wsdat,maxSys,1,sys); 40.81 !Why not print the system here at the end ??? ELSE ! Do not combine wave systems WRITE(103,*) 'Calling printFinalSys (NOT combining)...' ! sys=printFinalSys(wsdat,maxSys,3,sys); CALL printFinalSys (wsdat,maxSys,maxI,maxJ,3,sys) END IF ! Remove small and low-energy systems from set, based on ! the parameters maxPts and maxHgt. ALLOCATE( sysOut(maxSys) ) sysOut = sys(:)%sysInd mS = maxSys ss = 1 WRITE(6,*) '*** WARNING: ', & 'In combineWaveSystems: Redundant sys(:)%hs ', & 'fields not deallocated' DO so = 1, mS s = sysOut(so) ! ss=find([sys.sysInd]==s); DO WHILE (ss.LE.mS) IF ( sys(ss)%sysInd.EQ.s ) EXIT ss = ss + 1 END DO IF (ss.GT.mS) WRITE(103,*) & '*** In combineWaveSystems: sysInd not found!' ! WRITE(103,*) 'so, ss',so,ss leng = LENGTH(REAL(sys(ss)%hs), & SIZE(sys(ss)%hs),REAL(9999.)) dev = STD(sys(ss)%hs(1:leng),leng) hsCmp = sys(ss)%hsMean + 2.*dev maxHgt = (-0.1*sys(ss)%tpMean + hsKnob)/3.28 IF ( (hsCmp.LT.maxHgt).OR.(sys(ss)%nPoints.LT.maxPts) ) THEN IF ( hsCmp',num2str(ngbIndex(ngb))]); WRITE(103,*) 'Calling findIJV3...' ! [indSys1,indSys2]=findIJV3(sys(ss),sys(ngbIndex(ngb))); WRITE(103,*) 'ss, ngb, ngbIndex(ngb)', & ss, ngb, ngbIndex(ngb) CALL findIJV3 (sys(ss),sys(ngbIndex(ngb)), & maxI,indSys1,indSys2) ! WRITE(103,*) 'indSys1 =',indSys1 ! WRITE(103,*) 'indSys2 =',indSys2 IF ((SIZE(indSys1)>0).AND.(SIZE(indSys2)>0)) THEN ! % disp(['For wave system ',num2str(s),' with wave system:',num2str(ngbIndex(ngb))]); WRITE(103,*) 'Comparing wave system ',s, & ' with wave system: ',ngbIndex(ngb) lsys = SIZE(indSys1) lsys2 = SIZE(indSys2) Tb = SUM(sys(ss)%tp(indSys1))/lsys deltaPerB = (-0.06*Tb+2+freqKnob)*1.5 deltaDirB = (-Tb+(25+10*dirKnob))*1.5 absDir = ABS( & mean_angleV2(sys(ss)%dir(indSys1),lsys) - & mean_angleV2(sys(ngbIndex(ngb))%dir(indSys2), & lsys2) ) IF (absDir.GT.180) absDir = 360.-absDir WRITE(103,*) 'Tb_val,deltaPerB =',ABS(Tb - & SUM(sys(ngbIndex(ngb))%tp(indSys2))/lsys2), & deltaPerB WRITE(103,*) 'absDir,deltaDirB =',absDir,deltaDirB IF ( ( ABS(Tb - & SUM(sys(ngbIndex(ngb))%tp(indSys2))/lsys2) & .LT.deltaPerB ).AND.(absDir.LT.deltaDirB) ) THEN indMatch = ngbIndex(ngb) matchSys = sys(indMatch)%sysInd WRITE(103,*) 'indMatch, matchSys = ', & indMatch, matchSys ELSE ! continue; CYCLE END IF !---------------------Test out ! leng = LENGTH(REAL(sys(indMatch)%hs), ! & SIZE(sys(indMatch)%hs),REAL(9999.)) ! WRITE(103,*) 'indMatch = ',indMatch ! WRITE(103,*) 'sys(indMatch)%hs = ', ! & sys(indMatch)%hs(1:leng) ! WRITE(103,*) 'sys(indMatch)%tp = ', ! & sys(indMatch)%tp(1:leng) ! WRITE(103,*) 'sys(indMatch)%dir = ', ! & sys(indMatch)%dir(1:leng) !---------------------Test out ! indMatch=find([sys.sysInd]==matchSys); ??? Already computed above ??? WRITE(103,*) 'Combine wave system ',s, & ' mean tp=',sys(ss)%tpMean,' and dir=', & sys(ss)%dirMean,' with wave system ', & matchSys,' mean tp=',sys(indMatch)%tpMean, & ' and dir=',sys(indMatch)%dirMean ! %justOut(length(justOut)+1)=s; ! end ! %Find out if we need to update system list later, by checking if it is ! %part of the ones that has already been combined ! % indA=find(justOut==matchSys); ! % if isempty(indA)%make sure we didn't already deal with this one ! % update=0; ! % ! % if ~isempty(indMatch) ! % % disp('need to update matched system list cause can be used later'); ! % update=1; ! % else ! % disp('indMatch is empty! how come???'); ! % end ! % end ! %replace the system number is wsdat (first find 'where' the ! %system to replace is) keep = 0 keepInd(:) = 0 WRITE(103,*) 'sysic(ss)%nPoints=',sys(ss)%nPoints DO hh = 1, sys(ss)%nPoints ii = sys(ss)%i(hh) jj = sys(ss)%j(hh) ! wsdat.par(ii,jj).sys(wsdat.par(ii,jj).sys==s)=matchSys; ind = 0 ind = FINDFIRST(REAL(wsdat%par(ii,jj)%sys), & SIZE(wsdat%par(ii,jj)%sys),REAL(s)) ! WRITE(103,*) 'wsdat%par(ii,jj)%sys =', ! & wsdat%par(ii,jj)%sys IF (ind.NE.0) THEN wsdat%par(ii,jj)%sys(ind)=matchSys END IF ! WRITE(103,*) 's, ind, matchSys =',s,ind,matchSys ! Remove the "-1" system from the set ! oneLess = wsdat.par(ii,jj).sys(wsdat.par(ii,jj).sys~=-1); ind2 = 1 oneLess(:) = 9999 ??? Streamline this ??? leng = LENGTH(REAL(wsdat%par(ii,jj)%sys), & SIZE(wsdat%par(ii,jj)%sys),REAL(9999)) DO ind = 1, leng IF ( wsdat%par(ii,jj)%sys(ind).NE.-1 ) THEN oneLess(ind2) = wsdat%par(ii,jj)%sys(ind) ind2 = ind2+1 END IF END DO ind2 = ind2-1 ! WRITE(103,*) 'ind2, oneLess(:) =',ind2, oneLess ! Combine any partitions assigned to the same systems ! Check for duplicates CALL UNIQUE(REAL(oneLess(1:ind2)),ind2, & uniarr,outsize) IF (ind2.GT.outsize) THEN ! There is at least one duplicate, so combine systems WRITE(103,*) 'Call combinePartitionsV2...' CALL combinePartitionsV2(wsdat%par(ii,jj)) WRITE(103,*) 'wsdat%par(ii,jj)%sys(:) = ', & wsdat%par(ii,jj)%sys(:) ! Update the combined partitions values into the system we are keeping. ! Since partitions have been combined we don't know if the index is the same ! replacedInd=find(wsdat.par(ii,jj).sys==matchSys); replacedInd = & FINDFIRST(REAL(wsdat%par(ii,jj)%sys(:)), & SIZE(wsdat%par(ii,jj)%sys(:)), & REAL(matchSys)) ! WRITE(103,*) 'matchSys, replacedInd = ', ! & matchSys,replacedInd ! hhForIndMatch=find(sys(indMatch).i==ii & sys(indMatch).j==jj); hhForIndMatch = 1 DO WHILE (hhForIndMatch.LE. & sys(indMatch)%nPoints) IF ( (sys(indMatch)%i(hhForIndMatch) & .EQ.ii).AND. & (sys(indMatch)%j(hhForIndMatch) & .EQ.jj) ) EXIT hhForIndMatch = hhForIndMatch + 1 END DO sys(indMatch)%hs(hhForIndMatch) = & wsdat%par(ii,jj)%hs(replacedInd) sys(indMatch)%tp(hhForIndMatch) = & wsdat%par(ii,jj)%tp(replacedInd) sys(indMatch)%dir(hhForIndMatch) = & wsdat%par(ii,jj)%dir(replacedInd) ELSE keep = keep+1 keepInd(keep) = hh END IF DEALLOCATE(uniarr) END DO leng = LENGTH(REAL(sys(indMatch)%hs), & SIZE(sys(indMatch)%hs),REAL(9999.)) ! WRITE(103,*) 'indMatch = ',indMatch ! WRITE(103,*) 'sys(indMatch)%hs = ', ! & sys(indMatch)%hs(1:leng) ! WRITE(103,*) 'sys(indMatch)%tp = ', ! & sys(indMatch)%tp(1:leng) ! WRITE(103,*) 'sys(indMatch)%dir = ', ! & sys(indMatch)%dir(1:leng) ! WRITE(103,*) 'keep, keepInd(1:keep) = ', ! & keep, keepInd(1:keep) ! Update system info ! ------------------ ! First need to find which points were common to both systems => ! keepInd since that means partitions have not been combined for those ! points as a result of the combination of those 2 systems => ! distinct points ! keepInd = keepInd(1:keep) lMatch = LENGTH(REAL(sys(indMatch)%hs), & SIZE(sys(indMatch)%hs),REAL(9999.)) tot = lMatch + keep WRITE(103,*) 'tot =',tot ! sys(indMatch).ngbr = setdiff(union(sys(indMatch).ngbr,sys(ss).ngbr) , [sys(indMatch).sysInd sys(ss).sysInd]); CALL UNION (REAL(sys(indMatch)%ngbr), & SIZE(sys(indMatch)%ngbr), & REAL(sys(ss)%ngbr), & SIZE(sys(ss)%ngbr), & allngbr,outsize) CALL SETDIFF(allngbr,SIZE(allngbr), & REAL((/sys(indMatch)%sysInd, & sys(ss)%sysInd/)), & SIZE((/sys(indMatch)%sysInd, & sys(ss)%sysInd/)),difarr,outsize) sys(indMatch)%ngbr(:) = 9999 sys(indMatch)%ngbr(1:outsize) = NINT(difarr) ! WRITE(103,*) 'sys(indMatch)%ngbr', ! & sys(indMatch)%ngbr DEALLOCATE(allngbr) DEALLOCATE(difarr) ! % disp(['old ht mean 1=',num2str(sys(indMatch).hsMean),' ... ! % for npts=',num2str(sys(indMatch).nPoints),' and mean 2 =', ... ! % num2str(sys(ss).hsMean),' for npts=',num2str(sys(ss).nPoints)]); leng = LENGTH(REAL(sys(indMatch)%i), & SIZE(sys(indMatch)%i),REAL(9999)) sys(indMatch)%hsMean = SUM((/ & sys(ss)%hs(keepInd(1:keep)), & sys(indMatch)%hs(1:leng) /))/tot ! % disp(['new mean=',num2str(sys(indMatch).hsMean)]); sys(indMatch)%tpMean = SUM((/ & sys(ss)%tp(keepInd(1:keep)), & sys(indMatch)%tp(1:leng) /))/tot sys(indMatch)%dirMean = & mean_angleV2((/ sys(ss)%dir(keepInd(1:keep)), & sys(indMatch)%dir(1:leng) /),tot) ! WRITE(103,*) 'sys(indMatch)%hsMean = ', ! & sys(indMatch)%hsMean ! WRITE(103,*) 'sys(indMatch)%tpMean = ', ! & sys(indMatch)%tpMean ! WRITE(103,*) 'sys(indMatch)%dirMean = ', ! & sys(indMatch)%dirMean !>0725 leng = LENGTH(REAL(sys(indMatch)%i), !>0725 & SIZE(sys(indMatch)%i),REAL(9999)) sys(indMatch)%i(1:(keep+leng))= & (/sys(ss)%i(keepInd(1:keep)), & sys(indMatch)%i(1:leng)/) sys(indMatch)%j(1:(keep+leng))= & (/sys(ss)%j(keepInd(1:keep)), & sys(indMatch)%j(1:leng)/) sys(indMatch)%lat(1:(keep+leng)) = & (/sys(ss)%lat(keepInd(1:keep)), & sys(indMatch)%lat(1:leng)/) sys(indMatch)%lon(1:(keep+leng)) = & (/sys(ss)%lon(keepInd(1:keep)), & sys(indMatch)%lon(1:leng)/) sys(indMatch)%dir(1:(keep+leng)) = & (/sys(ss)%dir(keepInd(1:keep)), & sys(indMatch)%dir(1:leng)/) sys(indMatch)%hs(1:(keep+leng)) = & (/sys(ss)%hs(keepInd(1:keep)), & sys(indMatch)%hs(1:leng)/) sys(indMatch)%tp(1:(keep+leng)) = & (/sys(ss)%tp(keepInd(1:keep)), & sys(indMatch)%tp(1:leng)/) sys(indMatch)%nPoints = & LENGTH(REAL(sys(indMatch)%i), & SIZE(sys(indMatch)%i),REAL(9999)) sys(ss)%nPoints = 0 ! WRITE(103,*) 'sys(indMatch)%nPoints = ', ! & sys(indMatch)%nPoints ! WRITE(103,*) 'sys(ss)%nPoints = ', ! & sys(ss)%nPoints ! WRITE(103,*) 'sys(indMatch)%hs = ', ! & sys(indMatch)%hs(1:sys(indMatch)%nPoints) ! WRITE(103,*) 'sys(indMatch)%tp = ', ! & sys(indMatch)%tp(1:sys(indMatch)%nPoints) ! WRITE(103,*) 'sys(indMatch)%dir = ', ! & sys(indMatch)%dir(1:sys(indMatch)%nPoints) ! Loop through wsdat to update neighbouring system values DO i = 1, maxI DO j = 1, maxJ ind = FINDFIRST(REAL(wsdat%par(i,j)%ngbrSys), & SIZE(wsdat%par(i,j)%ngbrSys),REAL(s)) IF (ind.NE.0) THEN wsdat%par(i,j)%ngbrSys(ind)=matchSys END IF ! wsdat.par(i,j).ngbrSys=unique(wsdat.par(i,j).ngbrSys); leng = LENGTH(REAL(wsdat%par(i,j)%ngbrSys), & SIZE(wsdat%par(i,j)%ngbrSys),REAL(9999)) ! WRITE(103,*) 'leng = ',leng IF (leng.GT.0) THEN CALL UNIQUE( & REAL(wsdat%par(i,j)%ngbrSys(1:leng)), & leng,uniarr,outsize) wsdat%par(i,j)%ngbrSys(:) = 9999 wsdat%par(i,j)%ngbrSys(1:outsize) = & NINT(uniarr) DEALLOCATE(uniarr) ELSE wsdat%par(i,j)%ngbrSys(:) = 9999 END IF ! WRITE(103,*) 'i,j, wsdat%par(i,j)%ngbrSys = ', ! & i,j, wsdat%par(i,j)%ngbrSys END DO END DO ! Update neigbors in sys structure DO nn = 1, maxSys ! nbr=find(sys(nn).ngbr==s); nbr = FINDFIRST(REAL(sys(nn)%ngbr), & SIZE(sys(nn)%ngbr),REAL(s)) ! WRITE(103,*) 'nbr = ',nbr ! WRITE(103,*) 'for sys=',sys(nn)%sysInd, ! & ' had ngbr',sys(nn)%ngbr IF (nbr.NE.0) THEN ! WRITE(103,*) 'update' sys(nn)%ngbr(nbr)=matchSys END IF ! sys(nn).ngbr=unique(sys(nn).ngbr); leng2 = LENGTH(REAL(sys(nn)%ngbr), & SIZE(sys(nn)%ngbr),REAL(9999)) CALL UNIQUE(REAL(sys(nn)%ngbr(1:leng2)), & leng2,uniarr,outsize) sys(nn)%ngbr(:) = 9999 sys(nn)%ngbr(1:outsize) = uniarr DEALLOCATE(uniarr) ! WRITE(103,*) 'has now ngbr: ', ! & sys(nn)%ngbr(1:outsize) END DO ! break; ! leng = LENGTH(REAL(sys(:)%sysInd), ! & SIZE(sys(:)%sysInd),REAL(9999)) ! DO ic = 1,leng ! WRITE(103,*) ic,': For wave System number:', ! & sys(ic)%sysInd,' nGuys= ',sys(ic)%nPoints, ! & 'hsMean=',sys(ic)%hsMean,' tpMean=', ! & sys(ic)%tpMean,' dirMean=',sys(ic)%dirMean ! END DO EXIT END IF DEALLOCATE(indSys1) DEALLOCATE(indSys2) END DO END IF END DO WRITE(103,*) 'End of big loop:' DO ic = 1,maxSys WRITE(103,*) ic,': For wave System number:', & sys(ic)%sysInd,' nGuys= ',sys(ic)%nPoints, & 'hsMean=',sys(ic)%hsMean,' tpMean=', & sys(ic)%tpMean,' dirMean=',sys(ic)%dirMean END DO ! Loop one more time through sys to get rid of the empty ones. Can't do it ! in the first loop because it would change the index ! *** NOTE: maxSys is computed incorrectly in combineSys, ! True ??? ! but is not used further on, so does not affect the results. ! True ??? icEnd = maxSys WRITE(103,*) 'icEnd =',icEnd ! WRITE(6,*) 'maxSys =',maxSys ic = 1 WRITE(6,*) '*** WARNING: ', & 'In combineSys: Redundant sys(:)%hs ', & 'fields not deallocated' DO WHILE (ic.LE.icEnd) ! WRITE(6,*) 'Inside loop: icEnd =',icEnd ! WRITE(6,*) 'Inside loop: maxSys =',maxSys IF (sys(ic)%nPoints.EQ.0) THEN WRITE(103,*) 'In combineSys: Removing system ',ic,'...' ! WRITE(6,*) 'Before: maxSys =',maxSys ! WRITE(6,*) 'ASSOCIATED =',ASSOCIATED(sys(59)%hs) ! WRITE(6,*) 'ASSOCIATED =',ASSOCIATED(sys(60)%hs) ! WRITE(6,*) 'ASSOCIATED =',ASSOCIATED(sys(61)%hs) ! WRITE(6,*) 'ASSOCIATED =',ASSOCIATED(sys(62)%hs) ! WRITE(6,*) 'SIZE(sys(59)%hs)=',SIZE(sys(59)%hs) ! WRITE(6,*) 'SIZE(sys(60)%hs)=',SIZE(sys(60)%hs) ! WRITE(6,*) 'SIZE(sys(61)%hs)=',SIZE(sys(61)%hs) ! WRITE(6,*) 'SIZE(sys(62)%hs)=',SIZE(sys(62)%hs) ! sys(ic)=[]; ! Remove system, and shift up indices to fill the gap ! Shift up entries, deleting the duplicate partition sys( ic:(icEnd-1) ) = sys( (ic+1):icEnd ) ! Deallocate record at the end ! sys(icEnd)%hs(:) = 9999. ! sys(icEnd)%tp(:) = 9999. ! sys(icEnd)%dir(:) = 9999. ! sys(icEnd)%i(:) = 9999 ! sys(icEnd)%j(:) = 9999 ! sys(icEnd)%lat(:) = 9999. ! sys(icEnd)%lon(:) = 9999. !AW091311 DEALLOCATE( sys(icEnd)%hs ) !AW091311 DEALLOCATE( sys(icEnd)%tp ) !AW091311 DEALLOCATE( sys(icEnd)%dir ) !AW091311 DEALLOCATE( sys(icEnd)%i ) !AW091311 DEALLOCATE( sys(icEnd)%j ) !AW091311 DEALLOCATE( sys(icEnd)%lat ) !AW091311 DEALLOCATE( sys(icEnd)%lon ) sys(icEnd)%sysInd = 0 !Change to 9999 ???? sys(icEnd)%hsMean = 9999. sys(icEnd)%tpMean = 9999. sys(icEnd)%dirMean = 9999. sys(icEnd)%nPoints = 9999 ! sys(icEnd)%ngbr(:) = 9999 !AW091311 NULLIFY( sys(icEnd) ) icEnd = icEnd-1 maxSys = maxSys-1 ! WRITE(6,*) 'After: maxSys =',maxSys ! WRITE(6,*) 'ASSOCIATED =',ASSOCIATED(sys(59)%hs) ! WRITE(6,*) 'ASSOCIATED =',ASSOCIATED(sys(60)%hs) ! WRITE(6,*) 'ASSOCIATED =',ASSOCIATED(sys(61)%hs) ! WRITE(6,*) 'ASSOCIATED =',ASSOCIATED(sys(62)%hs) ! WRITE(6,*) 'SIZE(sys(59)%hs)=',SIZE(sys(59)%hs) ! WRITE(6,*) 'SIZE(sys(60)%hs)=',SIZE(sys(60)%hs) ! WRITE(6,*) 'SIZE(sys(61)%hs)=',SIZE(sys(61)%hs) ! WRITE(6,*) 'SIZE(sys(62)%hs)=',SIZE(sys(62)%hs) ELSE ic = ic+1 END IF END DO WRITE(103,*) 'Filtered result:' WRITE(103,*) 'maxSys =',maxSys ! WRITE(6,*) 'SIZE(sys) =',SIZE(sys) ! WRITE(6,*) 'ASSOCIATED =',ASSOCIATED(sys(1)%hs) ! WRITE(6,*) 'ASSOCIATED =',ASSOCIATED(sys(2)%hs) ! WRITE(6,*) 'ASSOCIATED =',ASSOCIATED(sys(3)%hs) ! WRITE(6,*) 'ASSOCIATED =',ASSOCIATED(sys(4)%hs) ! WRITE(6,*) SIZE(sys(1)%hs) ! WRITE(6,*) SIZE(sys(2)%hs) ! WRITE(6,*) SIZE(sys(3)%hs) ! WRITE(6,*) SIZE(sys(4)%hs) ! DO ic = 1,maxSys ! WRITE(103,*) 'sys(ic)%hs =',sys(ic)%hs ! END DO DO ic = 1,maxSys WRITE(103,*) ic,': For wave System number:', & sys(ic)%sysInd,' nGuys= ',sys(ic)%nPoints, & 'hsMean=',sys(ic)%hsMean,' tpMean=', & sys(ic)%tpMean,' dirMean=',sys(ic)%dirMean END DO DEALLOCATE(sysOrdered) DEALLOCATE(sysSortedInd) DEALLOCATE(sysOut) RETURN END SUBROUTINE combineSys !*********************************************************************** ! SUBROUTINE combinePartitionsV2 (dat) ! * !*********************************************************************** ! USE PARTDAT USE MATFUNC ! IMPLICIT NONE ! ! 0. Authors ! ! 40.81: Andre van der Westhuysen and Jeff Hanson ! ! 1. Updates ! ! 40.81, Jul. 11: Origination, based on Matlab code by Jeff Hanson ! and Eve-Marie Devaliere ! ! 2. Purpose ! ! Combine two partitions that have been assigned to the same system ! ! 3. Method ! ! Of all the partitions associated with a certain common system, ! add all the Hs values to the partition with the largest Hs, ! and delete the rest. NOTE that the tp and dir values of this ! maximum partition is not adjusted! ! ! 4. Argument variables ! ! dat TYPE(param) in/out Input data structure (partitions set) ! to combine ! TYPE(param) :: dat 40.PAR INTENT (IN OUT) dat ! ! 5. Local variables ! TYPE duplicate 40.PAR INTEGER :: val INTEGER :: ndup INTEGER :: ind(50) END TYPE duplicate TYPE(duplicate) :: dup(100) 40.PAR LOGICAL :: found INTEGER :: nsys, ndup, p, pp, maxInd, npart, s, ss, ppp REAL :: temp ! ! 6. Subroutines used ! ! - ! ! 7. Subroutines calling ! ! findSys ! ! 8. Error messages ! ! - ! ! 9. Remarks ! ! 10. Structure ! ! ! 11. Source text ! ! Find indices in dat%sys(:) of all partition associated with ! the same wave system, and store them in the data structure ! dup(1:nsys). Here nsys is the number of systems for which duplicates ! were found, and dup(s)%ndup the number of partitions assigned ! to the same system s. nsys = 0 dup(:)%ndup = 0 dup(:)%val = 9999 DO s = 1,100 dup(s)%ind(:) = 0 END DO npart = LENGTH(REAL(dat%ipart),SIZE(dat%ipart),REAL(0)) DO p = 1, npart-1 found = .FALSE. IF (ANY(dat%sys(p).EQ.dup(:)%val)) CYCLE !found = .TRUE. ! WRITE(103,*) 'In combinePartitionsV2: p,nsys,found=',p,nsys,found DO pp = (p+1), npart IF (dat%sys(p).EQ.dat%sys(pp)) THEN ! First value IF (.NOT.found) THEN nsys=nsys+1 dup(nsys)%val = dat%sys(p) dup(nsys)%ndup = 1 dup(nsys)%ind(dup(nsys)%ndup) = p found = .TRUE. END IF ! Subsequent duplicates IF (.NOT.ANY(pp.EQ.dup(nsys)%ind(:))) THEN ! WRITE(103,*) 'Found! at nsys, p,pp=',nsys, p, pp dup(nsys)%ndup = dup(nsys)%ndup+1 dup(nsys)%ind(dup(nsys)%ndup) = pp END IF END IF END DO END DO ! DO s = 1, nsys ! WRITE(103,*) 's, dup(s)%ind(:) =',s,dup(s)%ind(:) ! END DO ! Now go through array of duplicates for each of n systems ! to add all the wave energy to the most energetic of the ! duplicates, and then remove the rest. maxInd = 0 temp = -9999. DO s = 1, nsys ! WRITE(103,*) 's, dup(s)%ind(:) =',s,dup(s)%ind(:) ! Find duplicate partition with the largest Hs (most energy) DO p = 1, dup(s)%ndup IF ( temp.LT.dat%hs(dup(s)%ind(p)) ) THEN temp = dat%hs(dup(s)%ind(p)) maxInd = p END IF END DO ! WRITE(103,*) 's, dup(s)%ind(maxInd) =',s,dup(s)%ind(maxInd) ! Add all energy (Hs) to this partition dat%hs(dup(s)%ind(maxInd)) = & SQRT( SUM(dat%hs(dup(s)%ind(1:dup(s)%ndup))**2) ) ! Remove duplicate partitions which did not have the maximum Hs, ! and shift up indices to fill the gap DO p = 1, dup(s)%ndup ! Find index to remove IF (p.NE.maxInd) THEN ! Shift up entries, deleting the duplicate partition ! REPLACE WITH CSHIFT(ARRAY, SHIFT, dim) !!! ??? dat%hs( dup(s)%ind(p):(npart-1) ) = & dat%hs( (dup(s)%ind(p)+1):npart) dat%tp( dup(s)%ind(p):(npart-1) ) = & dat%tp( (dup(s)%ind(p)+1):npart) dat%dir( dup(s)%ind(p):(npart-1) ) = & dat%dir( (dup(s)%ind(p)+1):npart) dat%sys( dup(s)%ind(p):(npart-1) ) = & dat%sys( (dup(s)%ind(p)+1):npart) dat%ipart( dup(s)%ind(p):(npart-1) ) = & dat%ipart( (dup(s)%ind(p)+1):npart) ! Shift up indices DO ss = 1, nsys DO ppp = 1, dup(ss)%ndup IF (dup(ss)%ind(ppp).GT.dup(s)%ind(p)) & dup(ss)%ind(ppp) = dup(ss)%ind(ppp)-1 END DO END DO ! Add blank to end dat%hs(npart) = 9999. dat%tp(npart) = 9999. dat%dir(npart) = 9999. dat%sys(npart) = 9999 dat%ipart(npart) = 0 END IF END DO ! WRITE(103,*) 'dat%ipart(1) =',dat%ipart(1) ! WRITE(103,*) 'dat%ipart(2) =',dat%ipart(2) ! WRITE(103,*) 'dat%ipart(3) =',dat%ipart(3) ! WRITE(103,*) 'dat%ipart(4) =',dat%ipart(4) ! WRITE(103,*) 'dat%ipart(5) =',dat%ipart(5) ! WRITE(103,*) 'dat%ipart(6) =',dat%ipart(6) ! WRITE(103,*) 'dat%hs(1) =',dat%hs(1) ! WRITE(103,*) 'dat%hs(2) =',dat%hs(2) ! WRITE(103,*) 'dat%hs(3) =',dat%hs(3) ! WRITE(103,*) 'dat%hs(4) =',dat%hs(4) ! WRITE(103,*) 'dat%hs(5) =',dat%hs(5) ! WRITE(103,*) 'dat%hs(6) =',dat%hs(6) ! WRITE(103,*) 'dat%tp(1) =',dat%tp(1) ! WRITE(103,*) 'dat%tp(2) =',dat%tp(2) ! WRITE(103,*) 'dat%tp(3) =',dat%tp(3) ! WRITE(103,*) 'dat%tp(4) =',dat%tp(4) ! WRITE(103,*) 'dat%tp(5) =',dat%tp(5) ! WRITE(103,*) 'dat%tp(6) =',dat%tp(6) ! WRITE(103,*) 'dat%dir(1) =',dat%dir(1) ! WRITE(103,*) 'dat%dir(2) =',dat%dir(2) ! WRITE(103,*) 'dat%dir(3) =',dat%dir(3) ! WRITE(103,*) 'dat%dir(4) =',dat%dir(4) ! WRITE(103,*) 'dat%dir(5) =',dat%dir(5) ! WRITE(103,*) 'dat%dir(6) =',dat%dir(6) ! WRITE(103,*) 'dat%sys(1) =',dat%sys(1) ! WRITE(103,*) 'dat%sys(2) =',dat%sys(2) ! WRITE(103,*) 'dat%sys(3) =',dat%sys(3) ! WRITE(103,*) 'dat%sys(4) =',dat%sys(4) ! WRITE(103,*) 'dat%sys(5) =',dat%sys(5) ! WRITE(103,*) 'dat%sys(6) =',dat%sys(6) END DO RETURN END SUBROUTINE combinePartitionsV2 !*********************************************************************** ! * SUBROUTINE checkPoint (lat,lon,dat,alreadyIn,i) ! * !*********************************************************************** ! USE PARTDAT ! IMPLICIT NONE ! ! ! 0. Authors ! ! 40.81: Andre van der Westhuysen and Jeff Hanson ! ! 1. Updates ! ! 40.81, Jul. 11: Origination, based on Matlab code by Jeff Hanson ! and Eve-Marie Devaliere ! ! 2. Purpose ! ! Find whether a lon/lat point has already been processed. If so, ! return i = the last partition read for this point. ! ! 6. Subroutines used ! ! - ! ! 7. Subroutines calling ! ! waveTracking_NWS_V2 ! ! 8. Error messages ! ! - ! ! 9. Remarks ! ! 10. Structure ! ! 11. Source text ! REAL :: lat, lon 40.PAR INTEGER :: alreadyIn, i 40.PAR TYPE(dat1d) :: dat 40.PAR !Only a single time level here INTENT (IN) lat,lon,dat INTENT (OUT) alreadyIn,i i=0 ! WRITE(103,*) 'In checkPoint,lon/lat(1,1)=',dat%lon1d(1),dat%lat1d(1) ! WRITE(103,*) 'In checkPoint,lon/lat(2,2)=',dat%lon1d(2),dat%lat1d(2) ! WRITE(103,*) 'In checkPoint,lon/lat(3,3)=',dat%lon1d(3),dat%lat1d(3) ! WRITE(103,*) 'In checkPoint,lon/lat =',lon,lat IF ( ANY(dat%lat1d.EQ.lat).AND.ANY(dat%lon1d.EQ.lon) ) THEN DO i = 1,SIZE(dat%lat1d) IF ( (dat%lat1d(i).EQ.lat).AND.(dat%lon1d(i).EQ.lon) ) THEN ! Point found alreadyIn = MAXVAL(dat%par1d(i)%ipart(:)) GOTO 300 END IF END DO ! Point not found alreadyIn=0; i=0; ELSE ! Point not found alreadyIn=0; END IF 300 CONTINUE RETURN END SUBROUTINE checkPoint !*********************************************************************** ! REAL FUNCTION mean_angleV2(ang,ll) ! * !*********************************************************************** ! IMPLICIT NONE ! ! ! 0. Authors ! ! 40.81: Andre van der Westhuysen and Jeff Hanson ! ! 1. Updates ! ! 40.81, Jul. 11: Origination, based on Matlab code by ! J. Hanson, 24 April 2008 ! ! 2. Purpose ! ! Compute the mean direction from array of directions ! ! 3. Method ! ! ang is a column vector of angles ! m_ang is the mean from a unit-vector average of ang ! Assumes clockwise rotation from North = 0. ! ! 4. Argument variables ! ! ang Real input Array of angles to average ! ll Int input Length of ang ! REAL :: ang(ll) INTEGER :: ll ! ! 5. Local variables ! ! u,v Real Arrays of u,v dir components to average ! um,vm Real Mean u,v dir components ! theta Real Mean direction relative to North ! REAL :: PI PARAMETER (PI = 3.1416) REAL :: u(ll), v(ll), vm, um, theta ! ! 6. Subroutines used ! ! - ! ! 7. Subroutines calling ! ! findSys ! ! 8. Error messages ! ! - ! ! 9. Remarks ! ! 10. Structure ! ! ! 11. Source text ! ! North and East components v(:) = COS(ang(:)*(PI/180.)) u(:) = SIN(ang(:)*(PI/180.)) vm = SUM(v)/ll um = SUM(u)/ll ! Compute mean magnitude and direction relative to North (from Upolar.m) theta = (ATAN2(um,vm))*(180/PI) ! Convert inputs to radians, the to the -pi to pi range ! (incorporated from original function xunwrapV2.m) ! Convert to radians theta = theta*(PI/180) theta = PI*((ABS(theta)/PI) - & 2*CEILING(((ABS(theta)/PI)-1)/2))*SIGN(1.,theta) ! Shift the points in the -pi to 0 range to the pi to 2pi range IF (theta.LT.0.) theta = theta + 2*PI ! Convert back to degrees and return value mean_angleV2 = theta*(180/PI) RETURN END FUNCTION mean_angleV2 !*********************************************************************** ! SUBROUTINE findIJV3 (a ,b ,maxI ,ind_A ,ind_B ) ! * !*********************************************************************** ! USE PARTDAT USE MATFUNC ! IMPLICIT NONE ! ! 0. Authors ! ! 40.81: Andre van der Westhuysen and Jeff Hanson ! ! 1. Updates ! ! 40.81, Jul. 11: Origination, based on Matlab code by Jeff Hanson ! and Eve-Marie Devaliere ! ! 2. Purpose ! ! Find a(i,j) indices of system "a" that lie over or along the fringes ! of system "b". ! ! 3. Method ! ! Find indices of "a" such that a%i matches any b%(i-1), b%i or b%(i+1) ! AND a%j matches any b%(j-1), b%j or b%(j+1). ! ! 4. Argument variables ! ! a, b Type(system) output Final set of tracked systems, for one time level ! maxI Int input Maximum indices of wave field ! TYPE(system) :: a, b ??? INTEGER :: maxI, outsize 40.PAR REAL, ALLOCATABLE :: pos(:) 40.PAR INTEGER, ALLOCATABLE :: ind_A(:), ind_B(:) 40.PAR INTENT (IN) a, b, maxI INTENT (OUT) ind_A, ind_B ! ! 5. Local variables ! ! posA, posB Int Arr Array of neighbours ! indA*, indB* Int Arr Array of indices for combining systems ! INTEGER, ALLOCATABLE :: posA(:), posB(:) REAL, ALLOCATABLE :: posB_MM(:), posB_MP(:) REAL, ALLOCATABLE :: posB_PM(:), posB_PP(:) REAL, ALLOCATABLE :: indA_MM(:), indB_MM(:) REAL, ALLOCATABLE :: indA_MP(:), indB_MP(:) REAL, ALLOCATABLE :: indA_PM(:), indB_PM(:) REAL, ALLOCATABLE :: indA_PP(:), indB_PP(:) REAL, ALLOCATABLE :: indA_P(:), indA_M(:) REAL, ALLOCATABLE :: indB_P(:), indB_M(:) REAL, ALLOCATABLE :: TEMP_A(:), TEMP_B(:) REAL, ALLOCATABLE :: TEMP_A2(:), TEMP_B2(:) REAL, ALLOCATABLE :: TEMP_A3(:), TEMP_B3(:) INTEGER :: ind, leng_ai, leng_aj, leng_bi, leng_bj ! ! 6. Subroutines used ! ! INTERSECT ! UNION ! ! 7. Subroutines calling ! ! combineSys ! ! 8. Error messages ! ! - ! ! 9. Remarks ! ! 10. Structure ! ! ! 11. Source text ! ! Start by limiting search to possible common region. ! Find min for i and j for data A and B leng_ai = LENGTH(REAL(a%i),SIZE(a%i),REAL(9999)) leng_aj = LENGTH(REAL(a%j),SIZE(a%j),REAL(9999)) leng_bi = LENGTH(REAL(b%i),SIZE(b%i),REAL(9999)) leng_bj = LENGTH(REAL(b%j),SIZE(b%j),REAL(9999)) WRITE(103,*) 'leng_ai,leng_aj,leng_bi,leng_bj,maxI', & leng_ai,leng_aj,leng_bi,leng_bj,maxI ALLOCATE( posA(leng_ai) ) ALLOCATE( posB(leng_bi) ) DO ind = 1, leng_ai posA(ind)=(a%j(ind)-1)*maxI+a%i(ind) END DO DO ind = 1, leng_bi posB(ind)=(b%j(ind)-1)*maxI+b%i(ind) END DO WRITE(103,*) 'size(posA), size(posB) =',SIZE(posA),SIZE(posB) ! WRITE(103,*) 'posA =',posA ! WRITE(103,*) 'posB =',posB ! *** IS THIS NECESSARY??? ! WRITE(103,*) 'Calling INTERSECT...' CALL INTERSECT(REAL(posA),SIZE(posA),REAL(posB),SIZE(posB), & pos,outsize,TEMP_A,TEMP_B) ??? Right? ! WRITE(103,*) 'TEMP_A =',TEMP_A ! WRITE(103,*) 'TEMP_B =',TEMP_B ! *** IS THIS NECESSARY??? ALLOCATE( posB_MM(leng_bi) ) ALLOCATE( posB_MP(leng_bi) ) ALLOCATE( posB_PM(leng_bi) ) ALLOCATE( posB_PP(leng_bi) ) ! WRITE(103,*) 'size(posB_MM,posB_MP,posB_PM,posB_PP) =', ! & size(posB_MM),size(posB_MP),size(posB_PM),size(posB_PP) DO ind = 1, leng_bi posB_MM(ind)=(b%j(ind)-2)*maxI+b%i(ind)-1 posB_MP(ind)=(b%j(ind)-2)*maxI+b%i(ind)+1 posB_PM(ind)=b%j(ind)*maxI+b%i(ind)-1 posB_PP(ind)=b%j(ind)*maxI+b%i(ind)+1 END DO ! WRITE(103,*) 'posB_MM =',posB_MM ! WRITE(103,*) 'posB_MP =',posB_MP ! WRITE(103,*) 'posB_PM =',posB_PM ! WRITE(103,*) 'posB_PP =',posB_PP CALL INTERSECT(REAL(posA),SIZE(posA),REAL(posB_MM),SIZE(posB_MM), & pos,outsize,indA_MM,indB_MM) CALL INTERSECT(REAL(posA),SIZE(posA),REAL(posB_MP),SIZE(posB_MP), & pos,outsize,indA_MP,indB_MP) CALL INTERSECT(REAL(posA),SIZE(posA),REAL(posB_PM),SIZE(posB_PM), & pos,outsize,indA_PM,indB_PM) CALL INTERSECT(REAL(posA),SIZE(posA),REAL(posB_PP),SIZE(posB_PP), & pos,outsize,indA_PP,indB_PP) ! WRITE(103,*) 'SIZE(indA_MM, indB_MM):',SIZE(indA_MM),SIZE(indB_MM) ! WRITE(103,*) 'SIZE(indA_MP, indB_MP):',SIZE(indA_MP),SIZE(indB_MP) ! WRITE(103,*) 'SIZE(indA_PM, indB_PM):',SIZE(indA_PM),SIZE(indB_PM) ! WRITE(103,*) 'SIZE(indA_PP, indB_PP):',SIZE(indA_PP),SIZE(indB_PP) ! WRITE(103,*) 'indA_MM:',indA_MM ! WRITE(103,*) 'indB_MM:',indB_MM ! WRITE(103,*) 'indA_MP:',indA_MP ! WRITE(103,*) 'indB_MP:',indB_MP ! WRITE(103,*) 'indA_PM:',indA_PM ! WRITE(103,*) 'indB_PM:',indB_PM ! WRITE(103,*) 'indA_PP:',indA_PP ! WRITE(103,*) 'indB_PP:',indB_PP ! indA_P=union(indA_PM,indA_PP); ! indA_M=union(indA_MM,indA_MP); ! ind_A=union(ind_A,union(indA_P,indA_M)); CALL UNION (REAL(indA_PM),SIZE(indA_PM),REAL(indA_PP), & SIZE(indA_PP),indA_P,outsize) CALL UNION (REAL(indA_MM),SIZE(indA_MM),REAL(indA_MP), & SIZE(indA_MP),indA_M,outsize) ! WRITE(103,*) 'indA_P:',indA_P ! WRITE(103,*) 'indA_M:',indA_M CALL UNION (REAL(indA_P),SIZE(indA_P),REAL(indA_M), ??? Removed ind_A here !!! & SIZE(indA_M),TEMP_A2,outsize) CALL UNION (REAL(TEMP_A),SIZE(TEMP_A),REAL(TEMP_A2), ??? Removed ind_A here !!! & SIZE(TEMP_A2),TEMP_A3,outsize) ALLOCATE( ind_A(SIZE(TEMP_A3)) ) ind_A = NINT(TEMP_A3) ! indB_P=union(indB_PM,indB_PP); ! indB_M=union(indB_MM,indB_MP); ! ind_B=union(ind_B,union(indB_P,indB_M)); CALL UNION (REAL(indB_PM),SIZE(indB_PM),REAL(indB_PP), & SIZE(indB_PP),indB_P,outsize) CALL UNION (REAL(indB_MM),SIZE(indB_MM),REAL(indB_MP), & SIZE(indB_MP),indB_M,outsize) ! WRITE(103,*) 'indB_P:',indB_P ! WRITE(103,*) 'indB_M:',indB_M CALL UNION (REAL(indB_P),SIZE(indB_P),REAL(indB_M), ??? Removed ind_B here !!! & SIZE(indB_M),TEMP_B2,outsize) CALL UNION (REAL(TEMP_B),SIZE(TEMP_B),REAL(TEMP_B2), ??? Removed ind_B here !!! & SIZE(TEMP_B2),TEMP_B3,outsize) ALLOCATE( ind_B(SIZE(TEMP_B3)) ) ind_B = NINT(TEMP_B3) DEALLOCATE(posA,posB) DEALLOCATE(posB_MM,posB_MP) DEALLOCATE(posB_PM,posB_PP) DEALLOCATE(indA_MM,indB_MM) DEALLOCATE(indA_MP,indB_MP) DEALLOCATE(indA_PM,indB_PM) DEALLOCATE(indA_PP,indB_PP) DEALLOCATE(indA_P,indA_M) DEALLOCATE(indB_P,indB_M) DEALLOCATE(TEMP_A,TEMP_B) DEALLOCATE(TEMP_A2,TEMP_B2) DEALLOCATE(TEMP_A3,TEMP_B3) RETURN END SUBROUTINE findIJV3 END MODULE PARTFUNC MODULE PARTMAIN CONTAINS !*********************************************************************** ! * SUBROUTINE waveTracking_NWS_V2 (intype ,tmax , & tcur , & filename ,dirKnob , & perKnob ,hsKnob , & wetPts ,seedLat , & seedLon ,dirTimeKnob, & tpTimeKnob ,paramFile , & sysA ,wsdat , & maxSys ,maxGroup ) ! * !*********************************************************************** ! USE PARTDAT USE MATFUNC USE PARTFUNC ! IMPLICIT NONE ! ! 0. Authors ! ! 40.81: Andre van der Westhuysen and Jeff Hanson ! ! 1. Updates ! ! 40.81, Jul. 11: Origination, based on Matlab code by Jeff Hanson ! and Eve-Marie Devaliere ! ! 2. Purpose ! ! Main subroutine of spatial and temporal tracking algorithm ! ! 3. Method ! ! (1) Read the partitionning output from the NWS, e.g.: ! filename = 'partRes.txt' ! dirKnob = 2. ! perKnob = 2. ! hsKnob = 2. ! wetPts = 0.10 ! seedLat = 0. ! seedLon = 0. ! dirTimeKnob = 1. ! tpTimeKnob = 2. ! paramFile = 'params' ! ! The file looks like: ! code llat llon ts hs tp dir wndSpd wndDir wndFce ! ! (2) Perform tracking in space by calling subroutine spiralTrackV3 ! (3) Perform tracking in time by caling subroutine timeTrackingV2 ! ! 4. Argument variables ! ! intype Int input For coupling: Type of input (0 = from memory; 1 = from file) ! tmax Int input For coupling: Value of maxTs to apply (1 or 2) ! tcur Int input For coupling: Index of current time step (1 or 2) ! filename Char input File name of locally partitioned data output ! dirKnob Real input Parameter in direction for combining fields in space ! perKnob Real input Parameter in period for combining fields in space ! hsKnob Real input Parameter in wave height for combining fields in space ! wetPts Real input Percentage of wet points (fraction) ! seedLat Real input Start Lat for tracking spiral (if =0 centre of field is used) ! seedLon Real input Start Lon for tracking spiral (if =0 centre of field is used) ! dirTimeKnob Real input Parameter in direction for combining fields in time ! tpTimeKnob Real input Parameter in period for combining fields in time ! paramFile Char input File name of partitioning parameters Is this used??? ! sys TYPE(timsys) output Final set of spatially and temporally tracked systems ! wsdat TYPE(dat2d) output Final version of 2D (gridded) partition data ! maxGroup Int output Maximum number of wave systems ("groups") tracked in time ! CHARACTER :: filename *32, paramFile *32 40.PAR REAL :: dirKnob, perKnob, hsKnob, wetPts, seedLat, 40.PAR & seedLon, dirTimeKnob, tpTimeKnob 40.PAR INTEGER :: maxGroup, intype, tmax, tcur 40.PAR INTEGER, ALLOCATABLE :: maxSys(:) 40.PAR TYPE(dat2d), POINTER :: wsdat(:) 40.PAR TYPE(timsys), POINTER :: sysA(:) 40.PAR ! Note: Variables wsdat, sysA and maxSys have IN/OUT intent so that they ! can be manipulated outside of this subroutine, e.g. re-indexing of ! systems and groups during the simulation. INTENT (IN) intype, tmax, tcur, filename, paramFile, dirKnob, & perKnob, hsKnob, wetPts, seedLat, seedLon, & dirTimeKnob, tpTimeKnob INTENT (OUT) maxGroup INTENT (IN OUT) wsdat, sysA, maxSys ! ! 5. Local variables ! ! dat TYPE(dat1d) Data structure for 1D input raw partition data ! TYPE(dat1d), POINTER :: dat(:) 40.PAR ! ! code Char Partitioning field code, from input file (not used) ! llat Real Latitude of partition point, from input file ! llon Real Longitude of partition point, from input file ! ts Real Time step of partition, from input file ! hs0 Real Wave height of partition, from input file ! tp0 Real Peak period of partition, from input file ! dir0 Real Mean direction of partition, from input file ! wndSpd0 Real Wind speed of partition, from input file ! wndDir0 Real Wind direction of partition, from input file ! wndFce0 Real Wind force of partition, from input file (not, used) ! tss Int. Time step counter ! t0 Int Index of first time step to compute ! CHARACTER :: code(100000) *12 INTEGER :: ts(100000), line REAL :: llat(100000),llon(100000),hs0(100000),tp0(100000), & dir0(100000),wndSpd0(100000),wndDir0(100000), & wndFce0(100000) INTEGER :: maxTS, t0, nout, maxI, maxJ REAL, ALLOCATABLE :: uniqueTim(:), uniqueLat(:), uniqueLon(:), & mlon(:,:), mlat(:,:) INTEGER :: ierr, i, j, k, alreadyIn, ok, tss, tsA, counter INTEGER :: ic, leng ! REMOVE !!! ! ! 6. Subroutines used ! ! UNIQUE ! checkPoint ! spiralTrackV3 ! timeTrackingV2 ! ! 7. Subroutines calling ! ! Run in stand-alone mode, or called by host numerical model: ! SWAN: subroutine SWANOUT1 40.PAR CHECK ??? ! ! 8. Error messages ! ! - ! ! 9. Remarks ! ! 10. Structure ! ! The structure of the tracking algorithm is the following ! (parentheses indicate minor subroutines and functions): ! ! +---- SUBROUTINE waveTracking_NWS_V2 main subroutine of spatial and temporal tracking algorithm ! | (CALL UNIQUE) removes duplicate reals from an vector ! | (CALL checkPoint) find whether a lon/lat point has already been processed ! | CALL spiralTrackV3 See below ! | CALL timeTrackingV2 See below ! | ! +---+--- SUBROUTINE spiralTrackV3 performs the spatial spiral tracking for a given time step ! | | CALL findWay find direction and no. steps in spatial search spiral ! | | CALL findNext find next point on spatial search spiral ! | | CALL findSys See below ! | | CALL combineWaveSystems See below ! | | ! | +------ SUBROUTINE findSys find all neighboring wave systems for given grid point ! | | (CALL UNIQUE) ! | | CALL combinePartitionsV2 combine two partitions that have been assigned to the same system ! | | ! | +---+-- SUBROUTINE combineWaveSystems combine wave systems, then remove small and low-energy systems ! | | CALL printFinalSys See below ! | | CALL combineSys See below ! | | ! | +----- SUBROUTINE printFinalSys output the final output systems for this time step ! | | (CALL UNIQUE) ! | | (CALL SETDIFF) returns elements in vector1 that are not in vector2 ! | | (CALL SORT) sorts the vector in ascending or descending order ! | | ! | +----- SUBROUTINE combineSys combine wave systems ! | (CALL SORT) ! | (CALL UNIQUE) ! | (CALL UNION) returns the union of vector1 and vector2 ! | (CALL SETDIFF) ! | CALL findIJV3 Find indices of system "a" that lie over or next to system "b" ! | CALL combinePartitionsV2 ! | ! +------- SUBROUTINE timeTrackingV2 performs the time tracking of all wave systems ! (CALL SORT) ! (CALL SETDIFF) ! ! 11. Source text ! ! Select input typr for raw partitioning data IF (intype.EQ.1) THEN ! Raw partitioning data is coming from an input file. ! Read file here, and set up 2d array wsdat with the data. t0 = 1 WRITE(103,*) 'Reading partitioning file...' WRITE(6,*) 'Max number of lines in "partRes" input file is', & ' = 100000!' OPEN(unit=101,file=filename,status='old') line = 1 DO WHILE (line.LT.100000) READ (101, *, END=10) code(line),llat(line),llon(line), & ts(line),hs0(line),tp0(line),dir0(line), & wndSpd0(line),wndDir0(line),wndFce0(line) ! WRITE(103,*) 'Line ',line,': ',code(line),llat(line), ! & llon(line),ts(line) line = line+1 ENDDO ! end of file 10 IERR = -1 line = line-1 CLOSE(101) ts=ts+1; ! Find unique time steps (and sort in ascending order) CALL UNIQUE (REAL(ts),SIZE(ts),uniqueTim,maxTs) ! Find unique lat and lon values (and sort in ascending order) CALL UNIQUE (llat(1:line),SIZE(llat(1:line)),uniqueLat,nout) CALL UNIQUE (llon(1:line),SIZE(llon(1:line)),uniqueLon,nout) ! Replicate the Matlab function: [mlat,mlon]=meshgrid(uniqueLat,uniqueLon) ! - Map is transposed (rotated by 90 deg), so that in Matlab: ! I (matrix row) represents Longitute ! J (matrix column) represents Latitude ! i.e. from this point onwards the indices (i,j) represents Cart. coordinates WRITE(103,*) 'Allocating arrays...' ALLOCATE( mlon(SIZE(uniqueLon),SIZE(uniqueLat)) ) ALLOCATE( mlat(SIZE(uniqueLon),SIZE(uniqueLat)) ) maxI = SIZE(uniqueLon) maxJ = SIZE(uniqueLat) DO I = 1,maxI DO J = 1,maxJ mlon(I,J) = uniqueLon(I) mlat(I,J) = uniqueLat(J) END Do END DO ! Allocate an dinitialize arrays WRITE(103,*) 'Setting up dat structure...' NULLIFY(dat) ALLOCATE(dat(maxTs)) ! WRITE(103,*) 'line = ',line DO tss = uniqueTim(1),uniqueTim(maxTs) ALLOCATE(dat(tss)%lat1d(line)) ALLOCATE(dat(tss)%lon1d(line)) ALLOCATE(dat(tss)%par1d(line)) ALLOCATE(dat(tss)%wnd1d(line)) DO counter = 1,line 40.PAR dat(tss)%par1d(counter)%ipart(:) = 0. dat(tss)%par1d(counter)%hs(:) = 9999. dat(tss)%par1d(counter)%tp(:) = 9999. dat(tss)%par1d(counter)%dir(:) = 9999. dat(tss)%wnd1d(counter)%wdir = 9999. dat(tss)%wnd1d(counter)%wspd = 9999. dat(tss)%lat1d(counter) = 9999. END DO END DO ! WRITE(103,*) 'SIZE(dat) = ',SIZE(dat) ! WRITE(103,*) 'SIZE(dat(1)%lat1d) = ',SIZE(dat(1)%lat1d) i=0 alreadyIn=0 ok=0 tss=ts(1) ! DO counter=1,SIZE(llat(1:line)) DO counter=1,line ! Discard negated partitions or no data partitions IF ( hs0(counter).LE.0. ) THEN CYCLE END IF IF ( (counter.EQ.1).OR.( tss.LT.ts(counter) ) ) THEN i=1 j=1 tss=ts(counter) ELSE tss=ts(counter) IF ( (ok.EQ.1).AND.(SIZE(dat).GE.tss) ) THEN ! WRITE(103,*) 'Call checkPoint...', ! & llon(counter),llat(counter) CALL checkPoint (llat(counter),llon(counter), & dat(tss),alreadyIn,i) ! WRITE(103,*) 'alreadyIn,i =',alreadyIn,i END IF IF ( alreadyIn.EQ.0 ) THEN j=1 IF ( (SIZE(dat).GE.tss) ) THEN 40.PAR TO ADD: .AND.(~isempty(dat(tss))) ) THEN i = LENGTH(dat(tss)%lat1d(:), & SIZE(dat(tss)%lat1d(:)),9999.) ! WRITE(103,*) 'Maxval i =',i ELSE i=0 END IF i=i+1 ELSE j=alreadyIn+1 END IF END IF dat(tss)%lat1d(i) = llat(counter) dat(tss)%lon1d(i) = llon(counter) dat(tss)%par1d(i)%hs(j) = hs0(counter) dat(tss)%par1d(i)%ipart(j) = j dat(tss)%par1d(i)%tp(j) = tp0(counter) dat(tss)%par1d(i)%dir(j) = dir0(counter) dat(tss)%wnd1d(i)%wdir = wndDir0(counter) dat(tss)%wnd1d(i)%wspd = wndSpd0(counter) ok=1 ! WRITE(103,*) 'i, lon/lat =',i,dat(tss)%lon1d(i), ! & dat(tss)%lat1d(i) END DO ! Allocate the wsdat structure WRITE(103,*) 'Setting up wsdat structure...' NULLIFY(wsdat) WRITE(103,*) 'Allocating wsdat...' ALLOCATE(wsdat(maxTS)) WRITE(103,*) 'SIZE(wsdat) = ',SIZE(wsdat) DO tsA = 1,maxTs ALLOCATE(wsdat(tsA)%lat(maxI,maxJ)) ALLOCATE(wsdat(tsA)%lon(maxI,maxJ)) ALLOCATE(wsdat(tsA)%par(maxI,maxJ)) ALLOCATE(wsdat(tsA)%wnd(maxI,maxJ)) ! WRITE(103,*) 'SIZE(wsdat(tsA)%lat,1) = ',SIZE(wsdat(tsA)%lat,1) ! WRITE(103,*) 'SIZE(wsdat(tsA)%lat,2) = ',SIZE(wsdat(tsA)%lat,2) j=1 i=1 DO WHILE ( (i.LE.maxI).AND.(j.LE.maxJ) ) wsdat(tsA)%lat(i,j)=mlat(i,j) wsdat(tsA)%lon(i,j)=mlon(i,j) wsdat(tsA)%maxi=maxI wsdat(tsA)%maxj=maxJ wsdat(tsA)%par(i,j)%hs(1:10)=9999. wsdat(tsA)%par(i,j)%tp(1:10)=9999. wsdat(tsA)%par(i,j)%dir(1:10)=9999. wsdat(tsA)%par(i,j)%ipart(1:10)=9999 wsdat(tsA)%par(i,j)%sys(1:10)=9999 40.PAR !Increase this array, or make allocatable wsdat(tsA)%par(i,j)%ngbrSys(1:50)=9999 wsdat(tsA)%wnd(i,j)%wdir=9999. wsdat(tsA)%wnd(i,j)%wspd=9999. wsdat(tsA)%par(i,j)%checked=-1 DO k = 1,SIZE(dat(tsA)%lat1d) IF ( (dat(tsA)%lat1d(k).EQ.mlat(i,j)).AND. & (dat(tsA)%lon1d(k).EQ.mlon(i,j)) ) THEN ! WRITE(103,*) 'Match found... k=',k wsdat(tsA)%lat(i,j) = dat(tsA)%lat1d(k) wsdat(tsA)%lon(i,j) = dat(tsA)%lon1d(k) wsdat(tsA)%par(i,j)%hs = dat(tsA)%par1d(k)%hs(:) wsdat(tsA)%par(i,j)%ipart=dat(tsA)%par1d(k)%ipart(:) wsdat(tsA)%par(i,j)%tp = dat(tsA)%par1d(k)%tp(:) wsdat(tsA)%par(i,j)%dir = dat(tsA)%par1d(k)%dir(:) wsdat(tsA)%wnd(i,j)%wdir = dat(tsA)%wnd1d(i)%wdir wsdat(tsA)%wnd(i,j)%wspd = dat(tsA)%wnd1d(i)%wspd wsdat(tsA)%par(i,j)%checked = 0 END IF END DO IF (i.EQ.maxI) THEN IF (j.EQ.maxJ) THEN GOTO 400 ELSE j=j+1 i=0 END IF END IF i=i+1 END DO 400 CONTINUE END DO ! Allocate the sysA structure WRITE(103,*) 'Allocating sysA...' NULLIFY( sysA ) ALLOCATE( sysA(maxTs) ) WRITE(103,*) 'SIZE(sysA) = ',SIZE(sysA) ! Allocate maxSys WRITE(103,*) 'Allocating maxSys...' ALLOCATE( maxSys(maxTs) ) WRITE(103,*) 'SIZE(maxSys) = ',SIZE(maxSys) ELSE ! Raw partitioning data from wave model memory, via the array wsdat. ! Set maxTs to the time step to compute: 1=first time step, 2=otherwise maxTs = tmax t0 = tcur ! Allocate the sysA structure WRITE(103,*) 'Allocating sysA...' NULLIFY( sysA ) ALLOCATE( sysA(1) ) ??? Change to sysA(2) !!! WRITE(103,*) 'SIZE(sysA) = ',SIZE(sysA) ! Allocate maxSys WRITE(103,*) 'Allocating maxSys...' ALLOCATE( maxSys(1) ) ??? Change to maxSys(2) !!! WRITE(103,*) 'SIZE(maxSys) = ',SIZE(maxSys) END IF DO tsA = t0,maxTs WRITE(103,*) 'Call spiralTrackV3, tsA=',tsA,'...' WRITE(103,*) 'wsdat(tsA)%lat(:,:) =' WRITE(103,*) wsdat(tsA)%lat(:,:) WRITE(103,*) 'wsdat(tsA)%lon(:,:) =' WRITE(103,*) wsdat(tsA)%lon(:,:) WRITE(103,*) 'wsdat(tsA)%par(:,:)%hs(1) =' WRITE(103,*) wsdat(tsA)%par(:,:)%hs(1) WRITE(103,*) 'wsdat(tsA)%par(:,:)%hs(2) =' WRITE(103,*) wsdat(tsA)%par(:,:)%hs(2) WRITE(103,*) 'wsdat(tsA)%par(:,:)%hs(3) =' WRITE(103,*) wsdat(tsA)%par(:,:)%hs(3) WRITE(103,*) 'wsdat(tsA)%par(:,:)%tp(1) =' WRITE(103,*) wsdat(tsA)%par(:,:)%tp(1) WRITE(103,*) 'wsdat(tsA)%par(:,:)%tp(2) =' WRITE(103,*) wsdat(tsA)%par(:,:)%tp(2) WRITE(103,*) 'wsdat(tsA)%par(:,:)%tp(3) =' WRITE(103,*) wsdat(tsA)%par(:,:)%tp(3) WRITE(103,*) 'wsdat(tsA)%par(:,:)%dir(1) =' WRITE(103,*) wsdat(tsA)%par(:,:)%dir(1) WRITE(103,*) 'wsdat(tsA)%par(:,:)%dir(2) =' WRITE(103,*) wsdat(tsA)%par(:,:)%dir(2) WRITE(103,*) 'wsdat(tsA)%par(:,:)%dir(3) =' WRITE(103,*) wsdat(tsA)%par(:,:)%dir(3) WRITE(103,*) 'wsdat(tsA)%par(:,:)%ipart(1) =' WRITE(103,*) wsdat(tsA)%par(:,:)%ipart(1) WRITE(103,*) 'wsdat(tsA)%par(:,:)%ipart(2) =' WRITE(103,*) wsdat(tsA)%par(:,:)%ipart(2) WRITE(103,*) 'wsdat(tsA)%par(:,:)%ipart(3) =' WRITE(103,*) wsdat(tsA)%par(:,:)%ipart(3) WRITE(103,*) 'wsdat(tsA)%wnd(:,:)%wspd =' WRITE(103,*) wsdat(tsA)%wnd(:,:)%wspd WRITE(103,*) 'wsdat(tsA)%wnd(:,:)%wdir =' WRITE(103,*) wsdat(tsA)%wnd(:,:)%wdir ! [wsdat{tsA},maxSys{tsA},sysA{tsA}]=spiralTrackV3(mdat(tsA),dirKnob,perKnob,wetPts,hsKnob,seedLat,seedLon); CALL spiralTrackV3 ( wsdat(tsA), dirKnob, perKnob, wetPts, & hsKnob, seedLat, seedLon, & maxSys(tsA), sysA(tsA)%sys ) WRITE(103,*) 'After spiralTrackV3: maxSys(tsA) =',maxSys(tsA) ! WRITE(103,*) 'maxSys(tsA) =',maxSys(tsA) ! WRITE(103,*) 'Just after spiralTrackV3... during tsA =',tsA ! WRITE(103,*) '1: ASSOCIATED =',ASSOCIATED(sysA(1)%sys(1)%hs) ! WRITE(103,*) '1: ASSOCIATED =',ASSOCIATED(sysA(1)%sys(2)%hs) ! WRITE(103,*) '1: ASSOCIATED =',ASSOCIATED(sysA(1)%sys(3)%hs) ! WRITE(103,*) '1: ASSOCIATED =',ASSOCIATED(sysA(1)%sys(4)%hs) ! WRITE(103,*) SIZE(sysA(1)%sys(1)%hs) ! WRITE(103,*) SIZE(sysA(1)%sys(2)%hs) ! WRITE(103,*) SIZE(sysA(1)%sys(3)%hs) ! WRITE(103,*) SIZE(sysA(1)%sys(4)%hs) ! DO ic = 1,maxSys(1) ! WRITE(103,*) 'sysA(1)%sys(ic)%hs =', ! & sysA(1)%sys(ic)%hs ! END DO ! WRITE(103,*) 'ASSOCIATED =',ASSOCIATED(sysA(tsA)%sys(1)%hs) ! WRITE(103,*) 'ASSOCIATED =',ASSOCIATED(sysA(tsA)%sys(2)%hs) ! WRITE(103,*) 'ASSOCIATED =',ASSOCIATED(sysA(tsA)%sys(3)%hs) ! WRITE(103,*) 'ASSOCIATED =',ASSOCIATED(sysA(tsA)%sys(4)%hs) ! WRITE(103,*) SIZE(sysA(tsA)%sys(1)%hs) ! WRITE(103,*) SIZE(sysA(tsA)%sys(2)%hs) ! WRITE(103,*) SIZE(sysA(tsA)%sys(3)%hs) ! WRITE(103,*) SIZE(sysA(tsA)%sys(4)%hs) ! DO ic = 1,maxSys(tsA) ! WRITE(103,*) 'sysA(tsA)%sys(ic)%hs =', ! & sysA(tsA)%sys(ic)%hs ! END DO ! DO ic = 1, maxSys(tsA) ! leng = LENGTH(REAL(sysA(tsA)%sys(ic)%ngbr), ! & SIZE(sysA(tsA)%sys(ic)%ngbr),REAL(9999)) ! WRITE(103,*) ic,': For wave System number:', ! & sysA(tsA)%sys(ic)%sysInd, ! & ' nGuys= ',sysA(tsA)%sys(ic)%nPoints,' hsMean=', ! & sysA(tsA)%sys(ic)%hsMean,' hsMax=', ! & MAXVAL( ! & sysA(tsA)%sys(ic)%hs(1:sysA(tsA)%sys(ic)%nPoints)), ! & ' tpMean=',sysA(tsA)%sys(ic)%tpMean,' dirMean=', ! & sysA(tsA)%sys(ic)%dirMean ! & ,' ngbr=',sysA(tsA)%sys(ic)%ngbr(1:leng) ! END DO END DO NULLIFY(dat) !AW rawWsdat=mdat; !AW rawFlag=3; !------ Temp Output --------------------------------------------------------- ! WRITE(103,*) 'In waveTracking_NWS_V2...' ! WRITE(103,*) 'ASSOCIATED =',ASSOCIATED(sysA(1)%sys(1)%hs) ! WRITE(103,*) 'ASSOCIATED =',ASSOCIATED(sysA(1)%sys(2)%hs) ! WRITE(103,*) 'ASSOCIATED =',ASSOCIATED(sysA(1)%sys(3)%hs) ! WRITE(103,*) 'ASSOCIATED =',ASSOCIATED(sysA(1)%sys(4)%hs) ! WRITE(103,*) SIZE(sysA(1)%sys(1)%hs) ! WRITE(103,*) SIZE(sysA(1)%sys(2)%hs) ! WRITE(103,*) SIZE(sysA(1)%sys(3)%hs) ! WRITE(103,*) SIZE(sysA(1)%sys(4)%hs) ! DO ic = 1,maxSys(1) ! WRITE(103,*) 'sysA(1)%sys(ic)%hs =', ! & sysA(1)%sys(ic)%hs ! END DO !-------------------------------------------------------------------------- WRITE(103,*) 'Call timeTrackingV2...' ! [maxGroup,sys]=timeTrackingV2(sysA,tpTimeKnob,dirTimeKnob,1); CALL timeTrackingV2 (sysA, maxSys, tpTimeKnob, dirTimeKnob, 1, & maxGroup) WRITE(103,*) 'After timeTrackingV2: SIZE(sysA(tsA)%sys) =', & SIZE(sysA(tsA)%sys) WRITE(103,*) 'After timeTrackingV2: maxSys(tsA) =',maxSys(tsA) !AW ts1=1 !AW ts2=maxTs !AW [outFolder,date0,nhrs,dateSuffix]=printNetcdfNPlots(sys,wsdat,maxGroup,netCDFfilename,paramFile); !AW plotTrackedSystem_NWS(netCDFfilename,sys,maxGroup,dirKnob,perKnob,wetPts,ts1,ts2,dirTimeKnob,tpTimeKnob,outFolder,date0,nhrs,dateSuffix); !! Intermediate output ! OPEN(unit=2,file='mdat_out.txt',status='new') !! Write lon/lat ! DO counter = 1,1 ! WRITE(2,*) 'ipart = ',counter ! DO J = 100,1,-1 ! WRITE(2,'(100(F7.2))') mdat(1)%lon(:,j) ! END DO ! END DO ! DO counter = 1,1 ! WRITE(2,*) 'ipart = ',counter ! DO J = 100,1,-1 ! WRITE(2,'(100(F7.2))') mdat(1)%lat(:,j) ! END DO ! END DO ! DO counter = 1,1 ! WRITE(2,*) 'mdat%par(:,:)%checked...' ! DO J = 100,1,-1 ! WRITE(2,'(100(I7))') mdat(1)%par(:,j)%checked ! END DO ! END DO ! !! Write Hs ! DO counter = 1,10 ! WRITE(2,*) 'ts, ipart, maxi, maxj = ', ! & 1,counter,mdat(1)%maxi,mdat(1)%maxj ! DO J = 100,1,-1 ! WRITE(2,'(100(F8.2))') mdat(1)%par(:,j)%hs(counter) ! END DO ! END DO ! DO counter = 1,10 ! WRITE(2,*) 'ts, ipart, maxi, maxj = ', ! & 2,counter,mdat(2)%maxi,mdat(2)%maxj ! DO J = 100,1,-1 ! WRITE(2,'(100(F8.2))') mdat(2)%par(:,j)%hs(counter) ! END DO ! END DO ! ! CLOSE(2) RETURN END SUBROUTINE waveTracking_NWS_V2 !*********************************************************************** ! SUBROUTINE spiralTrackV3 (wsdat ,dirKnob ,freqKnob ,wetPts , & hsKnob ,seedLat ,seedLon , & maxSys ,sys ) ! * !*********************************************************************** ! USE PARTDAT USE MATFUNC USE PARTFUNC ! IMPLICIT NONE ! ! 0. Authors ! ! 40.81: Andre van der Westhuysen and Jeff Hanson ! ! 1. Updates ! ! 40.81, Jul. 11: Origination, based on Matlab code by Jeff Hanson ! and Eve-Marie Devaliere ! ! 2. Purpose ! ! Performs the spatial spiral tracking for a given time step ! ! 3. Method ! ! - ! ! 4. Argument variables ! ! dirKnob Real input Parameter in direction for combining fields in space ! freqKnob Real input Parameter in period for combining fields in space !Why the shift to frequency here ??? ! wetPts Real input Percentage of wet points (fraction) ! hsKnob Real input Parameter in wave height for combining fields in space ! seedLat Real input Start Lat for tracking spiral (if =0 centre of field is used) ! seedLon Real input Start Lon for tracking spiral (if =0 centre of field is used) ! wsdat Real arr output Input 2d (gridded) data structure to be spiral tracked ! maxSys Int output Maximum number of partition systems ! sys Type(system) output Final set of tracked systems, for one time level ! TYPE(dat2d) :: wsdat 40.PAR REAL :: dirKnob,freqKnob,wetPts,hsKnob,seedLat,seedLon 40.PAR INTEGER :: maxSys 40.PAR TYPE(system), POINTER :: sys(:) 40.PAR INTENT (IN) dirKnob,freqKnob,wetPts,hsKnob,seedLat,seedLon INTENT (IN OUT) wsdat INTENT (OUT) maxSys,sys ! ! 5. Local variables ! ! ngbrExt Int How far do we want the neigbor to be condsidered ! combine Int Toggle (1=combine wave systems; 0=do not combine) ! maxI,MaxJ Int Dimensions of the 2d (gridded) data wsdat ! deltaLat Real Delta in kilometers between 2 pts (in latitude) ! LOGICAL :: first 40.PAR CHARACTER :: way *1 40.PAR INTEGER :: ngbrExt, combine, maxI, maxJ, i, j, oldJ 40.PAR INTEGER :: horizStepCount, vertStepCount, checkCount, sc, 40.PAR & maxPts, landPts, horizBorder, vertBorder, m, k, 40.PAR & stepCount 40.PAR REAL :: deltaLat, minLat, maxLat, minLon, maxLon 40.PAR ! ! 6. Subroutines used ! ! findWay ! findNext ! findSys ! combineWaveSystems ! ! 7. Subroutines calling ! ! waveTracking_NWS_V2 ! ! 8. Error messages ! ! - ! ! 9. Remarks ! ! 10. Structure ! ! - ! ! 11. Source text ! ! Routine starts by identifying the starting point. Choose the 'center' of the domain ! Set the search distance for neighbors: ! 1: 1 row and column out, i.e. the 8 neighbors around the current point ! 2: 2 rows and columns out... etc. ngbrExt=1 combine=1 WRITE(103,*) 'In spiralTrackV3: combine = ',combine ! viewSpiral=0 ! display=1 maxI = wsdat%maxi maxJ = wsdat%maxj IF ( (seedLat.EQ.0).OR.(seedLon.EQ.0) ) THEN i=NINT(REAL(maxI)/2.) j=NINT(REAL(maxJ)/2.) WRITE(103,*) 'In spiralTrackV3, i=NINT(maxI/2.) =',i WRITE(103,*) 'In spiralTrackV3, j=NINT(maxJ/2.) =',j ELSE i=1 j=1 DO WHILE ( (wsdat%lat(1,j).LT.seedLat).AND.(j.LT.wsdat%maxj) ) 40.PAR !Improve with SWAN's indice identification j=j+1 END DO DO WHILE ( (wsdat%lon(1,i).LT.seedLon).AND.(i.LT.wsdat%maxi) ) i=i+1 END DO END IF ! In case center point is land point... IF (wsdat%par(i,j)%checked.EQ.-1) THEN oldJ=j DO WHILE (wsdat%par(i,j)%checked.EQ.-1) j=j+1 IF (j.EQ.maxJ) THEN j=oldJ i=i+1 oldJ=oldJ+1 END IF END DO END IF ! Compute distance in km between 2 grid points (at equator) deltaLat=(wsdat%lat(i,j)-wsdat%lat(i,j-1))*111.18 ! Starts the spiral ! Intitiate variables horizStepCount=0 vertStepCount=0 way='R' first=.TRUE. checkCount=1 maxSys=0 landPts=0 minLat=MINVAL(wsdat%lat) maxLat=MAXVAL(wsdat%lat) minLon=MINVAL(wsdat%lon) maxLon=MAXVAL(wsdat%lon) ! rawWsdat=wsdat horizBorder=0 vertBorder=0 DO WHILE (checkCount.LE.(maxI*maxJ-1) ) ! From the direction (way) we were going before, find which direction we ! are going now and how many 'step' we need to take ! WRITE(103,*) 'Call findWay...' CALL findWay(way, horizStepCount, vertStepCount, & vertBorder, horizBorder, stepCount) ! WRITE(103,*) 'way,stepCount =',way,stepCount IF (first) THEN m=0 DO k=1,LENGTH(wsdat%par(i,j)%hs, & SIZE(wsdat%par(i,j)%hs),9999.) IF ( (wsdat%par(i,j)%hs(k).EQ.0.).AND. & (wsdat%par(i,j)%tp(k).EQ.0.) ) THEN wsdat%par(i,j)%sys(k)=-1 ELSE m=m+1 wsdat%par(i,j)%sys(k)=m END IF END DO wsdat%par(i,j)%checked=1 checkCount=checkCount+1 first=.FALSE. END IF DO sc = 1, stepCount ! WRITE(103,*) 'Call findNext...' CALL findNext (i,j,maxI,maxJ,way,vertBorder,horizBorder) ! WRITE(103,*) i,j IF ( wsdat%par(i,j)%checked.EQ.-1 ) THEN ! Land point is one of our grid points, so we need to update counter checkCount=checkCount+1 landPts=landPts+1 ! So that we don't count the land points twice.... wsdat%par(i,j)%checked=-2 ELSE IF ( wsdat%par(i,j)%checked.EQ.0 ) THEN ! Hasn't been checked yet and is not land point checkCount=checkCount+1 ! WRITE(103,*) 'Call findSys...' CALL findSys(i, j, wsdat, maxSys, ngbrExt, maxI, maxJ, & freqKnob, dirKnob) END IF END DO END DO ! wetPts% of wet points maxPts=NINT(wetPts*(maxI*maxJ-1)) WRITE(103,*) 'Call combineWaveSystems...' ! [sys,wsdat]=combineWaveSystems(mdat,maxSys,maxPts); CALL combineWaveSystems(wsdat,maxSys,maxPts,maxI,maxJ, & freqKnob,dirKnob,hsKnob,combine,sys) ! WRITE(103,*) 'Just after combineWaveSystems, end of spiralTrackV3' ! WRITE(103,*) 'SIZE(sys) =',SIZE(sys) RETURN END SUBROUTINE spiralTrackV3 !*********************************************************************** ! SUBROUTINE timeTrackingV2 (sysA ,maxSys ,tpTimeKnob , & dirTimeKnob,ts1 ,maxGroup ) ! * !*********************************************************************** ! USE PARTDAT USE MATFUNC ! IMPLICIT NONE ! ! 0. Authors ! ! 40.81: Andre van der Westhuysen and Jeff Hanson ! ! 1. Updates ! ! 40.81, Jul. 11: Origination, based on Matlab code by Jeff Hanson ! and Eve-Marie Devaliere ! ! 2. Purpose ! ! Performs the time tracking of the systems identified within ! the subroutine spiralTrackV3. ! ! 3. Method ! ! - ! ! 4. Argument variables ! ! Note: perKnob, dirKnob in Matlab version replaced by tpTimeKnob, dirTimeKnob! ! ! sysA TYPE(timsys) in/out Final set of spatially and temporally tracked systems ! dirTimeKnob Real input Parameter in direction for combining fields in time ! tpTimeKnob Real input Parameter in period for combining fields in time ! ts1 Int input Time step to which default grp values are associated ! maxSys Int arr input Total number of systems per time level ! maxGroup Int output Maximum number of wave systems ("groups") tracked in time ! TYPE(timsys), POINTER :: sysA(:) 40.PAR INTEGER, ALLOCATABLE :: maxSys(:) 40.PAR REAL :: dirTimeKnob, tpTimeKnob 40.PAR INTEGER :: ts1, maxGroup 40.PAR INTENT (IN) maxSys, tpTimeKnob, dirTimeKnob, ts1 INTENT (IN OUT) sysA INTENT (OUT) maxGroup ! ! 5. Local variables ! ! ic Int Counter for wave systems ! INTEGER :: leng, l, i, ii, j, k, idir, numSys, & counter, new, DIFSIZE, tpMinInd, dirMinInd, used, ok REAL :: Tb, deltaPer, deltaDir, tpMinVal, dirMinVal, & dirForTpMin, tpForDirMin REAL, ALLOCATABLE :: sysOrdered(:), DIFARR(:), TEMP(:), dirs(:) INTEGER, ALLOCATABLE :: indSorted(:), alreadyUsed(:), allInd(:) INTEGER, ALLOCATABLE :: ind(:), ind2(:) ! ! 6. Subroutines used ! ! SORT ! SETDIFF ! ! 7. Subroutines calling ! ! waveTracking_NWS_V2 ! ! 8. Error messages ! ! - ! ! 9. Remarks ! ! 10. Structure ! ! - ! ! 11. Source text ! ! Associate default grp value to time step 1 WRITE(103,*) 'TIME TRACKING' ! WRITE(6,*) 'ASSOCIATED(sysA) =',ASSOCIATED(sysA) ! WRITE(6,*) 'SIZE(sysA) =',SIZE(sysA) ! WRITE(6,*) 'nPoints =',sysA(1)%sys(1)%nPoints ! WRITE(6,*) 'nPoints =',sysA(1)%sys(2)%nPoints ! WRITE(6,*) 'nPoints =',sysA(1)%sys(3)%nPoints ! WRITE(6,*) 'nPoints =',sysA(1)%sys(4)%nPoints ! WRITE(103,*) 'ASSOCIATED =',ASSOCIATED(sysA(1)%sys(1)%hs) ! WRITE(103,*) 'ASSOCIATED =',ASSOCIATED(sysA(1)%sys(2)%hs) ! WRITE(103,*) 'ASSOCIATED =',ASSOCIATED(sysA(1)%sys(3)%hs) ! WRITE(103,*) 'ASSOCIATED =',ASSOCIATED(sysA(1)%sys(4)%hs) ! WRITE(103,*) SIZE(sysA(1)%sys(1)%hs) ! WRITE(103,*) SIZE(sysA(1)%sys(2)%hs) ! WRITE(103,*) SIZE(sysA(1)%sys(3)%hs) ! WRITE(103,*) SIZE(sysA(1)%sys(4)%hs) ! DO i = 1,maxSys(1) ! WRITE(103,*) 'sysA(1)%sys(i)%hs =', ! & sysA(1)%sys(i)%hs ! END DO WRITE(103,*) 'Inside timeTrackingV2: SIZE(sysA(1)%sys) =', & SIZE(sysA(1)%sys) WRITE(103,*) 'Inside timeTrackingV2: maxSys(1) =',maxSys(1) WRITE(103,*) 'ts1 = ',ts1 ! Set up the group number array for the first time level to be tracked ALLOCATE( sysOrdered(maxSys(ts1)) ) ALLOCATE( indSorted(maxSys(ts1)) ) CALL SORT (REAL(sysA(ts1)%sys(1:maxSys(ts1))%nPoints), & maxSys(ts1),sysOrdered,indSorted,'D') sysA(ts1)%sys(1:maxSys(ts1)) = sysA(ts1)%sys(indSorted) DEALLOCATE(sysOrdered) DEALLOCATE(indSorted) ! Initialize groups ! Optimize ??? sysA(ts1)%sys(:)%grp = 9999 ! Optimize ??? DO i = 1, maxSys(ts1) sysA(ts1)%sys(i)%grp = i END DO maxGroup = maxSys(ts1) i = ts1 ! Loop over all time levels to track systems in time WRITE(103,*) 'Number of time levels = ',SIZE(sysA) DO i = (ts1+1), SIZE(sysA) WRITE(103,*) 'TS = ',i ! Initialize groups ! Optimize ??? sysA(i)%sys(:)%grp = 9999 ! Optimize ??? counter = 0 leng = LENGTH(REAL(sysA(i-1)%sys(:)%grp), & SIZE(sysA(i-1)%sys(:)%grp),REAL(9999)) ALLOCATE( alreadyUsed(leng) ) ALLOCATE( allInd(leng) ) alreadyUsed(:) = 0 allInd(:) = sysA(i-1)%sys(1:leng)%grp DO j = 1, maxSys(i) WRITE(103,*) 'System no. ',j,' of ',maxSys(i) new = 0 sysA(i)%sys(j)%grp = 9999 ??? Now redundant ??? ! Compute deltas Tb = sysA(i)%sys(j)%tpMean deltaPer = -0.06*Tb + 2 + tpTimeKnob deltaDir = (-Tb + (25 + 10*dirTimeKnob)) CALL SETDIFF(REAL( (/allInd, 0/) ),SIZE( (/allInd, 0/) ), & REAL(alreadyUsed),SIZE(alreadyUsed),DIFARR,DIFSIZE) ALLOCATE( ind(DIFSIZE) ) ind(:) = NINT(DIFARR) DEALLOCATE(DIFARR) ALLOCATE( ind2(SIZE(ind)) ) DO ii = 1, SIZE(ind) ind2(ii) = FINDFIRST(REAL(allInd),SIZE(allInd), & REAL(ind(ii))) END DO IF (SIZE(ind).GT.0) THEN IF (sysA(i)%sys(j)%tpMean.NE.9999.) THEN ALLOCATE( TEMP(SIZE(ind2)) ) TEMP = ABS( sysA(i)%sys(j)%tpMean - & sysA(i-1)%sys(ind2(:))%tpMean ) WRITE(103,*) 'tpMinVal list =', TEMP tpMinVal = MINVAL(TEMP) tpMinInd = FINDFIRST(TEMP,SIZE(TEMP),tpMinVal) DEALLOCATE(TEMP) ALLOCATE( dirs(SIZE(ind)) ) dirs(:)=ABS( sysA(i)%sys(j)%dirMean - & sysA(i-1)%sys(ind2(:))%dirMean ) ! Deal with wrap around DO idir = 1, SIZE(dirs) IF (dirs(idir).GE.180.) dirs(idir)=dirs(idir)-180. END DO WRITE(103,*) 'dirMinVal list =', dirs dirMinVal = MINVAL(dirs) dirMinInd = FINDFIRST(dirs,SIZE(dirs),dirMinVal) DEALLOCATE(dirs) WRITE(103,*) 'deltaPer,deltaDir =',deltaPer,deltaDir WRITE(103,*) 'tpMinVal,dirMinVal =',tpMinVal,dirMinVal WRITE(103,*) 'tpMinInd,dirMinInd =',tpMinInd,dirMinInd IF ( (tpMinInd.EQ.dirMinInd).AND. & (tpMinVal.LE.deltaPer).AND. & (dirMinVal.LE.deltaDir) ) THEN ! Success counter = counter+1 sysA(i)%sys(j)%grp = & sysA(i-1)%sys(ind2(dirMinInd))%grp alreadyUsed(counter) = sysA(i)%sys(j)%grp WRITE(103,*) 'Case 1: matched this ts (',i, & ') sys ',sysA(i)%sys(j)%sysInd,' (tp=', & sysA(i)%sys(j)%tpMean,' dir=', & sysA(i)%sys(j)%dirMean,') with sys ', & sysA(i-1)%sys(ind2(dirMinInd))%sysInd, ??? Is this "ind" or "ind2" ??? & ' (tp=',sysA(i-1)%sys(ind2(dirMinInd))%tpMean, & ' and dir=', & sysA(i-1)%sys(ind2(dirMinInd))%dirMean, & ') in grp=', & sysA(i-1)%sys(ind2(dirMinInd))%grp WRITE(103,*) 'Added ',ind2(dirMinInd), & ' in array *alreadyUsed*' ELSE dirForTpMin = ABS( sysA(i)%sys(j)%dirMean - & sysA(i-1)%sys(ind2(tpMinInd))%dirMean ) IF (dirForTpMin.GT.180.) & dirForTpMin = 360.-dirForTpMin IF ( (tpMinVal.LE.deltaPer).AND. & (dirForTpMin.LE.deltaDir) ) THEN ! Success counter = counter+1 sysA(i)%sys(j)%grp = & sysA(i-1)%sys(ind2(tpMinInd))%grp alreadyUsed(counter)=sysA(i)%sys(j)%grp WRITE(103,*) 'Case 2: matched this ts (',i, & ') sys ',sysA(i)%sys(j)%sysInd,' (tp=', & sysA(i)%sys(j)%tpMean,' dir=', & sysA(i)%sys(j)%dirMean,') with sys ', & sysA(i-1)%sys(ind2(tpMinInd))%sysInd, ??? Is this "ind" or "ind2" ??? & ' (tp=', & sysA(i-1)%sys(ind2(tpMinInd))%tpMean, & ' and dir=', & sysA(i-1)%sys(ind2(tpMinInd))%dirMean, & ') in grp=', & sysA(i-1)%sys(ind2(tpMinInd))%grp WRITE(103,*) 'Added ',ind2(tpMinInd), & ' in array *alreadyUsed*' ELSE tpForDirMin = ABS( sysA(i)%sys(j)%tpMean - & sysA(i-1)%sys(ind2(dirMinInd))%tpMean ) IF ((tpForDirMin.LE.deltaPer) .AND. & (dirMinVal.LE.deltaDir)) THEN ! Success counter = counter+1 sysA(i)%sys(j)%grp = & sysA(i-1)%sys(ind2(dirMinInd))%grp alreadyUsed(counter) = & sysA(i)%sys(j)%grp WRITE(103,*) 'Case 3: matched this ts (',i, & ') sys ',sysA(i)%sys(j)%sysInd, & ' (tp=',sysA(i)%sys(j)%tpMean, & ' dir=',sysA(i)%sys(j)%dirMean, & ') with sys ', & sysA(i-1)%sys(ind2(dirMinInd))%sysInd, ??? Is this "ind" or "ind2" ??? & ' (tp=', & sysA(i-1)%sys(ind2(dirMinInd))%tpMean, & ' and dir=', & sysA(i-1)%sys(ind2(dirMinInd))%dirMean, & ') in grp=', & sysA(i-1)%sys(ind2(dirMinInd))%grp WRITE(103,*) 'Added ',ind2(tpMinInd), & ' in array *alreadyUsed*' ELSE new = 1 END IF END IF END IF IF (new.EQ.1) THEN used = 0 DO k = 1, maxGroup ok = 1 ! Make sure it hasn't been used yet IF ((i.GT.2).AND. & (.NOT.ANY(alreadyUsed(:).EQ.k))) THEN DO l = 1, maxSys(i-1) IF (sysA(i-1)%sys(l)%grp.EQ.k) & ok = 0 END DO IF (ok.EQ.1) THEN sysA(i)%sys(j)%grp = k counter = counter+1; alreadyUsed(counter) = k used = 1 ! break; EXIT END IF END IF END DO IF (used.EQ.0) THEN maxGroup = maxGroup+1 sysA(i)%sys(j)%grp = maxGroup counter = counter+1 alreadyUsed(counter) = maxGroup END IF ! % disp(['NO GRP MATCH case 2 for ts (',num2str(i), ... ! % ') sys ',num2str(sys{i}(j).sysInd),' tp=', ... ! % num2str(sys{i}(j).tpMean),' dir=', ... ! % num2str(sys{i}(j).dirMean),' get grp num=', ... ! % num2str(sys{i}(j).grp)]); ! % ' and tp ts bef=',num2str([sys{i-1}.tpMean]), ... ! % ' and dir ts bef=',num2str([sys{i-1}.dirMean])]); WRITE(103,*) 'NO GRP MATCH case 2' END IF END IF ELSE WRITE(103,*) 'EMPTY IND' used = 0 DO k = 1, maxGroup ok = 1 ! Make sure it hasn't been used yet IF ((i.GT.2) .AND. & (.NOT.ANY(alreadyUsed(:).EQ.k)) ) THEN DO l = 1, maxSys(i-1) IF (sysA(i-1)%sys(l)%grp.EQ.k) & ok = 0 END DO IF (ok.EQ.1) THEN sysA(i)%sys(j)%grp = k used = 1 counter = counter+1 alreadyUsed(counter) = k ! break; EXIT END IF END IF END DO IF (used.EQ.0) THEN maxGroup = maxGroup+1 sysA(i)%sys(j)%grp = maxGroup counter = counter+1 alreadyUsed(counter) = maxGroup END IF ! % disp(['NO GRP MATCH case 3 for ts (',num2str(i),') sys ', ... ! % num2str(sys{i}(j).sysInd),' tp=',num2str(sys{i}(j).tpMean), ... ! % ' dir=',num2str(sys{i}(j).dirMean),' get grp num=', ... ! % num2str(sys{i}(j).grp)]); WRITE(103,*) 'NO GRP MATCH case 3' END IF ! clear ind ind2; DEALLOCATE(ind) DEALLOCATE(ind2) END DO ! %for j=1:length([sys{i}(:).sysInd]) ! % disp(['sys ',num2str(sys{i}(j).sysInd),' (tp=',num2str(sys{i}(j).tpMean), ... ! % ' dir=',num2str(sys{i}(j).dirMean),' grp=',num2str(sys{i}(j).grp),')']); ! %end DEALLOCATE(alreadyUsed) DEALLOCATE(allInd) END DO WRITE(103,*) 'End of timeTrackingV2: SIZE(sysA(1)%sys) =', & SIZE(sysA(1)%sys) WRITE(103,*) 'End of timeTrackingV2: maxSys(1) =',maxSys(1) WRITE(103,*) 'End of timeTrackingV2: sysA(1)%sys(1)%grp =', & sysA(1)%sys(1)%grp WRITE(103,*) 'End of timeTrackingV2: sysA(1)%sys(2)%grp =', & sysA(1)%sys(2)%grp RETURN END SUBROUTINE timeTrackingV2 END MODULE PARTMAIN