PROGRAM cluster2 !------------------------------------------------------------------------------- ! Program to cluster HYSPLIT trajectory files using a faster method ! This code, CLUSTER2, is distributed under the GNU General Public ! license https://www.gnu.org/licenses/gpl-3.0.en.html !------------------------------------------------------------------------------- ! Compilation requires the conformal mapping conversion library: MAP="libcmapf.a" ! www.arl.noaa.gov/wp_arl/wp-content/uploads/utilities/cmap/cmapf.v1_0.tar.gz ! OR ! www.iamg.org/documents/oldftp/VOL23/v23-1-5.tar.Z !------------------------------------------------------------------------------- ! Use the following statements to compile the code with gfortran: ! gfortran -ocluster2 -ffree-form cluster2.f !------------------------------------------------------------------------------- ! Author: roland.draxler@meteozone.com ! 27 Jun 2025 - Initial version !------------------------------------------------------------------------------- IMPLICIT none INTEGER :: narg ! number of command line arguments INTEGER :: iargc ! command line argument function INTEGER :: numt ! number of trajectories input INTEGER :: cnum ! cluster number to output INTEGER :: ngpt=100 ! number of trajectory grid pts INTEGER :: maxc=360 ! maximum number of clusters CHARACTER(256) :: label ! label field of data record CHARACTER(256) :: tfile ! trajectory file name CHARACTER(256) :: finp, fout ! input and output file names LOGICAL :: diag ! when true diagnostic output on LOGICAL :: ftest ! file exist test result INTEGER :: i,j,k ! index values INTEGER :: kret,kpts ! return flag, points along traj REAL :: lat1,lon1 ! minimum lat-lon REAL :: lat2,lon2 ! maximum lat-lon REAL :: olat,olon ! trajectory origin location REAL :: dlat ! latitude difference (max-origin) REAL :: tlat,tlon ! trajectory lat-lon REAL :: x0,y0 ! lower left corner position REAL :: xp,yp ! trajectory position grid units REAL :: xmax,ymax ! maximum y-position from origin REAL :: gsize ! conformal grid resolution (km) REAL :: parmap(9) ! conformal map parameters INTEGER(2), ALLOCATABLE :: yy(:),mm(:),dd(:),hh(:) CHARACTER(80), ALLOCATABLE :: tname(:) ! trajectory file names REAL, ALLOCATABLE :: tpti(:,:),tptj(:,:) ! traj points in grid units INTEGER, ALLOCATABLE :: clus(:) ! cluster assignment REAL, ALLOCATABLE :: cpti(:,:),cptj(:,:) ! clus points in grid units INTEGER, ALLOCATABLE :: numc(:) ! number of traj per cluster !-------------------- INTERFACE SUBROUTINE TRAJREAD(yr,mo,da,hr,kpts,olat,olon,lat1,lat2,lon1,lon2) IMPLICIT none INTEGER(2), INTENT(OUT) :: yr,mo,da,hr INTEGER, INTENT(INOUT) :: kpts REAL, INTENT(OUT) :: olat,olon REAL, INTENT(INOUT) :: lat1,lat2,lon1,lon2 END SUBROUTINE trajread SUBROUTINE MKGRID(ngpt,gsize,clat,clon,lat1,lon1,parmap) IMPLICIT none INTEGER, INTENT(IN) :: ngpt ! number of grid points REAL, INTENT(IN) :: gsize ! grid spacing in km REAL, INTENT(IN) :: clat,clon ! center of the grid REAL, INTENT(IN) :: lat1,lon1 ! lower left corner at grid (1,1) REAL, INTENT(OUT) :: parmap(9) ! conformal map (cmapf) parameters END SUBROUTINE mkgrid SUBROUTINE TRAJGRID(kt,parmap,tpti,tptj) IMPLICIT none INTEGER, INTENT(IN) :: kt ! trajectory number REAL, INTENT(IN) :: parmap(9) ! mapping parameters REAL, INTENT(INOUT) :: tpti(:,:),tptj(:,:) ! trajectory array END SUBROUTINE trajgrid SUBROUTINE SETCLUS (diag,parmap,xp,yp,ymax,kpts,numt,clus, & tpti,tptj,numc,maxc,cpti,cptj) IMPLICIT NONE LOGICAL, INTENT(IN) :: diag ! verbose output flag REAL, INTENT(IN) :: parmap(9) ! conformal map parameters REAL, INTENT(IN) :: xp,yp ! trajectory origin point REAL, INTENT(IN) :: ymax ! maximum y grid point INTEGER, INTENT(IN) :: kpts ! max points along the traj INTEGER, INTENT(IN) :: numt ! number of trajectories INTEGER, INTENT(INOUT) :: clus(:) ! cluster assign to each traj REAL, INTENT(IN) :: tpti(:,:),tptj(:,:) ! trajectory position array INTEGER, INTENT(INOUT) :: numc(:) ! number of traj per cluster INTEGER, INTENT(IN) :: maxc ! maximum number of clusters REAL, INTENT(INOUT) :: cpti(:,:),cptj(:,:) ! cluster position array END SUBROUTINE setclus SUBROUTINE MRGCLUS (diag,cnum,kpts,numt,clus,numc,maxc, & cpti,cptj,tpti,tptj) IMPLICIT NONE LOGICAL, INTENT(IN) :: diag ! verbose output flag INTEGER, INTENT(IN) :: cnum ! final desired cluster number INTEGER, INTENT(IN) :: kpts ! number of points along traj INTEGER, INTENT(IN) :: numt ! number of input traj INTEGER, INTENT(INOUT) :: clus(:) ! cluster assignment per traj INTEGER, INTENT(INOUT) :: numc(:) ! number of traj per cluster INTEGER, INTENT(IN) :: maxc ! max number of initial clusters REAL, INTENT(INOUT) :: cpti(:,:),cptj(:,:) ! cluster mean trajectory REAL, INTENT(IN) :: tpti(:,:),tptj(:,:) ! individual trajectories END SUBROUTINE mrgclus SUBROUTINE PRTCLUS (numt,maxc,clus,numc,yy,mm,dd,hh,tname) IMPLICIT NONE INTEGER, INTENT(IN) :: numt ! number of trajectories input INTEGER, INTENT(IN) :: maxc ! maximum number of clusters INTEGER, INTENT(IN) :: clus(:) ! cluster assignment per traj INTEGER, INTENT(IN) :: numc(:) ! number of traj per cluster INTEGER(2), INTENT(IN) :: yy(:),mm(:),dd(:),hh(:) CHARACTER(80), INTENT(IN) :: tname(:) ! trajectory file names END SUBROUTINE prtclus SUBROUTINE AVGCLUS (parmap,maxc,kpts,numc,cpti,cptj) IMPLICIT NONE REAL, INTENT(IN) :: parmap(9) INTEGER, INTENT(IN) :: maxc ! maximum number of clusters INTEGER, INTENT(IN) :: kpts ! number of points along traj INTEGER, INTENT(IN) :: numc(:) ! number of traj per cluster REAL, INTENT(IN) :: cpti(:,:),cptj(:,:) END SUBROUTINE avgclus END INTERFACE !--------------------- narg=IARGC() IF(NARG.EQ.0)THEN WRITE(*,*)'USAGE: cluster2 -[options (default)]' WRITE(*,*)' -h[Help information]' WRITE(*,*)' -i[Input file name containing trajectory file names (INFILE)]' WRITE(*,*)' -o[Output file name of the cluster results (CLUSLIST)]' WRITE(*,*)' -n[Number of clusters to output (6)]' WRITE(*,*)' -v[Verbose output]' STOP END IF ! command line defaults diag = .false. finp = 'INFILE' fout = 'CLUSLIST' cnum = 6 ! go through each argument DO WHILE (narg.GT.0) CALL GETARG(NARG,LABEL) SELECT CASE (LABEL(1:2)) CASE ('-h','-H') WRITE(*,*)'Program to cluster HYSPLIT formatted trajectory output files.' WRITE(*,*)'Each output file should contain only one trajectory. The' WRITE(*,*)'computation is similar to the clustering program distributed' WRITE(*,*)'with HYSPLIT using spatial variance. However, it uses a first-' WRITE(*,*)'guess solution to permit a faster convergence. The final number' WRITE(*,*)'of clusters desired should be selected on the command line (-n)' WRITE(*,*)'and the input files assigned to each cluster will be written' WRITE(*,*)'to a file named CLUSTER_{num}, where {num} equals the cluster' WRITE(*,*)'number. This file format is identical to the one produced by the' WRITE(*,*)'HYSPLIT cluster program, eliminating the need for clusplot and' WRITE(*,*)'cluslist. The remaining programs in the HYSPLIT distribution:' WRITE(*,*)'clusmem, trajmean, merglist, and trajplot, can be used to' WRITE(*,*)'create the final graphic. This should be scripted if possible.' WRITE(*,*)' ' STOP CASE ('-i','-I') finp=TRIM(label(3:)) CASE ('-o','-O') fout=TRIM(label(3:)) CASE ('-n','-N') k=LEN_TRIM(LABEL) READ(label(3:k),'(I8)')cnum CASE ('-v','-V') diag=.true. END SELECT narg=narg-1 END DO ! check arguments IF(finp.EQ.''.OR.fout.EQ.'')THEN WRITE(*,'(2A)')'Both input and output files need to be defined!' STOP END IF ! open input and output data files INQUIRE(FILE=finp,EXIST=ftest) IF(ftest)THEN IF(diag) WRITE(*,'(2A)')'Started processing: ',TRIM(finp) OPEN(10,FILE=finp) ELSE WRITE(*,'(2A)')'File not found:',finp STOP END IF ! count the number of input files to process numt=0 kret=0 DO WHILE (kret.EQ.0) READ(10,'(A)',IOSTAT=kret) numt=numt+1 END DO numt=numt-1 IF(diag)THEN WRITE(*,'(A,I5)')'Number of trajectories: ',numt END IF ALLOCATE (tname(numt)) ALLOCATE (yy(numt),mm(numt),dd(numt),hh(numt)) REWIND(10) ! read data for each trajectory file ! determine domain and resolution lat1=90.0 lon1=180.0 lat2=-90.0 lon2=-180.0 kpts=0 DO k=1,numt READ(10,'(A)') tfile tname(k)=TRIM(tfile) OPEN(20,FILE=tfile) ! read the trajectory file CALL TRAJREAD(yy(k),mm(k),dd(k),hh(k), & kpts,olat,olon,lat1,lat2,lon1,lon2) CLOSE(20) END DO kpts=kpts-1 ! first end point is origin point -> skip lat1=lat1-1.0 lon1=lon1-1.0 lat2=lat2+1.0 lon2=lon2+1.0 IF(diag)THEN WRITE(*,'(A,I5)')'Trajectory duration: ',kpts WRITE(*,'(A,2F6.1)')'Trajectory origin: ',olat,olon WRITE(*,'(A,4F10.3)')'Lat,Lon Limits: ',lat1,lat2,lon1,lon2 END IF REWIND(10) ! create trajectory grid mapping parameters ALLOCATE (tpti(numt,kpts),tptj(numt,kpts)) ALLOCATE (clus(numt)) IF(diag)THEN WRITE(*,'(A,2I6)')'Allocated traj endpoint array (numt,kpts): ',numt,kpts END IF ! create cluster array ALLOCATE (cpti(maxc,kpts),cptj(maxc,kpts)) ALLOCATE (numc(maxc)) IF(diag)THEN WRITE(*,'(A,2I6)')'Allocated cluster array (maxc,kpts): ',maxc,kpts END IF ! define the grid to map trajectory lat-lon to conformal x-y dlat=111.0*(lat2-lat1+1) ! vertical extent of domain in km gsize=dlat/ngpt ! average grid size CALL MKGRID(ngpt,gsize,olat,olon,lat1,lon1,parmap) CALL CXY2LL(parmap,1.0,1.0,tlat,tlon) CALL CLL2XY(parmap,lat2,lon2,x0,y0) CALL CLL2XY(parmap,olat,olon,xp,yp) CALL CLL2XY(parmap,lat2,olon,xmax,ymax) IF(diag)THEN WRITE(*,'(A,I5)')'Number of grid points: ',ngpt WRITE(*,'(A,F10.3)')'Grid size (km): ',gsize WRITE(*,'(A,2F8.1)')'Lower left corner (lat,lon): ',tlat,tlon WRITE(*,'(A,2F8.1)')'Upper right corner (x,y): ',x0,y0 WRITE(*,'(A,2F8.1)')'Origin point (i,j): ',xp,yp WRITE(*,'(A,2F8.1)')'Scale point (i,j): ',xmax,ymax END IF ! load trajectory data into array and convert to grid units DO k=1,numt READ(10,'(A)') tfile OPEN(20,FILE=tfile) CALL TRAJGRID(k,parmap,tpti,tptj) CLOSE(20) END DO CLOSE(10) IF(diag)THEN WRITE(*,'(A,6F6.1)')'Trajectories converted to grid units' END IF ! fill cluster array with the first guess and assign trajectories ! to their MIN(spv) initial cluster CALL SETCLUS (diag,parmap,xp,yp,ymax,kpts,numt,clus, & tpti,tptj,numc,maxc,cpti,cptj) IF(diag)THEN WRITE(*,'(A)')'Number of trajectories in each initial cluster' WRITE(*,'(A)')'Cluster number is in degrees from north ...' DO k=1,90 WRITE(*,'(2I6,3(10X,2I6))') k,numc(k), k+90,numc(k+90), & k+180,numc(k+180), k+270,numc(k+270) END DO END IF ! merge clusters together CALL MRGCLUS (diag,cnum,kpts,numt,clus,numc,maxc,cpti,cptj,tpti,tptj) ! print cluster results WRITE(LABEL,'(A,A1,I1)') TRIM(fout),'_',cnum OPEN(30,file=LABEL) CALL PRTCLUS (numt,maxc,clus,numc,yy,mm,dd,hh,tname) CLOSE(30) ! output cluster mean trajectories to individual files IF(diag) CALL AVGCLUS (parmap,maxc,kpts,numc,cpti,cptj) END PROGRAM cluster2 !------------------------------------------------------------------------------- ! reads a simple trajectory file with a single source location ! for all trajectories listed in file INFILE ! to determine how much array space to allocate SUBROUTINE TRAJREAD(yr,mo,da,hr,kpts,olat,olon,lat1,lat2,lon1,lon2) IMPLICIT none INTEGER(2),INTENT(OUT) :: yr,mo,da,hr ! traj start time INTEGER, INTENT(INOUT) :: kpts ! traj duration REAL, INTENT(OUT) :: olat,olon ! traj origin point REAL, INTENT(INOUT) :: lat1,lat2,lon1,lon2 ! min-max lat and lon ! trajectory file variables; not all used in this application CHARACTER(256) :: label CHARACTER(8) :: model, dirctn, motion, ldiag(2) INTEGER :: k, knum, kmet, kfmt, kret, kdiag, npts INTEGER :: fh,tn,mg,tyr,tmo,tda,thr,tmn,tfh REAL :: tlat,tlon,olvl REAL :: age,thgt,tval(2) ! number of meteorology files and file format (none/0 = old; 1=new) READ(20,'(2I6)') kmet,kfmt DO k=1,kmet READ(20,'(A8,5I6)') model,yr,mo,da,hr,fh END DO ! number of trajectories, direction, and vertical motion READ(20,'(I6,2(1X,A8))')knum,dirctn,motion IF(knum.GT.1)THEN WRITE(*,*)'Only one trajectory per file permitted: ',knum STOP END IF DO k=1,knum READ(20,'(4I6,2F9.3,F8.1)') yr,mo,da,hr,olat,olon,olvl END DO ! number of diagnostic variables READ(20,'(I6,10(1X,A8))')kdiag,ldiag IF(kdiag.GT.2)THEN WRITE(*,*)'No more than 2 diagnostic variables: ',kdiag STOP END IF kret=0 npts=0 DO WHILE (kret.EQ.0) ! traj#,meteo#,tyr,tmo,tda,thr,tmn,tfh,age READ(20,'(A)',IOSTAT=kret) label IF(kret.EQ.0)THEN npts=npts+1 READ(label,'(8I6,F8.1,2F9.3,5(1X,F8.1))') & tn,mg,tyr,tmo,tda,thr,tmn,tfh,age, & tlat,tlon,thgt,tval(1),tval(2) ! determine the lat-lon range over all trajectories lat1=MIN(lat1,tlat) lon1=MIN(lon1,tlon) lat2=MAX(lat2,tlat) lon2=MAX(lon2,tlon) END IF END DO ! determine the maximum number of points along the trajectory kpts=MAX(npts,kpts) END SUBROUTINE trajread !------------------------------------------------------------------------------- ! reads a simple trajectory file with a single source location ! and updates the trajectory counter grid, see TRAJREAD for more ! detail about the format of the input trajectory file; not all ! data are used in the clustering application SUBROUTINE TRAJGRID(kt,parmap,tpti,tptj) IMPLICIT none INTEGER, INTENT(IN) :: kt ! trajectory number REAL, INTENT(IN) :: parmap(9) ! mapping parameters REAL, INTENT(INOUT) :: tpti(:,:),tptj(:,:) ! trajectory array CHARACTER(256) :: label CHARACTER(8) :: model, dirctn, motion, ldiag(2) INTEGER :: k, kd, knum, kmet, kfmt, kret, kdiag INTEGER :: yr,mo,da,hr,fh REAL :: olat,olon,olvl INTEGER :: tn,mg,tyr,tmo,tda,thr,tmn,tfh REAL :: tlat,tlon REAL :: age,thgt,tval(2) ! number of meteorology files READ(20,'(2I6)') kmet,kfmt DO k=1,kmet READ(20,'(A8,5I6)') model,yr,mo,da,hr,fh END DO ! number of trajectories, direction, and vertical motion READ(20,'(I6,2(1X,A8))')knum,dirctn,motion DO k=1,knum READ(20,'(4I6,2F9.3,F8.1)') yr,mo,da,hr,olat,olon,olvl END DO ! number of diagnostic variables READ(20,'(I6,10(1X,A8))')kdiag,ldiag kret=0 DO WHILE (kret.EQ.0) ! traj#,meteo#,tyr,tmo,tda,thr,tmn,tfh,age READ(20,'(A)',IOSTAT=kret) label IF(kret.EQ.0)THEN READ(label,'(8I6,F8.1,2F9.3,5(1X,F8.1))') & tn,mg,tyr,tmo,tda,thr,tmn,tfh,age, & tlat,tlon,thgt,tval(1),tval(2) kd=INT(ABS(age)) ! skip the first output line (hour=0) which is the source location IF(kd.GT.0)CALL CLL2XY(parmap,tlat,tlon,tpti(kt,kd),tptj(kt,kd)) END IF END DO END SUBROUTINE trajgrid !------------------------------------------------------------------------------- ! sets the values for the output grid, including the map transformation ! parameters (parmap) based on the latitude of the trajectory origin point SUBROUTINE MKGRID(ngpt,gsize,clat,clon,lat1,lon1,parmap) IMPLICIT none INTEGER, INTENT(IN) :: ngpt ! number of grid points REAL, INTENT(IN) :: gsize ! grid spacing in km REAL, INTENT(IN) :: clat,clon ! center of the grid REAL, INTENT(IN) :: lat1,lon1 ! lower left corner at grid (1,1) REAL, INTENT(OUT) :: parmap(9) ! conformal map (cmapf) parameters REAL :: grids(12) ! user readable map definition ! defines a polar sterographic projection IF(ABS(clat).GT.60.0)THEN grids(1)=SIGN(90.0,clat) ! set the pole position at +90 or -90 grids(2)=clon ! pole longtitude (+180 from cut) grids(3)=clat ! reference lat/lon (at which grid size specified) grids(4)=clon ! reference longitude and grid alignment grids(7)=grids(1) ! tangent latitude ! defines a mercator projection ELSEIF(ABS(clat).LT.30.0)THEN grids(1)=0.0 grids(2)=clon ! pole lat/lon axis through pole grids(3)=clat ! reference lat grids(4)=clon ! reference lon grids(7)=0.0 ! tangent latitude ! defines a lambert conformal projection ELSE grids(1)=clat grids(2)=clon ! pole lat/lon axis through pole grids(3)=clat ! reference lat grids(4)=clon ! reference lon grids(7)=clat ! tangent latitude END IF grids(6)=0.0 ! grid orientation grids(5)=GSIZE ! delta=x grid size in km grids(8)=1.0+ngpt/2.0 ! synch point in x,y coordinates grids(9)=1.0+ngpt/2.0 ! ... grids(10)=clat ! synch point in lat/lon coordinates grids(11)=clon ! ... grids(12)=0.0 ! variable reserved for future use ! define the tangent latitude and reference longitude CALL STLMBR(parmap,grids(7),grids(4)) ! define the grid by a one-point specification CALL STCM1P(parmap,grids(8),grids(9),grids(10),grids(11), & grids(3),grids(4),grids(5),grids(6)) END SUBROUTINE mkgrid !------------------------------------------------------------------------------- ! creates the first guess clusters and assigns trajectories accordingly SUBROUTINE SETCLUS (diag,parmap,xp0,yp0,ymax,kpts,numt,clus, & tpti,tptj,numc,maxc,cpti,cptj) IMPLICIT NONE LOGICAL, INTENT(IN) :: diag ! verbose output flag REAL, INTENT(IN) :: parmap(9) ! conformal map parameters REAL, INTENT(IN) :: xp0,yp0 ! trajectory origin point REAL, INTENT(IN) :: ymax ! maximum y grid point INTEGER, INTENT(IN) :: kpts ! max points along the traj INTEGER, INTENT(IN) :: numt ! number of trajectories INTEGER, INTENT(INOUT) :: clus(:) ! cluster assign to each traj REAL, INTENT(IN) :: tpti(:,:),tptj(:,:) ! trajectory position array INTEGER, INTENT(INOUT) :: numc(:) ! number of traj per cluster INTEGER, INTENT(IN) :: maxc ! maximum number of clusters REAL, INTENT(INOUT) :: cpti(:,:),cptj(:,:) ! cluster position array INTEGER :: i,j,k REAL :: ang REAL :: tlen,tlen0 REAL :: xm,ym REAL :: frac,dist,vmin REAL, ALLOCATABLE :: xp(:),yp(:) ALLOCATE (xp(kpts),yp(kpts)) ! create first guess trajectories with maximum length of tlen0 ! in trajectory grid units tlen0=ymax-yp0+1.0 i=0 ! compute a first-guess trajectory (cluster) every degree ! resulting in a total of 360 possible clusters DO k=1,360 i=i+1 ! every other first-guess trajectory is half length tlen=tlen0 IF(MOD(i,2).EQ.0) tlen=tlen0/2.0 ! compute first-guess endpoint position ! in units relative to origin center ang=FLOAT(k)/57.3 ym=COS(ang)*tlen xm=SIN(ang)*tlen ! interpolate intermediate trajectory points DO j=1,kpts frac=FLOAT(j)/FLOAT(kpts) cpti(i,j)=xp0+xm*frac cptj(i,j)=yp0+ym*frac END DO ! number of trajectories within this cluster, the first guess ! does not count, because the trajectories have not yet been ! assigned to their nearest first-guess cluster numc(i)=0 END DO IF(i.NE.maxc)THEN WRITE(*,'(A)')'Number of clusters inconsistent: ',i,maxc STOP END IF ! assign each trajectory to a first-guess cluster DO k=1,numt ! for each trajectory vmin=1.0E+10 DO j=1,maxc ! check each cluster dist=0.0 DO i=1,kpts ! variance along trajectory dist=dist+(tpti(k,i)-cpti(j,i))**2+(tptj(k,i)-cptj(j,i))**2 END DO ! save the trajectory with the smallest spatial variance IF(dist.LT.vmin)THEN vmin=dist clus(k)=j ! assign trajectory k to cluster j END IF END DO j=clus(k) ! final cluster assignment for k trajectory numc(j)=numc(j)+1 ! increment trajectory count for cluster j END DO IF(diag) & WRITE(*,'(I5,A,I3,A)')SUM(numc),' trajectories initially assigned to ',maxc,' clusters' ! compute a cluster mean trajectory DO j=1,maxc ! check each cluster xp=0.0 yp=0.0 DO k=1,numt ! when k trajectory assigned to j cluster IF(clus(k).EQ.j)THEN ! sum positions along the trajectory DO i=1,kpts xp(i)=xp(i)+tpti(k,i) yp(i)=yp(i)+tptj(k,i) END DO END IF END DO ! compute average cpti(j,:)=xp(:)/numc(j) cptj(j,:)=yp(:)/numc(j) END DO DEALLOCATE (xp,yp) END SUBROUTINE setclus !---------------------------------------------------------------------- ! merge the clusters SUBROUTINE MRGCLUS (diag,cnum,kpts,numt,clus,numc,maxc, & cpti,cptj,tpti,tptj) IMPLICIT NONE LOGICAL, INTENT(IN) :: diag ! verbose output flag INTEGER, INTENT(IN) :: cnum ! final desired cluster number INTEGER, INTENT(IN) :: kpts ! number of points along traj INTEGER, INTENT(IN) :: numt ! number of input traj INTEGER, INTENT(INOUT) :: clus(:) ! cluster assignment per traj INTEGER, INTENT(INOUT) :: numc(:) ! number of traj per cluster INTEGER, INTENT(IN) :: maxc ! max number of initial clusters REAL, INTENT(INOUT) :: cpti(:,:),cptj(:,:) ! cluster mean trajectory REAL, INTENT(IN) :: tpti(:,:),tptj(:,:) ! individual trajectories INTEGER :: i,j,k INTEGER :: remc,jmin,npass REAL :: frac,dist,vmin REAL, ALLOCATABLE :: xp(:),yp(:) ALLOCATE (xp(kpts),yp(kpts)) IF(diag) WRITE(*,'(A,I2)')'Start merging clusters ... final goal: ',cnum npass=1 remc=maxc ! number of remaining clusters after merge DO WHILE (remc.GT.cnum) IF(diag) WRITE(*,'(A,I2)')'Pass: ',npass ! for each cluster check every other cluster to find ! the cluster with the minimum spatial variance DO i=1,maxc ! for each cluster IF(numc(i).GT.0)THEN ! that contains trajectories vmin=1.0E+10 jmin=0 ! check every other cluster that contains trajectories DO j=1,maxc IF(j.NE.i.AND.numc(j).GT.0)THEN dist=0.0 DO k=1,kpts ! sum the variance along trajectory dist=dist+(cpti(j,k)-cpti(i,k))**2+(cptj(j,k)-cptj(i,k))**2 END DO IF(dist.LT.vmin)THEN vmin=dist jmin=j END IF END IF END DO ! move trajectories in minimum cluster to i cluster IF(jmin.NE.0)THEN numc(i)=numc(i)+numc(jmin) numc(jmin)=0 ELSE ! this should never happen WRITE(*,'(A)')'Cluster merge failed to find minimum variance' STOP END IF ! reassign indiviudual trajectories to the new cluster DO k=1,numt IF(clus(k).EQ.jmin) clus(k)=i END DO ! compute a new cluster mean trajectory xp=0.0 yp=0.0 DO j=1,numt ! when j trajectory assigned to i cluster IF(clus(j).EQ.i)THEN ! sum positions along the trajectory DO k=1,kpts xp(k)=xp(k)+tpti(j,k) yp(k)=yp(k)+tptj(j,k) END DO END IF END DO ! compute average cpti(i,:)=xp(:)/numc(i) cptj(i,:)=yp(:)/numc(i) ! count how many clusters remain remc=0 DO k=1,maxc IF(numc(k).GT.0) remc=remc+1 END DO IF(diag)THEN WRITE(*,'(I5,A,I3,A)')SUM(numc),' trajectories assigned to ',remc,' clusters' WRITE(*,'(A)')'Number of trajectories in each cluster: ' DO k=1,90 WRITE(*,'(2I6,3(10X,2I6))') k,numc(k), k+90,numc(k+90), & k+180,numc(k+180), k+270,numc(k+270) END DO END IF IF(remc.EQ.cnum) THEN IF(diag) WRITE(*,'(A,I2,A)')'Cluster merging complete for ',remc,' clusters' DEALLOCATE (xp,yp) RETURN END IF END IF ! contains trajectories test END DO ! cluster loop ! another pass is required npass=npass+1 IF(diag) WRITE(*,'(2(A,I2))')'Cluster merging continues: ',remc,' < ',cnum END DO ! multiple pass loop END SUBROUTINE MRGCLUS !---------------------------------------------------------------------- ! print the clusters SUBROUTINE PRTCLUS (numt,maxc,clus,numc,yy,mm,dd,hh,tname) IMPLICIT NONE INTEGER, INTENT(IN) :: numt ! total number of trajectories INTEGER, INTENT(IN) :: maxc ! maximum number of clusters INTEGER, INTENT(IN) :: clus(:) ! cluster assignment per traj INTEGER, INTENT(IN) :: numc(:) ! number of traj per cluster INTEGER(2), INTENT(IN) :: yy(:),mm(:),dd(:),hh(:) CHARACTER(80), INTENT(IN) :: tname(:) ! trajectory file names INTEGER :: i,j,k j=0 ! check each cluster number DO k=1,maxc IF(numc(k).GT.0) THEN ! define a sequential cluster number because ! non-zero content k clusters may not be sequential j=j+1 ! find all trajectories associated with cluster k DO i=1,numt ! trajectory i is tagged with cluster k IF(clus(i).EQ.k)THEN WRITE(30,'(I5,I6,I5,3I3,I6,1X,A)') & j,numc(k),yy(i),mm(i),dd(i),hh(i),i,TRIM(tname(i)) END IF END DO END IF END DO END SUBROUTINE prtclus !---------------------------------------------------------------------- ! output the cluster mean trajectory to individal files but only ! when verbose output is turned on SUBROUTINE AVGCLUS (parmap,maxc,kpts,numc,cpti,cptj) IMPLICIT NONE REAL, INTENT(IN) :: parmap(9) INTEGER, INTENT(IN) :: maxc ! maximum number of clusters INTEGER, INTENT(IN) :: kpts ! number of points along traj INTEGER, INTENT(IN) :: numc(:) ! number of traj per cluster REAL, INTENT(IN) :: cpti(:,:),cptj(:,:) CHARACTER(80) :: label REAL :: tlat,tlon INTEGER :: i,j,k j=0 ! check each cluster number DO k=1,maxc IF(numc(k).GT.0) THEN ! define a sequential cluster number because ! non-zero content k clusters may not be sequential j=j+1 WRITE(label,'(A1,I1,A)')'C',j,'mean.txt' OPEN(40,FILE=label) ! convert all points along the cluster mean trajectory DO i=1,kpts CALL CXY2LL(parmap,cpti(k,i),cptj(k,i),tlat,tlon) WRITE(40,'(I2,2X,2F9.3)') i,tlat,tlon END DO CLOSE(40) END IF END DO END SUBROUTINE avgclus