WAVEWATCH III  beta 0.0.1
ww3_systrk.F90
Go to the documentation of this file.
1 
8 !
9 !/ ------------------------------------------------------------------- /
28 !
29 PROGRAM ww3_systrk
30  !/
31  !/ +-----------------------------------+
32  !/ | WAVEWATCH III NOAA/NCEP |
33  !/ | A. J. van der Westhuysen |
34  !/ | Jeff Hanson |
35  !/ | Eve-Marie Devaliere |
36  !/ | FORTRAN 95 |
37  !/ | Last update : 16-Jan-2017 |
38  !/ +-----------------------------------+
39  !/
40  !/ 03-Feb-2012 : Origination, based on Matlab code ( version 4.05 )
41  !/ by Jeff Hanson & Eve-Marie Devaliere
42  !/ 04-Jan-2013 : Inclusion in trunk ( version 4.08 )
43  !/ 29-Nov-2013 : Remove DOC control characters,
44  !/ update MPI! to MPI/! (H.L. Tolman). ( version 4.15 )
45  !/ 11-Feb-2014 : Add NetCDF output option. Both NetCDF-3 and
46  !/ NetCDF-4 are available. (B. Li). ( version 4.18 )
47  !/ 26-Sep-2016 : Optimization updates (A. van der Westhuysen)
48  !/ ( version 5.15 )
49  !/ 20-Sep-2016 : Add support for unformatted partition file.
50  !/ (S.Zieger BoM, Australia) ( version 5.16 )
51  !/ 20-Dec-2016 : Optimized search algorithms and
52  !/ set functions. (S.Zieger) ( version 5.16 )
53  !/
54  !/ Copyright 2009-2013 National Weather Service (NWS),
55  !/ National Oceanic and Atmospheric Administration. All rights
56  !/ reserved. WAVEWATCH III is a trademark of the NWS.
57  !/ No unauthorized use without permission.
58  !/
59  USE w3strkmd
60  USE w3timemd, ONLY: tdiff
61  IMPLICIT NONE
62 #ifdef W3_MPI
63 
64  include "mpif.h"
65 #endif
66  !
67  ! 1. Purpose :
68  !
69  ! Perform spatial and temporal tracking of wave systems, based
70  ! on spectral partition (bulletin) output.
71  !
72  ! 2. Method :
73  !
74  ! This is a controller program. It reads the input parameter file
75  ! ww3_systrk.inp and calls subroutine waveTracking_NWS_V2 to
76  ! perform the actual tracking procedure. Write output (fields and
77  ! point output).
78  !
79  ! 3. Parameters :
80  !
81  LOGICAL :: testout
82  parameter(testout = .false.)
83  CHARACTER :: filename*80, paramfile*32
84  REAL :: dirknob, perknob, hsknob, wetpts, seedlat, &
85  seedlon, dirtimeknob, tptimeknob, tint
86  REAL :: lonout(100), latout(100) !Increase dimension?
87  INTEGER :: maxgroup, ntint, noutp
88  INTEGER :: clkdt0(8),clkdt1(8)
89  REAL :: clkfel
90  TYPE(dat2d), POINTER :: wsdat(:)
91  TYPE(timsys), POINTER :: sysa(:)
92  INTEGER, POINTER :: maxsys(:)
93  !
94  ! Local parameters.
95  ! ----------------------------------------------------------------
96  ! intype Int input Type of input (0 = from memory; 1 = from file)
97  ! tmax Int input Value of maxTs to apply (1 or 2, used for model coupling)
98  ! tcur Int input Index of current time step (1 or 2, used for model coupling)
99  ! ulimGroup Int input Upper limit of number of wave systems to output
100  !
101  LOGICAL :: file_exists
102  CHARACTER :: inpstr*72
103  INTEGER :: intype, tmax, tcur, maxi, maxj
104  INTEGER :: it, igrp, sysmatch, ind, ip
105  INTEGER :: i, j, leng, ulimgroup
106  REAL, ALLOCATABLE :: dum(:,:)
107 #ifdef W3_TRKNC
108  REAL, ALLOCATABLE :: dum2nc(:,:,:,:)
109  REAL, ALLOCATABLE :: hsprt_nc(:,:,:)
110  REAL, ALLOCATABLE :: tpprt_nc(:,:,:)
111  REAL, ALLOCATABLE :: dirprt_nc(:,:,:)
112  REAL, ALLOCATABLE :: longitude_nc(:),latitude_nc(:)
113  REAL, ALLOCATABLE :: lonprt_nc(:),latprt_nc(:)
114 #endif
115  INTEGER ntime_nc
116  INTEGER :: outputtype
117  LOGICAL :: outputcheck1
118  DOUBLE PRECISION :: date1, date2, tstart, tend
119  REAL :: dlon, dlat, lonprt, latprt
120  REAL :: dt
121  REAL :: minlon, maxlon, minlat, maxlat
122  INTEGER :: mxcwt, mycwt
123 #ifdef W3_MPI
124  INTEGER :: rank, nproc, ierr
125  CHARACTER :: rankstr*4
126 #endif
127 
128  ! For point output (bilinear interpolation)
129  REAL :: hsprt(10),tpprt(10),dirprt(10)
130  REAL :: bl_hsprt(10),br_hsprt(10),tr_hsprt(10),tl_hsprt(10), &
131  bl_tpprt(10),br_tpprt(10),tr_tpprt(10),tl_tpprt(10), &
132  bl_dirprt(10),br_dirprt(10),tr_dirprt(10),tl_dirprt(10)
133  REAL :: bl_dirx,br_dirx,tr_dirx,tl_dirx, &
134  bl_diry,br_diry,tr_diry,tl_diry
135  REAL :: bl_lonprt,br_lonprt,tr_lonprt,tl_lonprt, &
136  bl_latprt,br_latprt,tr_latprt,tl_latprt
137  REAL :: t, u, bl_w, br_w, tr_w, tl_w
138  REAL :: pi
139  parameter(pi = 3.1416)
140  !
141  ! 4. Subroutines used :
142  !
143  ! waveTracking_NWS_V2
144  !
145  ! 5. Called by :
146  !
147  ! None, stand-alone program.
148  !
149  ! 6. Error messages :
150  !
151  ! 7. Remarks :
152  !
153  ! 8. Structure :
154  !
155  ! Calls subroutine waveTracking_NWS_V2 in trackmd.95 - see that
156  ! file for structure.
157  !
158  ! 9. Switches :
159  !
160  ! !/SHRD Switch for shared / distributed memory architecture.
161  ! !/MPI Id.
162  !
163  ! 10. Source code :
164  !
165  !/ ------------------------------------------------------------------- /
166  !
167 #ifdef W3_MPI
168  ! Start of parallel region
169  CALL mpi_init(ierr)
170 
171  CALL mpi_comm_rank(mpi_comm_world, rank, ierr)
172  CALL mpi_comm_size(mpi_comm_world, nproc, ierr)
173 
174 #endif
175  ! Open log file
176 #ifdef W3_MPI
177  WRITE(rankstr,'(i4.4)') rank
178  OPEN(unit=20,file='sys_log'//rankstr//'.ww3',status='unknown')
179 #endif
180 #ifdef W3_SHRD
181  OPEN(unit=20,file='sys_log.ww3',status='unknown')
182 #endif
183 
184  ! Print code version
185 #ifdef W3_MPI
186  IF (rank.EQ.0) THEN
187 #endif
188  WRITE(6,900)
189 #ifdef W3_MPI
190  END IF
191 #endif
192  WRITE(20,900)
193 900 FORMAT (/15x,' *** WAVEWATCH III Wave system tracking *** '/ &
194  15x,'==============================================='/)
195 
196  ! Since this program reads the raw partitioning input from file,
197  ! we set intype=1 or 2, and tmax and tcur to dummy values (not used).
198  intype = 2
199  ! intype = 1
200  IF (intype.EQ.1) WRITE(6,*) &
201  '*** WARNING: partRes format input used!'
202  tmax = 0
203  tcur = 0
204 
205  ! Read input parameter file
206 #ifdef W3_MPI
207  IF (rank.EQ.0) THEN
208 #endif
209  INQUIRE(file='ww3_systrk.inp', exist=file_exists)
210  IF (.NOT.file_exists) THEN
211  WRITE(20,2000)
212  WRITE(6,2000)
213  CALL abort
214  END IF
215  OPEN(unit=10,file='ww3_systrk.inp',status='old')
216 
217  READ(10,'(A72)') inpstr
218  DO WHILE (inpstr(1:1).EQ.'$')
219  READ(10,'(A72)') inpstr
220  END DO
221  backspace(10)
222  READ(10,*) filename
223 
224  READ(10,'(A72)') inpstr
225  DO WHILE (inpstr(1:1).EQ.'$')
226  READ(10,'(A72)') inpstr
227  END DO
228  backspace(10)
229  READ(10,*) date1, date2, dt, ntint
230  tstart = date1 + date2/1000000
231 
232  READ(10,'(A72)') inpstr
233  DO WHILE (inpstr(1:1).EQ.'$')
234  READ(10,'(A72)') inpstr
235  END DO
236  backspace(10)
237  READ(10,*) outputtype
238 
239  !Check for correct outputType option:
240  IF (outputtype.EQ.1) THEN
241  !ASCII output
242  ELSEIF (outputtype.EQ.3) THEN
243  !NetCDF 3 - requrires !/TRKNC switch
244  outputcheck1 = .true.
245 #ifdef W3_TRKNC
246  outputcheck1 = .false.
247 #endif
248  IF(outputcheck1) THEN
249  WRITE(6,993)
250  stop
251  END IF
252  ELSEIF (outputtype.EQ.4) THEN
253  !NetCDF 4 - requrires !/TRKNC switch
254  outputcheck1 = .true.
255 #ifdef W3_TRKNC
256  outputcheck1 = .false.
257 #endif
258  IF(outputcheck1) THEN
259  WRITE(6,994)
260  stop
261  END IF
262  ELSE
263  !Not a valid outputType
264  WRITE(6,995) outputtype
265  stop
266  ENDIF
267 
268  READ(10,'(A72)') inpstr
269  DO WHILE (inpstr(1:1).EQ.'$')
270  READ(10,'(A72)') inpstr
271  END DO
272  backspace(10)
273  READ(10,*) minlon, maxlon, mxcwt
274 
275  READ(10,'(A72)') inpstr
276  DO WHILE (inpstr(1:1).EQ.'$')
277  READ(10,'(A72)') inpstr
278  END DO
279  backspace(10)
280  READ(10,*) minlat, maxlat, mycwt
281 
282  READ(10,'(A72)') inpstr
283  DO WHILE (inpstr(1:1).EQ.'$')
284  READ(10,'(A72)') inpstr
285  END DO
286  backspace(10)
287  READ(10,*) dirknob, perknob, hsknob, wetpts, &
288  dirtimeknob, tptimeknob
289 
290  READ(10,'(A72)') inpstr
291  DO WHILE (inpstr(1:1).EQ.'$')
292  READ(10,'(A72)') inpstr
293  END DO
294  backspace(10)
295  READ(10,*) seedlat, seedlon
296 
297  READ(10,'(A72)') inpstr
298  DO WHILE (inpstr(1:1).EQ.'$')
299  READ(10,'(A72)') inpstr
300  END DO
301  backspace(10)
302  noutp = 1
303  lonout(:) = 9999.
304  latout(:) = 9999.
305  DO WHILE (.true.)
306  READ(10,*) lonout(noutp),latout(noutp)
307  IF ((lonout(noutp).EQ.0.).AND.(latout(noutp).EQ.0.)) EXIT
308  noutp = noutp + 1
309  END DO
310  noutp = noutp - 1
311 
312  CLOSE(10)
313 
314  WRITE(20,*) 'Raw partition file = ',filename
315  WRITE(20,'(A,F15.6)') 'Start time = ',tstart
316  WRITE(20,*) 'dt = ',dt
317  WRITE(20,*) 'No. time levels = ',ntint
318  WRITE(20,'(A,2F7.2)') 'Domain limits: Longitude =',minlon, maxlon
319  WRITE(20,'(A,2F7.2)') ' Latitude =',minlat, maxlat
320  WRITE(20,*) 'No. increments: Long, Lat =',mxcwt, mycwt
321  WRITE(20,*) 'dirKnob, perKnob, hsKnob, wetPts, &
322  dirTimeKnob, tpTimeKnob, seedLat, seedLon ='
323  WRITE(20,'(8F6.2)') dirknob, perknob, hsknob, wetpts, &
324  dirtimeknob, tptimeknob, seedlat, seedlon
325  WRITE(20,*) 'No. output points =',noutp
326  DO i = 1,noutp
327  WRITE(20,*) lonout(i), latout(i)
328  END DO
329 
330  INQUIRE(file=filename, exist=file_exists)
331  IF (.NOT.file_exists) THEN
332  WRITE(20,2200) filename
333  WRITE(6,2200) filename
334  CALL exit(1)
335  END IF
336 
337 
338  CALL date_and_time ( values=clkdt0 )
339 
340 #ifdef W3_MPI
341  END IF
342 #endif
343 
344 #ifdef W3_MPI
345  ! MPI communication block
346  CALL mpi_bcast(filename,80,mpi_character,0,mpi_comm_world,ierr)
347  CALL mpi_bcast(tstart,1,mpi_double_precision,0, &
348  mpi_comm_world,ierr)
349  CALL mpi_bcast(tend,1,mpi_double_precision,0, &
350  mpi_comm_world,ierr)
351  CALL mpi_bcast(dt,1,mpi_real,0,mpi_comm_world,ierr)
352  CALL mpi_bcast(ntint,1,mpi_integer,0,mpi_comm_world,ierr)
353  CALL mpi_bcast(minlon,1,mpi_real,0,mpi_comm_world,ierr)
354  CALL mpi_bcast(maxlon,1,mpi_real,0,mpi_comm_world,ierr)
355  CALL mpi_bcast(minlat,1,mpi_real,0,mpi_comm_world,ierr)
356  CALL mpi_bcast(maxlat,1,mpi_real,0,mpi_comm_world,ierr)
357  CALL mpi_bcast(mxcwt,1,mpi_integer,0,mpi_comm_world,ierr)
358  CALL mpi_bcast(mycwt,1,mpi_integer,0,mpi_comm_world,ierr)
359  CALL mpi_bcast(dirknob,1,mpi_real,0,mpi_comm_world,ierr)
360  CALL mpi_bcast(perknob,1,mpi_real,0,mpi_comm_world,ierr)
361  CALL mpi_bcast(hsknob,1,mpi_real,0,mpi_comm_world,ierr)
362  CALL mpi_bcast(wetpts,1,mpi_real,0,mpi_comm_world,ierr)
363  CALL mpi_bcast(dirtimeknob,1,mpi_real,0,mpi_comm_world,ierr)
364  CALL mpi_bcast(tptimeknob,1,mpi_real,0,mpi_comm_world,ierr)
365  CALL mpi_bcast(seedlon,1,mpi_real,0,mpi_comm_world,ierr)
366  CALL mpi_bcast(seedlat,1,mpi_real,0,mpi_comm_world,ierr)
367  CALL mpi_bcast(noutp,1,mpi_integer,0,mpi_comm_world,ierr)
368  CALL mpi_bcast(lonout,100,mpi_real,0,mpi_comm_world,ierr)
369  CALL mpi_bcast(latout,100,mpi_real,0,mpi_comm_world,ierr)
370 #endif
371 
372 #ifdef W3_MPI
373  CALL mpi_barrier(mpi_comm_world,ierr)
374 #endif
375 
376  CALL wavetracking_nws_v2 (intype ,tmax , &
377  tcur ,filename , &
378  tstart ,tend , &
379  dt ,ntint , &
380  minlon ,maxlon , &
381  minlat ,maxlat , &
382  mxcwt ,mycwt , &
383  dirknob , &
384  perknob ,hsknob , &
385  wetpts ,seedlat , &
386  seedlon ,dirtimeknob, &
387  tptimeknob ,paramfile , &
388  sysa ,wsdat , &
389  maxsys ,maxgroup )
390 
391 #ifdef W3_MPI
392  IF (rank.EQ.0) THEN
393 #endif
394 
395  CALL date_and_time ( values=clkdt1 )
396  clkfel = tdiff( clkdt0,clkdt1 )
397  WRITE (6,998) clkfel
398  WRITE (6,*) 'Final system output...'
399 
400  ! Set upper limit for wave systems to output (limited by AWIPS display)
401  ulimgroup = 9
402 
403  !-----Output systems as plain text----------------------------------------
404 
405  maxi = SIZE(wsdat(1)%lon,1)
406  maxj = SIZE(wsdat(1)%lon,2)
407  dlon = wsdat(1)%lon(2,2)-wsdat(1)%lon(1,1)
408  dlat = wsdat(1)%lat(2,2)-wsdat(1)%lat(1,1)
409  WRITE(20,*) 'dlon, dlat =',dlon,dlat
410 
411  !-----Final SYSTEM output: Coordinates
412  OPEN(unit=21,file='sys_coord.ww3', status='unknown')
413 
414  WRITE(21,'(I6,69X,A)') maxj,'Number of rows'
415  WRITE(21,'(I6,69X,A)') maxi,'Number of cols'
416 #ifdef W3_TRKNC
417  ALLOCATE( longitude_nc(maxi) )
418  ALLOCATE( latitude_nc(maxj) )
419 #endif
420 
421  WRITE(21,*) 'Longitude ='
422  DO j = maxj,1,-1
423  DO i = 1,maxi
424  WRITE(21,'(F7.2)',advance='NO') wsdat(1)%lon(i,j)
425 #ifdef W3_TRKNC
426  longitude_nc(i)=wsdat(1)%lon(i,1)
427 #endif
428  END DO
429  WRITE(21,'(A)',advance='YES') ''
430  END DO
431 
432  WRITE(21,*) 'Latitude = '
433  DO j = maxj,1,-1
434  DO i = 1,maxi
435  WRITE(21,'(F7.2)',advance='NO') wsdat(1)%lat(i,j)
436 #ifdef W3_TRKNC
437  latitude_nc(j)=wsdat(1)%lat(1,j)
438 #endif
439  END DO
440  WRITE(21,'(A)',advance='YES') ''
441  END DO
442 
443  CLOSE(21)
444 
445  !-----Final SYSTEM output: hs
446  IF(outputtype == 1) THEN
447  OPEN(unit=22,file='sys_hs.ww3', status='unknown')
448 
449  WRITE(22,'(I6,69X,A)') maxj,'Number of rows'
450  WRITE(22,'(I6,69X,A)') maxi,'Number of cols'
451  ENDIF
452 
453  ntime_nc=SIZE(sysa)
454  ALLOCATE( dum(maxi,maxj) )
455 #ifdef W3_TRKNC
456  IF(outputtype == 3 .OR. outputtype == 4) THEN
457  ALLOCATE( dum2nc(maxi,maxj,maxgroup,ntime_nc) )
458  ENDIF
459 #endif
460 
461  DO it = 1,SIZE(sysa)
462  ! Loop through identified groups, limiting the output in file to ulimGroup
463  IF(outputtype == 1) THEN
464  WRITE(22,'(F15.6,60x,A)') wsdat(it)%date,'Time'
465  WRITE(22,'(I6,69x,A)') min(ulimgroup,maxgroup), &
466  'Tot number of systems'
467  ENDIF
468  DO igrp = 1,min(ulimgroup,maxgroup)
469  dum(1:maxi,1:maxj) = 9999.00
470  ! Find system with this group tag
471  sysmatch = 1
472  DO WHILE (sysmatch.LE.maxsys(it))
473  IF (sysa(it)%sys(sysmatch)%grp.EQ.igrp) EXIT
474  sysmatch = sysmatch+1
475  END DO
476  IF (sysmatch.LE.maxsys(it)) THEN
477  ! Match found: fill the output matrix with this data
478  leng = sysa(it)%sys(sysmatch)%nPoints
479  DO ind = 1, leng
480  dum(sysa(it)%sys(sysmatch)%i(ind), &
481  sysa(it)%sys(sysmatch)%j(ind)) = &
482  sysa(it)%sys(sysmatch)%hs(ind)
483  END DO
484  ELSE
485  leng = 0
486  END IF
487 
488  IF(outputtype == 1) THEN
489  WRITE(22,'(I6,69x,A)') igrp,'System number'
490  WRITE(22,'(I6,69x,A)') leng,'Number of points in system'
491 
492  DO j = maxj,1,-1
493  DO i = 1,maxi
494  WRITE(22,'(F8.2)',advance='NO') dum(i,j)
495  END DO
496  WRITE(22,'(A)',advance='YES') ''
497  END DO
498  ELSE
499 #ifdef W3_TRKNC
500  DO j = maxj,1,-1
501  DO i = 1,maxi
502  dum2nc(i,j,igrp,it)=dum(i,j)
503  END DO
504  END DO
505 #endif
506  ENDIF
507 
508  END DO
509  END DO
510 
511 #ifdef W3_TRKNC
512  IF(outputtype == 3 .OR. outputtype == 4 ) THEN
513  call t2netcdf(longitude_nc,latitude_nc,dum2nc,maxi,maxj,&
514  maxgroup,date1,date2,dt,ntime_nc,1,outputtype)
515  ENDIF
516 #endif
517 
518  IF(outputtype.EQ.1) CLOSE(22)
519 
520  !-----Final SYSTEM output: tp
521  IF(outputtype == 1) THEN
522  OPEN(unit=23,file='sys_tp.ww3',status='unknown')
523 
524  WRITE(23,'(I6,69X,A)') maxj,'Number of rows'
525  WRITE(23,'(I6,69X,A)') maxi,'Number of cols'
526  ENDIF
527 
528  DO it = 1,SIZE(sysa)
529  ! Loop through identified groups, limiting the output in file to ulimGroup
530  IF(outputtype == 1) THEN
531  WRITE(23,'(F15.6,60x,A)') wsdat(it)%date,'Time'
532  WRITE(23,'(I6,69X,A)') min(ulimgroup,maxgroup), &
533  'Tot number of systems'
534  ENDIF
535  DO igrp = 1,min(ulimgroup,maxgroup)
536  dum(1:maxi,1:maxj) = 9999.00
537  ! Find system with this group tag
538  sysmatch = 1
539  DO WHILE (sysmatch.LE.maxsys(it))
540  IF (sysa(it)%sys(sysmatch)%grp.EQ.igrp) EXIT
541  sysmatch = sysmatch+1
542  END DO
543  IF (sysmatch.LE.maxsys(it)) THEN
544  ! Match found: fill the output matrix with this data
545  leng = sysa(it)%sys(sysmatch)%nPoints
546  DO ind = 1, leng
547  dum(sysa(it)%sys(sysmatch)%i(ind), &
548  sysa(it)%sys(sysmatch)%j(ind)) = &
549  sysa(it)%sys(sysmatch)%tp(ind)
550  END DO
551  ELSE
552  leng = 0
553  END IF
554 
555  IF(outputtype == 1) THEN
556  WRITE(23,'(I6,69X,A)') igrp,'System number'
557  WRITE(23,'(I6,69X,A)') leng,'Number of points in system'
558  DO j = maxj,1,-1
559  DO i = 1,maxi
560  WRITE(23,'(F8.2)',advance='NO') dum(i,j)
561  END DO
562  WRITE(23,'(A)',advance='YES') ''
563  END DO
564  ELSE
565 
566 #ifdef W3_TRKNC
567  DO j = maxj,1,-1
568  DO i = 1,maxi
569  dum2nc(i,j,igrp,it)=dum(i,j)
570  END DO
571  END DO
572 #endif
573  ENDIF
574 
575  END DO
576  END DO
577 
578 #ifdef W3_TRKNC
579  IF(outputtype.EQ.3 .OR. outputtype.EQ. 4 ) THEN
580  call t2netcdf(longitude_nc,latitude_nc,dum2nc,maxi,maxj,&
581  maxgroup,date1,date2,dt,ntime_nc,2,outputtype)
582  ENDIF
583 #endif
584 
585  IF(outputtype.EQ.1) CLOSE(23)
586 
587  !-----Final SYSTEM output: dir
588  IF(outputtype == 1) THEN
589  OPEN(unit=24,file='sys_dir.ww3',status='unknown')
590 
591  WRITE(24,'(I6,69X,A)') maxj,'Number of rows'
592  WRITE(24,'(I6,69X,A)') maxi,'Number of cols'
593  ENDIF
594 
595  DO it = 1,SIZE(sysa)
596  ! Loop through identified groups, limiting the output in file to
597  ! ulimGroup
598  IF(outputtype == 1) THEN
599  WRITE(24,'(F15.6,60x,A)') wsdat(it)%date,'Time'
600  WRITE(24,'(I6,69X,A)') min(ulimgroup,maxgroup), &
601  'Tot number of systems'
602  ENDIF
603  DO igrp = 1,min(ulimgroup,maxgroup)
604  dum(1:maxi,1:maxj) = 9999.00
605  ! Find system with this group tag
606  sysmatch = 1
607  DO WHILE (sysmatch.LE.maxsys(it))
608  IF (sysa(it)%sys(sysmatch)%grp.EQ.igrp) EXIT
609  sysmatch = sysmatch+1
610  END DO
611  IF (sysmatch.LE.maxsys(it)) THEN
612  ! Match found: fill the output matrix with this data
613  leng = sysa(it)%sys(sysmatch)%nPoints
614  DO ind = 1, leng
615  dum(sysa(it)%sys(sysmatch)%i(ind), &
616  sysa(it)%sys(sysmatch)%j(ind)) = &
617  sysa(it)%sys(sysmatch)%dir(ind)
618  END DO
619  ELSE
620  leng = 0
621  END IF
622 
623  IF(outputtype == 1) THEN
624  WRITE(24,'(I6,69X,A)') igrp,'System number'
625  WRITE(24,'(I6,69X,A)') leng,'Number of points in system'
626  DO j = maxj,1,-1
627  DO i = 1,maxi
628  WRITE(24,'(F8.2)',advance='NO') dum(i,j)
629  END DO
630  WRITE(24,'(A)',advance='YES') ''
631  END DO
632  ELSE
633 #ifdef W3_TRKNC
634  DO j = maxj,1,-1
635  DO i = 1,maxi
636  dum2nc(i,j,igrp,it)=dum(i,j)
637  END DO
638  END DO
639 #endif
640  END IF
641 
642  END DO
643  END DO
644 
645 #ifdef W3_TRKNC
646  IF(outputtype.EQ.3 .OR. outputtype.EQ.4 ) THEN
647  call t2netcdf(longitude_nc,latitude_nc,dum2nc,maxi,maxj,&
648  maxgroup,date1,date2,dt,ntime_nc,3,outputtype)
649  ENDIF
650 #endif
651  IF(outputtype.EQ.1) CLOSE(24)
652 
653  !-----Final SYSTEM output: dspr
654  IF(outputtype == 1) THEN
655  OPEN(unit=25,file='sys_dspr.ww3',status='unknown')
656 
657  WRITE(25,'(I6,69X,A)') maxj,'Number of rows'
658  WRITE(25,'(I6,69X,A)') maxi,'Number of cols'
659  ENDIF
660 
661  DO it = 1,SIZE(sysa)
662  ! Loop through identified groups, limiting the output in file to ulimGroup
663  IF(outputtype == 1) THEN
664  WRITE(25,'(F15.6,60x,A)') wsdat(it)%date,'Time'
665  WRITE(25,'(I6,69X,A)') min(ulimgroup,maxgroup), &
666  'Tot number of systems'
667  ENDIF
668  DO igrp = 1,min(ulimgroup,maxgroup)
669  dum(1:maxi,1:maxj) = 9999.00
670  ! Find system with this group tag
671  sysmatch = 1
672  DO WHILE (sysmatch.LE.maxsys(it))
673  IF (sysa(it)%sys(sysmatch)%grp.EQ.igrp) EXIT
674  sysmatch = sysmatch+1
675  END DO
676  IF (sysmatch.LE.maxsys(it)) THEN
677  ! Match found: fill the output matrix with this data
678  leng = sysa(it)%sys(sysmatch)%nPoints
679  DO ind = 1, leng
680  dum(sysa(it)%sys(sysmatch)%i(ind), &
681  sysa(it)%sys(sysmatch)%j(ind)) = &
682  sysa(it)%sys(sysmatch)%dspr(ind)
683  END DO
684  ELSE
685  leng = 0
686  END IF
687 
688  IF(outputtype == 1) THEN
689  WRITE(25,'(I6,69X,A)') igrp,'System number'
690  WRITE(25,'(I6,69X,A)') leng,'Number of points in system'
691  DO j = maxj,1,-1
692  DO i = 1,maxi
693  WRITE(25,'(F8.2)',advance='NO') dum(i,j)
694  END DO
695  WRITE(25,'(A)',advance='YES') ''
696  END DO
697  ELSE
698 #ifdef W3_TRKNC
699  DO j = maxj,1,-1
700  DO i = 1,maxi
701  dum2nc(i,j,igrp,it)=dum(i,j)
702  END DO
703  END DO
704 #endif
705  ENDIF
706 
707  END DO
708  END DO
709 
710 #ifdef W3_TRKNC
711  IF(outputtype.EQ.3 .OR. outputtype.EQ.4 ) THEN
712  call t2netcdf(longitude_nc,latitude_nc,dum2nc,maxi,maxj,&
713  maxgroup,date1,date2,dt,ntime_nc,4,outputtype)
714  ENDIF
715 #endif
716  IF(outputtype.EQ.1) CLOSE(25)
717 
718  IF (ALLOCATED(dum)) DEALLOCATE(dum)
719 #ifdef W3_TRKNC
720  IF (ALLOCATED(dum2nc)) DEALLOCATE(dum2nc)
721 #endif
722 
723 #ifdef W3_TRKNC
724  IF(outputtype.EQ.3.OR.outputtype.EQ.4) THEN
725  ALLOCATE( hsprt_nc(10,noutp,ntime_nc) )
726  ALLOCATE( tpprt_nc(10,noutp,ntime_nc) )
727  ALLOCATE( dirprt_nc(10,noutp,ntime_nc) )
728  ALLOCATE( lonprt_nc(noutp) )
729  ALLOCATE( latprt_nc(noutp) )
730  ENDIF
731 #endif
732 
733  !-----Final SYSTEM output: point output
734  IF(outputtype == 1) THEN
735  OPEN(unit=26,file='sys_pnt.ww3',status='unknown')
736  WRITE(26,'(A)') '%'
737  WRITE(26,'(A)') '%'
738  WRITE(26,'(A)') '% WW3 Wave tracking point output'
739  WRITE(26,'(A)') '%'
740  WRITE(26,'(10A)') '% Xp Yp ', &
741  'HsSY01 HsSY02 HsSY03 HsSY04 ', &
742  'HsSY05 HsSY06 HsSY07 HsSY08 ', &
743  'HsSY09 HsSY10 ', &
744  'TpSY01 TpSY02 TpSY03 TpSY04 ', &
745  'TpSY05 TpSY06 TpSY07 TpSY08 ', &
746  'TpSY09 TpSY10 ', &
747  'DrSY01 DrSY02 DrSY03 DrSY04 ', &
748  'DrSY05 DrSY06 DrSY07 DrSY08 ', &
749  'DrSY09 DrSY10'
750  WRITE(26,'(10A)') '% [degr] [degr] ', &
751  '[m] [m] [m] [m] ', &
752  '[m] [m] [m] [m] ', &
753  '[m] [m] ', &
754  '[sec] [sec] [sec] [sec] ', &
755  '[sec] [sec] [sec] [sec] ', &
756  '[sec] [sec] ', &
757  '[degr] [degr] [degr] [degr] ', &
758  '[degr] [degr] [degr] [degr] ', &
759  '[degr] [degr]'
760  WRITE(26,'(A)') '%'
761  ENDIF
762 
763  DO it = 1,SIZE(sysa)
764  IF(outputtype == 1) THEN
765  WRITE(26,'(A,F15.6)') 'Time : ',wsdat(it)%date
766  ENDIF
767 
768  DO ip = 1,noutp
769  hsprt(1:10) = 999.9999
770  tpprt(1:10) = 999.9999
771  dirprt(1:10) = 999.9999
772  lonprt = 999.9999
773  latprt = 999.9999
774  bl_hsprt(1:10) = 999.9999
775  bl_tpprt(1:10) = 999.9999
776  bl_dirprt(1:10) = 999.9999
777  br_hsprt(1:10) = 999.9999
778  br_tpprt(1:10) = 999.9999
779  br_dirprt(1:10) = 999.9999
780  tl_hsprt(1:10) = 999.9999
781  tl_tpprt(1:10) = 999.9999
782  tl_dirprt(1:10) = 999.9999
783  tr_hsprt(1:10) = 999.9999
784  tr_tpprt(1:10) = 999.9999
785  tr_dirprt(1:10) = 999.9999
786  bl_lonprt = 999.9999
787  bl_latprt = 999.9999
788  br_lonprt = 999.9999
789  br_latprt = 999.9999
790  tl_lonprt = 999.9999
791  tl_latprt = 999.9999
792  tr_lonprt = 999.9999
793  tr_latprt = 999.9999
794  bl_w = 999
795  br_w = 999
796  tr_w = 999
797  tl_w = 999
798 
799  DO j = 1, (maxj-1)
800  DO i = 1, (maxi-1)
801  IF ( ( ((lonout(ip).GE. &
802  wsdat(1)%lon(i,j)).AND. &
803  (lonout(ip).LT. &
804  wsdat(1)%lon(i+1,j))).OR. &
805  ((lonout(ip).GT. &
806  wsdat(1)%lon(i,j)).AND. &
807  (lonout(ip).LE. &
808  wsdat(1)%lon(i+1,j))) ).AND. &
809  ( ((latout(ip).GE. &
810  wsdat(1)%lat(i,j)).AND. &
811  (latout(ip).LT. &
812  wsdat(1)%lat(i,j+1))).OR. &
813  ((latout(ip).GT. &
814  wsdat(1)%lat(i,j)).AND. &
815  (latout(ip).LE. &
816  wsdat(1)%lat(i,j+1))) ) ) &
817  THEN
818  bl_lonprt = wsdat(1)%lon(i,j)
819  bl_latprt = wsdat(1)%lat(i,j)
820  br_lonprt = wsdat(1)%lon(i+1,j)
821  br_latprt = wsdat(1)%lat(i+1,j)
822  tl_lonprt = wsdat(1)%lon(i,j+1)
823  tl_latprt = wsdat(1)%lat(i,j+1)
824  tr_lonprt = wsdat(1)%lon(i+1,j+1)
825  tr_latprt = wsdat(1)%lat(i+1,j+1)
826  ! Compute weights for this point
827  t = (lonout(ip)-bl_lonprt)/(br_lonprt-bl_lonprt)
828  u = (latout(ip)-bl_latprt)/(tl_latprt-bl_latprt)
829  bl_w = (1-t)*(1-u)
830  br_w = t*(1-u)
831  tr_w = t*u
832  tl_w = (1-t)*u
833  ! Compute output values using weights
834  lonprt = bl_w*bl_lonprt + br_w*br_lonprt + &
835  tl_w*tl_lonprt + tr_w*tr_lonprt
836  latprt = bl_w*bl_latprt + br_w*br_latprt + &
837  tl_w*tl_latprt + tr_w*tr_latprt
838  END IF
839  END DO
840  END DO
841  ! Loop through identified groups, limiting the output in file to 10
842  DO igrp = 1,min(10,maxgroup)
843  ! Find system with this group tag
844  sysmatch = 1
845  DO WHILE (sysmatch.LE.maxsys(it))
846  IF (sysa(it)%sys(sysmatch)%grp.EQ.igrp) EXIT
847  sysmatch = sysmatch+1
848  END DO
849  IF (sysmatch.LE.maxsys(it)) THEN
850  ! Match found: fill the output matrix with this data
851  leng = sysa(it)%sys(sysmatch)%nPoints
852  DO ind = 1, leng
853  ! Write output point data with bilinear interpolation
854  IF ( (sysa(it)%sys(sysmatch)%lon(ind).EQ.&
855  bl_lonprt).AND.&
856  (sysa(it)%sys(sysmatch)%lat(ind).EQ.&
857  bl_latprt) ) THEN
858  bl_hsprt(igrp) = sysa(it)%sys(sysmatch)%hs(ind)
859  bl_tpprt(igrp) = sysa(it)%sys(sysmatch)%tp(ind)
860  bl_dirprt(igrp) = sysa(it)%sys(sysmatch)%dir(ind)
861  ELSE IF ( (sysa(it)%sys(sysmatch)%lon(ind).EQ.&
862  br_lonprt).AND.&
863  (sysa(it)%sys(sysmatch)%lat(ind).EQ.&
864  br_latprt)) THEN
865  br_hsprt(igrp) = sysa(it)%sys(sysmatch)%hs(ind)
866  br_tpprt(igrp) = sysa(it)%sys(sysmatch)%tp(ind)
867  br_dirprt(igrp) = sysa(it)%sys(sysmatch)%dir(ind)
868  ELSE IF ( (sysa(it)%sys(sysmatch)%lon(ind).EQ.&
869  tl_lonprt).AND.&
870  (sysa(it)%sys(sysmatch)%lat(ind).EQ.&
871  tl_latprt)) THEN
872  tl_hsprt(igrp) = sysa(it)%sys(sysmatch)%hs(ind)
873  tl_tpprt(igrp) = sysa(it)%sys(sysmatch)%tp(ind)
874  tl_dirprt(igrp) = sysa(it)%sys(sysmatch)%dir(ind)
875  ELSE IF ( (sysa(it)%sys(sysmatch)%lon(ind).EQ.&
876  tr_lonprt).AND.&
877  (sysa(it)%sys(sysmatch)%lat(ind).EQ.&
878  tr_latprt)) THEN
879  tr_hsprt(igrp) = sysa(it)%sys(sysmatch)%hs(ind)
880  tr_tpprt(igrp) = sysa(it)%sys(sysmatch)%tp(ind)
881  tr_dirprt(igrp) = sysa(it)%sys(sysmatch)%dir(ind)
882  END IF
883  END DO
884  ! Compute output value using weights
885  ! (only if output point is surrounded by valid points)
886  IF ( (bl_hsprt(igrp).NE.999.9999).AND. &
887  (br_hsprt(igrp).NE.999.9999).AND. &
888  (tl_hsprt(igrp).NE.999.9999).AND. &
889  (tr_hsprt(igrp).NE.999.9999) ) THEN
890  hsprt(igrp) = bl_w * bl_hsprt(igrp) + &
891  br_w * br_hsprt(igrp) + &
892  tl_w * tl_hsprt(igrp) + &
893  tr_w * tr_hsprt(igrp)
894  tpprt(igrp) = bl_w * bl_tpprt(igrp) + &
895  br_w * br_tpprt(igrp) + &
896  tl_w * tl_tpprt(igrp) + &
897  tr_w * tr_tpprt(igrp)
898  bl_dirx = cos((270-bl_dirprt(igrp))*pi/180.)
899  br_dirx = cos((270-br_dirprt(igrp))*pi/180.)
900  tr_dirx = cos((270-tr_dirprt(igrp))*pi/180.)
901  tl_dirx = cos((270-tl_dirprt(igrp))*pi/180.)
902  bl_diry = sin((270-bl_dirprt(igrp))*pi/180.)
903  br_diry = sin((270-br_dirprt(igrp))*pi/180.)
904  tr_diry = sin((270-tr_dirprt(igrp))*pi/180.)
905  tl_diry = sin((270-tl_dirprt(igrp))*pi/180.)
906  dirprt(igrp)=270 - 180./pi* &
907  atan2(bl_w*bl_diry+br_w*br_diry+ &
908  tl_w*tl_diry+tr_w*tr_diry, &
909  bl_w*bl_dirx+br_w*br_dirx+ &
910  tl_w*tl_dirx+tr_w*tr_dirx)
911  IF (dirprt(igrp).GT.360.) THEN
912  dirprt(igrp) = dirprt(igrp) - 360.
913  END IF
914  ELSE
915  hsprt(igrp) = 999.9999
916  tpprt(igrp) = 999.9999
917  dirprt(igrp) = 999.9999
918  END IF
919  END IF
920  END DO
921  IF(outputtype == 1) THEN
922  WRITE(26,'(32F14.4)') lonprt,latprt, &
923  hsprt(1:10),tpprt(1:10),dirprt(1:10)
924  ENDIF
925 #ifdef W3_TRKNC
926  IF(outputtype.EQ.3.OR.outputtype.EQ.4) THEN
927  lonprt_nc(ip)=lonprt
928  latprt_nc(ip)=latprt
929  do igrp=1,10
930  hsprt_nc(igrp,ip,it)=hsprt(igrp)
931  tpprt_nc(igrp,ip,it)=tpprt(igrp)
932  dirprt_nc(igrp,ip,it)=dirprt(igrp)
933  enddo
934  ENDIF
935 #endif
936 
937  END DO
938  END DO
939 #ifdef W3_TRKNC
940  IF(outputtype.EQ.3.OR.outputtype.EQ.4) THEN
941  call pt2netcdf(lonprt_nc,latprt_nc,hsprt_nc,tpprt_nc, &
942  dirprt_nc,noutp,date1,date2,dt,ntime_nc,outputtype)
943  ENDIF
944 #endif
945 
946  IF(outputtype.EQ.1) CLOSE(26)
947 
948  !-----Final SYSTEM output: point output (Nearest neighbor, as a double check)
949  IF (testout) THEN
950  OPEN(unit=28,file='sys_pnt_nn.ww3',status='unknown')
951  WRITE(28,'(A)') '%'
952  WRITE(28,'(A)') '%'
953  WRITE(28,'(A)') '% WW3 Wave tracking point output'
954  WRITE(28,'(A)') '%'
955  WRITE(28,'(10A)') '% Xp Yp ', &
956  'HsSY01 HsSY02 HsSY03 HsSY04 ', &
957  'HsSY05 HsSY06 HsSY07 HsSY08 ', &
958  'HsSY09 HsSY10 ', &
959  'TpSY01 TpSY02 TpSY03 TpSY04 ', &
960  'TpSY05 TpSY06 TpSY07 TpSY08 ', &
961  'TpSY09 TpSY10 ', &
962  'DrSY01 DrSY02 DrSY03 DrSY04 ', &
963  'DrSY05 DrSY06 DrSY07 DrSY08 ', &
964  'DrSY09 DrSY10'
965  WRITE(28,'(10A)') '% [degr] [degr] ', &
966  '[m] [m] [m] [m] ', &
967  '[m] [m] [m] [m] ', &
968  '[m] [m] ', &
969  '[sec] [sec] [sec] [sec] ', &
970  '[sec] [sec] [sec] [sec] ', &
971  '[sec] [sec] ', &
972  '[degr] [degr] [degr] [degr] ', &
973  '[degr] [degr] [degr] [degr] ', &
974  '[degr] [degr]'
975  WRITE(28,'(A)') '%'
976 
977  DO it = 1,SIZE(sysa)
978  WRITE(28,'(A,F15.6)') 'Time : ',wsdat(it)%date
979 
980  DO ip = 1,noutp
981  hsprt(1:10) = 999.9999
982  tpprt(1:10) = 999.9999
983  dirprt(1:10) = 999.9999
984  lonprt = 999.9999
985  latprt = 999.9999
986 
987  DO j = 1, maxj
988  DO i = 1, maxi
989  ! Write nearest nearbor output (no bilinear interpolation)
990  IF ( (lonout(ip).GE. &
991  (wsdat(1)%lon(i,j)-dlon/2)).AND. &
992  (lonout(ip).LT. &
993  (wsdat(1)%lon(i,j)+dlon/2)).AND. &
994  (latout(ip).GE. &
995  (wsdat(1)%lat(i,j)-dlat/2)).AND. &
996  (latout(ip).LT. &
997  (wsdat(1)%lat(i,j)+dlat/2)) ) &
998  THEN
999  lonprt = wsdat(1)%lon(i,j)
1000  latprt = wsdat(1)%lat(i,j)
1001  END IF
1002  END DO
1003  END DO
1004  ! Loop through identified groups, limiting the output in file to 10
1005  DO igrp = 1,min(10,maxgroup)
1006  ! Find system with this group tag
1007  sysmatch = 1
1008  DO WHILE (sysmatch.LE.maxsys(it))
1009  IF (sysa(it)%sys(sysmatch)%grp.EQ.igrp) EXIT
1010  sysmatch = sysmatch+1
1011  END DO
1012  IF (sysmatch.LE.maxsys(it)) THEN
1013  ! Match found: fill the output matrix with this data
1014  leng = sysa(it)%sys(sysmatch)%nPoints
1015  DO ind = 1, leng
1016  ! Write nearest nearbor output (no bilinear interpolation)
1017  IF ( (lonout(ip).GE. &
1018  (sysa(it)%sys(sysmatch)%lon(ind)-dlon/2)).AND. &
1019  (lonout(ip).LT. &
1020  (sysa(it)%sys(sysmatch)%lon(ind)+dlon/2)).AND. &
1021  (latout(ip).GE. &
1022  (sysa(it)%sys(sysmatch)%lat(ind)-dlat/2)).AND. &
1023  (latout(ip).LT. &
1024  (sysa(it)%sys(sysmatch)%lat(ind)+dlat/2)) ) &
1025  THEN
1026  hsprt(igrp) = sysa(it)%sys(sysmatch)%hs(ind)
1027  tpprt(igrp) = sysa(it)%sys(sysmatch)%tp(ind)
1028  dirprt(igrp) = sysa(it)%sys(sysmatch)%dir(ind)
1029  END IF
1030  END DO
1031  END IF
1032  END DO
1033  WRITE(28,'(32F14.4)') lonprt,latprt, &
1034  hsprt(1:10),tpprt(1:10),dirprt(1:10)
1035  END DO
1036  END DO
1037 
1038  CLOSE(28)
1039  END IF
1040 
1041  !-------------------------------------------------------------------------
1042 
1043  WRITE(20,*) 'In ww3_systrk: Deallocating wsdat ...'
1044  DO it=1,size(wsdat)
1045  IF (ASSOCIATED(wsdat(it)%lat)) DEALLOCATE(wsdat(it)%lat)
1046  IF (ASSOCIATED(wsdat(it)%lon)) DEALLOCATE(wsdat(it)%lon)
1047  IF (ASSOCIATED(wsdat(it)%par)) DEALLOCATE(wsdat(it)%par)
1048  IF (ASSOCIATED(wsdat(it)%wnd)) DEALLOCATE(wsdat(it)%wnd)
1049  END DO
1050  IF (ASSOCIATED(wsdat)) DEALLOCATE(wsdat)
1051  WRITE(20,*) ' Deallocating sysA ...'
1052  DO i=1,size(sysa)
1053  DO j=1,size(sysa(i)%sys)
1054  IF (ASSOCIATED(sysa(i)%sys(j)%i)) DEALLOCATE(sysa(i)%sys(j)%i)
1055  IF (ASSOCIATED(sysa(i)%sys(j)%j)) DEALLOCATE(sysa(i)%sys(j)%j)
1056  IF (ASSOCIATED(sysa(i)%sys(j)%lon)) &
1057  DEALLOCATE(sysa(i)%sys(j)%lon)
1058  IF (ASSOCIATED(sysa(i)%sys(j)%lat)) &
1059  DEALLOCATE(sysa(i)%sys(j)%lat)
1060  IF (ASSOCIATED(sysa(i)%sys(j)%hs)) &
1061  DEALLOCATE(sysa(i)%sys(j)%hs)
1062  IF (ASSOCIATED(sysa(i)%sys(j)%tp)) &
1063  DEALLOCATE(sysa(i)%sys(j)%tp)
1064  IF (ASSOCIATED(sysa(i)%sys(j)%dir)) &
1065  DEALLOCATE(sysa(i)%sys(j)%dir)
1066  IF (ASSOCIATED(sysa(i)%sys(j)%dspr)) &
1067  DEALLOCATE(sysa(i)%sys(j)%dspr)
1068  END DO
1069  END DO
1070  IF (ASSOCIATED(sysa)) DEALLOCATE(sysa)
1071  WRITE(20,*) ' Deallocating maxSys ...'
1072  IF (ASSOCIATED(maxsys)) DEALLOCATE(maxsys)
1073  CLOSE(20)
1074 
1075  WRITE(6,*) '... ww3_systrk completed successfully.'
1076 
1077  WRITE(6,999)
1078 
1079 #ifdef W3_MPI
1080  END IF !/IF (rank.EQ.0)
1081 #endif
1082 
1083 #ifdef W3_MPI
1084  CALL mpi_finalize(ierr)
1085  ! End of parallel region
1086 #endif
1087 
1088 998 FORMAT ( ' ... finished. Elapsed time : ',f10.2,' s')
1089 993 FORMAT (/' *** WAVEWATCH III ERROR IN WW3_SYSTRK : '/ &
1090  ' OutputType=3 needs TRKNC switch ')
1091 994 FORMAT (/' *** WAVEWATCH III ERROR IN WW3_SYSTRK : '/ &
1092  ' OutputType=4 needs TRKNC switch ')
1093 995 FORMAT (/' *** WAVEWATCH III ERROR IN WW3_SYSTRK : '/ &
1094  ' OutputType,',i3,'not valid. Options: 1,3,4')
1095 
1096 999 FORMAT (/15x,'End of program '/ &
1097  15x,'==============================================='/ &
1098  15x,' *** WAVEWATCH III Wave system tracking *** ')
1099 
1100 2000 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/ &
1101  ' ERROR IN OPENING INPUT FILE')
1102 2200 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/ &
1103  ' ERROR IN OPENING PARTITION FILE : ',a)
1104 
1105 END PROGRAM ww3_systrk
1106 !
1107 #ifdef W3_TRKNC
1108 
1127 subroutine t2netcdf(lons,lats,data_in,nlons,nlats,nsys,date1,date2,&
1128  dt,ntime,ivar, outputType)
1129  USE w3timemd
1130  use netcdf
1131  implicit none
1132  character (len = 15) :: file_name
1133  integer, parameter :: ndims = 4
1134  integer, parameter :: deflate = 1
1135  integer :: outputType, ncid, oldMode
1136  integer :: nlons,nlats,nsys,rec,ntime,ivar
1137  double precision :: date1,date2,timenc
1138  real :: data_in(nlons, nlats, nsys,ntime)
1139  real :: lats(nlats), lons(nlons),dt
1140  double precision :: times(ntime)
1141  integer :: iyc,imc,idc,ihc,iminc,isc,Jday,Jday0
1142  integer :: iret
1143 #endif
1144  !
1145 #ifdef W3_TRKNC
1146  integer :: lon_varid, lat_varid, rec_varid
1147  character (len = *), parameter :: lsys_name = "system_index"
1148  character (len = *), parameter :: lat_name = "latitude"
1149  character (len = *), parameter :: lon_name = "longitude"
1150  character (len = *), parameter :: time_name = "time"
1151  integer :: sys_dimid, lon_dimid, lat_dimid, rec_dimid
1152  integer :: start(ndims), count(ndims)
1153 #endif
1154  !
1155 #ifdef W3_TRKNC
1156  character (len = *), parameter :: var1_name="hs"
1157  character (len = *), parameter :: var2_name="tp"
1158  character (len = *), parameter :: var3_name="dir"
1159  character (len = *), parameter :: var4_name="dspr"
1160  integer :: var1_varid, var2_varid, var3_varid,var4_varid
1161  integer :: dimids(ndims)
1162 #endif
1163  !
1164 #ifdef W3_TRKNC
1165  character (len = *), parameter :: units = "units"
1166  character (len = *), parameter :: var1_units = "m"
1167  character (len = *), parameter :: var2_units = "s"
1168  character (len = *), parameter :: var3_units = "degrees"
1169  character (len = *), parameter :: var4_units = "degrees"
1170  character (len = *), parameter :: lat_units = "degrees_north"
1171  character (len = *), parameter :: lon_units = "degrees_east"
1172  iyc=date1/10000
1173  imc=(date1-iyc*10000)/100
1174  idc=int(date1-dble(iyc*10000)-dble(imc*100))
1175  ihc=date2/10000
1176  iminc=(date2-ihc*10000)/100
1177  isc=date2-ihc*10000-100*iminc
1178  timenc=dble(julday(idc,imc,iyc))+(dble(ihc)+(dble(iminc)+ &
1179  (dble(isc)/60.0d0))/60.0d0)/24.0d0
1180  jday0=julday(1,1,1990)
1181  timenc=timenc-jday0
1182  do rec=1,ntime
1183  times(rec)=timenc+dble( (rec-1)*dt)/3600.0d0/24.0d0
1184  enddo
1185  if( ivar == 1) then
1186  file_name = "sys_hs.ww3.nc"
1187  else if( ivar == 2) then
1188  file_name = "sys_tp.ww3.nc"
1189  else if( ivar == 3) then
1190  file_name = "sys_dir.ww3.nc"
1191  else
1192  file_name = "sys_dspr.ww3.nc"
1193  endif
1194 #endif
1195  !
1196 #ifdef W3_TRKNC
1197  ! create the netcdf file.
1198  if (outputtype.EQ.3) then
1199  call check( nf90_create(file_name, nf90_clobber, ncid) )
1200  endif
1201  if(outputtype.EQ.4) call check( nf90_create(file_name,nf90_netcdf4,ncid))
1202  call check ( nf90_set_fill(ncid,nf90_nofill,oldmode) )
1203  call check( nf90_def_dim(ncid, lsys_name, nsys, sys_dimid) )
1204  call check( nf90_def_dim(ncid, lat_name, nlats, lat_dimid) )
1205  call check( nf90_def_dim(ncid, lon_name, nlons, lon_dimid) )
1206  call check( nf90_def_dim(ncid, time_name, ntime, rec_dimid) )
1207  call check( nf90_def_var(ncid, lat_name, nf90_real, lat_dimid,lat_varid))
1208  if(outputtype.EQ.4) call check( nf90_def_var_deflate(ncid,lat_varid,1,1,deflate) )
1209  call check( nf90_def_var(ncid, lon_name, nf90_real, lon_dimid,lon_varid))
1210  if(outputtype.EQ.4) call check( nf90_def_var_deflate(ncid,lon_varid,1,1,deflate) )
1211  call check( nf90_def_var(ncid,time_name,nf90_double,rec_dimid,rec_varid))
1212  if(outputtype.EQ.4) call check( nf90_def_var_deflate(ncid,rec_varid,1,1,deflate) )
1213 #endif
1214  !
1215 #ifdef W3_TRKNC
1216  call check( nf90_put_att(ncid, lat_varid, units, lat_units) )
1217  call check( nf90_put_att(ncid, lat_varid, 'long_name', 'latitude') )
1218  call check( nf90_put_att(ncid, lat_varid, 'standard_name', 'latitude') )
1219  call check( nf90_put_att(ncid, lat_varid, 'axis','Y'))
1220  call check( nf90_put_att(ncid, lon_varid, units, lon_units) )
1221  call check( nf90_put_att(ncid, lon_varid, 'long_name', 'longitude') )
1222  call check( nf90_put_att(ncid, lon_varid, 'standard_name', 'longitude') )
1223  call check( nf90_put_att(ncid, lon_varid, 'axis','X'))
1224  call check(nf90_put_att(ncid,rec_varid,units,&
1225  'days since 1990-01-01 00:00:00'))
1226  call check(nf90_put_att(ncid,rec_varid,'long_name','julian day (UT)'))
1227  call check( nf90_put_att(ncid, rec_varid,'standard_name', 'time') )
1228  call check( nf90_put_att(ncid, rec_varid, 'conventions',&
1229  'relative julian day with decimal part (as part of the day)' ) )
1230  call check( nf90_put_att(ncid, rec_varid, 'axis','T'))
1231 #endif
1232  !
1233 #ifdef W3_TRKNC
1234  dimids = (/ lon_dimid, lat_dimid, sys_dimid, rec_dimid /)
1235  if( ivar == 1) then
1236  call check( nf90_def_var(ncid, var1_name, nf90_real, dimids,var1_varid) )
1237  if(outputtype.EQ.4) call check( nf90_def_var_deflate(ncid,var1_varid,1,1,deflate) )
1238  call check( nf90_put_att(ncid, var1_varid, units, var1_units) )
1239  call check( nf90_put_att(ncid, var1_varid,'long_name','significant_wave_height') )
1240  call check( nf90_put_att(ncid, var1_varid,'missing_value','9999.00'))
1241  else if( ivar == 2) then
1242  call check( nf90_def_var(ncid, var2_name, nf90_real, dimids, var2_varid) )
1243  if(outputtype.EQ.4) call check( nf90_def_var_deflate(ncid,var2_varid,1,1,deflate) )
1244  call check( nf90_put_att(ncid, var2_varid, units, var2_units) )
1245  call check( nf90_put_att(ncid, var2_varid,'long_name','peak_period') )
1246  call check( nf90_put_att(ncid, var2_varid,'missing_value','9999.00') )
1247  else if ( ivar ==3 ) then
1248  call check( nf90_def_var(ncid, var3_name, nf90_real, dimids, var3_varid) )
1249  if(outputtype.EQ.4) call check( nf90_def_var_deflate(ncid,var3_varid,1,1,deflate) )
1250  call check( nf90_put_att(ncid, var3_varid, units, var3_units) )
1251  call check( nf90_put_att(ncid, var3_varid,'long_name','peak_direction') )
1252  call check( nf90_put_att(ncid, var3_varid,'missing_value','9999.00') )
1253  else
1254  call check( nf90_def_var(ncid, var4_name, nf90_real, dimids, var4_varid) )
1255  if(outputtype.EQ.4) call check( nf90_def_var_deflate(ncid,var4_varid,1,1,deflate) )
1256  call check( nf90_put_att(ncid, var4_varid, units, var4_units) )
1257  call check( nf90_put_att(ncid,var4_varid,'long_name','directional_spread') )
1258  call check( nf90_put_att(ncid, var4_varid,'missing_value','9999.00') )
1259  endif
1260  call check( nf90_enddef(ncid) )
1261 #endif
1262  !
1263 #ifdef W3_TRKNC
1264  call check( nf90_put_var(ncid, lat_varid, lats) )
1265  call check( nf90_put_var(ncid, lon_varid, lons) )
1266  call check( nf90_put_var(ncid, rec_varid, times) )
1267 #endif
1268  !
1269 #ifdef W3_TRKNC
1270  count = (/ nlons, nlats, nsys, ntime /)
1271  start = (/ 1, 1, 1, 1 /)
1272  if( ivar == 1) then
1273  call check( nf90_put_var(ncid, var1_varid, data_in, start = start, &
1274  count = count) )
1275  else if( ivar == 2) then
1276  call check( nf90_put_var(ncid, var2_varid, data_in, start = start, &
1277  count = count) )
1278  else if( ivar == 3) then
1279  call check( nf90_put_var(ncid, var3_varid, data_in, start = start, &
1280  count = count) )
1281  else
1282  call check( nf90_put_var(ncid, var4_varid, data_in, start = start, &
1283  count = count) )
1284  endif
1285  call check( nf90_close(ncid) )
1286 end subroutine t2netcdf
1287 #endif
1288 !
1289 #ifdef W3_TRKNC
1290 
1298 subroutine check(status)
1299  use netcdf
1300  integer, intent ( in) :: status
1301  if(status /= nf90_noerr) then
1302  write(6,996)
1303 996 FORMAT (/' *** WAVEWATCH III ERROR IN WW3_SYSTRK:'/ &
1304  'netCDF error:')
1305  print *, trim(nf90_strerror(status))
1306  stop "Stopped in netcdf output part"
1307  endif
1308 end subroutine check
1309 #endif
1310 !
1311 #ifdef W3_TRKNC
1312 
1330 subroutine pt2netcdf(longitude,latitude,hs,tp,&
1331  dir,npoints,date1,date2,dt,ntime,outputType)
1332  USE w3timemd
1333  use netcdf
1334  implicit none
1335  integer :: ntime,npoints,outputType
1336  integer, parameter :: deflate = 1
1337  integer :: iret, oldMode
1338  integer :: ncid
1339  integer :: system_index_dim
1340  integer :: point_dim,rec_dim
1341  integer :: nsys
1342  integer :: start(3), count(3)
1343  parameter(nsys = 10)
1344  integer :: latitude_id
1345  integer :: longitude_id
1346  integer :: time_id
1347  integer :: hs_id
1348  integer :: tp_id
1349  integer :: dir_id
1350  integer :: time_rank
1351  integer :: hs_rank
1352  integer :: tp_rank
1353  integer :: dir_rank
1354  parameter(time_rank = 1)
1355  parameter(hs_rank = 3)
1356  parameter(tp_rank = 3)
1357  parameter(dir_rank = 3)
1358 #endif
1359  !
1360 #ifdef W3_TRKNC
1361  integer :: hs_dims(hs_rank)
1362  integer :: tp_dims(tp_rank)
1363  integer :: dir_dims(dir_rank)
1364  real :: latitude(npoints),dt
1365  real :: longitude(npoints)
1366  real :: hs(nsys, npoints, ntime)
1367  real :: tp(nsys, npoints, ntime)
1368  real :: dir(nsys, npoints, ntime)
1369  integer :: iyc,imc,idc,ihc,iminc,isc,Jday,Jday0,rec
1370  double precision date1,date2,timenc
1371  double precision times(ntime)
1372 #endif
1373  !
1374 #ifdef W3_TRKNC
1375  iyc=date1/10000
1376  imc=(date1-iyc*10000)/100
1377  idc=int(date1-dble(iyc*10000)-dble(imc*100))
1378  ihc=date2/10000
1379  iminc=(date2-ihc*10000)/100
1380  isc=date2-ihc*10000-100*iminc
1381  timenc=dble(julday(idc,imc,iyc))+(dble(ihc)+(dble(iminc)+ &
1382  (dble(isc)/60.0d0))/60.0d0)/24.0d0
1383  jday0=julday(1,1,1990)
1384  timenc=timenc-jday0
1385  do rec=1,ntime
1386  times(rec)=timenc+dble( (rec-1)*dt)/3600.0d0/24.0d0
1387  enddo
1388 #endif
1389  !
1390 #ifdef W3_TRKNC
1391  if(outputtype.EQ.3) then
1392  iret = nf90_create('sys_pnt.ww3.nc', nf90_clobber, ncid)
1393  endif
1394  if (outputtype.EQ.4) iret = nf90_create('sys_pnt.ww3.nc',nf90_netcdf4, ncid)
1395  call check(iret)
1396  iret = nf90_set_fill(ncid,nf90_nofill,oldmode)
1397  call check(iret)
1398  ! define dimensions
1399  iret = nf90_def_dim(ncid, 'system_index', nsys, system_index_dim)
1400  call check(iret)
1401  iret = nf90_def_dim(ncid, 'point', npoints, point_dim)
1402  call check(iret)
1403  iret = nf90_def_dim(ncid, 'time', ntime, rec_dim)
1404  call check(iret)
1405  ! define variables
1406  iret = nf90_def_var(ncid, 'latitude', nf90_real, point_dim, &
1407  latitude_id)
1408  call check(iret)
1409  if (outputtype.EQ.4) call check( nf90_def_var_deflate(ncid,latitude_id,1,1,deflate))
1410  iret = nf90_def_var(ncid, 'longitude', nf90_real, point_dim, &
1411  longitude_id)
1412  call check(iret)
1413  if (outputtype.EQ.4) call check( nf90_def_var_deflate(ncid,longitude_id,1,1,deflate))
1414  iret = nf90_def_var(ncid, 'time', nf90_double, rec_dim, &
1415  time_id)
1416  call check(iret)
1417  if (outputtype.EQ.4) call check( nf90_def_var_deflate(ncid,time_id,1,1,deflate) )
1418  hs_dims(3) = rec_dim
1419  hs_dims(2) = point_dim
1420  hs_dims(1) = system_index_dim
1421  iret = nf90_def_var(ncid, 'hs', nf90_real, &
1422  hs_dims, hs_id)
1423  call check(iret)
1424  if (outputtype.EQ.4) call check( nf90_def_var_deflate(ncid,hs_id,1,1,deflate))
1425  tp_dims(3) = rec_dim
1426  tp_dims(2) = point_dim
1427  tp_dims(1) = system_index_dim
1428  iret = nf90_def_var(ncid, 'tp', nf90_real, &
1429  tp_dims, tp_id)
1430  call check(iret)
1431  if (outputtype.EQ.4) call check( nf90_def_var_deflate(ncid,tp_id,1,1,deflate))
1432  dir_dims(3) = rec_dim
1433  dir_dims(2) = point_dim
1434  dir_dims(1) = system_index_dim
1435  iret = nf90_def_var(ncid, 'dir', nf90_real, &
1436  dir_dims, dir_id)
1437  call check(iret)
1438  if (outputtype.EQ.4) call check( nf90_def_var_deflate(ncid,dir_id,1,1,deflate))
1439  ! assign attributes
1440  iret = nf90_put_att(ncid, latitude_id, 'units', 'degrees_north')
1441  call check(iret)
1442  iret = nf90_put_att(ncid, latitude_id, 'long_name', 'latitude')
1443  call check(iret)
1444  iret = nf90_put_att(ncid, latitude_id, 'standard_name', 'latitude')
1445  call check(iret)
1446  iret = nf90_put_att(ncid, latitude_id, 'axis', 'Y')
1447  call check(iret)
1448  iret = nf90_put_att(ncid, longitude_id, 'units', 'degrees_east')
1449  call check(iret)
1450  iret = nf90_put_att(ncid, longitude_id,'long_name','longitude')
1451  call check(iret)
1452  iret = nf90_put_att(ncid, longitude_id,'standard_name','longitude')
1453  call check(iret)
1454  iret = nf90_put_att(ncid, longitude_id, 'axis', 'X')
1455  call check(iret)
1456  iret = nf90_put_att(ncid, time_id, 'units', &
1457  'days since 1990-01-01 00:00:00')
1458  call check(iret)
1459  iret = nf90_put_att(ncid, time_id, 'long_name','julian day(UT)')
1460  call check(iret)
1461  iret = nf90_put_att(ncid, time_id, 'standard_name','time')
1462  call check(iret)
1463  iret = nf90_put_att(ncid, time_id, 'conventions', &
1464  'relative julian day with decimal part (as part of the day)')
1465  call check(iret)
1466  iret = nf90_put_att(ncid, time_id, 'axis', 'T')
1467  call check(iret)
1468  iret = nf90_put_att(ncid, hs_id, 'units', 'm')
1469  call check(iret)
1470  iret = nf90_put_att(ncid, hs_id,'long_name','significant_wave_height')
1471  call check(iret)
1472  iret = nf90_put_att(ncid, hs_id, 'missing_value', &
1473  '999.9999')
1474  call check(iret)
1475  iret = nf90_put_att(ncid, tp_id, 'units', 's')
1476  call check(iret)
1477  iret = nf90_put_att(ncid, tp_id,'long_name','peak_period')
1478  call check(iret)
1479  iret = nf90_put_att(ncid, tp_id, 'missing_value', &
1480  '999.9999')
1481  call check(iret)
1482  iret = nf90_put_att(ncid, dir_id, 'units', 'degrees')
1483  call check(iret)
1484  iret = nf90_put_att(ncid, dir_id,'long_name','peak_direction')
1485  call check(iret)
1486  iret = nf90_put_att(ncid, dir_id, 'missing_value',&
1487  '999.9999')
1488  call check(iret)
1489  ! leave define mode
1490  iret = nf90_enddef(ncid)
1491  call check(iret)
1492  iret = nf90_put_var(ncid, latitude_id, latitude)
1493  call check(iret)
1494 #endif
1495  !
1496 #ifdef W3_TRKNC
1497  iret = nf90_put_var(ncid, longitude_id, longitude)
1498  call check(iret)
1499 #endif
1500  !
1501 #ifdef W3_TRKNC
1502  iret = nf90_put_var(ncid, time_id, times)
1503  call check(iret)
1504 #endif
1505  !
1506 #ifdef W3_TRKNC
1507  start = (/ 1, 1, 1 /)
1508  count = (/ nsys,npoints,ntime /)
1509 #endif
1510  !
1511 #ifdef W3_TRKNC
1512  iret = nf90_put_var(ncid, hs_id, hs,&
1513  start = start, count = count )
1514  call check(iret)
1515  iret = nf90_put_var(ncid, tp_id, tp, &
1516  start = start, count = count )
1517  call check(iret)
1518  iret = nf90_put_var(ncid, dir_id, dir,&
1519  start = start, count = count )
1520  call check(iret)
1521  iret = nf90_close(ncid)
1522  call check(iret)
1523  return
1524 end subroutine pt2netcdf
1525 #endif
constants::pi
real, parameter pi
PI Value of Pi.
Definition: constants.F90:71
include
cmake src_list cmake include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/check_switches.cmake) check_switches("$
Definition: CMakeLists.txt:15
ww3_systrk
program ww3_systrk
Perform spatial and temporal tracking of wave systems, based on spectral partition (bulletin) output.
Definition: ww3_systrk.F90:29
w3timemd::tdiff
real function tdiff(T1, T2)
Definition: w3timemd.F90:576
w3strkmd::timsys
Definition: w3strkmd.F90:174
w3strkmd::dat2d
Definition: w3strkmd.F90:96
t2netcdf
subroutine t2netcdf(lons, lats, data_in, nlons, nlats, nsys, date1, date2, dt, ntime, ivar, outputType)
N/A.
Definition: ww3_systrk.F90:1129
pt2netcdf
subroutine pt2netcdf(longitude, latitude, hs, tp, dir, npoints, date1, date2, dt, ntime, outputType)
N/A.
Definition: ww3_systrk.F90:1332
file
file(STRINGS ${CMAKE_BINARY_DIR}/switch switch_strings) separate_arguments(switches UNIX_COMMAND $
Definition: CMakeLists.txt:3
w3strkmd
Definition: w3strkmd.F90:3
w3strkmd::wavetracking_nws_v2
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)
Definition: w3strkmd.F90:298
w3timemd::julday
integer function julday(id, mm, iyyy)
Definition: w3timemd.F90:756
w3timemd
Definition: w3timemd.F90:3
check
subroutine check(status)
N/A.
Definition: ww3_systrk.F90:1299