#include "w3macros.h" !/ ------------------------------------------------------------------- / MODULE W3STRKMD !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 03-Mar-2016 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ 29-Nov-2013 : Remove DOC control characters, !/ update MPI! to MPI/! (H. L. Tolman). ( version 4.15 ) !/ 26-Sep-2016 : Optimization updates (A. van der Westhuysen) !/ ( version 5.15 ) !/ 03-Mar-2016 : Optimization updates for INTERSECT, !/ UNION, UNIQUE, SORT, SETDIFF, FINDIJ !/ (S. Zieger, BoM Australia) ( version 5.16 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! Module containing data structures and subroutines for spatial and ! temporal tracking (part of wave partitioning). ! ! 2. Method : ! ! Read raw partitioning data. ! Perform tracking in space. ! Perform tracking in time. ! ! 3. Variables and types : ! ! 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. ! ! Name Type Description ! ---------------------------------------------------------------- ! 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 mean direction partitions ! dspr Real arr array of mean directional spread (one-sided) of partitions ! wf Real arr array of wind fraction ! 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) REAL :: tp(10) REAL :: dir(10) REAL :: dspr(10) ! REAL :: wf(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 ! ! 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*8 :: date 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 ! wfVal Real arr array of wind fraction values of matched systems ! TYPE mtchsys INTEGER :: sysVal(50) REAL :: tpVal(50) REAL :: dirVal(50) REAL :: hsVal(50) ! REAL :: wfVal(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) ! dspr Real arr mean directional spread field assoc with wave system (in 1d array) ! wf Real arr wind fraction 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 ! wfMean Real spatial mean wind fraction 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(:) REAL, POINTER :: dspr(:) ! REAL, POINTER :: wf(:) INTEGER, POINTER :: i(:) INTEGER, POINTER :: j(:) INTEGER, POINTER :: indx(:) REAL, POINTER :: lat(:) REAL, POINTER :: lon(:) INTEGER :: sysInd REAL :: hsMean REAL :: tpMean REAL :: dirMean ! REAL :: wfMean 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 ! ! sysmemory Der. type Structure to store key characteristics of systems over multiple ! time levels. Used during the time tracking routine. ! TYPE sysmemory INTEGER :: grp INTEGER :: nPoints INTEGER, POINTER :: indx(:) INTEGER :: updated INTEGER :: length REAL :: lonMean REAL :: latMean REAL :: tpMean REAL :: dirMean END TYPE sysmemory ! ! 4. Subroutines and functions used : ! ! a. Main subroutines for spatial/temporal tracking: ! ! 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 ! ! b. Auxiliary subroutines and functions for tracking: ! ! 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 ! func. mean_angleV2 compute the mean direction from array of directions ! findIJV4 Find indices of system "a" that lie over or next to system "b" ! ! c. Simple data manipulation (based on Matlab intrinsic functions): ! ! 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 ! ! 5. Called by : ! ! WW3_SYSTRK (main program) ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. 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 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 findIJV4 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) ! ! 9. Switches : ! ! !/SHRD Switch for shared / distributed memory architecture. ! !/MPI Id. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ CONTAINS !/ ------------------------------------------------------------------- / SUBROUTINE waveTracking_NWS_V2 (intype ,tmax , & tcur ,filename , & tstart ,tend , & dt ,ntint , & minlon ,maxlon , & minlat ,maxlat , & mxcwt ,mycwt , & dirKnob , & perKnob ,hsKnob , & wetPts ,seedLat , & seedLon ,dirTimeKnob, & tpTimeKnob ,paramFile , & sysA ,wsdat , & maxSys ,maxGroup ) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 4-Jan-2013 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE !/MPI !/MPI INCLUDE "mpif.h" ! ! 1. Purpose : ! ! Main subroutine of spatial and temporal tracking algorithm ! ! 2. Method ! ! (1) Read the raw partitioning output from one of two file formats: ! (a) "partRes" format of IFP-SWAN (intype=1), or ! (b) WW3 spectral bulletin format (intype=2). ! If intype=0, the partition data is read from memory (not activated yet). ! ! (2) Perform tracking in space by calling subroutine spiralTrackV3 ! (3) Perform tracking in time by calling subroutine timeTrackingV2 ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! intype Int input For coupling: Type of input (0 = memory; 1 = partRes file; 2 = WW3 part 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 ! tstart Char input Start time in raw partition file (if used) ! tend Char input End time in raw partition file (if used) ! minlon Real input Lower lon boundary of domain to be processed ! maxlon Real input Upper lon boundary of domain to be processed ! minlat Real input Lower lat boundary of domain to be processed ! maxlat Real input Upper lat boundary of domain to be processed ! 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 purging fields ! wetPts Real input Percentage of wet points for purging fields (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*50, paramFile*32 REAL :: dirKnob, perKnob, hsKnob, wetPts, seedLat, & seedLon, dirTimeKnob, tpTimeKnob REAL*8 :: tstart, tend INTEGER :: maxGroup, intype, tmax, tcur, ntint INTEGER, POINTER :: maxSys(:) TYPE(dat2d), POINTER :: wsdat(:) TYPE(timsys), POINTER :: sysA(:), sysAA(:) INTEGER :: NumConsSys, iConsSys REAL :: dt REAL :: minlon, maxlon, minlat, maxlat INTEGER :: mxcwt, mycwt ! 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, & minlon, maxlon, minlat, maxlat, & hsKnob, wetPts, seedLat, seedLon, & dirKnob, perKnob, dirTimeKnob, tpTimeKnob INTENT (OUT) maxGroup ! INTENT (IN OUT) wsdat, sysA, maxSys ! ! Local variables ! ---------------------------------------------------------------- ! 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 ! dspr0 Real Mean directional spread of partition, from input file ! wf0 Real wind fraction of partition, from input file (removed) ! 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; removed) ! tss Int. Time step counter ! t0 Int Index of first time step to compute ! LOGICAL :: file_exists, FLFORM, LOOP LOGICAL :: testout PARAMETER (testout = .FALSE.) CHARACTER :: dummy*10, dummyc*12 CHARACTER(LEN=10) :: VERPRT CHARACTER(LEN=35) :: IDSTR CHARACTER(LEN=78) :: headln1 CHARACTER(LEN=51) :: headln2 INTEGER :: line INTEGER, ALLOCATABLE :: ts(:), tmp_i4(:) REAL, ALLOCATABLE :: llat(:),llon(:),hs0(:), & tp0(:),dir0(:),dspr0(:),& wndSpd0(:),wndDir0(:) REAL*8, ALLOCATABLE :: date0(:),tmp_r8(:) INTEGER :: maxTs, t0, nout1, nout2, maxI, maxJ REAL, ALLOCATABLE :: mlon(:,:), mlat(:,:), tmp_r4(:) REAL, POINTER :: uniqueTim(:),uniqueLatraw(:),uniqueLonraw(:), & uniqueLat(:),uniqueLon(:) INTEGER :: ioerr,ierr, i, j, k, l, alreadyIn, ok, tss, tsA INTEGER :: maxPart, DATETIME(2) INTEGER :: tstep, iline, numpart, skipln, readln, filesize REAL :: x,y,wnd,wnddir REAL :: invar1, invar2, invar3, invar4 REAL :: invar5, invar6, invar7 REAL, ALLOCATABLE :: phs(:),ptp(:),pdir(:),pspr(:),pwf(:) ! current partition values REAL*8 :: date1, date2, ttest, ttemp INTEGER :: ic, leng, maxpartout ! Remove? REAL :: dx INTEGER :: latind1, latind2, lonind1, lonind2 REAL :: lonext, latext !/MPI INTEGER :: rank, irank, nproc, EXTENT, DOMSIZE, tag1, tag2 !/MPI/! INTEGER :: MPI_INT_DOMARR, MPI_REAL_DOMARR !/MPI INTEGER :: MPI_STATUS(MPI_STATUS_SIZE) !/MPI INTEGER :: REQ(16) !/MPI/! INTEGER :: ISTAT(MPI_STATUS_SIZE,16) !/MPI REAL :: COMMARR1(44) !/MPI INTEGER :: COMMARR2(11) ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! UNIQUE ! spiralTrackV3 ! timeTrackingV2 ! ! 5. Subroutines calling ! ! WW3_SYSTRK (main program) ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! See above ! ! 9. Switches : ! ! None defined yet. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/MPI CALL MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr) !/MPI CALL MPI_COMM_SIZE(MPI_COMM_WORLD, nproc, ierr) NULLIFY( sysA ) NULLIFY( maxSys ) ! Select input type for raw partitioning data IF ((intype.EQ.1).OR.(intype.EQ.2)) THEN ! Raw partitioning data is coming from an input file. ! Read file here, and set up 2d array wsdat with the data. t0 = 1 IF (intype.EQ.1) THEN !/MPI IF (rank.EQ.0) THEN ! Read partRes format file WRITE(20,*) 'Reading partRes partitioning file...' filesize = 7500000 ALLOCATE(ts(filesize)) ALLOCATE(llat(filesize)) ALLOCATE(llon(filesize)) ALLOCATE(hs0(filesize)) ALLOCATE(tp0(filesize)) ALLOCATE(dir0(filesize)) ALLOCATE(dspr0(filesize)) ! ALLOCATE(wf0(filesize)) ALLOCATE(wndSpd0(filesize)) ALLOCATE(wndDir0(filesize)) ALLOCATE(date0(filesize)) WRITE(20,*) '*** Max number of lines read from "partRes" ', & 'input file is = ',filesize,'!' WRITE(6,*) 'Reading partRes file...' INQUIRE(FILE=filename, EXIST=file_exists) IF (.NOT.file_exists) THEN WRITE(20,2001) WRITE(6,2001) STOP 1 END IF OPEN(unit=11,file=filename,status='old') line = 1 DO WHILE (.TRUE.) READ (11, *, END=113) dummyc,llat(line),llon(line), & ts(line),hs0(line),tp0(line),dir0(line), & wndSpd0(line),wndDir0(line),invar7 !partRes file does not contain the dspr variable dspr0(line) = 9999. ! wf0(line) = 9999. line = line+1 ENDDO 113 IERR = -1 CLOSE(11) line = line-1 WRITE(6,*) '... finished' ! DEALLOCATE(date0) !/MPI END IF ELSE IF (intype.EQ.2) THEN !/MPI IF (rank.EQ.0) THEN ! Read WW3 Spectral Partition format file ! Query input file to determine required array sizes INQUIRE(FILE=filename, EXIST=file_exists) IF (.NOT.file_exists) THEN WRITE(20,2001) WRITE(6,2001) STOP 1 END IF !/ ------------------------------------------------- !/ Test unformatted read !/ ------------------------------------------------- OPEN(UNIT=11,FILE=FILENAME,FORM='UNFORMATTED',STATUS='OLD',ACCESS='STREAM') READ(11,ERR=802,IOSTAT=IOERR) I CLOSE(11) !/ --- First four-byte integer could possibly be byte-swapped, ! if ww3_shel was compiled on a different architecture. --- K = SWAPI4(I) FLFORM = .NOT.(I.EQ.(LEN(IDSTR)+LEN(VERPRT)).OR.& K.EQ.(LEN(IDSTR)+LEN(VERPRT)) ) ! ======== COUNT LOOP =========== IF (FLFORM) THEN ! Input file in formatted ASCII WRITE(6,*) 'Reading formatted ASCII file...' OPEN(unit=11,file=filename,status='old') READ(11,'(78A)') headln1 IDSTR = headln1(1:LEN(IDSTR)) READ(11,'(78A)') headln1 READ(11,'(51A)') headln2 ELSE IF (K.EQ.(LEN(IDSTR)+LEN(VERPRT))) THEN !/ --- Stop here. The file appears to be endian encoded ! different from the native machine format. And, the ! compiler option will override support for FORTRAN ! convert statements convert='big_endian' or ! convert='little_endian'. --- WRITE(20,1200) WRITE(6,1200) STOP 1 ELSE ! Input file in unformatted binary WRITE(6,*) 'Reading binary formatted file...' OPEN(unit=11,file=filename,form='UNFORMATTED', & status='OLD') ENDIF REWIND(11) READ(11,ERR=802,IOSTAT=IOERR) IDSTR,VERPRT READ(11,ERR=802,IOSTAT=IOERR) headln1 READ(11,ERR=802,IOSTAT=IOERR) headln2 END IF !/ IF (IDSTR(1:9).ne.'WAVEWATCH') THEN CLOSE(11) WRITE(20,1300) WRITE(6,1300) STOP 1 ENDIF !/ ------------------------------------------------- !/ Skip to start time !/ ------------------------------------------------- skipln = 3 ttest = 0 DO WHILE (ttest.LT.tstart) IF (FLFORM) THEN READ (11,1000,ERR=802,END=112) date1,date2,x,y, & numpart,wnd,wnddir,invar6,invar7 !/del write(*,*) '0:',x,y,numpart skipln = skipln+1 ELSE READ (11,ERR=802,IOSTAT=IOERR) DATETIME,x,y, & dummy,numpart,invar1,wnd,wnddir, & invar5,invar6 ! write(*,*) '0:',DATETIME,numpart date1=dble(DATETIME(1)) date2=dble(DATETIME(2)) END IF ttest = date1 + date2*1.0E-6 IF (FLFORM) THEN DO line = 1,numpart+1 READ(11,1010,END=111,ERR=802,IOSTAT=IOERR) & invar1,invar2,invar3,invar4 ! write(*,*) '0+:',line,numpart+1,invar1,invar2,invar3,invar4 skipln = skipln+1 END DO ELSE DO line = 1,numpart+1 READ (11,ERR=802,IOSTAT=IOERR) iline,invar1, & invar2,invar3,invar4,invar5,invar6 ! write(*,*) '0+:',line,iline,invar1,invar2,invar3,invar4,invar5,invar6 END DO END IF END DO skipln = skipln-numpart-1-1 !/ ------------------------------------------------- ! Read file for ntint time levels !/ ------------------------------------------------- readln = numpart tstep = 1 ttemp = tstart maxPart = numpart DO WHILE (tstep.LE.ntint) IF (readln.GT.0) THEN IF (FLFORM) THEN READ (11,1000,ERR=802,END=111) date1,date2,x,y, & numpart,wnd,wnddir,invar6,invar7 ELSE READ (11,END=111,ERR=802,IOSTAT=IOERR) DATETIME, & x,y,dummy,numpart,wnd,wnddir,invar5,invar6,invar7 ! write(*,*) '1:',numpart,x,y date1=dble(DATETIME(1)) date2=dble(DATETIME(2)) END IF maxPart = MAX(maxPart,numpart) END IF ttest = date1 + date2*1.E-6 IF (ttest.GT.ttemp) THEN tstep = tstep+1 ttemp = ttest IF (tstep.GT.ntint) EXIT END IF IF (FLFORM) THEN DO line = 1,numpart+1 READ (11,1010,END=111,ERR=802,IOSTAT=IOERR) & invar1,invar2,invar3,invar4 !/del write(*,'(A,2I6,4F7.2)') '1+:',line,numpart+1,invar1,invar2,invar3,invar4 readln = readln+1 END DO ELSE DO line = 1,numpart+1 READ (11,END=111,ERR=802,IOSTAT=IOERR) iline,invar1,& invar2,invar3,invar4,invar5,invar6 readln = readln+1 END DO END IF ENDDO 111 CONTINUE CLOSE(11) ! ===== END COUNT LOOP ===== ! ===== START READ LOOP ===== ALLOCATE(ts(readln)) ALLOCATE(llat(readln)) ALLOCATE(llon(readln)) ALLOCATE(hs0(readln)) ALLOCATE(tp0(readln)) ALLOCATE(dir0(readln)) ALLOCATE(dspr0(readln)) ! ALLOCATE(wf0(readln)) ALLOCATE(wndSpd0(readln)) ALLOCATE(wndDir0(readln)) ALLOCATE(date0(readln)) ts(1:readln) = -1 llat(1:readln) = 9999. llon(1:readln) = 9999. hs0(1:readln) = 9999. tp0(1:readln) = 9999. dir0(1:readln) = 9999. dspr0(1:readln) = 9999. IF (FLFORM) THEN OPEN(unit=11,file=filename,status='old') ELSE OPEN(unit=11,file=filename,status='old', & form='unformatted') END IF line = 1 tstep = 1 !/ ------------------------------------------------- !/ Skip to start time !/ ------------------------------------------------- IF (FLFORM) THEN DO i = 1,skipln READ (11, *) END DO ELSE ! --- Repeat from above since access='DIRECT' ! does not support fseek and ftell. --- READ(11,END=112,ERR=802,IOSTAT=IOERR) IDSTR,VERPRT READ(11,END=112,ERR=802,IOSTAT=IOERR) headln1 READ(11,END=112,ERR=802,IOSTAT=IOERR) headln2 !/ --- allocate buffer for all partition parameters !/ for a single grid point --- IF (.NOT.ALLOCATED(PHS)) ALLOCATE(PHS(maxPart)) IF (.NOT.ALLOCATED(PTP)) ALLOCATE(PTP(maxPart)) IF (.NOT.ALLOCATED(PDIR)) ALLOCATE(PDIR(maxPart)) IF (.NOT.ALLOCATED(PSPR)) ALLOCATE(PSPR(maxPart)) IF (.NOT.ALLOCATED(PWF)) ALLOCATE(PWF(maxPart)) ttest = 0 DO WHILE (ttest.LT.tstart) READ (11,END=112,ERR=802,IOSTAT=IOERR) DATETIME, & invar1,invar2,dummy,numpart,invar3, & invar4,invar5,invar6,invar7 date1=dble(DATETIME(1)) date2=dble(DATETIME(2)) ttest = date1 + date2*1.0E-6 !/ --- reset buffer --- PHS(:) = 0. PTP(:) = 0. PDIR(:) = 0. PSPR(:) = 0. PWF(:) = 0. !/ --- fill buffer with partition data --- READ (11,END=112,ERR=802,IOSTAT=IOERR) iline,invar1, & invar2,invar3,invar4,invar5,invar6 DO i = 1,numpart READ (11,END=112,ERR=802,IOSTAT=IOERR) iline, & phs(i),ptp(i),invar3,pdir(i),pspr(i),pwf(i) END DO END DO !/ --- move buffer content to data array --- DO i=1,numpart hs0(line) = phs(i) tp0(line) = ptp(i) dir0(line) = pdir(i) dspr0(line) = pspr(i) date0(line) = date1 + date2*1.0E-6 ts(line) = tstep llat(line) = x llon(line) = y wndSpd0(line) = wnd wndDir0(line) = wnddir line = line + 1 END DO END IF !/ ------------------------------------------------- ! Read file for ntint time levels !/ ------------------------------------------------- ttemp = tstart DO WHILE (line.LE.readln) IF (FLFORM) THEN READ (11,1000,END=112) date1,date2,x,y,numpart, & wnd,wnddir,invar6,invar7 ELSE READ (11,ERR=802,IOSTAT=IOERR) DATETIME,x,y, & dummy,numpart,wnd,wnddir,invar5,invar6,invar7 date1=dble(DATETIME(1)) date2=dble(DATETIME(2)) END IF ttest = date1 + date2*1.0E-6 IF (ttest.GT.ttemp) THEN tstep = tstep+1 ttemp = ttest IF (tstep.GT.ntint) EXIT END IF IF (FLFORM) THEN READ (11,1010,END=112) invar1,invar2,invar3,invar4 ! Skip total integral parameters DO i = 1,numpart IF (line.LE.readln) THEN READ (11,1010,END=112) hs0(line),tp0(line), & dir0(line),dspr0(line) date0(line) = ttest ts(line) = tstep llat(line) = x llon(line) = y wndSpd0(line) = wnd wndDir0(line) = wnddir line = line+1 END IF END DO ELSE READ (11,ERR=802,IOSTAT=IOERR) k,invar1,invar2, & invar3,invar4,invar5 DO i = 1,numpart IF (line.LE.readln) THEN READ (11,END=112,ERR=802,IOSTAT=IOERR) k, & hs0(line),tp0(line),invar3,dir0(line), & dspr0(line) date0(line) = ttest ts(line) = tstep llat(line) = x llon(line) = y wndSpd0(line) = wnd wndDir0(line) = wnddir line = line+1 END IF END DO END IF END DO 110 IERR = -1 CLOSE(11) 112 CONTINUE IF (line.EQ.1) THEN WRITE(20,2002) WRITE(6,2002) STOP 1 END IF CLOSE(11) ! ===== READ LOOP FINISHED ===== LINE=LINE-1 WRITE(6,*) '... finished' IF (ttest.LT.tstart) THEN WRITE(20,2003) TSTART WRITE(6,2003) TSTART STOP 1 END IF IF (ALLOCATED(PHS)) DEALLOCATE(PHS) IF (ALLOCATED(PTP)) DEALLOCATE(PTP) IF (ALLOCATED(PDIR)) DEALLOCATE(PDIR) IF (ALLOCATED(PSPR)) DEALLOCATE(PSPR) IF (ALLOCATED(PWF)) DEALLOCATE(PWF) !/MPI END IF END IF !/MPI IF (rank.EQ.0) THEN ! Find unique time steps (and sort in ascending order) CALL UNIQUE(REAL(ts(1:line)),line,uniqueTim,maxTs) ! Find unique lat and lon values (and sort in ascending order) CALL UNIQUE(llat(1:line),SIZE(llat(1:line)),uniqueLatraw,nout1) CALL UNIQUE(llon(1:line),SIZE(llon(1:line)),uniqueLonraw,nout2) !--042916----------------------- ! ! Redefine uniqueLatraw and uniqueLonrawto based on domain extent WRITE(20,*) 'uniqueLatraw(:) =', uniqueLatraw(:) WRITE(20,*) 'uniqueLonraw(:) =', uniqueLonraw(:) WRITE(20,*) 'No. increments: Longitude, Latitue =', mxcwt, mycwt DEALLOCATE(uniqueLatraw) DEALLOCATE(uniqueLonraw) ALLOCATE(uniqueLatraw(mycwt+1)) ALLOCATE(uniqueLonraw(mxcwt+1)) DO i = 1,(mycwt+1) uniqueLatraw(i) = minlat + & (REAL(i)-1)/REAL(mycwt)*(maxlat-minlat) END DO DO i = 1,(mxcwt+1) uniqueLonraw(i) = minlon + & (REAL(i)-1)/REAL(mxcwt)*(maxlon-minlon) END DO WRITE(20,*) 'uniqueLatraw(:) =', uniqueLatraw(:) WRITE(20,*) 'uniqueLonraw(:) =', uniqueLonraw(:) ! !--042916----------------------- ! Filter out lats and lons outside of domain of interest DO latind1 = 1,SIZE(uniqueLatraw) IF (uniqueLatraw(latind1).GE.minlat) EXIT END DO DO latind2 = SIZE(uniqueLatraw),1,-1 IF (uniqueLatraw(latind2).LE.maxlat) EXIT END DO DO lonind1 = 1,SIZE(uniqueLonraw) IF (uniqueLonraw(lonind1).GE.minlon) EXIT END DO DO lonind2 = SIZE(uniqueLonraw),1,-1 IF (uniqueLonraw(lonind2).LE.maxlon) EXIT END DO WRITE(20,*) 'latind1, latind2, lonind1, lonind2 =', & latind1, latind2, lonind1, lonind2 IF ((latind1.GE.latind2).OR.(lonind1.GE.lonind2)) THEN WRITE(20,1400) WRITE(6,1400) STOP 1 END IF NULLIFY(uniqueLat) NULLIFY(uniqueLon) ALLOCATE(uniqueLat(latind2-latind1+1)) ALLOCATE(uniqueLon(lonind2-lonind1+1)) uniqueLat = uniqueLatraw(latind1:latind2) uniqueLon = uniqueLonraw(lonind1:lonind2) WRITE(20,*) 'In waveTracking_NWS_V2: Longitude range =', & uniqueLon(1), uniqueLon(SIZE(uniqueLon)) WRITE(20,*) ' Latitude range =', & uniqueLat(1), uniqueLat(SIZE(uniqueLat)) ! Map is transposed (rotated by 90 deg), so that: ! I (matrix row) represents Longitute ! J (matrix column) represents Latitude ! i.e. from this point onwards the indices (i,j) represents Cart. coordinates 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 !/MPI END IF !/MPI CALL MPI_BCAST(maxI,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR) !/MPI CALL MPI_BCAST(maxJ,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR) !/MPI CALL MPI_BCAST(maxTs,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR) ! Allocate the wsdat structure !/MPI IF (rank.EQ.0) THEN WRITE(20,*) 'Allocating wsdat...' !/MPI END IF NULLIFY(wsdat) ALLOCATE(wsdat(maxTs)) !/MPI IF (rank.EQ.0) THEN WRITE(20,*) 'SIZE(wsdat) = ',SIZE(wsdat) !/MPI END IF ! Allocate and initialize the wsdat array !/MPI IF (rank.EQ.0) THEN 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)) DO j = 1,maxJ DO i = 1,maxI 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)%dspr(1:10)=9999. ! wsdat(tsA)%par(i,j)%wf(1:10)=9999. wsdat(tsA)%par(i,j)%ipart(1:10)=0 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 END DO END DO END DO ! Assign to each line in partition file an entry in wsdat ! At each time step each point contains all numpart partitions. ! Only store the first 10 partitions. l = 1 DO WHILE (l.LE.line) DO j = 1,maxJ DO i = 1,maxI !>042916 IF ( (llat(l).EQ.mlat(i,j)).AND. & !>042916 (llon(l).EQ.mlon(i,j)) ) THEN IF ( (ABS(llat(l)-mlat(i,j)).LT.1.E-2).AND. & (ABS(llon(l)-mlon(i,j)).LT.1.E-2) ) THEN ! WRITE(20,*) 'MATCHED! ',l,& ! llat(l),mlat(i,j),ABS(llat(l)-mlat(i,j)),& ! llon(l),mlon(i,j),ABS(llon(l)-mlon(i,j)) wsdat(ts(l))%lat(i,j) = llat(l) wsdat(ts(l))%lon(i,j) = llon(l) ! --- Find ALL partition values associated with ! lat(i,j) and lon(i,j). Keep list index l ! fixed and recycle iline as variable index. --- iline = l k = 1 DO WHILE ( & ABS(wsdat(ts(l))%lat(i,j)-llat(iline)).LT.1.E-3 & .AND.ABS(wsdat(ts(l))%lon(i,j)-llon(iline)).LT.1.E-3 ) IF (k.LE.10) THEN wsdat(ts(iline))%par(i,j)%ipart(k) = k wsdat(ts(iline))%par(i,j)%hs(k) = hs0(iline) wsdat(ts(iline))%par(i,j)%tp(k) = tp0(iline) wsdat(ts(iline))%par(i,j)%dir(k) = dir0(iline) wsdat(ts(iline))%par(i,j)%dspr(k) = dspr0(iline) ! wsdat(ts(k))%par(i,j)%wf(k) = wf0(l) IF (k.EQ.1) THEN wsdat(ts(iline))%date = date0(iline) wsdat(ts(iline))%wnd(i,j)%wdir = wndDir0(iline) wsdat(ts(iline))%wnd(i,j)%wspd = wndSpd0(iline) wsdat(ts(iline))%par(i,j)%checked = 0 END IF END IF k = k + 1 iline = iline + 1 if (iline.GT.line) EXIT END DO ! --- Account for increment at the end of loop (400 CONTINUE) ! and go one element back in list because of increment. --- l = iline-1 GOTO 400 END IF END DO END DO 400 CONTINUE IF (l+1.le.line) THEN IF (ts(l).LT.ts(l+1)) THEN K = line-l ! --- With each time step completed, deallocate processed 1:l ! elements from 1d array. Create a temporary array size of ! (l+1:line) with k elements and reallocate original array. --- IF (ALLOCATED(tmp_i4)) DEALLOCATE(tmp_i4) ! --- REALLOCATE(integer arrays) --- ALLOCATE(tmp_i4(k)) tmp_i4(1:k) = ts((l+1):line) DEALLOCATE(ts) ALLOCATE(ts(k)) ts(1:k) = tmp_i4(1:k) DEALLOCATE(tmp_i4) ! --- REALLOCATE(double precision arrays) --- IF (ALLOCATED(tmp_r8)) DEALLOCATE(tmp_r8) ALLOCATE(tmp_r8(k)) tmp_r8(1:k) = date0((l+1):line) DEALLOCATE(date0) ALLOCATE(date0(k)) date0(1:k) = tmp_r8(1:k) DEALLOCATE(tmp_r8) ! --- REALLOCATE(single precision arrays) --- IF (ALLOCATED(tmp_r4)) DEALLOCATE(tmp_r4) ALLOCATE(tmp_r4(k)) tmp_r4(1:k) = llat((l+1):line) DEALLOCATE(llat) ALLOCATE(llat(k)) llat(1:k) = tmp_r4(1:k) tmp_r4(1:k) = llon((l+1):line) DEALLOCATE(llon) ALLOCATE(llon(k)) llon(1:k) = tmp_r4(1:k) tmp_r4(1:k) = hs0((l+1):line) DEALLOCATE(hs0) ALLOCATE(hs0(k)) hs0(1:k) = tmp_r4(1:k) tmp_r4(1:k) = tp0((l+1):line) DEALLOCATE(tp0) ALLOCATE(tp0(k)) tp0(1:k) = tmp_r4(1:k) tmp_r4(1:k) = dir0((l+1):line) DEALLOCATE(dir0) ALLOCATE(dir0(k)) dir0(1:k) = tmp_r4(1:k) tmp_r4(1:k) = dspr0((l+1):line) DEALLOCATE(dspr0) ALLOCATE(dspr0(k)) dspr0(1:k) = tmp_r4(1:k) tmp_r4(1:k) = wndSpd0((l+1):line) DEALLOCATE(wndSpd0) ALLOCATE(wndSpd0(k)) wndSpd0(1:k) = tmp_r4(1:k) tmp_r4(1:k) = wndDir0((l+1):line) DEALLOCATE(wndDir0) ALLOCATE(wndDir0(k)) wndDir0(1:k) = tmp_r4(1:k) DEALLOCATE(tmp_r4) line = k l = 0 END IF END IF l = l + 1 END DO ! IF (ALLOCATED(ts)) DEALLOCATE(ts) IF (ALLOCATED(llat)) DEALLOCATE(llat) IF (ALLOCATED(llon)) DEALLOCATE(llon) IF (ALLOCATED(mlat)) DEALLOCATE(mlat) IF (ALLOCATED(mlon)) DEALLOCATE(mlon) IF (ALLOCATED(date0)) DEALLOCATE(date0) IF (ALLOCATED(hs0)) DEALLOCATE(hs0) IF (ALLOCATED(tp0)) DEALLOCATE(tp0) IF (ALLOCATED(dir0)) DEALLOCATE(dir0) IF (ALLOCATED(dspr0)) DEALLOCATE(dspr0) ! IF (ALLOCATED(wf0)) DEALLOCATE(wf0) IF (ALLOCATED(wndSpd0)) DEALLOCATE(wndSpd0) IF (ALLOCATED(wndDir0)) DEALLOCATE(wndDir0) !/MPI END IF !/MPI/! Communicate the wsdat entries from rank=0 to other ranks !/MPI DO tsA = t0,maxTs !/MPI irank = MOD((tsA-t0),MIN(nproc,maxTS)) !/MPI/! WRITE(20,*) 'Rank,irank=',rank,irank !/MPI IF (irank.NE.0) THEN !/MPI/! WRITE(20,*) 'Communicating for Rank,irank=',rank,irank !/MPI !/MPI IF (rank.EQ.irank) THEN !/MPI ALLOCATE(wsdat(tsA)%lat(maxI,maxJ)) !/MPI ALLOCATE(wsdat(tsA)%lon(maxI,maxJ)) !/MPI ALLOCATE(wsdat(tsA)%par(maxI,maxJ)) !/MPI ALLOCATE(wsdat(tsA)%wnd(maxI,maxJ)) !/MPI !/MPI DO j = 1,maxJ !/MPI DO i = 1,maxI !/MPI wsdat(tsA)%maxi=maxI !/MPI wsdat(tsA)%maxj=maxJ !/MPI wsdat(tsA)%par(i,j)%hs(1:10)=9999. !/MPI wsdat(tsA)%par(i,j)%tp(1:10)=9999. !/MPI wsdat(tsA)%par(i,j)%dir(1:10)=9999. !/MPI wsdat(tsA)%par(i,j)%dspr(1:10)=9999. !/MPI wsdat(tsA)%par(i,j)%ipart(1:10)=0 !/MPI wsdat(tsA)%par(i,j)%sys(1:10)=9999 ! 40.PAR Increase this array, or make allocatable !/MPI wsdat(tsA)%par(i,j)%ngbrSys(1:50)=9999 !/MPI wsdat(tsA)%wnd(i,j)%wdir=9999. !/MPI wsdat(tsA)%wnd(i,j)%wspd=9999. !/MPI wsdat(tsA)%par(i,j)%checked=-1 !/MPI END DO !/MPI END DO !/MPI END IF !/MPI !/MPI DO j = 1,maxJ !/MPI DO i = 1,maxI !/MPI tag1 = ((j-1)*maxI+i)*10 !/MPI !/MPI IF (rank.EQ.0) THEN !/MPI/! WRITE(6,*) '>> Sending: rank,irank,tag1=', & !/MPI/! rank,irank,(tag1+1) !/MPI COMMARR1 = (/wsdat(tsA)%par(i,j)%hs(:), & !/MPI wsdat(tsA)%par(i,j)%tp(:), & !/MPI wsdat(tsA)%par(i,j)%dir(:), & !/MPI wsdat(tsA)%par(i,j)%dspr(:), & !/MPI wsdat(tsA)%wnd(i,j)%wdir, & !/MPI wsdat(tsA)%wnd(i,j)%wspd, & !/MPI wsdat(tsA)%lat(i,j), & !/MPI wsdat(tsA)%lon(i,j)/) !/MPI CALL MPI_SEND(COMMARR1,44,MPI_REAL,irank, & !/MPI (tag1+1),MPI_COMM_WORLD,IERR) !/MPI END IF !/MPI IF (rank.EQ.irank) THEN !/MPI/! WRITE(6,*) '<< Receiving: rank,irank,tag1=', & !/MPI/! rank,irank,(tag1+1) !/MPI CALL MPI_RECV(COMMARR1,44,MPI_REAL,0,(tag1+1), & !/MPI MPI_COMM_WORLD,MPI_STATUS,IERR) !/MPI wsdat(tsA)%par(i,j)%hs = COMMARR1(1:10) !/MPI wsdat(tsA)%par(i,j)%tp = COMMARR1(11:20) !/MPI wsdat(tsA)%par(i,j)%dir = COMMARR1(21:30) !/MPI wsdat(tsA)%par(i,j)%dspr = COMMARR1(31:40) !/MPI wsdat(tsA)%wnd(i,j)%wdir = COMMARR1(41) !/MPI wsdat(tsA)%wnd(i,j)%wspd = COMMARR1(42) !/MPI wsdat(tsA)%lat(i,j) = COMMARR1(43) !/MPI wsdat(tsA)%lon(i,j) = COMMARR1(44) !/MPI END IF !/MPI !/MPI IF (rank.EQ.0) THEN !/MPI CALL MPI_SEND(wsdat(tsA)%date,1, & !/MPI MPI_DOUBLE_PRECISION,irank, & !/MPI (tag1+2),MPI_COMM_WORLD,IERR) !/MPI END IF !/MPI IF (rank.EQ.irank) THEN !/MPI CALL MPI_RECV(wsdat(tsA)%date,1, & !/MPI MPI_DOUBLE_PRECISION,0,(tag1+2), & !/MPI MPI_COMM_WORLD,MPI_STATUS,IERR) !/MPI END IF !/MPI !/MPI IF (rank.EQ.0) THEN !/MPI/! WRITE(6,*) '>> Sending: rank,irank,tag1=', & !/MPI/! rank,irank,(tag1+3) !/MPI COMMARR2 = (/wsdat(tsA)%par(i,j)%ipart(:), & !/MPI wsdat(tsA)%par(i,j)%checked/) !/MPI CALL MPI_SEND(COMMARR2,11, & !/MPI MPI_INTEGER,irank,(tag1+3),MPI_COMM_WORLD,IERR) !/MPI END IF !/MPI IF (rank.EQ.irank) THEN !/MPI/! WRITE(6,*) '<< Receiving: rank,irank,tag1=', & !/MPI/! rank,irank,(tag1+3) !/MPI CALL MPI_RECV(COMMARR2,11, & !/MPI MPI_INTEGER,0,(tag1+3), & !/MPI MPI_COMM_WORLD,MPI_STATUS,IERR) !/MPI wsdat(tsA)%par(i,j)%ipart(:) = COMMARR2(1:10) !/MPI wsdat(tsA)%par(i,j)%checked = COMMARR2(11) !/MPI END IF !/MPI !/MPI END DO !/MPI END DO !/MPI END IF !/MPI END DO !/MPI !/MPI CALL MPI_Barrier(MPI_COMM_WORLD,IERR) !/MPI IF (rank.EQ.0) THEN ! ----*** Test Output *** -------------------------------------------------- IF (testout) THEN !-----RAW PARTITION output: Coordinates OPEN(unit=31,file='PART_COORD.OUT',status='unknown') WRITE(31,*) 'Longitude =' DO j = maxJ,1,-1 DO i = 1,maxI WRITE(31,'(F7.2)',ADVANCE='NO') wsdat(1)%lon(i,j) END DO WRITE(31,'(A)',ADVANCE='YES') '' END DO WRITE(31,*) 'Latitude = ' DO j = maxJ,1,-1 DO i = 1,maxI WRITE(31,'(F7.2)',ADVANCE='NO') wsdat(1)%lat(i,j) END DO WRITE(31,'(A)',ADVANCE='YES') '' END DO CLOSE(31) !-----RAW PARTITION output: hs OPEN(unit=32, file='PART_HSIGN.OUT', & status='unknown') maxpartout = 5 DO tsA = 1,SIZE(wsdat) WRITE(32,'(I4,71x,A)') tsA,'Time step' WRITE(32,'(I4,71x,A)') maxpartout,'Tot number of raw partitions' DO k = 1,maxpartout WRITE(32,'(I4,71x,A)') k,'System number' WRITE(32,'(I4,71x,A)') 9999,'Number of points in system' DO j = maxJ,1,-1 DO i = 1,maxI WRITE(32,'(F8.2)',ADVANCE='NO') wsdat(tsA)%par(i,j)%hs(k) END DO WRITE(32,'(A)',ADVANCE='YES') '' END DO END DO END DO CLOSE(32) !-----RAW PARTITION output: tp ! OPEN(unit=33,recl=2147483646, file='PART_TP.OUT', & ! status='unknown') OPEN(unit=33, file='PART_TP.OUT', & status='unknown') DO tsA = 1,SIZE(wsdat) WRITE(33,'(I4,71x,A)') tsA,'Time step' WRITE(33,'(I4,71x,A)') maxpartout,'Tot number of raw partitions' DO k = 1,maxpartout WRITE(33,'(I4,71x,A)') k,'System number' WRITE(33,'(I4,71x,A)') 9999,'Number of points in system' DO j = maxJ,1,-1 DO i = 1,maxI WRITE(33,'(F8.2)',ADVANCE='NO') wsdat(tsA)%par(i,j)%tp(k) END DO WRITE(33,'(A)',ADVANCE='YES') '' END DO END DO END DO CLOSE(33) !-----RAW PARTITION output: dir OPEN(unit=34, file='PART_DIR.OUT', & status='unknown') DO tsA = 1,SIZE(wsdat) WRITE(34,'(I4,71x,A)') tsA,'Time step' WRITE(34,'(I4,71x,A)') maxpartout,'Tot number of raw partitions' DO k = 1,maxpartout WRITE(34,'(I4,71x,A)') k,'System number' WRITE(34,'(I4,71x,A)') 9999,'Number of points in system' DO j = maxJ,1,-1 DO i = 1,maxI WRITE(34,'(F8.2)',ADVANCE='NO') wsdat(tsA)%par(i,j)%dir(k) END DO WRITE(34,'(A)',ADVANCE='YES') '' END DO END DO END DO CLOSE(34) !-----RAW PARTITION output: dspr OPEN(unit=35, file='PART_DSPR.OUT', & status='unknown') DO tsA = 1,SIZE(wsdat) WRITE(35,'(I4,71x,A)') tsA,'Time step' WRITE(35,'(I4,71x,A)') maxpartout,'Tot number of raw partitions' DO k = 1,maxpartout WRITE(35,'(I4,71x,A)') k,'System number' WRITE(35,'(I4,71x,A)') 9999,'Number of points in system' DO j = maxJ,1,-1 DO i = 1,maxI WRITE(35,'(F8.2)',ADVANCE='NO') & wsdat(tsA)%par(i,j)%dspr(k) END DO WRITE(35,'(A)',ADVANCE='YES') '' END DO END DO END DO CLOSE(35) END IF !/MPI END IF ! ------------------------------------------------------------------------ ! Allocate the sysA structure !/MPI IF (rank.EQ.0) THEN WRITE(20,*) 'Allocating sysA...' !/MPI END IF ALLOCATE( sysA(maxTs) ) !/MPI IF (rank.EQ.0) THEN WRITE(20,*) 'SIZE(sysA) = ',SIZE(sysA) WRITE(6,1020) ' Number of time levels being processed:',SIZE(sysA) 1020 FORMAT(A,I4) !/MPI END IF ! Allocate maxSys ALLOCATE( maxSys(maxTs) ) 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 ALLOCATE( sysA(1) ) !Change to sysA(2)? ! Allocate maxSys ALLOCATE( maxSys(1) ) !Change to maxSys(2)? END IF ! Big loop over all time levels !/MPI IF (rank.EQ.0) THEN WRITE(6,*) 'Performing spatial tracking...' !/MPI END IF !/MPI/! WRITE(20,*) 'rank,t0,maxTs,nproc =',rank,t0,maxTs,nproc !/MPI DO tsA = (t0+rank),maxTs,MIN(nproc,maxTS) !/MPI/! WRITE(20,*) 'Computing: Rank, tsA =',rank,tsA !/SHRD DO tsA = t0,maxTs WRITE(20,*) 'Call spiralTrackV3, tsA=',tsA,'...' CALL spiralTrackV3 ( wsdat(tsA), dirKnob, perKnob, wetPts, & hsKnob, seedLat, seedLon, & maxSys(tsA), sysA(tsA)%sys ) WRITE(20,*) '*** SIZE(sysA(1:tsA)%sys) at end of time step', & tsA,':' WRITE(20,*) SIZE(sysA(tsA)%sys) END DO !/MPI CALL MPI_Barrier(MPI_COMM_WORLD,IERR) !/MPI !/MPI/!! Define communicator for array of integers in structure "system" !/MPI/! DOMSIZE = maxI*maxJ !/MPI/! WRITE(20,*) 'Rank',rank,'DOMSIZE =',DOMSIZE !/MPI/! CALL MPI_TYPE_CONTIGUOUS(DOMSIZE,MPI_INTEGER,MPI_INT_DOMARR,IERR) !/MPI/! CALL MPI_TYPE_COMMIT(MPI_INT_DOMARR,IERR) !/MPI/! CALL MPI_TYPE_EXTENT(MPI_INT_DOMARR,EXTENT,IERR) !/MPI/! WRITE(20,*) 'Rank',rank,'has set up communicator MPI_INT_DOMARR, & !/MPI/! size =',EXTENT !/MPI !/MPI/!! Define communicator for array of reals in structure "system" !/MPI/! CALL MPI_TYPE_CONTIGUOUS(DOMSIZE,MPI_REAL,MPI_REAL_DOMARR,IERR) !/MPI/! CALL MPI_TYPE_COMMIT(MPI_REAL_DOMARR,IERR) !/MPI/! CALL MPI_TYPE_EXTENT(MPI_REAL_DOMARR,EXTENT,IERR) !/MPI/! WRITE(20,*) 'Rank',rank,'has set up communicator MPI_REAL_DOMARR, & !/MPI/! size =',EXTENT !/MPI !/MPI/! Communicate results back to rank 0 !/MPI DO tsA = t0,maxTs !/MPI irank = MOD((tsA-t0),MIN(nproc,maxTS)) !/MPI/! WRITE(20,*) 'Rank,irank=',rank,irank !/MPI IF (irank.NE.0) THEN !/MPI/! WRITE(20,*) 'Communicating for Rank,irank=',rank,irank !/MPI !/MPI/! Send maxSys(tsA) at each time level to rank 0 !/MPI tag1 = tsA !/MPI IF (rank.EQ.irank) THEN !/MPI/! Send results from current rank to rank 0 (blocking) !/MPI/! WRITE(20,*) '>> Sending: rank,tsA,tag1=',rank,tsA,tag1 !/MPI CALL MPI_SEND(maxSys(tsA),1,MPI_INTEGER,0,tag1, & !/MPI MPI_COMM_WORLD,IERR) !/MPI/! WRITE(20,*) 'Rank, IERR=',rank,IERR !/MPI END IF !/MPI IF (rank.EQ.0) THEN !/MPI/! WRITE(20,*) '<< Receiving: rank,tsA,tag1=',rank,tsA,tag1 !/MPI CALL MPI_RECV(maxSys(tsA),1,MPI_INTEGER, & !/MPI irank,tag1,MPI_COMM_WORLD,MPI_STATUS,IERR) !/MPI/! Allocate structure at this time level !/MPI ALLOCATE( sysA(tsA)%sys(maxSys(tsA)) ) !/MPI DO ic = 1,maxSys(tsA) !/MPI NULLIFY( sysA(tsA)%sys(ic)%i ) !/MPI NULLIFY( sysA(tsA)%sys(ic)%j ) !/MPI NULLIFY( sysA(tsA)%sys(ic)%lon ) !/MPI NULLIFY( sysA(tsA)%sys(ic)%lat ) !/MPI NULLIFY( sysA(tsA)%sys(ic)%hs ) !/MPI NULLIFY( sysA(tsA)%sys(ic)%tp ) !/MPI NULLIFY( sysA(tsA)%sys(ic)%dir) !/MPI NULLIFY( sysA(tsA)%sys(ic)%dspr) !/MPI ALLOCATE( sysA(tsA)%sys(ic)%i(maxI*maxJ) ) !/MPI ALLOCATE( sysA(tsA)%sys(ic)%j(maxI*maxJ) ) !/MPI ALLOCATE( sysA(tsA)%sys(ic)%lon(maxI*maxJ) ) !/MPI ALLOCATE( sysA(tsA)%sys(ic)%lat(maxI*maxJ) ) !/MPI ALLOCATE( sysA(tsA)%sys(ic)%hs(maxI*maxJ) ) !/MPI ALLOCATE( sysA(tsA)%sys(ic)%tp(maxI*maxJ) ) !/MPI ALLOCATE( sysA(tsA)%sys(ic)%dir(maxI*maxJ) ) !/MPI ALLOCATE( sysA(tsA)%sys(ic)%dspr(maxI*maxJ) ) !/MPI sysA(tsA)%sys(ic)%i(:) = 9999 !/MPI sysA(tsA)%sys(ic)%j(:) = 9999 !/MPI sysA(tsA)%sys(ic)%lon(:) = 9999. !/MPI sysA(tsA)%sys(ic)%lat(:) = 9999. !/MPI sysA(tsA)%sys(ic)%hs(:) = 9999. !/MPI sysA(tsA)%sys(ic)%tp(:) = 9999. !/MPI sysA(tsA)%sys(ic)%dir(:) = 9999. !/MPI sysA(tsA)%sys(ic)%dspr(:) = 9999. !/MPI sysA(tsA)%sys(ic)%hsMean = 9999. !/MPI sysA(tsA)%sys(ic)%tpMean = 9999. !/MPI sysA(tsA)%sys(ic)%dirMean = 9999. !/MPI sysA(tsA)%sys(ic)%sysInd = 9999 !/MPI sysA(tsA)%sys(ic)%nPoints = 9999 !/MPI sysA(tsA)%sys(ic)%grp = 9999 !/MPI END DO !/MPI END IF !/MPI !/MPI/! Send data fields at each (tsA,ic) combination !/MPI IF ((rank.EQ.0).OR.(rank.EQ.irank)) THEN !/MPI DO ic = 1, maxSys(tsA) !/MPI/! Construct a unique tag for each message !/MPI tag2 = tsA*10000 + ic*100 !/MPI DOMSIZE = maxI*maxJ !/MPI !/MPI IF (rank.EQ.irank) THEN !/MPI/! WRITE(20,*) '>> Sending: rank,irank,tag2=', & !/MPI/! rank,irank,(tag2+1) !/MPI CALL MPI_SEND(sysA(tsA)%sys(ic)%i(:),DOMSIZE, & !/MPI MPI_INTEGER,0,(tag2+1),MPI_COMM_WORLD,REQ(1),IERR) !/MPI END IF !/MPI IF (rank.EQ.0) THEN !/MPI/! WRITE(20,*) '<< Receiving: rank,irank,tag2=', & !/MPI/! rank,irank,(tag2+1) !/MPI CALL MPI_RECV(sysA(tsA)%sys(ic)%i(:),DOMSIZE, & !/MPI MPI_INTEGER,irank,(tag2+1), & !/MPI MPI_COMM_WORLD,MPI_STATUS,REQ(2),IERR) !/MPI END IF !/MPI/! CALL MPI_WAITALL(2,REQ,ISTAT,IERR) !/MPI !/MPI IF (rank.EQ.irank) THEN !/MPI/! WRITE(20,*) '>> Sending: rank,irank,tag2=', & !/MPI/! rank,irank,(tag2+2) !/MPI CALL MPI_SEND(sysA(tsA)%sys(ic)%j(:),DOMSIZE, & !/MPI MPI_INTEGER,0,(tag2+2),MPI_COMM_WORLD,REQ(1),IERR) !/MPI END IF !/MPI IF (rank.EQ.0) THEN !/MPI/! WRITE(20,*) '<< Receiving: rank,irank,tag2=', & !/MPI/! rank,irank,(tag2+2) !/MPI CALL MPI_RECV(sysA(tsA)%sys(ic)%j(:),DOMSIZE, & !/MPI MPI_INTEGER,irank,(tag2+2), & !/MPI MPI_COMM_WORLD,MPI_STATUS,REQ(2),IERR) !/MPI END IF !/MPI/! CALL MPI_WAITALL(2,REQ,ISTAT,IERR) !/MPI !/MPI IF (rank.EQ.irank) THEN !/MPI/! WRITE(20,*) '>> Sending: rank,tag2=',rank,(tag2+3) !/MPI CALL MPI_SEND(sysA(tsA)%sys(ic)%lon(:),DOMSIZE, & !/MPI MPI_REAL,0,(tag2+3),MPI_COMM_WORLD,REQ(1),IERR) !/MPI END IF !/MPI IF (rank.EQ.0) THEN !/MPI/! WRITE(20,*) '<< Receiving: rank,tag2=',rank,(tag2+3) !/MPI CALL MPI_RECV(sysA(tsA)%sys(ic)%lon(:),DOMSIZE, & !/MPI MPI_REAL,irank,(tag2+3), & !/MPI MPI_COMM_WORLD,MPI_STATUS,REQ(2),IERR) !/MPI END IF !/MPI/! CALL MPI_WAITALL(2,REQ,ISTAT,IERR) !/MPI !/MPI IF (rank.EQ.irank) THEN !/MPI/! WRITE(20,*) '>> Sending: rank,tag2=',rank,(tag2+4) !/MPI CALL MPI_SEND(sysA(tsA)%sys(ic)%lat(:),DOMSIZE, & !/MPI MPI_REAL,0,(tag2+4),MPI_COMM_WORLD,REQ(1),IERR) !/MPI END IF !/MPI IF (rank.EQ.0) THEN !/MPI/! WRITE(20,*) '<< Receiving: rank,tag2=',rank,(tag2+4) !/MPI CALL MPI_RECV(sysA(tsA)%sys(ic)%lat(:),DOMSIZE, & !/MPI MPI_REAL,irank,(tag2+4), & !/MPI MPI_COMM_WORLD,MPI_STATUS,REQ(2),IERR) !/MPI END IF !/MPI/! CALL MPI_WAITALL(2,REQ,ISTAT,IERR) !/MPI !/MPI IF (rank.EQ.irank) THEN !/MPI/! WRITE(20,*) '>> Sending: rank,tag2=',rank,(tag2+5) !/MPI CALL MPI_SEND(sysA(tsA)%sys(ic)%hs(:),DOMSIZE, & !/MPI MPI_REAL,0,(tag2+5),MPI_COMM_WORLD,REQ(1),IERR) !/MPI END IF !/MPI IF (rank.EQ.0) THEN !/MPI/! WRITE(20,*) '<< Receiving: rank,tag2=',rank,(tag2+5) !/MPI CALL MPI_RECV(sysA(tsA)%sys(ic)%hs(:),DOMSIZE, & !/MPI MPI_REAL,irank,(tag2+5), & !/MPI MPI_COMM_WORLD,MPI_STATUS,REQ(2),IERR) !/MPI END IF !/MPI/! CALL MPI_WAITALL(2,REQ,ISTAT,IERR) !/MPI !/MPI IF (rank.EQ.irank) THEN !/MPI/! WRITE(20,*) '>> Sending: rank,tag2=',rank,(tag2+6) !/MPI CALL MPI_SEND(sysA(tsA)%sys(ic)%tp(:),DOMSIZE, & !/MPI MPI_REAL,0,(tag2+6),MPI_COMM_WORLD,REQ(1),IERR) !/MPI END IF !/MPI IF (rank.EQ.0) THEN !/MPI/! WRITE(20,*) '<< Receiving: rank,tag2=',rank,(tag2+6) !/MPI CALL MPI_RECV(sysA(tsA)%sys(ic)%tp(:),DOMSIZE, & !/MPI MPI_REAL,irank,(tag2+6), & !/MPI MPI_COMM_WORLD,MPI_STATUS,REQ(2),IERR) !/MPI END IF !/MPI/! CALL MPI_WAITALL(2,REQ,ISTAT,IERR) !/MPI !/MPI IF (rank.EQ.irank) THEN !/MPI/! WRITE(20,*) '>> Sending: rank,tag2=',rank,(tag2+7) !/MPI CALL MPI_SEND(sysA(tsA)%sys(ic)%dir(:),DOMSIZE, & !/MPI MPI_REAL,0,(tag2+7),MPI_COMM_WORLD,REQ(1),IERR) !/MPI END IF !/MPI IF (rank.EQ.0) THEN !/MPI/! WRITE(20,*) '<< Receiving: rank,tag2=',rank,(tag2+7) !/MPI CALL MPI_RECV(sysA(tsA)%sys(ic)%dir(:),DOMSIZE, & !/MPI MPI_REAL,irank,(tag2+7), & !/MPI MPI_COMM_WORLD,MPI_STATUS,REQ(2),IERR) !/MPI END IF !/MPI/! CALL MPI_WAITALL(2,REQ,ISTAT,IERR) !/MPI !/MPI IF (rank.EQ.irank) THEN !/MPI/! WRITE(20,*) '>> Sending: rank,tag2=',rank,(tag2+8) !/MPI CALL MPI_SEND(sysA(tsA)%sys(ic)%dspr(:),DOMSIZE, & !/MPI MPI_REAL,0,(tag2+8),MPI_COMM_WORLD,REQ(1),IERR) !/MPI END IF !/MPI IF (rank.EQ.0) THEN !/MPI/! WRITE(20,*) '<< Receiving: rank,tag2=',rank,(tag2+8) !/MPI CALL MPI_RECV(sysA(tsA)%sys(ic)%dspr(:),DOMSIZE, & !/MPI MPI_REAL,irank,(tag2+8), & !/MPI MPI_COMM_WORLD,MPI_STATUS,REQ(2),IERR) !/MPI END IF !/MPI/! CALL MPI_WAITALL(2,REQ,ISTAT,IERR) !/MPI !/MPI IF (rank.EQ.irank) THEN !/MPI/! WRITE(20,*) '>> Sending: rank,irank,tag2=', & !/MPI/! rank,irank,(tag2+9) !/MPI CALL MPI_SEND(sysA(tsA)%sys(ic)%hsMean,1,MPI_REAL, & !/MPI 0,(tag2+9),MPI_COMM_WORLD,IERR) !/MPI END IF !/MPI IF (rank.EQ.0) THEN !/MPI/! WRITE(20,*) '<< Receiving: rank,irank,tag2=', & !/MPI/! rank,irank,(tag2+9) !/MPI CALL MPI_RECV(sysA(tsA)%sys(ic)%hsMean,1,MPI_REAL, & !/MPI irank,(tag2+9),MPI_COMM_WORLD,MPI_STATUS,IERR) !/MPI END IF !/MPI !/MPI IF (rank.EQ.irank) THEN !/MPI/! WRITE(20,*) '>> Sending: rank,irank,tag2=', & !/MPI/! rank,irank,(tag2+10) !/MPI CALL MPI_SEND(sysA(tsA)%sys(ic)%tpMean,1,MPI_REAL, & !/MPI 0,(tag2+10),MPI_COMM_WORLD,IERR) !/MPI END IF !/MPI IF (rank.EQ.0) THEN !/MPI/! WRITE(20,*) '<< Receiving: rank,irank,tag2=', & !/MPI/! rank,irank,(tag2+10) !/MPI CALL MPI_RECV(sysA(tsA)%sys(ic)%tpMean,1,MPI_REAL, & !/MPI irank,(tag2+10),MPI_COMM_WORLD,MPI_STATUS,IERR) !/MPI END IF !/MPI !/MPI IF (rank.EQ.irank) THEN !/MPI/! WRITE(20,*) '>> Sending: rank,irank,tag2=', & !/MPI/! rank,irank,(tag2+11) !/MPI CALL MPI_SEND(sysA(tsA)%sys(ic)%dirMean,1,MPI_REAL, & !/MPI 0,(tag2+11),MPI_COMM_WORLD,IERR) !/MPI END IF !/MPI IF (rank.EQ.0) THEN !/MPI/! WRITE(20,*) '<< Receiving: rank,irank,tag2=', & !/MPI/! rank,irank,(tag2+11) !/MPI CALL MPI_RECV(sysA(tsA)%sys(ic)%dirMean,1,MPI_REAL, & !/MPI irank,(tag2+11),MPI_COMM_WORLD,MPI_STATUS,IERR) !/MPI END IF !/MPI !/MPI IF (rank.EQ.irank) THEN !/MPI/! WRITE(20,*) '>> Sending: rank,irank,tag2=', & !/MPI/! rank,irank,(tag2+12) !/MPI CALL MPI_SEND(sysA(tsA)%sys(ic)%sysInd,1,MPI_INTEGER,& !/MPI 0,(tag2+12),MPI_COMM_WORLD,IERR) !/MPI END IF !/MPI IF (rank.EQ.0) THEN !/MPI/! WRITE(20,*) '<< Receiving: rank,irank,tag2=', & !/MPI/! rank,irank,(tag2+12) !/MPI CALL MPI_RECV(sysA(tsA)%sys(ic)%sysInd,1,MPI_INTEGER,& !/MPI irank,(tag2+12),MPI_COMM_WORLD,MPI_STATUS,IERR) !/MPI END IF !/MPI !/MPI IF (rank.EQ.irank) THEN !/MPI/! WRITE(20,*) '>> Sending: rank,irank,tag2=', & !/MPI/! rank,irank,(tag2+13) !/MPI CALL MPI_SEND(sysA(tsA)%sys(ic)%nPoints,1,MPI_INTEGER,& !/MPI 0,(tag2+13),MPI_COMM_WORLD,IERR) !/MPI END IF !/MPI IF (rank.EQ.0) THEN !/MPI/! WRITE(20,*) '<< Receiving: rank,irank,tag2=', & !/MPI/! rank,irank,(tag2+13) !/MPI CALL MPI_RECV(sysA(tsA)%sys(ic)%nPoints,1,MPI_INTEGER,& !/MPI irank,(tag2+13),MPI_COMM_WORLD,MPI_STATUS,IERR) !/MPI END IF !/MPI !/MPI IF (rank.EQ.irank) THEN !/MPI/! WRITE(20,*) '>> Sending: rank,irank,tag2=', & !/MPI/! rank,irank,(tag2+14) !/MPI CALL MPI_SEND(sysA(tsA)%sys(ic)%grp,1,MPI_INTEGER,& !/MPI 0,(tag2+14),MPI_COMM_WORLD,IERR) !/MPI END IF !/MPI IF (rank.EQ.0) THEN !/MPI/! WRITE(20,*) '<< Receiving: rank,irank,tag2=', & !/MPI/! rank,irank,(tag2+14) !/MPI CALL MPI_RECV(sysA(tsA)%sys(ic)%grp,1,MPI_INTEGER,& !/MPI irank,(tag2+14),MPI_COMM_WORLD,MPI_STATUS,IERR) !/MPI END IF !/MPI END DO !/MPI END IF !/MPI END IF !/MPI END DO !/MPI !/MPI CALL MPI_Barrier(MPI_COMM_WORLD,IERR) !/MPI !/MPI/! CALL MPI_TYPE_FREE(MPI_INT_DOMARR,IERR) !/MPI/! CALL MPI_TYPE_FREE(MPI_REAL_DOMARR,IERR) !/MPI !/MPI IF (rank.EQ.0) THEN WRITE(6,*) 'Performing temporal tracking...' WRITE(20,*) 'Calling timeTrackingV2...' lonext = wsdat(1)%lon(maxI,1)-wsdat(1)%lon(1,1) latext = wsdat(1)%lat(1,maxJ)-wsdat(1)%lat(1,1) CALL timeTrackingV2 (sysA, maxSys, tpTimeKnob, dirTimeKnob, 1, & maxGroup, dt, lonext, latext, maxI, maxJ) ! !/MPI END IF ! RETURN ! 802 CONTINUE WRITE (6,990) IOERR STOP 1 990 FORMAT (/' *** WAVEWATCH III ERROR IN W3STRKMD : '/ & ' ERROR IN READING FROM PARTITION FILE'/ & ' IOSTAT =',I5/) 1000 FORMAT (F9.0,F7.0,F8.3,F8.3,14X,I3,7X,F5.1,F6.1,F5.1,F6.1) 1010 FORMAT (3X,F8.2,F8.2,8X,F9.2,F9.2) 1200 FORMAT (/' *** WAVEWATCH III ERROR IN W3STRKMD : '/ & ' ERROR IN READING PARTITION FILE '/ & ' INCOMPATIBLE ENDIANESS'/ ) 1300 FORMAT (/' *** WAVEWATCH III ERROR IN W3STRKMD : '/ & ' ERROR IN READING PARTITION FILE '/ & ' EXPECTED IDSTR "WAVEWATCH III PARTITIONED DATA FILE"'/ ) 1400 FORMAT (/' *** WAVEWATCH III ERROR IN W3STRKMD : '/ & ' ERROR IN FINDING DOMAIN TO PROCESS - '/ & ' SPECIFIED LAT/LON LIMITS WITHIN DOMAIN '/ & ' OF RAW PARTITION FILE?'/ ) 2001 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/ & ' ERROR IN OPENING INPUT FILE'/ ) 2002 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/ & ' PREMATURE END OF INPUT FILE'/ ) 2003 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/ & ' PREMATURE END OF PARTITION FILE - '/ & ' TSTART=',F13.4/ ) ! ! END SUBROUTINE waveTracking_NWS_V2 !/ End of waveTracking_NWS_V2 ---------------------------------------- / !/ !/ ------------------------------------------------------------------- / SUBROUTINE spiralTrackV3 (wsdat ,dirKnob ,perKnob ,wetPts , & hsKnob ,seedLat ,seedLon , & maxSys ,sys ) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 4-Jan-2013 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! Performs the spatial spiral tracking for a given time step ! ! 2. Method ! ! Index convention on grid: ! ! j ! ^ ! |+(1,maxJ) +(maxI,maxJ) ! | ! | ! | ! | ! | ! | ! |(1,1) +(maxI,1) ! +----------------------> i ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! dirKnob Real input Parameter in direction for combining fields in space ! perKnob Real input Parameter in period for combining fields in space ! wetPts Real input Percentage of wet points for purging fields (fraction) ! hsKnob Real input Parameter in wave height for purging fields ! 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 REAL :: dirKnob,perKnob,wetPts,hsKnob,seedLat,seedLon INTEGER :: maxSys TYPE(system), POINTER :: sys(:) INTENT (IN) wetPts,dirKnob,perKnob,hsKnob,seedLat,seedLon INTENT (IN OUT) wsdat ! INTENT (OUT) maxSys,sys ! ! Local variables ! ---------------------------------------------------------------- ! ngbrExt Int How far do we want the neighbour to be considered ! 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 CHARACTER :: way *1 INTEGER :: ngbrExt, combine, maxI, maxJ, i, j, oldJ INTEGER :: horizStepCount, vertStepCount, checkCount, sc, & maxPts, landPts, horizBorder, vertBorder, m, k, & stepCount REAL :: deltaLat, minLat, maxLat, minLon, maxLon ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! findWay ! findNext ! findSys ! combineWaveSystems ! ! 5. Subroutines calling ! ! waveTracking_NWS_V2 ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! - ! ! 9. Switches : ! ! None defined yet. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / ! 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(20,*) 'In spiralTrackV3: combine = ',combine 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(20,*) 'In spiralTrackV3, i=NINT(maxI/2.) =',i WRITE(20,*) '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(i,1).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) horizBorder=0 vertBorder=0 DO WHILE (checkCount.LE.(maxI*maxJ-3) ) ! From the direction (way) we were going before, find which direction we ! are going now and how many 'step' we need to take CALL findWay(way, horizStepCount, vertStepCount, & vertBorder, horizBorder, 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 CALL findNext (i,j,maxI,maxJ,way,vertBorder,horizBorder) 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 CALL findSys(i, j, wsdat, maxSys, ngbrExt, maxI, maxJ, & perKnob, dirKnob, hsKnob) END IF END DO END DO ! wetPts% of wet points maxPts=NINT(wetPts*(maxI*maxJ-1)) ! WRITE(20,*) 'Call combineWaveSystems...' CALL combineWaveSystems(wsdat,maxSys,maxPts,maxI,maxJ, & perKnob,dirKnob,hsKnob,combine,sys) RETURN END SUBROUTINE spiralTrackV3 !/ End of spiralTrackV3 ---------------------------------------------- / !/ !/ ------------------------------------------------------------------- / SUBROUTINE timeTrackingV2 (sysA ,maxSys ,tpTimeKnob , & dirTimeKnob,ts0 ,maxGroup , & dt ,lonext ,latext , & maxI ,maxJ ) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 4-Jan-2013 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! Performs the time tracking of the systems identified within ! the subroutine spiralTrackV3. ! ! 2. Method ! ! - ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! 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 ! ts0 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 ! lonext Real input Longitudinal extent of domain ! latext Real input Latitudinal extent of domain ! maxI, maxJ Int input Maximum indices of wave field ! TYPE(timsys), POINTER :: sysA(:) INTEGER, POINTER :: maxSys(:) REAL :: dirTimeKnob, tpTimeKnob INTEGER :: ts0, maxGroup REAL :: dt REAL :: lonext, latext INTEGER :: maxI, maxJ INTENT (IN) tpTimeKnob, dirTimeKnob, ts0, maxI, maxJ ! INTENT (IN OUT) sysA INTENT (OUT) maxGroup ! ! Local variables ! ---------------------------------------------------------------- ! ic Int Counter for wave systems ! ts1 Int Adjusted initial time step in case ts0 has only empty systems ! LOGICAL :: file_exists CHARACTER :: dummy*23 TYPE(sysmemory) :: sysMem(50) !!! 50 memory spaces should be enough Check!!! INTEGER :: leng, l, i, ii, j, k, kk, idir, numSys, & counter, new, DIFSIZE, tpMinInd, dirMinInd, used, ok REAL :: Tb, deltaPer, deltaDir, tpMinVal, dirMinVal, & dirForTpMin, tpForDirMin REAL, ALLOCATABLE :: sysOrdered(:), TEMP(:), dirs(:) REAL, POINTER :: DIFARR(:) INTEGER, ALLOCATABLE :: indSorted(:), alreadyUsed(:), allInd(:) INTEGER, ALLOCATABLE :: ind(:), ind2(:) INTEGER :: ts1 REAL, ALLOCATABLE :: GOF(:,:), GOFMinVal(:), GOFMinInd(:), & Tbsysmem(:), deltaDirsysmem(:), & deltaPersysmem(:),m1sysmem(:),m2sysmem(:) REAL :: m1, m2 REAL :: lonmean, latmean, dmndiag INTEGER :: npnts, npnts2 REAL, ALLOCATABLE :: mnlonlist(:), mnlatlist(:), mndist(:) REAL, POINTER :: dummy1(:),dummy2(:),dummy3(:) INTEGER, ALLOCATABLE :: olsize(:) REAL :: TEMP1, TEMP2 INTEGER :: iii, jj, ll, idup ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! SORT ! SETDIFF ! ! 5. Subroutines calling ! ! waveTracking_NWS_V2 ! ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! - ! ! 9. Switches : ! ! None defined yet. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / ! Associate default grp value to time step 1 WRITE(20,*) 'TIME TRACKING' WRITE(20,*) 'Inside timeTrackingV2: SIZE(sysA(1)%sys) =', & SIZE(sysA(1)%sys) WRITE(20,*) 'Inside timeTrackingV2: maxSys(1) =',maxSys(1) WRITE(20,*) 'ts0 = ',ts0 ts1 = ts0 ! Skip initial time steps with empty systems (e.g. when starting from rest) DO i = ts1, SIZE(sysA) IF (SIZE(sysA(ts1)%sys).EQ.0) ts1 = ts1+1 ! No non-empty systems found IF (ts1.GT.SIZE(sysA)) THEN maxGroup = 0 GOTO 2000 END IF END DO WRITE(20,*) 'TS = ',ts1 IF (SIZE(sysA(ts1)%sys).GT.0) THEN ! Initialize system memory groups sysA(ts1)%sys(:)%grp = 9999 sysMem(:)%grp = 9999 sysMem(:)%nPoints = 0 sysMem(:)%lonMean = 9999. sysMem(:)%latMean = 9999. sysMem(:)%tpMean = 9999. sysMem(:)%dirMean = 9999. sysMem(:)%updated = -9999 sysMem(:)%length = 0 DO iii = 1,50 ALLOCATE(sysMem(iii)%indx(maxI*maxJ)) sysMem(iii)%indx = 9999 END DO INQUIRE(FILE="sys_restart.ww3", EXIST=file_exists) IF (file_exists) THEN ! Use groups from wave tracking hotfile WRITE(20,*) '*** Using group memory hotfile' OPEN(unit=12,file='sys_restart.ww3',status='old') READ(12,'(A23,I10)') dummy,maxGroup WRITE(20,*) 'Reading ',maxGroup,' systems' DO k = 1,maxGroup READ(12,'(A23,I10)') dummy,sysMem(k)%grp READ(12,'(A23,I10)') dummy,sysMem(k)%nPoints READ(12,'(A23,F10.4)') dummy,sysMem(k)%lonMean READ(12,'(A23,F10.4)') dummy,sysMem(k)%latMean READ(12,'(A23,F10.3)') dummy,sysMem(k)%tpMean READ(12,'(A23,F10.3)') dummy,sysMem(k)%dirMean READ(12,'(A23,I10)') dummy,sysMem(k)%updated READ(12,'(A23,I10)') dummy,sysMem(k)%length DO j = maxJ,1,-1 READ(12,*) (sysMem(k)%indx((j-1)*maxI+i), i = 1,maxI) END DO !Reset update counter sysMem(k)%updated = 0 END DO CLOSE(12) ts1 = ts1-1 ELSE ! 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) IF (ALLOCATED(sysOrdered)) DEALLOCATE(sysOrdered) IF (ALLOCATED(indSorted)) DEALLOCATE(indSorted) ! Set the initial long-term system memory DO i = 1, maxSys(ts1) sysA(ts1)%sys(i)%grp = i ! Set initial values of long-term system memory sysMem(i)%grp = i sysMem(i)%nPoints = sysA(ts1)%sys(i)%nPoints sysMem(i)%lonMean = & SUM(sysA(ts1)%sys(i)%lon(1:sysMem(i)%nPoints))/& sysMem(i)%nPoints sysMem(i)%latMean = & SUM(sysA(ts1)%sys(i)%lat(1:sysMem(i)%nPoints))/& sysMem(i)%nPoints !070512----------- Weight averages with Hm0 --------------------- TEMP1 = 0. TEMP2 = 0. DO iii = 1,sysMem(i)%nPoints TEMP1 = TEMP1 + & (sysA(ts1)%sys(i)%hs(iii)**2)*sysA(ts1)%sys(i)%lon(iii) TEMP2 = TEMP2 + & (sysA(ts1)%sys(i)%hs(iii)**2)*sysA(ts1)%sys(i)%lat(iii) END DO sysMem(i)%lonMean = TEMP1/& MAX(SUM(sysA(ts1)%sys(i)%hs(1:sysMem(i)%nPoints)**2),& 0.001) sysMem(i)%latMean = TEMP2/& MAX(SUM(sysA(ts1)%sys(i)%hs(1:sysMem(i)%nPoints)**2),& 0.001) !070512----------- Weight averages with Hm0 --------------------- sysMem(i)%tpMean = sysA(ts1)%sys(i)%tpMean sysMem(i)%dirMean = sysA(ts1)%sys(i)%dirMean sysMem(i)%updated = ts1 sysMem(i)%length = 1 !071012----------- Grid point indexing -------------------------- DO iii = 1,sysMem(i)%nPoints sysMem(i)%indx(iii) = (sysA(ts1)%sys(i)%j(iii)-1)*maxI +& sysA(ts1)%sys(i)%i(iii) END DO !071012----------- Grid point indexing -------------------------- END DO maxGroup = maxSys(ts1) ! i = ts1 END IF !******** Test output *********************** DO i = 1, maxGroup WRITE(20,*) 'sysMem(',i,')%grp =',sysMem(i)%grp WRITE(20,*) 'sysMem(',i,')%nPoints =',sysMem(i)%nPoints WRITE(20,*) 'sysMem(',i,')%lonMean =',sysMem(i)%lonMean WRITE(20,*) 'sysMem(',i,')%latMean =',sysMem(i)%latMean WRITE(20,*) 'sysMem(',i,')%tpMean =',sysMem(i)%tpMean WRITE(20,*) 'sysMem(',i,')%dirMean =',sysMem(i)%dirMean WRITE(20,*) 'sysMem(',i,')%updated =',sysMem(i)%updated WRITE(20,*) 'sysMem(',i,')%length =',sysMem(i)%length END DO !******************************************** END IF ! Loop over all time levels to track systems in time WRITE(20,*) 'Number of time levels = ',SIZE(sysA) DO i = (ts1+1), SIZE(sysA) WRITE(20,*) 'TS = ',i IF (SIZE(sysA(i)%sys).GT.0) THEN ! *** Added: 02/29/12 ************************************* ! Sort groups, so that larger systems get associated first ALLOCATE( sysOrdered(maxSys(i)) ) ALLOCATE( indSorted(maxSys(i)) ) CALL SORT (REAL(sysA(i)%sys(1:maxSys(i))%nPoints), & maxSys(i),sysOrdered,indSorted,'D') sysA(i)%sys(1:maxSys(i)) = sysA(i)%sys(indSorted) IF (ALLOCATED(sysOrdered)) DEALLOCATE(sysOrdered) IF (ALLOCATED(indSorted)) DEALLOCATE(indSorted) ! *** Added: 02/29/12 ************************************* ! Initialize groups ! Optimize? sysA(i)%sys(:)%grp = 9999 ! Optimize? counter = 0 leng = LENGTH(REAL(sysMem(:)%grp), & SIZE(sysMem(:)%grp),REAL(9999)) ALLOCATE( alreadyUsed(leng+10) ) !Make space for 10 new potential entries. Improve!!! WRITE(20,*) 'sysMem(1:leng)%grp =', & sysMem(1:leng)%grp ALLOCATE( allInd(leng) ) alreadyUsed(:) = 0 allInd(:) = sysMem(1:leng)%grp !071212-----GoF 2D------------------------------- ALLOCATE( ind(SIZE(allInd)) ) ind(:) = allInd ALLOCATE( ind2(SIZE(ind)) ) DO ii = 1, SIZE(ind) ind2(ii) = FINDFIRST(REAL(allInd),SIZE(allInd), & REAL(ind(ii))) END DO ! Define 2D array for evaluating best fit for systems ALLOCATE( GOF(maxSys(i),maxGroup) ) ALLOCATE( GOFMinVal(maxGroup) ) ALLOCATE( GOFMinInd(maxGroup) ) ALLOCATE( Tbsysmem(maxGroup) ) ALLOCATE( deltaDirsysmem(maxGroup) ) ALLOCATE( deltaPersysmem(maxGroup) ) ALLOCATE( m1sysmem(maxGroup) ) ALLOCATE( m2sysmem(maxGroup) ) !071212-----GoF 2D------------------------------- DO j = 1, maxSys(i) npnts = sysA(i)%sys(j)%nPoints lonmean = SUM(sysA(i)%sys(j)%lon(1:npnts))/npnts latmean = SUM(sysA(i)%sys(j)%lat(1:npnts))/npnts !070512----------- Weight averages with Hm0 --------------------- TEMP1 = 0. TEMP2 = 0. DO iii = 1,npnts TEMP1 = TEMP1 + & (sysA(i)%sys(j)%hs(iii)**2)*sysA(i)%sys(j)%lon(iii) TEMP2 = TEMP2 + & (sysA(i)%sys(j)%hs(iii)**2)*sysA(i)%sys(j)%lat(iii) END DO lonmean=TEMP1/MAX(SUM(sysA(i)%sys(j)%hs(1:npnts)**2),0.001) latmean=TEMP2/MAX(SUM(sysA(i)%sys(j)%hs(1:npnts)**2),0.001) !070512----------- Weight averages with Hm0 --------------------- !071012----------- Grid point indexing -------------------------- ALLOCATE(sysA(i)%sys(j)%indx(maxI*maxJ)) sysA(i)%sys(j)%indx = 9999 DO iii = 1,sysA(i)%sys(j)%nPoints sysA(i)%sys(j)%indx(iii) = & (sysA(i)%sys(j)%j(iii)-1)*maxI + & sysA(i)%sys(j)%i(iii) END DO !071012----------- Grid point indexing -------------------------- WRITE(20,*) 'System no. ',j,' of ',maxSys(i) WRITE(20,*) 'Size =', npnts WRITE(20,*) 'lonMean =', lonmean WRITE(20,*) 'latMean =', latmean WRITE(20,*) 'tpMean =', sysA(i)%sys(j)%tpMean WRITE(20,*) 'dirMean =', sysA(i)%sys(j)%dirMean sysA(i)%sys(j)%grp = 9999 !Now redundant? ! Compute deltas Tbsysmem = sysMem(1:maxGroup)%tpMean WRITE(20,*) 'Tbsysmem(:) = ', Tbsysmem(:) ! Compute deltas the same way as for field combining - they should ! be of the same degree of strictness as the latter, otherwise ! the time combining will lose track! !3stddev m1 = -3.645*Tb + 63.211 !3stddev m1 = MAX(m1,10.) !3stddev m2 = -0.346*Tb + 3.686 !3stddev m2 = MAX(m2,0.6) !1stddev m1 = -2.219*Tb + 35.734 !1stddev m1 = MAX(m1,5.) !1stddev m2 = -0.226*Tb + 2.213 !1stddev m2 = MAX(m2,0.35) !071412 m1 = -5.071*Tb + 90.688 !071412 m1 = MAX(m1,16.) !071412 m2 = -0.467*Tb + 5.161 !071412 m2 = MAX(m2,1.0) !071412 deltaDir = (m1*1. + dirTimeKnob)*1. !071412 deltaPer = (m2*1. + tpTimeKnob)*1. DO ii = 1,SIZE(ind2) m1sysmem(ii) = MAX((-3.645*Tbsysmem(ii)+63.211),10.) m2sysmem(ii) = MAX((-0.346*Tbsysmem(ii)+3.686),0.6) END DO deltaDirsysmem = m1sysmem(:)*1. + dirTimeKnob deltaPersysmem = m2sysmem(:)*1. + tpTimeKnob WRITE(20,*) 'deltaDirsysmem(:) = ',deltaDirsysmem WRITE(20,*) 'deltaPersysmem(:) = ',deltaPersysmem ! Criterion 1: Mean period ALLOCATE( TEMP(SIZE(ind2)) ) TEMP = ABS( sysA(i)%sys(j)%tpMean - & sysMem(ind2(:))%tpMean ) WRITE(20,*) 'tpMean list =', & sysMem(ind2(:))%tpMean WRITE(20,*) 'tpMinVal list =', TEMP tpMinVal = MINVAL(TEMP) tpMinInd = FINDFIRST(TEMP,SIZE(TEMP),tpMinVal) ! Criterion 2: Mean direction ALLOCATE( dirs(SIZE(ind2)) ) dirs(:)=ABS( sysA(i)%sys(j)%dirMean - & sysMem(ind2(:))%dirMean ) ! Deal with wrap around DO idir = 1, SIZE(dirs) IF (dirs(idir).GE.180.) dirs(idir)=360-dirs(idir) END DO WRITE(20,*) 'dirMean list =', & sysMem(ind2(:))%dirMean WRITE(20,*) 'dirMinVal list =', dirs ! Criterion 3: Size WRITE(20,*) 'Size list =', & sysMem(ind2(:))%nPoints ! Criterion 4: Distance between systems ALLOCATE (mnlonlist(SIZE(ind2))) ALLOCATE (mnlatlist(SIZE(ind2))) ALLOCATE (mndist(SIZE(ind2))) DO ii = 1,SIZE(ind2) mnlonlist(ii) = sysMem(ind2(ii))%lonMean mnlatlist(ii) = sysMem(ind2(ii))%latMean mndist(ii) = SQRT((lonmean-mnlonlist(ii))**2 + & (latmean-mnlatlist(ii))**2) END DO dmndiag = SQRT(lonext**2+latext**2) WRITE(20,*) 'Distance list =',mndist(:) WRITE(20,*) 'Domain diagonal =',dmndiag ! Criterion 5: Overlap of systems ALLOCATE (olsize(SIZE(ind2))) DO ii = 1,SIZE(ind2) IF (sysMem(ind2(ii))%nPoints.GT.0) THEN CALL INTERSECT(REAL(sysA(i)%sys(j)%indx(1:npnts)),npnts, & REAL(sysMem(ind2(ii))%indx(1:sysMem(ind2(ii))%nPoints)),& sysMem(ind2(ii))%nPoints,dummy1,olsize(ii),dummy2,dummy3) ELSE olsize(ii) = 0 END IF END DO GOF(j,1:SIZE(ind2)) = (TEMP/deltaPersysmem(:))**2 + & (dirs/deltaDirsysmem(:))**2 + & ! (4*mndist(:)/dmndiag)**2 ( (REAL(olsize(:)) - & REAL(sysMem(ind2(:))%nPoints) )/& (0.50*MAX(REAL(sysMem(ind2(:))%nPoints),0.001)) )**2 ! Remove GoF entries which exceed predifined tolerances DO ii = 1,SIZE(ind2) WRITE(20,*) 'Testing: ii,olsize(ii),size,frac =',& ii,olsize(ii),sysMem(ind2(ii))%nPoints,& REAL(olsize(ii))/& MAX(REAL(sysMem(ind2(ii))%nPoints),0.001) IF ( REAL(olsize(ii)).LT.& 0.50*REAL(sysMem(ind2(ii))%nPoints) ) THEN GOF(j,ii) = 9999. END IF IF ( (TEMP(ii).GT.deltaPersysmem(ii)).OR.& (dirs(ii).GT.deltaDirsysmem(ii)) ) THEN GOF(j,ii) = 9999. END IF END DO WRITE(20,*) 'GOF(j,:) =',GOF(j,:) IF (ALLOCATED(TEMP)) DEALLOCATE(TEMP) IF (ALLOCATED(dirs)) DEALLOCATE(dirs) IF (ALLOCATED(mnlonlist)) DEALLOCATE(mnlonlist) IF (ALLOCATED(mnlatlist)) DEALLOCATE(mnlatlist) IF (ALLOCATED(mndist)) DEALLOCATE(mndist) IF (ALLOCATED(olsize)) DEALLOCATE(olsize) !071212-----------GoF 2D------------- END DO IF (ALLOCATED(Tbsysmem)) DEALLOCATE(Tbsysmem) IF (ALLOCATED(deltaDirsysmem)) DEALLOCATE(deltaDirsysmem) IF (ALLOCATED(deltaPersysmem)) DEALLOCATE(deltaPersysmem) IF (ALLOCATED(m1sysmem)) DEALLOCATE(m1sysmem) IF (ALLOCATED(m2sysmem)) DEALLOCATE(m2sysmem) WRITE(20,*) 'GoF3:' DO jj = 1,maxSys(i) WRITE(20,*) GOF(jj,:) END DO ! Find minima in GoF DO k = 1,maxGroup GOFMinVal(k) = MINVAL(GOF(:,k)) GOFMinInd(k) = FINDFIRST(GOF(:,k),SIZE(GOF,1),GOFMinVal(k)) IF (GOFMinVal(k).EQ.9999) THEN GOFMinInd(k) = 0 END IF END DO IF (ALLOCATED(GOF)) DEALLOCATE(GOF) DO j = 1, maxSys(i) new = 0 ! Look up sysMem match for this current system. If no match ! is found, the index value 0 is returned. tpMinInd = 0 TEMP1 = 9999. DO jj = 1, SIZE(GOFMinInd) IF (GOFMinInd(jj).EQ.j) THEN IF (GOFMinVal(jj).LT.TEMP1) THEN tpMinInd = jj TEMP1 = GOFMinVal(jj) END IF END IF END DO dirMinInd = tpMinInd WRITE(20,*) 'System, GOFMinInd: ',j,tpMinInd IF (tpMinInd.NE.0) THEN ! Success !071212-----------GoF 2D------------- counter = counter+1 sysA(i)%sys(j)%grp = & sysMem(ind2(dirMinInd))%grp alreadyUsed(counter) = sysA(i)%sys(j)%grp WRITE(20,*) '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 grp ', & sysMem(ind2(dirMinInd))%grp WRITE(20,*) 'Added ',alreadyUsed(counter), & ' in array *alreadyUsed*' ELSE new = 1 END IF IF (new.EQ.1) THEN used = 0 DO k = 1, maxGroup ok = 1 WRITE(20,*) 'maxGroup,k,ok,used =', & maxGroup,k,ok,used ! Make sure it hasn't been used yet (at current time level) IF ((i.GT.2).AND. & (.NOT.ANY(alreadyUsed(:).EQ.k))) THEN ! Make sure it hasn't been used yet (at previous time level) DO l = 1, maxGroup ! If last update of system was more that *6* time steps ! ago, system can be released (TO CALIBRATE) IF ( (sysMem(l)%grp.EQ.k).AND. & ((i-sysMem(l)%updated).LT.6) ) ok = 0 WRITE(20,*) 'l, ok = ',l,ok END DO IF (ok.EQ.1) THEN sysA(i)%sys(j)%grp = k counter = counter+1; alreadyUsed(counter) = k used = 1 WRITE(20,*) 'k,used,counter =', & k,used,counter EXIT END IF END IF END DO IF (used.EQ.0) THEN maxGroup = maxGroup+1 sysA(i)%sys(j)%grp = maxGroup ! Increase sysMem by one slot sysMem(maxGroup)%grp = maxGroup counter = counter+1 alreadyUsed(counter) = maxGroup END IF WRITE(20,*) 'counter,maxGroup,sysA(i)%sys(j)%grp =',& counter,maxGroup,sysA(i)%sys(j)%grp WRITE(20,*) 'NO GRP MATCH case 2' END IF END DO IF (ALLOCATED(ind)) DEALLOCATE(ind) !071212 Shifted IF (ALLOCATED(ind2)) DEALLOCATE(ind2) !071212 Shifted IF (ALLOCATED(GOFMinVal)) DEALLOCATE(GOFMinVal) IF (ALLOCATED(GOFMinInd)) DEALLOCATE(GOFMinInd) IF (ALLOCATED(alreadyUsed)) DEALLOCATE(alreadyUsed) IF (ALLOCATED(allInd)) DEALLOCATE(allInd) ! Update sysMem DO k = 1, maxGroup DO kk = 1, maxSys(i) IF (sysA(i)%sys(kk)%grp.EQ.sysMem(k)%grp) THEN sysMem(k)%nPoints = sysA(i)%sys(kk)%nPoints sysMem(k)%lonMean = & SUM(sysA(i)%sys(kk)%lon(1:sysMem(k)%nPoints))/& sysMem(k)%nPoints sysMem(k)%latMean = & SUM(sysA(i)%sys(kk)%lat(1:sysMem(k)%nPoints))/& sysMem(k)%nPoints !070512----------- Weight averages with Hm0 --------------------- TEMP1 = 0. TEMP2 = 0. DO iii = 1,sysMem(k)%nPoints TEMP1 = TEMP1 + & (sysA(i)%sys(kk)%hs(iii)**2)*sysA(i)%sys(kk)%lon(iii) TEMP2 = TEMP2 + & (sysA(i)%sys(kk)%hs(iii)**2)*sysA(i)%sys(kk)%lat(iii) END DO sysMem(k)%lonMean = TEMP1/& MAX(SUM(sysA(i)%sys(kk)%hs(1:sysMem(k)%nPoints)**2),& 0.001) sysMem(k)%latMean = TEMP2/& MAX(SUM(sysA(i)%sys(kk)%hs(1:sysMem(k)%nPoints)**2),& 0.001) !070512----------- Weight averages with Hm0 --------------------- sysMem(k)%tpMean = sysA(i)%sys(kk)%tpMean sysMem(k)%dirMean = sysA(i)%sys(kk)%dirMean !071012----------- Grid point indexing -------------------------- sysMem(k)%indx(:) = 9999 DO iii = 1,sysMem(k)%nPoints sysMem(k)%indx(iii) = & (sysA(i)%sys(kk)%j(iii)-1)*maxI + & sysA(i)%sys(kk)%i(iii) END DO !071012----------- Grid point indexing -------------------------- sysMem(k)%updated = i sysMem(k)%length = sysMem(k)%length + 1 END IF END DO !Test for expired groups IF ((i-sysMem(k)%updated).GE.6) THEN sysMem(k)%nPoints = 0 sysMem(k)%lonMean = 9999. sysMem(k)%latMean = 9999. sysMem(k)%tpMean = 9999. sysMem(k)%dirMean = 9999. sysMem(k)%indx(:) = 9999 sysMem(k)%updated = -9999 sysMem(k)%length = 0 END IF END DO !083012 !Filter out duplicates groups that can develop DO l = 1, maxGroup DO ll = (l+1), maxGroup deltaDir = MAX((-3.645*sysMem(l)%tpMean+63.211),10.)*1. deltaPer = MAX((-0.346*sysMem(l)%tpMean+3.686),0.6)*1. IF ( (ABS(sysMem(l)%tpMean-sysMem(ll)%tpMean).LT.& deltaPer).AND. & (ABS(sysMem(l)%dirMean-sysMem(ll)%dirMean).LT.& deltaDir).AND. & (sysMem(l)%updated.NE.sysMem(ll)%updated).AND. & (sysMem(ll)%nPoints.NE.0) ) THEN !Find the more recent entry, and delete from group IF (sysMem(ll)%length.LT.sysMem(l)%length) THEN idup = ll WRITE(20,*) 'Deleting memgroup ',ll, & '(updated',sysMem(ll)%updated,', length', & sysMem(ll)%length,'), duplicate of memgroup', & l,'(updated',sysMem(l)%updated,', length', & sysMem(l)%length,'):' ELSE idup = l WRITE(20,*) 'Deleting memgroup ',l, & '(updated',sysMem(l)%updated,', length', & sysMem(l)%length,'), duplicate of memgroup', & ll,'(updated',sysMem(ll)%updated,', length', & sysMem(ll)%length,'):' END IF WRITE(20,*) 'deltaPer, diff Per:',deltaPer,& ABS(sysMem(l)%tpMean-sysMem(ll)%tpMean) WRITE(20,*) 'deltaDir, diff Dir:',deltaDir,& ABS(sysMem(l)%dirMean-sysMem(ll)%dirMean) sysMem(idup)%nPoints = 0 sysMem(idup)%lonMean = 9999. sysMem(idup)%latMean = 9999. sysMem(idup)%tpMean = 9999. sysMem(idup)%dirMean = 9999. sysMem(idup)%indx(:) = 9999 sysMem(idup)%updated = -9999 sysMem(idup)%length = 0 END IF END DO END DO ELSE WRITE(20,*) '*** No systems at this time level. ', & 'No. systems =',SIZE(sysA(i)%sys) !Test for expired groups DO k = 1, maxGroup IF ((i-sysMem(k)%updated).GE.6) THEN sysMem(k)%nPoints = 0 sysMem(k)%lonMean = 9999. sysMem(k)%latMean = 9999. sysMem(k)%tpMean = 9999. sysMem(k)%dirMean = 9999. sysMem(k)%indx(:) = 9999 sysMem(k)%updated = -9999 sysMem(k)%length = 0 END IF END DO END IF ! ******** Test output *********************** DO k = 1, maxGroup WRITE(20,*) 'sysMem(',k,')%grp =',sysMem(k)%grp WRITE(20,*) 'sysMem(',k,')%nPoints =',sysMem(k)%nPoints WRITE(20,*) 'sysMem(',k,')%lonMean =',sysMem(k)%lonMean WRITE(20,*) 'sysMem(',k,')%latMean =',sysMem(k)%latMean WRITE(20,*) 'sysMem(',k,')%tpMean =',sysMem(k)%tpMean WRITE(20,*) 'sysMem(',k,')%dirMean =',sysMem(k)%dirMean WRITE(20,*) 'sysMem(',k,')%updated =',sysMem(k)%updated WRITE(20,*) 'sysMem(',k,')%length =',sysMem(k)%length END DO ! ******************************************** END DO ! Write hotfile of wave groups OPEN(unit=27,file='sys_restart1.ww3',status='unknown') WRITE(27,'(A23,I10)') 'maxGroup =',maxGroup DO k = 1, maxGroup WRITE(27,'(A8,I3,A12,I10)') 'sysMem( ',k, & ' )%grp =',sysMem(k)%grp WRITE(27,'(A8,I3,A12,I10)') 'sysMem( ',k, & ' )%nPoints =',sysMem(k)%nPoints WRITE(27,'(A8,I3,A12,F10.4)') 'sysMem( ',k, & ' )%lonMean =',sysMem(k)%lonMean WRITE(27,'(A8,I3,A12,F10.4)') 'sysMem( ',k, & ' )%latMean =',sysMem(k)%latMean WRITE(27,'(A8,I3,A12,F10.3)') 'sysMem( ',k, & ' )%tpMean =',sysMem(k)%tpMean WRITE(27,'(A8,I3,A12,F10.3)') 'sysMem( ',k, & ' )%dirMean =',sysMem(k)%dirMean WRITE(27,'(A8,I3,A12,I10)') 'sysMem( ',k, & ' )%updated =',sysMem(k)%updated WRITE(27,'(A8,I3,A12,I10)') 'sysMem( ',k, & ' )%length =',sysMem(k)%length DO j = maxJ,1,-1 DO i = 1,maxI WRITE(27,'(I8)',ADVANCE='NO') sysMem(k)%indx((j-1)*maxI+i) END DO WRITE(27,'(A)',ADVANCE='YES') '' END DO END DO CLOSE(27) 2000 CONTINUE RETURN END SUBROUTINE timeTrackingV2 !/ End of timeTrackingV2 --------------------------------------------- / !/ !/ ------------------------------------------------------------------- / SUBROUTINE findWay (way ,horizStepCount,vertStepCount , & vertBorder ,horizBorder ,stepCount ) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 4-Jan-2013 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! From the direction (way) we were going before, find which direction we ! are going now and how many 'steps' we need to take ! ! 2. Method ! ! - ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! 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 INTEGER :: horizStepCount, vertStepCount, & vertBorder, horizBorder, stepCount INTENT (IN) vertBorder, horizBorder INTENT (OUT) stepCount INTENT (IN OUT) way ! ! Local variables ! ---------------------------------------------------------------- ! - ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! - ! ! 5. Subroutines calling ! ! spiralTrackV3 ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! See above ! ! 9. Switches : ! ! None defined yet. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / 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(20,*) 'In spaTack:findWay should NOT go here!' END SELECT RETURN END SUBROUTINE findWay !/ End of findWay ---------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / SUBROUTINE findNext (i ,j ,maxI ,maxJ , & way ,vertBorder ,horizBorder ) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | H. L. Tolman | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 4-Jan-2013 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! Find next point on spatial search spiral ! ! 2. Method ! ! - ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! 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 INTEGER :: i, j, maxI, maxJ, vertBorder, horizBorder INTENT (IN) maxI, maxJ, way INTENT (IN OUT) i, j INTENT (OUT) vertBorder, horizBorder ! ! Local variables ! ---------------------------------------------------------------- ! - ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! - ! ! 5. Subroutines calling ! ! spiralTrackV3 ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! - ! ! 9. Switches : ! ! None defined yet. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / 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 RETURN END SUBROUTINE findNext !/ End of findNext --------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / SUBROUTINE findSys (i ,j ,wsdat ,maxSys , & ngbrExt ,maxI ,maxJ ,perKnob , & dirKnob ,hsKnob ) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 4-Jan-2013 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! Find all wave systems that neighbour the grid point (i,j), and ! match these with the systems at (i,j). ! ! 2. Method ! ! For the given point (i,j), find all wave systems at neighbouring grid ! points within the reach specified by ngbrExt. ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! 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 INTEGER :: i, j, maxI, maxJ, ngbrExt, maxSys REAL :: perKnob ,dirKnob, hsKnob INTENT (IN) i, j, maxI, maxJ, ngbrExt, perKnob ,dirKnob INTENT (IN OUT) wsdat, maxSys ! ! Local variables ! ---------------------------------------------------------------- ! tmpsys TYPE(system) Temporary instance of the wave system variable ! nngbr Int Number of neighbours found ! TYPE(system), ALLOCATABLE :: tmpsys(:) 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, leng INTEGER :: allFullSys(50) REAL, POINTER :: realarr(:) INTEGER, ALLOCATABLE :: allSys(:) REAL :: hsAll(50),tpAll(50),dirAll(50),GOF(50) REAL :: absDir,absPer,absHs,T,& deltaPer,deltaDir,deltaHs,temp REAL :: dx, m1, m2 REAL :: GOFMinVal INTEGER :: GOFMinInd ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! UNIQUE ! combinePartitionsV2 ! ! 5. Subroutines calling ! ! spiralTrackV3 ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! - ! ! 9. Switches : ! ! None defined yet. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / NULLIFY(realarr) ! WRITE(20,*) '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) !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) leng = LENGTH(REAL(ngbr(l)%par%sys), & SIZE(ngbr(l)%par%sys),REAL(9999)) endCount = startCount+leng-1 allFullSys(startCount:endCount) = ngbr(l)%par%sys(1:leng) startCount=endCount+1 l=l+1 END DO IF (endCount.EQ.0) WRITE(20,*) '***1.Calling UNIQUE w. len=0!' CALL UNIQUE (REAL(allFullSys),endCount,realarr,nout) !Can one do this? ALLOCATE(allSys(nout)) allSys = INT(realarr) !Can one do this? IF (ASSOCIATED(realarr)) DEALLOCATE(realarr) maxS = MAXVAL(allSys) IF (maxSys.LT.maxS) THEN 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. ! wfAll(:) = 0. n=1 countAll=0 DO WHILE (n.LE.nngbr) ! Calculate mean of common neighbor wave system ! for every neigbor wave system found = .FALSE. DO ind = 1, SIZE(ngbr(n)%par%sys) !Optimize this? IF ( ngbr(n)%par%sys(ind).EQ.allSys(s) ) THEN !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) ! wfAll(countAll)=ngbr(n)%par%wf(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) ! tmpsys(s)%wfMean = SUM(wfAll(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).LT.hsKnob).OR. & (wsdat%par(i,j)%tp(p).EQ.0.) ) THEN wsdat%par(i,j)%sys(p)=-1 CYCLE END IF ind=0 !Replaced 'index' by 'ind' match%sysVal(:) = 9999 match%tpVal(:) = 9999. match%dirVal(:) = 9999. ! match%wfVal(:) = 9999. ! Cycle through the neighbouring systems identified above DO s=1,SIZE(allSys) absHs = ABS(wsdat%par(i,j)%hs(p)-tmpsys(s)%hsMean) absPer = ABS(wsdat%par(i,j)%tp(p)-tmpsys(s)%tpMean) absDir = ABS(wsdat%par(i,j)%dir(p)-tmpsys(s)%dirMean) ! absWf = ABS(wsdat%par(i,j)%wf(p)-tmpsys(s)%wfMean) IF (absDir.GT.180) THEN absDir = 360 - absDir IF (absDir.LT.0) THEN WRITE(20,*) '*** WARNING: absDir negative!' WRITE(20,*) 'wsdat%par(i,j)%dir(p) =', & wsdat%par(i,j)%dir(p) WRITE(20,*) 'tmpsys(s)%dirMean) =', & tmpsys(s)%dirMean END IF END IF ! Calculate delta dir and freq as a function of the partition ! dir and freq T = tmpsys(s)%tpMean dx = 0.5*( (wsdat%lon(2,1)-wsdat%lon(1,1)) + & (wsdat%lat(1,2)-wsdat%lat(1,1)) ) m1 = -3.645*T + 63.211 m1 = MAX(m1,10.) m2 = -0.346*T + 3.686 m2 = MAX(m2,0.6) !1stddev m1 = -2.219*T + 35.734 !1stddev m1 = MAX(m1,5.) !1stddev m2 = -0.226*T + 2.213 !1stddev m2 = MAX(m2,0.35) !5stddev m1 = -5.071*T + 90.688 !5stddev m1 = MAX(m1,16.) !5stddev m2 = -0.467*T + 5.161 !5stddev m2 = MAX(m2,1.0) deltaDir = m1*dx + dirKnob deltaPer = m2*dx + perKnob deltaHs = 0.25*tmpsys(s)%hsMean IF ((absPer.LT.deltaPer).AND.(absDir.LT.deltaDir)) THEN ind=ind+1 match%sysVal(ind) = allSys(s) match%tpVal(ind) = absPer match%dirVal(ind) = absDir match%hsVal(ind) = absHs ! match%wfVal(ind) = absWf END IF END DO IF (ind.GT.0) THEN IF (ind.EQ.1) THEN wsdat%par(i,j)%sys(p) = match%sysVal(1) ELSE ! Take the closest match, using GoF function GOF(:) = 9999. GOF(1:ind) = (match%tpVal(1:ind)/deltaPer)**2 + & (match%dirVal(1:ind)/deltaDir)**2 + & (match%hsVal(1:ind)/deltaHs)**2 GOFMinVal = MINVAL(GOF(1:ind)) GOFMinInd = FINDFIRST(GOF(1:ind),ind,GOFMinVal) wsdat%par(i,j)%sys(p) = match%sysVal(GOFMinInd) !The index of the system is swapped - the remaining info stays the same! END IF END IF 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)) DO p = 1, (npart-1) !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 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 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(tmpsys)) DEALLOCATE(tmpsys) RETURN END SUBROUTINE findSys !/ End of findSys ---------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / SUBROUTINE combineWaveSystems (wsdat ,maxSys ,maxPts , & maxI ,maxJ ,perKnob , & dirKnob ,hsKnob ,combine , & sys ) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 4-Jan-2013 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! Combine wave systems. Then remove small and low-energy systems from set, ! based on the parameters maxPts and maxHgt. ! ! 2. Method ! ! - ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! 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 TYPE(system), POINTER :: sys(:), systemp(:) INTEGER :: maxSys, maxPts, maxI, maxJ, combine REAL :: perKnob ,dirKnob, hsKnob INTENT (IN) maxPts, maxI, maxJ, hsKnob, combine INTENT (IN OUT) wsdat, maxSys !In the Matlab code maxSys is only input ??? ! INTENT (OUT) sys ! ! Local variables ! ---------------------------------------------------------------- ! nSys Int Number of wave systems (for checking iterative combining loop) ! LOGICAL :: found INTEGER, ALLOCATABLE :: sysOut(:) INTEGER, ALLOCATABLE :: actSysInd(:) INTEGER :: iter, ok, nSys, mS, s, so, ss, ind, leng, & iw, jw, iloop INTEGER :: actSys REAL :: dev, hsCmp, maxHgt, temp(5) ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! printFinalSys ! combineSys ! ! 5. Subroutines calling ! ! spiralTrackV3 ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! - ! ! 9. Switches : ! ! None defined yet. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !012912 WRITE(20,*) 'maxSys,maxPts,maxI,maxJ,hsKnob,combine =', & !012912 maxSys,maxPts,maxI,maxJ,hsKnob,combine ! Set up initial index array of active systems IF (.NOT.ALLOCATED(actSysInd)) ALLOCATE( actSysInd(maxSys) ) actSysInd(1:maxSys) = (/ (ind, ind = 1, maxSys) /) !opt WRITE(20,*) 'actSysInd =',actSysInd IF (combine.EQ.1) THEN ! Combine wave systems WRITE(20,*) 'Calling printFinalSys...' CALL printFinalSys (wsdat,maxSys,actSysInd,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 IF (ALLOCATED(actSysInd)) THEN nSys = SIZE(actSysInd) ELSE nSys = maxSys END IF WRITE(20,'(A,A,I3,A,I5,A)') 'Calling combineSys for ', & 'iteration',iter,' (maxSys =',nSys,').' !opt WRITE(20,*) 'SIZE(sys)=',SIZE(sys) CALL combineSys (wsdat,sys,maxSys,maxI,maxJ, & actSysInd,perKnob,dirKnob) ! No of systems after combining !opt WRITE(20,*) 'maxSys,nSys,SIZE(actSysInd) =', & !opt maxSys,nSys,SIZE(actSysInd) ! IF (maxSys.EQ.nSys) ok = 1 IF (SIZE(actSysInd).EQ.nSys) ok = 1 END DO ELSE ! Do not combine wave systems CALL printFinalSys (wsdat,maxSys,actSysInd,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(1:maxSys)%sysInd ! mS = maxSys mS = SIZE(actSysInd) ss = 1 WRITE(20,*) 'Filtering the set of',mS,'systems on size and mag.' DO so = 1, mS s = actSysInd(so) !opt NOTE: if we deallocate the individual records without !opt compressing sys, then s and sysInd will remain the same ss = s leng = LENGTH(sys(ss)%hs,SIZE(sys(ss)%hs),9999.) dev = STD(sys(ss)%hs(1:leng),leng) hsCmp = sys(ss)%hsMean + 2.*dev maxHgt = hsKnob IF ( (hsCmp.LT.maxHgt).OR.(sys(ss)%nPoints.LT.maxPts) ) THEN ! Remove system, and shift up indices to fill the gap DO ind = 1, maxSys ! Find index to remove IF (ind.EQ.ss) THEN ! Shift up entries, deleting the duplicate partition ! REPLACE WITH CSHIFT(ARRAY, SHIFT, dim)? ! IF (ind.LT.maxSys) & ! sys( ind:(maxSys-1) ) = sys( (ind+1):maxSys ) IF (ind.LE.maxSys) THEN ! Since we use pointers, we have to copy each index and ! field individually. Otherwise memory corruption occurs. DO iloop = ind,ind sys(iloop)%sysInd = 9999 sys(iloop)%nPoints = 0 sys(iloop)%grp = 9999 DEALLOCATE( sys(iloop)%hs ) DEALLOCATE( sys(iloop)%tp ) DEALLOCATE( sys(iloop)%dir ) DEALLOCATE( sys(iloop)%dspr ) ! DEALLOCATE( sys(iloop)%wf ) DEALLOCATE( sys(iloop)%i ) DEALLOCATE( sys(iloop)%j ) DEALLOCATE( sys(iloop)%lat ) DEALLOCATE( sys(iloop)%lon ) ! DEALLOCATE( sys(iloop)%hsMean ) ! DEALLOCATE( sys(iloop)%tpMean ) ! DEALLOCATE( sys(iloop)%dirMean ) ! DEALLOCATE( sys(iloop)%ngbr ) END DO END IF END IF END DO ! Update wsdat as well DO iw = 1, maxI DO jw = 1, maxJ leng = LENGTH(REAL(wsdat%par(iw,jw)%sys), & SIZE(wsdat%par(iw,jw)%sys),REAL(9999)) ind = 1 found = .FALSE. ! Identify system index (there are no duplicate ! systems at this point. DO WHILE (ind.LE.leng) IF ( wsdat%par(iw,jw)%sys(ind).EQ.s ) THEN found = .TRUE. EXIT END IF ind = ind + 1 END DO IF (found) THEN ! Blank out used record wsdat%par(iw,jw)%sys(ind) = 9999 wsdat%par(iw,jw)%ipart(ind) = 9999 END IF END DO END DO END IF END DO ! Compile array index of active systems in sys actSys = 0 DO so = 1,maxSys IF (sys(so)%nPoints>0) actSys = actSys + 1 END DO IF (ALLOCATED(actSysInd)) DEALLOCATE(actSysInd) ALLOCATE( actSysInd(actSys) ) actSys = 0 DO so = 1,maxSys IF (sys(so)%nPoints>0) THEN actSys = actSys + 1 actSysInd(actSys) = sys(so)%sysInd END IF END DO !opt WRITE(20,*) 'actSysInd =',actSysInd DO so = 1,SIZE(actSysInd) s = actSysInd(so) !opt WRITE(20,*) 'sys(',s,')%sysInd =',sys(s)%sysInd END DO CALL printFinalSys (wsdat,maxSys,actSysInd,maxI,maxJ,1,sys) !opt WRITE(20,*) 'actSysInd =',actSysInd !opt DO so = 1,maxSys !opt WRITE(20,*) 'sys(',so,')%sysInd =',sys(so)%sysInd, & !opt ', sys(',so,')%nPoints =',sys(so)%nPoints !opt END DO RETURN END SUBROUTINE combineWaveSystems !/ End of combineWaveSystems ----------------------------------------- / !/ !/ ------------------------------------------------------------------- / SUBROUTINE printFinalSys (wsdat ,maxSys ,actSysInd , & maxI ,maxJ ,flag ,sys ) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 4-Jan-2013 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! Output (print) the final output systems for this time step. ! ! 2. Method ! ! - ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! wsdat Type(dat2d) input Combined data structure ! maxI, maxJ Int input Maximum indices of wave field ! maxSys Int input Maximum number of systems identified ! flag Int input Flag for printing system ! sys Type(system) output Final set of tracked systems, for one time level ! TYPE(dat2d) :: wsdat TYPE(system), POINTER :: sys(:) INTEGER :: maxSys, maxI, maxJ, flag INTEGER, ALLOCATABLE :: actSysInd(:) INTENT (IN) wsdat, actSysInd, maxI, maxJ, flag INTENT (OUT) maxSys ! INTENT (IN OUT) sys ! ! Local variables ! ---------------------------------------------------------------- ! ic Int Counter for wave systems ! INTEGER :: ic, nGuys, startInd, endInd, i, j, ind, leng, leng2 INTEGER :: UNISIZE, DIFSIZE REAL, ALLOCATABLE :: sysOrdered(:) REAL, POINTER :: UNIARR(:), DIFARR(:) INTEGER, ALLOCATABLE :: ngbrSysAll(:), sysSortedInd(:) REAL :: TEMP(2), TEMP1, TEMP2 ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! UNIQUE ! SETDIFF ! SORT ! ! 5. Subroutines calling ! ! combineWaveSystems ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! - ! ! 9. Switches : ! ! None defined yet. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / ! Initialize sys structure IF (flag.NE.2) THEN ! Allocate data structure with the final wave systems WRITE(20,*) 'In printFinalSys...' maxSys = SIZE(actSysInd) NULLIFY(sys) ALLOCATE( sys(maxSys) ) WRITE(20,*) 'Allocated sys okay, SIZE(sys) =',SIZE(sys) ALLOCATE( ngbrSysAll(50*maxI*maxJ) ) !Large enough? DO ic = 1, maxSys NULLIFY( sys(ic)%hs ) NULLIFY( sys(ic)%tp ) NULLIFY( sys(ic)%dir ) NULLIFY( sys(ic)%dspr ) ! NULLIFY( sys(ic)%wf ) NULLIFY( sys(ic)%i ) NULLIFY( sys(ic)%j ) NULLIFY( sys(ic)%lat ) NULLIFY( sys(ic)%lon ) ALLOCATE( sys(ic)%hs(maxI*maxJ) ) ALLOCATE( sys(ic)%tp(maxI*maxJ) ) ALLOCATE( sys(ic)%dir(maxI*maxJ) ) ALLOCATE( sys(ic)%dspr(maxI*maxJ) ) ! ALLOCATE( sys(ic)%wf(maxI*maxJ) ) ALLOCATE( sys(ic)%i(maxI*maxJ) ) ALLOCATE( sys(ic)%j(maxI*maxJ) ) ALLOCATE( sys(ic)%lat(maxI*maxJ) ) ALLOCATE( sys(ic)%lon(maxI*maxJ) ) sys(ic)%hs(:) = 9999. !Optimize this further? sys(ic)%tp(:) = 9999. sys(ic)%dir(:) = 9999. sys(ic)%dspr(:) = 9999. ! sys(ic)%wf(:) = 9999. sys(ic)%i(:) = 9999 sys(ic)%j(:) = 9999 sys(ic)%lat(:) = 9999. sys(ic)%lon(:) = 9999. sys(ic)%sysInd = 9999 sys(ic)%hsMean = 9999. sys(ic)%tpMean = 9999. sys(ic)%dirMean = 9999. sys(ic)%nPoints = 0 sys(ic)%ngbr(:) = 9999 sys(ic)%grp = 9999 ngbrSysAll(:) = 0 startInd=1 nGuys=0 DO i = 1, maxI DO j = 1, maxJ ! ind=wsdat.par(i,j).sys==ic; DO ind = 1, SIZE(wsdat%par(i,j)%sys) !40.81 !Optimize this? IF (wsdat%par(i,j)%sys(ind).EQ.actSysInd(ic)) & THEN nGuys=nGuys+1 sys(ic)%hs(nGuys)=wsdat%par(i,j)%hs(ind) sys(ic)%tp(nGuys)=wsdat%par(i,j)%tp(ind) sys(ic)%dir(nGuys)=wsdat%par(i,j)%dir(ind) sys(ic)%dspr(nGuys)=wsdat%par(i,j)%dspr(ind) ! sys(ic)%wf(nGuys)=wsdat%par(i,j)%wf(ind) sys(ic)%i(nGuys)=i sys(ic)%j(nGuys)=j sys(ic)%lat(nGuys)=wsdat%lat(i,j) sys(ic)%lon(nGuys)=wsdat%lon(i,j) leng = LENGTH(REAL(wsdat%par(i,j)%ngbrSys), & SIZE(wsdat%par(i,j)%ngbrSys),REAL(9999)) endInd = startInd + leng-1 ngbrSysAll(startInd:endInd) = & wsdat%par(i,j)%ngbrSys(1:leng) startInd=endInd+1 END IF END DO END DO END DO ! if ~isempty(sys) IF (nGuys.GT.0) THEN sys(ic)%sysInd=ic sys(ic)%hsMean = SUM(sys(ic)%hs(1:nGuys))/nGuys sys(ic)%tpMean = SUM(sys(ic)%tp(1:nGuys))/nGuys ! sys(ic)%dirMean=mean_angle_single(sys(ic).dir) 40.81 Replaced with two-argument mean_angleV2 sys(ic)%dirMean = & mean_angleV2(sys(ic)%dir(1:nGuys),nGuys) !070512----------- Weight averages with Hm0 --------------------- TEMP1 = 0. TEMP2 = 0. DO i = 1,nGuys TEMP1 = TEMP1 + (sys(ic)%hs(i)**2)*sys(ic)%hs(i) TEMP2 = TEMP2 + (sys(ic)%hs(i)**2)*sys(ic)%tp(i) END DO sys(ic)%hsMean = & TEMP1/MAX(SUM(sys(ic)%hs(1:nGuys)**2),0.001) sys(ic)%tpMean = & TEMP2/MAX(SUM(sys(ic)%hs(1:nGuys)**2),0.001) sys(ic)%dirMean = mean_angleV3(sys(ic)%dir(1:nGuys), & sys(ic)%hs(1:nGuys),nGuys) !070512----------- Weight averages with Hm0 --------------------- sys(ic)%nPoints = nGuys IF (endInd.GT.0) THEN CALL UNIQUE(REAL(ngbrSysAll(1:endInd)),endInd, & UNIARR,UNISIZE) TEMP = (/REAL(sys(ic)%sysInd),REAL(sys(ic)%sysInd)/) CALL SETDIFF(REAL(UNIARR),UNISIZE, & TEMP,2,DIFARR,DIFSIZE) DIFSIZE = MIN(DIFSIZE,SIZE(sys(ic)%ngbr)) sys(ic)%ngbr(1:DIFSIZE) = NINT(DIFARR(1:DIFSIZE)) IF (ASSOCIATED(UNIARR)) DEALLOCATE(UNIARR) IF (ASSOCIATED(DIFARR)) DEALLOCATE(DIFARR) END IF ELSE CYCLE END IF END DO IF (ALLOCATED(ngbrSysAll)) DEALLOCATE(ngbrSysAll) END IF ! Print the sorted field to the screen leng = LENGTH(REAL(sys(:)%nPoints), & SIZE(sys(:)%nPoints),REAL(9999)) ALLOCATE( sysOrdered(leng) ) ALLOCATE( sysSortedInd(leng) ) CALL SORT (REAL(sys(:)%nPoints),leng, & sysOrdered,sysSortedInd,'D') leng = LENGTH(REAL(sysOrdered), & SIZE(sysOrdered),REAL(0)) DO ic = 1, leng leng2 = LENGTH(REAL(sys(sysSortedInd(ic))%ngbr), & SIZE(sys(sysSortedInd(ic))%ngbr),REAL(9999)) END DO IF (ALLOCATED(sysOrdered)) DEALLOCATE(sysOrdered) IF (ALLOCATED(sysSortedInd)) DEALLOCATE(sysSortedInd) RETURN END SUBROUTINE printFinalSys !/ End of printFinalSys ---------------------------------------------- / !/ !/ ------------------------------------------------------------------- / SUBROUTINE combineSys (wsdat ,sys ,maxSys ,maxI , & maxJ ,actSysInd,perKnob ,dirKnob ) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 4-Jan-2013 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! Combine wave systems ! ! 2. Method ! ! - ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! wsdat Type(dat2d) input Combined data structure ! maxI, maxJ Int input Maximum indices of wave field ! sys Type(system) output Final set of tracked systems, for one time level ! maxSys Int input Number of systems ! dirKnob Real input Parameter in direction for combining fields in space ! perKnob Real input Parameter in period for combining fields in space ! TYPE(dat2d) :: wsdat !40.PAR TYPE(system), POINTER :: sys(:) !40.PAR INTEGER :: maxSys, maxI, maxJ !40.PAR INTEGER, ALLOCATABLE :: actSysInd(:) REAL :: perKnob ,dirKnob REAL :: dx, m1, m2 INTENT (IN) maxI, maxJ, perKnob, dirKnob !40.PAR ! INTENT (IN OUT) wsdat, sys, maxSys !40.PAR ! ! Local variables ! ---------------------------------------------------------------- ! ngbIndex Int Arr Array of neighbours ! INTEGER, ALLOCATABLE :: sysSortedInd(:), sysOut(:) INTEGER, POINTER :: indSys1(:), indSys2(:) REAL, ALLOCATABLE :: sysOrdered(:), rounded(:) REAL, POINTER :: uniarr(:), difarr(:), allngbr(:) INTEGER :: leng, leng2, s, ss, so, ngb, lsys, lsys2, hh, i, j, & ii, jj, ind, ind2, nn, nbr, icEnd,ic,iii,iloop INTEGER :: myngbr, indMatch, matchSys, keep, replacedInd, & hhForIndMatch, lMatch, tot, outsize INTEGER :: ngbIndex(10000), keepInd(maxI*maxJ), oneLess(1000) !Array large enough? ! REAL :: Tb,deltaPerB,deltaDirB,absDir,absPer,absHs,absWf REAL :: Tb,deltaPerB,deltaDirB,deltaHsB,absDir,absPer,absHs LOGICAL :: file_exists INTEGER :: MASK(maxI,maxJ) REAL :: lonmean, latmean, DIST !061512 ----------------------------------------------- LOGICAL :: ZIPMATCH INTEGER :: counter, count2, izp, izp2, in, jn, icnt, ngbrExt REAL :: T, ngb_tp, ngb_dir REAL :: ngbmatch(maxI*maxJ) TYPE(neighbr) :: ngbr(50) !061512 ----------------------------------------------- REAL :: TEMP1, TEMP2 INTEGER :: actSys ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! SORT ! findIJV4 ! UNIQUE ! combinePartitionsV2 ! UNION ! SETDIFF ! ! 5. Subroutines calling ! ! combineWaveSystems ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! - ! ! 9. Switches : ! ! None defined yet. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / ! Initialize pointer (first use) NULLIFY(indSys1) NULLIFY(indSys2) ! Flag to combine systems on a point-by-point basis along boundary, ! instead of using mean values. ZIPMATCH = .FALSE. ngbrExt = 1 ! Combine systems on the basis of tpMean ALLOCATE( sysOrdered(maxSys) ) ALLOCATE( sysSortedInd(maxSys) ) ALLOCATE( sysOut(maxSys) ) ALLOCATE( rounded(maxSys) ) ! Sort in descending Tp: the following improves the iterative combining in ! the special case that the wave period is constant over the domain, but ! tpMean is not because of truncation errors at very high decimals. rounded = REAL(INT(sys(1:maxSys)%tpMean*1.E4))*1.E-4 CALL SORT(rounded,maxSys,sysOrdered,sysSortedInd,'D') sysOut=sys(sysSortedInd)%sysInd IF (ALLOCATED(rounded)) DEALLOCATE(rounded) !051612 --- Land mask addition MASK(:,:) = 0 INQUIRE(FILE="sys_mask.ww3", EXIST=file_exists) IF (file_exists) THEN WRITE(20,*) '*** Using land mask' OPEN(unit=13,file='sys_mask.ww3',status='old') DO j = maxJ,1,-1 READ(13,*) (MASK(i,j), i=1,maxI) END DO CLOSE(13) END IF !051612 --- Land mask addition !opt WRITE(20,*) 'SIZE(sysOut)=',SIZE(sysOut) DO so = 1, SIZE(sysOut) ! WRITE(20,*) 'so =',so s = sysOut(so) ss = FINDFIRST(REAL(sys(:)%sysInd),SIZE(sys(:)%sysInd), & REAL(s)) !opt WRITE(20,*) 's,ss=',s,ss ngbIndex(:) = 0 ii = 1 leng = LENGTH(REAL(sys(ss)%ngbr),SIZE(sys(ss)%ngbr), & REAL(9999)) ! Identify the indices of all the systems that neighbour the current system s, ! store in ngbIndex(:) DO ngb = 1, leng IF ( sys(ss)%ngbr(ngb).NE.s ) THEN myngbr = 1 DO WHILE (myngbr.LE.SIZE(sysOut)) IF (sys(myngbr)%sysInd.EQ.sys(ss)%ngbr(ngb)) THEN ngbIndex(ii) = myngbr ii = ii+1 IF (ii.GT.1000) & WRITE(20,*) '*** WARNING: ngbIndex(:) exceeded!' END IF myngbr = myngbr+1 END DO END IF END DO ii = ii-1 !opt WRITE(20,*) so,'. sys =',s,', Tp =',sys(s)%tpMean, & !opt ', size=',sys(s)%nPoints,', #neighbours=',ii IF ( ii.GT.0 ) THEN DO ngb = 1, ii ! We first need to find the (i,j) points that are either common ! to both these systems, or at the boundary of the two systems. Here ! sys 1 will carry the 'ss' index and sys 2 the ngbIndex(ngb) index. CALL findIJV4 (sys(ss),sys(ngbIndex(ngb)), & maxI,maxJ,indSys1,indSys2) IF ((SIZE(indSys1)>10).AND.(SIZE(indSys2)>10).AND. & (sys(ss)%nPoints.GT.sys(ngbIndex(ngb))%nPoints)) & THEN lsys = SIZE(indSys1) lsys2 = SIZE(indSys2) !061512---------------Add zipper compare IF (ZIPMATCH) THEN ! Omit small systems to save time IF ((sys(ss)%nPoints.LT.5).OR. & (sys(ngbIndex(ngb))%nPoints.LT.5)) THEN CYCLE END IF dx=0.5*((wsdat%lon(2,1)-wsdat%lon(1,1)) + & (wsdat%lat(1,2)-wsdat%lat(1,1))) ngbmatch(:)=0. DO izp = 1,lsys ! Find neighbors of this point counter=0 DO in=(sys(ss)%i(indSys1(izp))-ngbrExt), & (sys(ss)%i(indSys1(izp))+ngbrExt) DO jn=(sys(ss)%j(indSys1(izp))-ngbrExt), & (sys(ss)%j(indSys1(izp))+ngbrExt) counter=counter+1 ngbr(counter)%i = in ngbr(counter)%j = jn END DO END DO ! Find these points in neighboring system ngb_tp = 0. ngb_dir = 0. count2 = 0 DO izp2 = 1,lsys2 DO icnt = 1,counter IF ((sys(ngbIndex(ngb))%i(indSys2(izp2)) & .EQ.ngbr(icnt)%i).AND. & (sys(ngbIndex(ngb))%j(indSys2(izp2)) & .EQ.ngbr(icnt)%j)) THEN count2 = count2+1 ngb_tp = ngb_tp + & sys(ngbIndex(ngb))%tp(indSys2(izp2)) ngb_dir = ngb_dir + & sys(ngbIndex(ngb))%dir(indSys2(izp2)) END IF END DO END DO IF (count2.GT.0) THEN absPer = ABS(sys(ss)%tp(indSys1(izp))-ngb_tp/count2) absDir = ABS(sys(ss)%dir(indSys1(izp))-ngb_dir/count2) T = sys(ss)%tp(indSys1(izp)) m1 = -3.645*T + 63.211 m1 = MAX(m1,10.) m2 = -0.346*T + 3.686 m2 = MAX(m2,0.6) deltaDirB = (m1*dx + dirKnob)*1. deltaPerB = (m2*dx + perKnob)*1. IF ( (absPer.LT.deltaPerB).AND. & (absDir.LT.deltaDirB) ) THEN ngbmatch(izp)=1. END IF END IF END DO ! If >80% of neighbors fall within criteria, system is matched IF ((SUM(ngbmatch(1:lsys))/lsys).GT.0.50) THEN indMatch = ngbIndex(ngb) matchSys = sys(indMatch)%sysInd ELSE CYCLE END IF ELSE !061512--------------------------------- Tb = MAX(SUM(sys(ss)%tp(indSys1))/lsys, & SUM(sys(ngbIndex(ngb))%tp(indSys2))/lsys2) ! deltaPerB = (-0.06*Tb+2+perKnob)*1.5 ! deltaDirB = (-Tb+(25+10*dirKnob))*1.5 ! deltaPerB = (-0.06*Tb+2+2)*1.5 ! deltaDirB = (-Tb+(25+10*2))*1.5 dx=0.5*((wsdat%lon(2,1)-wsdat%lon(1,1)) + & (wsdat%lat(1,2)-wsdat%lat(1,1))) m1 = -3.523*Tb + 64.081 m1 = MAX(m1,10.) m2 = -0.337*Tb + 3.732 m2 = MAX(m2,0.6) !1stddev m1 = -2.219*Tb + 35.734 !1stddev m1 = MAX(m1,5.) !1stddev m2 = -0.226*Tb + 2.213 !1stddev m2 = MAX(m2,0.35) !5stddev m1 = -5.071*Tb + 90.688 !5stddev m1 = MAX(m1,16.) !5stddev m2 = -0.467*Tb + 5.161 !5stddev m2 = MAX(m2,1.0) deltaDirB = (m1*1. + dirKnob)*1. deltaPerB = (m2*1. + perKnob)*1. deltaHsB = 0.50*SUM(sys(ss)%hs(indSys1))/lsys ! deltaHsB = 0.25*SUM(sys(ss)%hs(indSys1))/lsys !051612 --- Land mask addition ! Option 1: If system centroid is near a land mask (e.g. 3 arc-deg), ! increase the tolerances IF (ANY(MASK.EQ.1)) THEN lonmean = SUM(sys(ss)%lon(indSys1))/lsys latmean = SUM(sys(ss)%lat(indSys1))/lsys DO j = 1,maxJ DO i = 1,maxI IF (MASK(i,j).EQ.1) THEN ! Land point found. Compute distance to system centroid DIST = SQRT((lonmean-wsdat%lon(i,j))**2 +& (latmean-wsdat%lat(i,j))**2) IF (DIST.LT.3.) THEN ! System assumed to be influenced by land, ! increase tolerances to deltaDirB=30,deltaPerB=3 ! deltaDirB = (m1*1. + 30)*1. ! deltaPerB = (m2*1. + 3)*1. deltaDirB = (m1*1. + 30)*1. deltaPerB = (m2*1. + 3)*1. !Remove dHs limitation from criteria deltaHsB = 9999. GOTO 500 END IF END IF END DO END DO END IF 500 CONTINUE !051612 --- Land mask addition absHs = ABS( SUM(sys(ss)%hs(indSys1))/lsys - & SUM(sys(ngbIndex(ngb))%hs(indSys2))/lsys2 ) absPer = ABS( SUM(sys(ss)%tp(indSys1))/lsys - & SUM(sys(ngbIndex(ngb))%tp(indSys2))/lsys2 ) absDir = ABS( & mean_angleV2(sys(ss)%dir(indSys1),lsys) - & mean_angleV2(sys(ngbIndex(ngb))%dir(indSys2), & lsys2) ) IF (absDir.GT.180) absDir = 360.-absDir ! absWf = ABS( SUM(sys(ss)%wf(indSys1))/lsys - & ! SUM(sys(ngbIndex(ngb))%wf(indSys2))/lsys2 ) IF ( (absPer.LT.deltaPerB).AND. & (absDir.LT.deltaDirB).AND. & (absHs.LT.deltaHsB) ) THEN indMatch = ngbIndex(ngb) matchSys = sys(indMatch)%sysInd !opt WRITE(20,*) '-> Matched sys',s, & !opt 'with neighbor sys',matchSys ELSE CYCLE END IF !061512--------------------------------- END IF !061512--------------------------------- keep = 0 keepInd(:) = 0 DO hh = 1, sys(ss)%nPoints ii = sys(ss)%i(hh) jj = sys(ss)%j(hh) ind = 0 ind = FINDFIRST(REAL(wsdat%par(ii,jj)%sys), & SIZE(wsdat%par(ii,jj)%sys),REAL(s)) !Shouldn't REAL(s) be matchSys... IF (ind.NE.0) THEN wsdat%par(ii,jj)%sys(ind)=matchSys !...and matchSys be s, (i.e. add the matching neigbour to the base?) END IF ! Remove the "-1" system from the set 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 ! Combine any partitions assigned to the same systems ! Check for duplicates IF (ind2.EQ.0) & WRITE(20,*) '***2.Calling UNIQUE w. len=0!' CALL UNIQUE(REAL(oneLess(1:ind2)),ind2, & uniarr,outsize) IF (ASSOCIATED(uniarr)) DEALLOCATE(uniarr) IF (ind2.GT.outsize) THEN ! There is at least one duplicate, so combine systems CALL combinePartitionsV2(wsdat%par(ii,jj)) ! 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 = & FINDFIRST(REAL(wsdat%par(ii,jj)%sys(:)), & SIZE(wsdat%par(ii,jj)%sys(:)), & REAL(matchSys)) 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) sys(indMatch)%dspr(hhForIndMatch) = & wsdat%par(ii,jj)%dspr(replacedInd) ! sys(indMatch)%wf(hhForIndMatch) = & ! wsdat%par(ii,jj)%wf(replacedInd) ELSE keep = keep+1 keepInd(keep) = hh END IF END DO leng = LENGTH(REAL(sys(indMatch)%hs), & SIZE(sys(indMatch)%hs),REAL(9999.)) ! 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 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 outsize = MIN(outsize,size(sys(indMatch)%ngbr)) sys(indMatch)%ngbr(1:outsize) = NINT(difarr(1:outsize)) IF (ASSOCIATED(allngbr)) DEALLOCATE(allngbr) IF (ASSOCIATED(difarr)) DEALLOCATE(difarr) 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 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) !070512----------- Weight averages with Hm0 --------------------- TEMP1 = 0. TEMP2 = 0. DO iii = 1,keep TEMP1 = TEMP1 + (sys(ss)%hs(keepInd(iii))**2)*& sys(ss)%hs(keepInd(iii)) TEMP2 = TEMP2 + (sys(ss)%hs(keepInd(iii))**2)*& sys(ss)%tp(keepInd(iii)) END DO DO iii = 1,leng TEMP1 = TEMP1 + (sys(indMatch)%hs(iii)**2)*& sys(indMatch)%hs(iii) TEMP2 = TEMP2 + (sys(indMatch)%hs(iii)**2)*& sys(indMatch)%tp(iii) END DO sys(indMatch)%hsMean = TEMP1/MAX(SUM((/ & sys(ss)%hs(keepInd(1:keep))**2, & sys(indMatch)%hs(1:leng)**2 /)),0.001) sys(indMatch)%tpMean = TEMP2/MAX(SUM((/ & sys(ss)%hs(keepInd(1:keep))**2, & sys(indMatch)%hs(1:leng)**2 /)),0.001) sys(indMatch)%dirMean = & mean_angleV3((/ sys(ss)%dir(keepInd(1:keep)), & sys(indMatch)%dir(1:leng) /), & (/ sys(ss)%hs(keepInd(1:keep)), & sys(indMatch)%hs(1:leng) /),tot) !070512----------- Weight averages with Hm0 --------------------- 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)%dspr(1:(keep+leng)) = & (/sys(ss)%dspr(keepInd(1:keep)), & sys(indMatch)%dspr(1:leng)/) ! sys(indMatch)%wf(1:(keep+leng)) = & ! (/sys(ss)%wf(keepInd(1:keep)), & ! sys(indMatch)%wf(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)) ! Clear array of system that has just been combined with another sys(ss)%nPoints = 0 sys(ss)%ngbr(:) = 9999 WRITE(20,*) 'Deallocating sys',s DEALLOCATE( sys(ss)%hs ) !opt DEALLOCATE( sys(ss)%tp ) !opt DEALLOCATE( sys(ss)%dir ) !opt DEALLOCATE( sys(ss)%dspr ) !opt ! DEALLOCATE( sys(ss)%wf ) !opt DEALLOCATE( sys(ss)%i ) !opt DEALLOCATE( sys(ss)%j ) !opt DEALLOCATE( sys(ss)%lat ) !opt DEALLOCATE( sys(ss)%lon ) !opt ! DEALLOCATE( sys(ss)%hsMean ) !opt ! DEALLOCATE( sys(ss)%tpMean ) !opt ! DEALLOCATE( sys(ss)%dirMean ) !opt ! 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 leng = LENGTH(REAL(wsdat%par(i,j)%ngbrSys), & SIZE(wsdat%par(i,j)%ngbrSys),REAL(9999)) 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) IF (ASSOCIATED(uniarr)) DEALLOCATE(uniarr) ELSE wsdat%par(i,j)%ngbrSys(:) = 9999 END IF END DO END DO ! Update neigbors in sys structure DO nn = 1, maxSys nbr = FINDFIRST(REAL(sys(nn)%ngbr), & SIZE(sys(nn)%ngbr),REAL(s)) IF (nbr.NE.0) THEN ! WRITE(20,*) 'update' sys(nn)%ngbr(nbr)=matchSys END IF leng2 = LENGTH(REAL(sys(nn)%ngbr), & SIZE(sys(nn)%ngbr),REAL(9999)) IF (leng2.GT.0) THEN CALL UNIQUE(REAL(sys(nn)%ngbr(1:leng2)), & leng2,uniarr,outsize) sys(nn)%ngbr(:) = 9999 sys(nn)%ngbr(1:outsize) = NINT(uniarr) IF (ASSOCIATED(uniarr)) DEALLOCATE(uniarr) ! WRITE(20,*) 'has now ngbr: ', & ! sys(nn)%ngbr(1:outsize) END IF END DO EXIT END IF IF (ASSOCIATED(indSys1)) DEALLOCATE(indSys1) IF (ASSOCIATED(indSys2)) DEALLOCATE(indSys2) END DO END IF END DO IF (ALLOCATED(sysOrdered)) DEALLOCATE(sysOrdered) IF (ALLOCATED(sysSortedInd)) DEALLOCATE(sysSortedInd) IF (ALLOCATED(sysOut)) DEALLOCATE(sysOut) ! Compile array index of active systems in sys actSys = 0 DO ic = 1,maxSys IF (sys(ic)%nPoints>0) actSys = actSys + 1 END DO IF (ALLOCATED(actSysInd)) DEALLOCATE(actSysInd) ALLOCATE( actSysInd(actSys) ) actSys = 0 DO ic = 1,maxSys IF (sys(ic)%nPoints>0) THEN actSys = actSys + 1 actSysInd(actSys) = sys(ic)%sysInd END IF END DO !opt WRITE(20,*) 'actSys =',actSys !opt WRITE(20,*) 'actSysInd =',actSysInd !opt DO ic = 1,SIZE(actSysInd) !opt s = actSysInd(ic) !opt WRITE(20,*) 'sys(',s,')%sysInd =',sys(s)%sysInd !opt END DO WRITE(20,*) 'Leaving combineSys...' RETURN END SUBROUTINE combineSys !/ End of combineSys ------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / SUBROUTINE combinePartitionsV2 (dat) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 4-Jan-2013 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! Combine two partitions that have been assigned to the same system ! ! 2. 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! ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! dat TYPE(param) in/out Input data structure (partitions set) ! to combine ! TYPE(param) :: dat INTENT (IN OUT) dat ! ! Local variables ! ---------------------------------------------------------------- TYPE duplicate 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 ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! - ! ! 5. Subroutines calling ! ! findSys ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! - ! ! 9. Switches : ! ! None defined yet. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / ! 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. 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 dup(nsys)%ndup = dup(nsys)%ndup+1 dup(nsys)%ind(dup(nsys)%ndup) = pp END IF END IF END DO 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 ! 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 ! 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%dspr( dup(s)%ind(p):(npart-1) ) = & dat%dspr( (dup(s)%ind(p)+1):npart) ! dat%wf( dup(s)%ind(p):(npart-1) ) = & ! dat%wf( (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%dspr(npart) = 9999. ! dat%wf(npart) = 9999. dat%sys(npart) = 9999 dat%ipart(npart) = 0 END IF END DO END DO RETURN END SUBROUTINE combinePartitionsV2 !/ End of combinePartitionsV2 ---------------------------------------- / !/ !/ ------------------------------------------------------------------- / REAL FUNCTION mean_angleV2(ang,ll) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 4-Jan-2013 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! Compute the mean direction from array of directions ! ! 2. 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. ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ang Real input Array of angles to average ! ll Int input Length of ang ! REAL :: ang(ll) INTEGER :: ll ! ! 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 ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! - ! ! 5. Subroutines calling ! ! findSys ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! - ! ! 9. Switches : ! ! None defined yet. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / ! 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 !/ End of mean_angleV2 ----------------------------------------------- / !/ !/ ------------------------------------------------------------------- / REAL FUNCTION mean_angleV3(ang,hsign,ll) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 4-Jan-2013 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! Compute the mean direction from array of directions, ! INCLUDING WEIGHTING WITH HMO ! ! 2. 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. ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ang Real input Array of angles to average ! ll Int input Length of ang ! REAL :: ang(ll), hsign(ll) REAL :: TEMP1, TEMP2 INTEGER :: ll ! ! 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 INTEGER :: i ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! - ! ! 5. Subroutines calling ! ! findSys ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! - ! ! 9. Switches : ! ! None defined yet. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / ! North and East components v(:) = COS(ang(:)*(PI/180.)) u(:) = SIN(ang(:)*(PI/180.)) TEMP1 = 0. TEMP2 = 0. DO i = 1,ll TEMP1 = TEMP1 + (hsign(i)**2)*v(i) TEMP2 = TEMP2 + (hsign(i)**2)*u(i) END DO vm = TEMP1/MAX(SUM(hsign**2),0.001) um = TEMP2/MAX(SUM(hsign**2),0.001) ! 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_angleV3 = theta*(180/PI) RETURN END FUNCTION mean_angleV3 !/ End of mean_angleV3 ----------------------------------------------- / !/ !/ ------------------------------------------------------------------- / SUBROUTINE UNIQUE (INARRAY,INSIZE,OUTARRAY,OUTSIZE) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 22-Dec-2016 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ 12-Dec-2016 : Change algorithm from N*N to N*log(N) !/ (S. Zieger BoM Australia) ( version 5.16 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! Returns the sorted elements that are unique in INARRAY. ! ! 2. Method ! ! 1. Sort input array with quicksort ! 2. Copy sequential-elements if the 'current' element ! is not equal to the 'previous' element in array. ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! INARRAY REAL ARR input Input array ! INSIZE INTEGER input Size of input array ! OUTARRAY REAL ARR output Output array (sorted) ! OUTSIZE INTEGER output Size of output array (number of unique elements) ! INTEGER, INTENT(IN) :: INSIZE INTEGER, INTENT(OUT) :: OUTSIZE REAL, INTENT(IN) :: INARRAY(INSIZE) REAL, POINTER :: OUTARRAY(:) ! ! Local variables ! ---------------------------------------------------------------- INTEGER :: I, K REAL :: ARRAY(INSIZE), TEMP(INSIZE) ! ! 4. Subroutines used : ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! QSORT Subr. Private Quicksort algorithm ! ! 5. Subroutines calling ! ! waveTracking_NWS_V2 ! findSys ! printFinalSys ! combineSys ! findIJV4 ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! - ! ! 9. Switches : ! ! None defined yet. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / K = 1 IF ( INSIZE.EQ.0 ) THEN WRITE(20,*) '*** In Subr. UNIQUE: Input array has length=0!' ELSE !/ --- Setup input arrays and temporary arrays. --- DO I=1,INSIZE ARRAY(I) = INARRAY(I) TEMP(I) = REAL(I) END DO !/ !/ --- Sort input arrays (use temporary array to store indices). --- CALL QSORT(ARRAY,TEMP,1,INSIZE) !/ !/ --- Reset temporary array. --- TEMP(:) = 9999. !/ !/ --- Initialise first values and array index. --- K = 1 TEMP(K) = ARRAY(K) K = K + 1 !/ --- Iterate over elements in array. --- DO I=2,INSIZE !/ --- Compare sequential array values ('previous' less than 'next') !/ and test against the last list element check in. --- IF ( ARRAY(I).GT.ARRAY(I-1) .AND. & ARRAY(I).GT.TEMP(K-1) ) THEN TEMP(K) = ARRAY(I) K = K + 1 END IF END DO !/ --- Allocate output array --- OUTSIZE = K - 1 ALLOCATE(OUTARRAY(OUTSIZE)) !/ --- Transfer output from temporary array to output array. --- IF ( OUTSIZE.GE.1 ) THEN DO I=1,OUTSIZE OUTARRAY(I) = TEMP(I) END DO END IF END IF !/ RETURN !/ END SUBROUTINE UNIQUE !/ End of UNIQUE ----------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / SUBROUTINE SORT (INARRAY,INSIZE,OUTARRAY,IY,DIRECTION) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 20-Dec-2016 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ 20-Dec-2016 : Add quicksort algorithm (S. Zieger) ( version 5.16 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! 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. ! ! 2. Method ! ! Sort algorithm based on quicksort. ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! 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 INTEGER :: IY(INSIZE) REAL :: INARRAY(INSIZE), OUTARRAY(INSIZE) INTENT (IN) INARRAY, INSIZE, DIRECTION INTENT (OUT) OUTARRAY, IY ! ! 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 :: I REAL :: IND(INSIZE) ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! - ! ! 5. Subroutines calling ! ! printFinalSys ! combineSys ! timeTrackingV2 ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! - ! ! 9. Switches : ! ! None defined yet. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / ! Sort OUTARRAY in as/decending order IF (INSIZE.EQ.0) THEN WRITE(20,*) '*** In Subr. SORT: Input array has length=0 !!!' ELSE DO I = 1, INSIZE OUTARRAY(I) = INARRAY(I) IND(I) = REAL(I) END DO IF (DIRECTION .EQ. 'A') THEN CALL QSORT(OUTARRAY,IND,1,INSIZE) ELSE IF (DIRECTION .EQ. 'D') THEN CALL QSORT_DESC(OUTARRAY,IND,1,INSIZE) END IF END IF ! !/ --- Cast index array to integer. --- DO I = 1, INSIZE IY(I) = INT(IND(I)) END DO ! RETURN ! END SUBROUTINE SORT !/ End of SORT ------------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / SUBROUTINE SETDIFF (INARRAY1, INSIZE1, INARRAY2, INSIZE2, & OUTARRAY, OUTSIZE) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 20-Dec-2016 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ 20-Dec-2016 : Add quicksort algorithm (S.Zieger) ( version 5.16 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! (i) Returns the elements in INARRAY1 that are not in INARRAY2. ! (ii) Sort the resulting array in ascending order. ! ! 2. Method ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! 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, POINTER :: OUTARRAY(:) INTENT (IN) INARRAY1, INSIZE1, INARRAY2, INSIZE2 INTENT (OUT) OUTSIZE ! ! Local variables ! ---------------------------------------------------------------- INTEGER :: I,J,K REAL :: TEMP(INSIZE1) REAL :: ARRAY1(INSIZE1),ARRAY2(INSIZE2) REAL :: ID1(INSIZE1),ID2(INSIZE2) LOGICAL :: LOOP ! ! 4. Subroutines used : ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! QSORT Subr. Private Quicksort algorithm ! ! 5. Subroutines calling ! ! printFinalSys ! combineSys ! timeTrackingV2 ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! - ! ! 9. Switches : ! ! None defined yet. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / IF ( (INSIZE1).EQ.0 ) THEN OUTSIZE = 0 ALLOCATE(OUTARRAY(OUTSIZE)) ELSE IF ( INSIZE2.EQ.0 ) THEN CALL UNIQUE(INARRAY1,INSIZE1,OUTARRAY,OUTSIZE) ELSE !/ --- Setup input arrays. --- DO I=1,INSIZE1 ARRAY1(I) = INARRAY1(I) ID1(I) = REAL(I) END DO DO I=1,INSIZE2 ARRAY2(I) = INARRAY2(I) ID2(I) = REAL(I) END DO !/ !/ --- Sort input arrays. --- CALL QSORT(ARRAY1,ID1,1,INSIZE1) CALL QSORT(ARRAY2,ID2,1,INSIZE2) !/ !/ --- Initialise indices. --- I = 1 J = 1 K = 1 !/ !/ --- Allocate and initialize temporary output --- TEMP(:) = 9999. !/ !/ --- Loop though both arrays by incrementing I,J. --- LOOP = .TRUE. DO WHILE ( LOOP ) !/ IF ( ARRAY1(I).LT.ARRAY2(J) .OR. & ARRAY1(I).GT.ARRAY2(INSIZE2) ) THEN !/ --- Populate output array. Check for dumplicates !/ in output array. --- IF ( K.EQ.1 ) THEN TEMP(K) = ARRAY1(I) K = K + 1 ELSE IF ( TEMP(K-1).LT.ARRAY1(I) ) THEN TEMP(K) = ARRAY1(I) K = K + 1 END IF I = I + 1 ELSE IF ( ARRAY2(J).LT.ARRAY1(I) ) THEN J = J + 1 ELSE I = I + 1 J = J + 1 END IF !/ --- Check for exit the loop. --- IF ( I.GT.INSIZE1 ) THEN LOOP = .FALSE. END IF !/ --- Make sure array pointer I,J are within array bounds. --- I = MIN(I,INSIZE1) J = MIN(J,INSIZE2) !/ END DO !/ !/ --- Allocate output array --- OUTSIZE = K-1 ALLOCATE(OUTARRAY(OUTSIZE)) !/ --- Transfer output from temporary array to output array. --- DO I=1,OUTSIZE OUTARRAY(I) = TEMP(I) END DO END IF !/ RETURN !/ END SUBROUTINE SETDIFF !/ End of SETDIFF ---------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / SUBROUTINE INTERSECT (INARRAY1 ,INSIZE1 ,INARRAY2 ,INSIZE2 , & OUTARRAY ,OUTSIZE ,IND1 ,IND2 ) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 20-Dec-2016 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ 20-Dec-2016 : Add count-histogram method based on !/ algorithm from Mirko Velic (BoM) !/ (S. Zieger BoM, Australia) ( version 5.16 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! (i) Returns the elements that are mutual in INARRAY1 and INARRAY2. ! (ii) Sort the resulting array in ascending order. ! ! 2. Method ! ! Sort with counting/histogram method with input array being ! cast as integer. ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! 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 ! intersects) ! INTEGER :: INSIZE1, INSIZE2, OUTSIZE REAL :: INARRAY1(INSIZE1), INARRAY2(INSIZE2) REAL, POINTER :: OUTARRAY(:) REAL, POINTER :: IND1(:), IND2(:) ! INTENT (IN) INARRAY1, INSIZE1, INARRAY2, INSIZE2 INTENT (OUT) OUTSIZE ! ! Local variables ! ---------------------------------------------------------------- ! VIDX1, VIDX2 - array(s) in which the value is represented by ! its index (i.e. histogram with frequency 1) ! N - data range and size of possible intersections. ! LOGICAL,ALLOCATABLE :: VIDX1(:),VIDX2(:) INTEGER,ALLOCATABLE :: IPOS1(:),IPOS2(:) ! INTEGER :: I, J INTEGER :: N, IMIN, IMAX INTEGER :: MINV1,MAXV1, MINV2, MAXV2 ! ! 4. Subroutines used : ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! ! 5. Subroutines calling ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! - ! ! 9. Switches : ! ! None defined yet. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / ! OUTSIZE = 0 !/ --- Calculate the range of the two sets. --- MINV1 = INT(MINVAL(INARRAY1)) MAXV1 = INT(MAXVAL(INARRAY1)) MINV2 = INT(MINVAL(INARRAY2)) MAXV2 = INT(MAXVAL(INARRAY2)) !/ --- Check if ranges overlap. --- IF ( MAXV1.LT.MINV2.OR.INSIZE1.EQ.0.OR.INSIZE2.EQ.0 ) THEN ALLOCATE(OUTARRAY(OUTSIZE)) ALLOCATE(IND1(OUTSIZE)) ALLOCATE(IND2(OUTSIZE)) ELSE !/ --- Calculate size of temporary output arrays. Allow !/ extra elements: ZERO, and make sure index is 1:N. --- IMIN = MIN(MINV1,MINV2)-1 IMAX = MAX(MAXV1,MAXV2)+1 N = IMAX-IMIN ALLOCATE(VIDX1(N),VIDX2(N)) ALLOCATE(IPOS1(N),IPOS2(N)) VIDX1(1:N) = .FALSE. VIDX2(1:N) = .FALSE. DO I=1,INSIZE1 J = INT(INARRAY1(I)-IMIN) VIDX1(J) = .TRUE. IPOS1(J) = I END DO DO I=1,INSIZE2 J = INT(INARRAY2(I)-IMIN) !/ --- Intersect arrays and check for !/ duplicate elements in array2. --- IF ( VIDX1(J).AND..NOT.VIDX2(J) ) THEN OUTSIZE = OUTSIZE + 1 VIDX2(J) = .TRUE. IPOS2(J) = I END IF END DO !/ --- Allocate output arrays. --- ALLOCATE(OUTARRAY(OUTSIZE)) ALLOCATE(IND1(OUTSIZE)) ALLOCATE(IND2(OUTSIZE)) !/ --- Transfer contents. --- I = 1 DO J=1,N IF ( VIDX1(J).AND.VIDX2(J).AND.I.LE.OUTSIZE ) THEN OUTARRAY(I) = INARRAY1(IPOS1(J)) IND1(I) = IPOS1(J) IND2(I) = IPOS2(J) I = I + 1 END IF END DO !/ --- Free memory. --- DEALLOCATE(VIDX1,VIDX2) DEALLOCATE(IPOS1,IPOS2) END IF !/ RETURN !/ END SUBROUTINE INTERSECT !/ End of INTERSECT -------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / SUBROUTINE UNION (INARRAY1, INSIZE1, INARRAY2, INSIZE2, & OUTARRAY, OUTSIZE) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 4-Jan-2013 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ 20-Dec-2016 : Add count-histogram method similarly !/ to INTERSECT (S. Zieger) ( version 5.16 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! (i) Returns the union of INARRAY1 and INARRAY2. ! (ii) Sort the resulting array in ascending order. ! ! 2. Method ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! INARRAY REAL ARR input Input array ! INSIZE INTEGER input Size of input array ! OUTARRAY REAL ARR output Output array (sorted) ! OUTSIZE INTEGER output Size of output array (number of ! unique elements) INTEGER :: INSIZE1, INSIZE2, OUTSIZE REAL :: INARRAY1(INSIZE1), INARRAY2(INSIZE2) REAL, POINTER :: OUTARRAY(:) ! INTENT (IN) INARRAY1, INSIZE1, INARRAY2, INSIZE2 INTENT (OUT) OUTSIZE ! ! Local variables ! ---------------------------------------------------------------- ! VIDX1, VIDX2 - array(s) in which the value is represented by ! its index (i.e. histogram with frequency 1) ! N - data range and size of possible intersections. ! LOGICAL,ALLOCATABLE :: VIDX1(:),VIDX2(:) INTEGER,ALLOCATABLE :: IPOS1(:),IPOS2(:) REAL,ALLOCATABLE :: TEMP(:) ! INTEGER :: I, J INTEGER :: N, IMIN, IMAX INTEGER :: MINV1,MAXV1, MINV2, MAXV2 ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! QSORT Subr. Private Quicksort algorithm ! ! 5. Subroutines calling ! ! combineSys ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! - ! ! 9. Switches : ! ! None defined yet. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/ --- Setup input arrays. --- IF ( (INSIZE1+INSIZE2).EQ.0 ) THEN OUTSIZE = 0 ALLOCATE(OUTARRAY(OUTSIZE)) ELSEIF ( INSIZE1.EQ.0 ) THEN OUTSIZE = INSIZE2 ALLOCATE(OUTARRAY(OUTSIZE)) ALLOCATE(TEMP(OUTSIZE)) DO I=1,OUTSIZE OUTARRAY(I) = INARRAY2(I) TEMP(I) = REAL(I) END DO CALL QSORT(OUTARRAY,TEMP,1,OUTSIZE) ELSEIF ( INSIZE2.EQ.0 ) THEN OUTSIZE = INSIZE1 ALLOCATE(OUTARRAY(OUTSIZE),TEMP(OUTSIZE)) DO I=1,OUTSIZE OUTARRAY(I) = INARRAY1(I) TEMP(I) = REAL(I) END DO CALL QSORT(OUTARRAY,TEMP,1,OUTSIZE) ELSE OUTSIZE = 0 !/ --- Calculate the range of the two sets. --- MINV1 = INT(MINVAL(INARRAY1)) MAXV1 = INT(MAXVAL(INARRAY1)) MINV2 = INT(MINVAL(INARRAY2)) MAXV2 = INT(MAXVAL(INARRAY2)) ! !/ --- Allow extra elementes: ZERO, and make sure index is 1:N. --- IMIN = MIN(MINV1,MINV2)-1 IMAX = MAX(MAXV1,MAXV2)+1 N = IMAX-IMIN ALLOCATE(VIDX1(N),VIDX2(N)) ALLOCATE(IPOS1(N),IPOS2(N)) VIDX1(1:N) = .FALSE. VIDX2(1:N) = .FALSE. IPOS1(1:N) = -9999 IPOS2(1:N) = -9999 !/ DO I=1,INSIZE1 J = INT(INARRAY1(I)-IMIN) IF ( .NOT.VIDX1(J) ) THEN OUTSIZE = OUTSIZE + 1 VIDX1(J) = .TRUE. IPOS1(J) = I END IF END DO DO I=1,INSIZE2 J = INT(INARRAY2(I)-IMIN) IF ( .NOT.VIDX1(J).AND..NOT.VIDX2(J) ) THEN OUTSIZE = OUTSIZE + 1 VIDX2(J) = .TRUE. IPOS2(J) = I END IF END DO ALLOCATE(OUTARRAY(OUTSIZE)) I = 1 DO J=1,N IF ( VIDX1(J).AND.I.LE.OUTSIZE ) THEN OUTARRAY(I) = INARRAY1(IPOS1(J)) I = I + 1 ELSEIF ( VIDX2(J).AND.I.LE.OUTSIZE ) THEN OUTARRAY(I) = INARRAY2(IPOS2(J)) I = I + 1 END IF END DO DEALLOCATE(VIDX1,VIDX2) DEALLOCATE(IPOS1,IPOS2) END IF !/ RETURN !/ END SUBROUTINE UNION !/ End of UNION ------------------------------------------------------ / !/ !/ ------------------------------------------------------------------- / INTEGER FUNCTION LENGTH(ARRAY,ARRSIZE,VAL) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 4-Jan-2013 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! 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. ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- INTEGER :: ARRSIZE REAL :: ARRAY(ARRSIZE) REAL :: VAL ! ! Local variables ! ---------------------------------------------------------------- REAL :: FIELD INTEGER :: I ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / IF (ARRSIZE.GT.0) THEN 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 ELSE LENGTH = 0 END IF RETURN END FUNCTION LENGTH !/ End of LENGTH ----------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / INTEGER FUNCTION FINDFIRST(ARRAY,ARRSIZE,VAL) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 4-Jan-2013 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! Fast algorithm to find the *first* index IND in ARRAY ! for which ARRAY(IND) = VAL. Use only when there are ! no duplicates in ARRAY! ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- INTEGER :: ARRSIZE REAL :: ARRAY(ARRSIZE) REAL :: VAL ! ! Local variables ! ---------------------------------------------------------------- INTEGER :: IND ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / IND = 1 DO WHILE (IND.LE.ARRSIZE) IF ( ARRAY(IND).EQ.VAL ) EXIT IND = IND + 1 END DO IF (IND.GT.ARRSIZE) THEN FINDFIRST = 0 ELSE FINDFIRST = IND ENDIF RETURN END FUNCTION FINDFIRST !/ End of FINDFIRST -------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / REAL FUNCTION STD(ARRAY,N) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 4-Jan-2013 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! Computes standard deviation. ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! 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 ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / 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 FUNCTION STD !/ End of STD -------------------------------------------------------- / !/ RECURSIVE SUBROUTINE QSORT(ARRAY,IDX,LO,HI) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | Stefan Zieger | !/ | FORTRAN 95 | !/ | Last update : 6-Sep-2016 | !/ +-----------------------------------+ !/ !/ 06-Sep-2016 : Origination, based on code by Mirko ( version 5.16 ) !/ Velic (BoM, Australia) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ ! ! 1. Purpose : ! ! Quicksort algorithm. ! ! 2. Method ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ARRAY REAL ARR in/out Input array ! IDX REAL ARR in/out Original indices of input array ! LO INTEGER input First element ! HI INTEGER input Last element ! IMPLICIT NONE !/ INTEGER, INTENT(IN) :: LO,HI REAL,INTENT(INOUT) :: ARRAY(:),IDX(:) !/ ! Local variables ! ---------------------------------------------------------------- LOGICAL :: LOOP INTEGER :: TOP, BOT REAL :: VAL, TMP ! ! 4. Subroutines used : ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! ! 5. Subroutines calling ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / ! !/ --- Check array size and bounds. --- IF ( SIZE(ARRAY).EQ. 0 ) THEN WRITE(6,199) CALL ABORT ELSE IF ( SIZE(ARRAY).NE.SIZE(IDX) ) THEN WRITE(6,201) CALL ABORT ELSE IF ( LBOUND(ARRAY,1).GT.LO ) THEN WRITE(6,203) CALL ABORT ELSE IF ( UBOUND(ARRAY,1).LT.HI ) THEN WRITE(6,205) CALL ABORT END IF ! TOP = LO BOT = HI VAL = ARRAY(INT((LO+HI)/2)) ! LOOP = .TRUE. DO WHILE ( LOOP ) DO WHILE ( ARRAY(TOP).LT.VAL ) TOP = TOP + 1 END DO DO WHILE ( VAL.LT.ARRAY(BOT) ) BOT = BOT - 1 END DO IF ( TOP.LT.BOT ) THEN !/ --- Swap values at indices TOP and BOT --- TMP = ARRAY(TOP) ARRAY(TOP) = ARRAY(BOT) ARRAY(BOT) = TMP !/ --- Swap index values at indices TOP and BOT --- TMP = IDX(TOP) IDX(TOP) = IDX(BOT) IDX(BOT) = TMP ! TOP = TOP + 1 BOT = BOT - 1 ELSE LOOP = .FALSE. END IF END DO !/ --- Recursive call quicksort --- IF (LO.LT.TOP-1) CALL QSORT(ARRAY,IDX,LO,TOP-1) IF (BOT+1.LT.HI) CALL QSORT(ARRAY,IDX,BOT+1,HI) ! RETURN !/ 199 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/ & ' QSORT ARRAY IS EMPTY' ) 201 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/ & ' QSORT ARRAY SIZE AND INDEX ARRAY SIZE MISMATCH' ) 203 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/ & ' QSORT ARRAY INDEX OUT OF LOWER BOUND' ) 205 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/ & ' QSORT ARRAY INDEX OUT OF UPPER BOUND' ) !/ END SUBROUTINE QSORT !/ ------------------------------------------------------------------- / !/ RECURSIVE SUBROUTINE QSORT_DESC(ARRAY,IDX,LO,HI) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | Stefan Zieger | !/ | FORTRAN 95 | !/ | Last update : 6-Sep-2016 | !/ +-----------------------------------+ !/ !/ 06-Sep-2016 : Origination, based on code by Mirko ( version 5.16 ) !/ Velic (BoM, AUstralia) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ ! ! 1. Purpose : ! ! Quicksort algorithm with descending sort order. ! ! 2. Method ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ARRAY REAL ARR in/out Input array ! LO INTEGER input First element ! HI INTEGER input Last element ! IMPLICIT NONE !/ INTEGER, INTENT(IN) :: LO,HI REAL,INTENT(INOUT) :: ARRAY(:),IDX(:) !/ ! Local variables ! ---------------------------------------------------------------- INTEGER :: TOP, BOT, I REAL :: VAL, TMP LOGICAL :: LOOP ! ! 4. Subroutines used : ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! ! 5. Subroutines calling ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / ! !/ --- Check array size and bounds. --- IF ( SIZE(ARRAY).EQ. 0 ) THEN WRITE(6,199) CALL ABORT ELSE IF ( SIZE(ARRAY).NE.SIZE(IDX) ) THEN WRITE(6,201) CALL ABORT ELSE IF ( LBOUND(ARRAY,1).GT.LO ) THEN WRITE(6,203) CALL ABORT ELSE IF ( UBOUND(ARRAY,1).LT.HI ) THEN WRITE(6,205) CALL ABORT END IF ! TOP = LO BOT = HI VAL = ARRAY(INT((LO+HI)/2)) ! LOOP = .TRUE. DO WHILE ( LOOP ) DO WHILE ( ARRAY(TOP).GT.VAL ) TOP = TOP + 1 END DO DO WHILE ( VAL.GT.ARRAY(BOT) ) BOT = BOT - 1 END DO IF ( TOP.LT.BOT ) THEN !/ --- Swap values at indices TOP and BOT --- TMP = ARRAY(TOP) ARRAY(TOP) = ARRAY(BOT) ARRAY(BOT) = TMP !/ --- Swap index values at indices TOP and BOT --- TMP = IDX(TOP) IDX(TOP) = IDX(BOT) IDX(BOT) = TMP ! TOP = TOP + 1 BOT = BOT - 1 ELSE LOOP = .FALSE. END IF END DO !/ --- Recursive call quicksort --- IF (LO.LT.TOP-1) CALL QSORT_DESC(ARRAY,IDX,LO,TOP-1) IF (BOT+1.LT.HI) CALL QSORT_DESC(ARRAY,IDX,BOT+1,HI) ! RETURN !/ 199 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/ & ' QSORT ARRAY IS EMPTY' ) 201 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/ & ' QSORT ARRAY SIZE AND INDEX ARRAY SIZE MISMATCH' ) 203 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/ & ' QSORT ARRAY INDEX OUT OF LOWER BOUND' ) 205 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/ & ' QSORT ARRAY INDEX OUT OF UPPER BOUND' ) !/ END SUBROUTINE QSORT_DESC !/ ------------------------------------------------------------------- / !/ FUNCTION SWAPI4(INT4) RESULT(INT4SWP) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | S. Zieger | !/ | FORTRAN 90 | !/ | Last update : 03-Jan-2017 | !/ +-----------------------------------+ !/ !/ 03-Jan-2017 : Origination ( version 5.16 ) !/ (S. Zieger) !/ ! 1. Purpose : ! ! Return a Byte-swapped integer (size of 4 bytes) ! ! 2. Source code : ! !/ ------------------------------------------------------------------- / IMPLICIT NONE INTEGER(KIND=4), INTENT(IN) :: INT4 INTEGER(KIND=4) :: INT4SWP !/ ! Local variables ! ---------------------------------------------------------------- INTEGER(KIND=1), DIMENSION(4) :: BYTEIN, BYTEOUT !/ BYTEIN = TRANSFER(INT4, BYTEIN) BYTEOUT = (/BYTEIN(4),BYTEIN(3),BYTEIN(2),BYTEIN(1)/) INT4SWP = TRANSFER(BYTEOUT, INT4SWP) !/ RETURN !/ END FUNCTION SWAPI4 !/ ------------------------------------------------------------------- / SUBROUTINE findIJV4 (a ,b ,maxI, maxJ ,indA ,indB ) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. J. van der Westhuysen | !/ | Jeff Hanson | !/ | Eve-Marie Devaliere | !/ | FORTRAN 95 | !/ | Last update : 03-Mar-2017 | !/ +-----------------------------------+ !/ !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 ) !/ by Jeff Hanson & Eve-Marie Devaliere !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 ) !/ 03-Mar-2017 : Calls to INTERSECT and UNION ( version 5.16 ) !/ replaced (S. Zieger, BoM, Australia) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ IMPLICIT NONE ! ! 1. Purpose : ! ! Find a(i,j) indices of system "a" that lie over or along the ! fringes of system "b". ! ! 2. Method ! ! (i) Use an index matrix to map locations of wave systems in B ! (ii) Avoid multiple use of INTERSECT and UNION as in findIJV3 ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! a, b Type(system) input Final set of tracked systems, for one time level ! maxI Int input Number rows indices of wave field ! maxJ Int input Number column indices of wave field ! indA*, indB* Int.A. output Pointer array of indices for combining systems ! TYPE(system) :: a, b INTEGER :: maxI, maxJ INTEGER, POINTER :: indA(:), indB(:) ! INTENT (IN) a, b, maxI,maxJ ! ! Local variables ! ---------------------------------------------------------------- ! posB Int Neighbour index ! posB_MM Int Neighbour index (-1,-1) ! posB_MP Int Neighbour index (-1,+1) ! posB_PM Int Neighbour index (+1,-1) ! posB_PP Int Neighbour index (+1,+1) ! tmpA*, tmpB* Int.A. Array of indices for combining ! systems ! INTEGER :: LENG_AI,LENG_BI INTEGER :: OUTA,OUTB,I,J,IND,OUTDUMB INTEGER :: POSB,POSB_MM,POSB_PM,POSB_MP,POSB_PP INTEGER :: IND_B2(maxI,maxJ) REAL,ALLOCATABLE :: TMPA(:),DUMA(:),TMPB(:) REAL,POINTER :: DUMB(:) LOGICAL :: FOUND ! ! 4. Subroutines used : ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! QSORT Subr. Private Quicksort algorithm ! UNIQUE Subr. Private Return sorted unique numbers of an array ! ! 5. Subroutines calling ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / ! NULLIFY(DUMB) ! IF (ASSOCIATED(INDA)) DEALLOCATE(INDA) IF (ASSOCIATED(INDB)) DEALLOCATE(INDB) ! LENG_AI = LENGTH(REAL(a%i),SIZE(a%i),REAL(9999)) LENG_BI = LENGTH(REAL(b%i),SIZE(b%i),REAL(9999)) ! ALLOCATE(TMPA(LENG_AI)) ALLOCATE(TMPB(5*LENG_AI)) ! TMPA(:) = 9999. TMPB(:) = 9999. ! OUTA = 0 OUTB = 0 IND_B2(:,:) = 0 ! DO IND=1,LENG_BI I = B%I(IND) J = B%J(IND) IF (IND_B2(I,J).EQ.0) IND_B2(I,J) = IND END DO ! DO IND=1,LENG_AI I = A%I(IND) J = A%J(IND) POSB = IND_B2(I,J) POSB_MM = 0 POSB_PP = 0 POSB_MP = 0 POSB_PM = 0 IF (I.GT.1.AND.J.GT.1) POSB_MM = IND_B2(i-1,j-1) IF (I.GT.1.AND.J.LT.MAXJ) POSB_MP = IND_B2(i-1,j+1) IF (I.LT.MAXI.AND.J.LT.MAXJ) POSB_PP = IND_B2(i+1,j+1) IF (I.LT.MAXI.AND.J.GT.1) POSB_PM = IND_B2(i+1,j-1) FOUND = .FALSE. IF (POSB.NE.0) THEN OUTB = OUTB + 1 TMPB(OUTB) = REAL(POSB) IF (.NOT.FOUND) THEN OUTA = OUTA + 1 TMPA(OUTA) = REAL(IND) FOUND = .TRUE. END IF END IF IF (POSB_MM.NE.0) THEN OUTB = OUTB + 1 TMPB(OUTB) = REAL(POSB_MM) IF (.NOT.FOUND) THEN OUTA = OUTA + 1 TMPA(OUTA) = REAL(IND) FOUND = .TRUE. END IF END IF IF (POSB_MP.NE.0) THEN OUTB = OUTB + 1 TMPB(OUTB) = REAL(POSB_MP) IF (.NOT.FOUND) THEN OUTA = OUTA + 1 TMPA(OUTA) = REAL(IND) FOUND = .TRUE. END IF END IF IF (POSB_PM.NE.0) THEN OUTB = OUTB + 1 TMPB(OUTB) = REAL(POSB_PM) IF (.NOT.FOUND) THEN OUTA = OUTA + 1 TMPA(OUTA) = REAL(IND) FOUND = .TRUE. END IF END IF IF (POSB_PP.NE.0) THEN OUTB = OUTB + 1 TMPB(OUTB) = REAL(POSB_PP) IF (.NOT.FOUND) THEN OUTA = OUTA + 1 TMPA(OUTA) = REAL(IND) FOUND = .TRUE. END IF END IF END DO ! !/ Compact indices for wave systems in B. !/ Check for empty arrays first. IF (OUTB.GT.0) THEN CALL UNIQUE(TMPB,OUTB,DUMB,OUTDUMB) OUTB = OUTDUMB END IF ALLOCATE(INDB(OUTB)) IF (OUTB.GT.0) INDB(1:OUTB) = INT(DUMB(1:OUTB)) IF (ASSOCIATED(DUMB)) DEALLOCATE(DUMB) ! !/ Allocate output array and transfer content !/ for wave systems in A. ALLOCATE(INDA(OUTA)) IF (OUTA.GT.0) THEN ALLOCATE(DUMA(OUTA)) DUMA(:) = 0 CALL QSORT(TMPA(1:OUTA),DUMA(1:OUTA),1,OUTA) IF (ALLOCATED(DUMA)) DEALLOCATE(DUMA) INDA(1:OUTA) = INT(TMPA(1:OUTA)) END IF !/ IF (ALLOCATED(TMPA)) DEALLOCATE(TMPA) IF (ALLOCATED(TMPB)) DEALLOCATE(TMPB) !/ RETURN !/ END SUBROUTINE findIJV4 !/ End of findIJV4 --------------------------------------------------- / !/ !/ ------------------------------------------------------------------- / !/ END MODULE W3STRKMD !/ !/ End of module W3STRKMD -------------------------------------------- / !/