!/ ------------------------------------------------------------------- /
PROGRAM WW3_SYSTRK
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | A. J. van der Westhuysen |
!/ | Jeff Hanson |
!/ | Eve-Marie Devaliere |
!/ | FORTRAN 95 |
!/ | Last update : 16-Jan-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 )
!/ 29-Nov-2013 : Remove DOC control characters,
!/ update MPI! to MPI/! (H.L. Tolman). ( version 4.15 )
!/ 11-Feb-2014 : Add NetCDF output option. Both NetCDF-3 and
!/ NetCDF-4 are available. (B. Li). ( version 4.18 )
!/ 26-Sep-2016 : Optimization updates (A. van der Westhuysen)
!/ ( version 5.15 )
!/ 20-Sep-2016 : Add support for unformatted partition file.
!/ (S.Zieger BoM, Australia) ( version 5.16 )
!/ 20-Dec-2016 : Optimized search algorithms and
!/ set functions. (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.
!/
USE W3STRKMD
!/F90 USE W3TIMEMD, ONLY: TDIFF
IMPLICIT NONE
!/MPI
!/MPI INCLUDE "mpif.h"
!
! 1. Purpose :
!
! Perform spatial and temporal tracking of wave systems, based
! on spectral partition (bulletin) output.
!
! 2. Method :
!
! This is a controller program. It reads the input parameter file
! ww3_systrk.inp and calls subroutine waveTracking_NWS_V2 to
! perform the actual tracking procedure. Write output (fields and
! point output).
!
! 3. Parameters :
!
LOGICAL :: testout
PARAMETER (testout = .FALSE.)
CHARACTER :: filename*80, paramFile*32
REAL :: dirKnob, perKnob, hsKnob, wetPts, seedLat, &
seedLon, dirTimeKnob, tpTimeKnob, tint
REAL :: lonout(100), latout(100) !Increase dimension?
INTEGER :: maxGroup, ntint, noutp
!/F90 INTEGER :: CLKDT0(8),CLKDT1(8)
!/F90 REAL :: CLKFEL
TYPE(dat2d), POINTER :: wsdat(:)
TYPE(timsys), POINTER :: sysA(:)
INTEGER, POINTER :: maxSys(:)
!
! Local parameters.
! ----------------------------------------------------------------
! intype Int input Type of input (0 = from memory; 1 = from file)
! tmax Int input Value of maxTs to apply (1 or 2, used for model coupling)
! tcur Int input Index of current time step (1 or 2, used for model coupling)
! ulimGroup Int input Upper limit of number of wave systems to output
!
LOGICAL :: file_exists
CHARACTER :: inpstr*72
INTEGER :: intype, tmax, tcur, maxI, maxJ
INTEGER :: it, igrp, sysmatch, ind, ip
INTEGER :: i, j, leng, ulimGroup
REAL, ALLOCATABLE :: dum(:,:)
!/TRKNC REAL, ALLOCATABLE :: dum2nc(:,:,:,:)
!/TRKNC REAL, ALLOCATABLE :: hsprt_nc(:,:,:)
!/TRKNC REAL, ALLOCATABLE :: tpprt_nc(:,:,:)
!/TRKNC REAL, ALLOCATABLE :: dirprt_nc(:,:,:)
!/TRKNC REAL, ALLOCATABLE :: longitude_nc(:),latitude_nc(:)
!/TRKNC REAL, ALLOCATABLE :: lonprt_nc(:),latprt_nc(:)
INTEGER NTIME_NC
INTEGER :: outputType
LOGICAL :: outputCheck1,outputCheck2
DOUBLE PRECISION :: date1, date2, tstart, tend
REAL :: dlon, dlat, lonprt, latprt
REAL :: dt
REAL :: minlon, maxlon, minlat, maxlat
INTEGER :: mxcwt, mycwt
!/MPI INTEGER :: rank, nproc, ierr
!/MPI CHARACTER :: rankstr*4
! For point output (bilinear interpolation)
REAL :: hsprt(10),tpprt(10),dirprt(10)
REAL :: BL_hsprt(10),BR_hsprt(10),TR_hsprt(10),TL_hsprt(10), &
BL_tpprt(10),BR_tpprt(10),TR_tpprt(10),TL_tpprt(10), &
BL_dirprt(10),BR_dirprt(10),TR_dirprt(10),TL_dirprt(10)
REAL :: BL_dirx,BR_dirx,TR_dirx,TL_dirx, &
BL_diry,BR_diry,TR_diry,TL_diry
REAL :: BL_lonprt,BR_lonprt,TR_lonprt,TL_lonprt, &
BL_latprt,BR_latprt,TR_latprt,TL_latprt
REAL :: t, u, BL_W, BR_W, TR_W, TL_W
REAL :: PI
PARAMETER (PI = 3.1416)
!
! 4. Subroutines used :
!
! waveTracking_NWS_V2
!
! 5. Called by :
!
! None, stand-alone program.
!
! 6. Error messages :
!
! 7. Remarks :
!
! 8. Structure :
!
! Calls subroutine waveTracking_NWS_V2 in trackmd.95 - see that
! file for structure.
!
! 9. Switches :
!
! !/SHRD Switch for shared / distributed memory architecture.
! !/MPI Id.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
!
!/MPI/! Start of parallel region
!/MPI CALL MPI_INIT(ierr)
!/MPI
!/MPI CALL MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)
!/MPI CALL MPI_COMM_SIZE(MPI_COMM_WORLD, nproc, ierr)
!/MPI
! Open log file
!/MPI WRITE(rankstr,'(i4.4)') rank
!/MPI OPEN(unit=20,file='sys_log'//rankstr//'.ww3',status='unknown')
!/SHRD OPEN(unit=20,file='sys_log.ww3',status='unknown')
! Print code version
!/MPI IF (rank.EQ.0) THEN
WRITE(6,900)
!/MPI END IF
WRITE(20,900)
900 FORMAT (/15X,' *** WAVEWATCH III Wave system tracking *** '/ &
15X,'==============================================='/)
! Since this program reads the raw partitioning input from file,
! we set intype=1 or 2, and tmax and tcur to dummy values (not used).
intype = 2
! intype = 1
IF (intype.EQ.1) WRITE(6,*) &
'*** WARNING: partRes format input used!'
tmax = 0
tcur = 0
! Read input parameter file
!/MPI IF (rank.EQ.0) THEN
INQUIRE(FILE='ww3_systrk.inp', EXIST=file_exists)
IF (.NOT.file_exists) THEN
WRITE(20,2000)
WRITE(6,2000)
CALL ABORT
END IF
OPEN(unit=10,file='ww3_systrk.inp',status='old')
READ(10,'(A72)') inpstr
DO WHILE (inpstr(1:1).EQ.'$')
READ(10,'(A72)') inpstr
END DO
BACKSPACE(10)
READ(10,*) filename
READ(10,'(A72)') inpstr
DO WHILE (inpstr(1:1).EQ.'$')
READ(10,'(A72)') inpstr
END DO
BACKSPACE(10)
READ(10,*) date1, date2, dt, ntint
tstart = date1 + date2/1000000
READ(10,'(A72)') inpstr
DO WHILE (inpstr(1:1).EQ.'$')
READ(10,'(A72)') inpstr
END DO
BACKSPACE(10)
READ(10,*) outputType
!Check for correct outputType option:
IF (outputType.EQ.1) THEN
!ASCII output
ELSEIF (outputType.EQ.3) THEN
!NetCDF 3 - requrires !/TRKNC switch
outputCheck1 = .TRUE.
!/TRKNC outputCheck1 = .FALSE.
IF(outputCheck1) THEN
WRITE(6,993)
STOP
END IF
ELSEIF (outputType.EQ.4) THEN
!NetCDF 4 - requrires !/TRKNC and !/NC4 switch
outputCheck1 = .TRUE.
outputCheck2 = .TRUE.
!/TRKNC outputCheck1 = .FALSE.
!/NC4 outputCheck2 = .FALSE.
IF(outputCheck1.OR.outputCheck2) THEN
WRITE(6,994)
STOP
END IF
ELSE
!Not a valid outputType
WRITE(6,995) outputType
STOP
ENDIF
READ(10,'(A72)') inpstr
DO WHILE (inpstr(1:1).EQ.'$')
READ(10,'(A72)') inpstr
END DO
BACKSPACE(10)
READ(10,*) minlon, maxlon, mxcwt
READ(10,'(A72)') inpstr
DO WHILE (inpstr(1:1).EQ.'$')
READ(10,'(A72)') inpstr
END DO
BACKSPACE(10)
READ(10,*) minlat, maxlat, mycwt
READ(10,'(A72)') inpstr
DO WHILE (inpstr(1:1).EQ.'$')
READ(10,'(A72)') inpstr
END DO
BACKSPACE(10)
READ(10,*) dirKnob, perKnob, hsKnob, wetPts, &
dirTimeKnob, tpTimeKnob
READ(10,'(A72)') inpstr
DO WHILE (inpstr(1:1).EQ.'$')
READ(10,'(A72)') inpstr
END DO
BACKSPACE(10)
READ(10,*) seedLat, seedLon
READ(10,'(A72)') inpstr
DO WHILE (inpstr(1:1).EQ.'$')
READ(10,'(A72)') inpstr
END DO
BACKSPACE(10)
noutp = 1
lonout(:) = 9999.
latout(:) = 9999.
DO WHILE (.TRUE.)
READ(10,*) lonout(noutp),latout(noutp)
IF ((lonout(noutp).EQ.0.).AND.(latout(noutp).EQ.0.)) EXIT
noutp = noutp + 1
END DO
noutp = noutp - 1
CLOSE(10)
WRITE(20,*) 'Raw partition file = ',filename
WRITE(20,'(A,F15.6)') 'Start time = ',tstart
WRITE(20,*) 'dt = ',dt
WRITE(20,*) 'No. time levels = ',ntint
WRITE(20,'(A,2F7.2)') 'Domain limits: Longitude =',minlon, maxlon
WRITE(20,'(A,2F7.2)') ' Latitude =',minlat, maxlat
WRITE(20,*) 'No. increments: Long, Lat =',mxcwt, mycwt
WRITE(20,*) 'dirKnob, perKnob, hsKnob, wetPts, &
dirTimeKnob, tpTimeKnob, seedLat, seedLon ='
WRITE(20,'(8F6.2)') dirKnob, perKnob, hsKnob, wetPts, &
dirTimeKnob, tpTimeKnob, seedLat, seedLon
WRITE(20,*) 'No. output points =',noutp
DO i = 1,noutp
WRITE(20,*) lonout(i), latout(i)
END DO
INQUIRE(FILE=filename, EXIST=file_exists)
IF (.NOT.file_exists) THEN
WRITE(20,2200) filename
WRITE(6,2200) filename
CALL EXIT(1)
END IF
!/F90 CALL DATE_AND_TIME ( VALUES=CLKDT0 )
!/MPI END IF
!/MPI/! MPI communication block
!/MPI CALL MPI_BCAST(filename,80,MPI_CHARACTER,0,MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(tstart,1,MPI_DOUBLE_PRECISION,0, &
!/MPI MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(tend,1,MPI_DOUBLE_PRECISION,0, &
!/MPI MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(dt,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(ntint,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(minlon,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(maxlon,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(minlat,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(maxlat,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(mxcwt,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(mycwt,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(dirKnob,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(perKnob,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(hsKnob,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(wetPts,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(dirTimeKnob,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(tpTimeKnob,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(seedLon,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(seedLat,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(noutp,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(lonout,100,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_BCAST(latout,100,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI CALL MPI_Barrier(MPI_COMM_WORLD,IERR)
CALL 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 )
!/MPI IF (rank.EQ.0) THEN
!/F90 CALL DATE_AND_TIME ( VALUES=CLKDT1 )
!/F90 CLKFEL = TDIFF ( CLKDT0,CLKDT1 )
!/F90 WRITE (6,998) CLKFEL
WRITE (6,*) 'Final system output...'
! Set upper limit for wave systems to output (limited by AWIPS display)
ulimGroup = 9
!-----Output systems as plain text----------------------------------------
maxI = SIZE(wsdat(1)%lon,1)
maxJ = SIZE(wsdat(1)%lon,2)
dlon = wsdat(1)%lon(2,2)-wsdat(1)%lon(1,1)
dlat = wsdat(1)%lat(2,2)-wsdat(1)%lat(1,1)
WRITE(20,*) 'dlon, dlat =',dlon,dlat
!-----Final SYSTEM output: Coordinates
OPEN(unit=21,file='sys_coord.ww3', status='unknown')
WRITE(21,'(I6,69X,A)') maxJ,'Number of rows'
WRITE(21,'(I6,69X,A)') maxI,'Number of cols'
!/TRKNC ALLOCATE( longitude_nc(maxI) )
!/TRKNC ALLOCATE( latitude_nc(maxJ) )
WRITE(21,*) 'Longitude ='
DO j = maxJ,1,-1
DO i = 1,maxI
WRITE(21,'(F7.2)',ADVANCE='NO') wsdat(1)%lon(i,j)
!/TRKNC longitude_nc(i)=wsdat(1)%lon(i,1)
END DO
WRITE(21,'(A)',ADVANCE='YES') ''
END DO
WRITE(21,*) 'Latitude = '
DO j = maxJ,1,-1
DO i = 1,maxI
WRITE(21,'(F7.2)',ADVANCE='NO') wsdat(1)%lat(i,j)
!/TRKNC latitude_nc(j)=wsdat(1)%lat(1,j)
END DO
WRITE(21,'(A)',ADVANCE='YES') ''
END DO
CLOSE(21)
!-----Final SYSTEM output: hs
IF(outputType == 1) THEN
OPEN(unit=22,file='sys_hs.ww3', status='unknown')
WRITE(22,'(I6,69X,A)') maxJ,'Number of rows'
WRITE(22,'(I6,69X,A)') maxI,'Number of cols'
ENDIF
NTIME_NC=SIZE(sysA)
ALLOCATE( dum(maxI,maxJ) )
!/TRKNC IF(outputType == 3 .OR. outputType == 4) THEN
!/TRKNC ALLOCATE( dum2nc(maxI,maxJ,maxGroup,NTIME_NC) )
!/TRKNC ENDIF
DO it = 1,SIZE(sysA)
! Loop through identified groups, limiting the output in file to ulimGroup
IF(outputType == 1) THEN
WRITE(22,'(F15.6,60x,A)') wsdat(it)%date,'Time'
WRITE(22,'(I6,69x,A)') MIN(ulimGroup,maxGroup), &
'Tot number of systems'
ENDIF
DO igrp = 1,MIN(ulimGroup,maxGroup)
dum(1:maxI,1:maxJ) = 9999.00
! Find system with this group tag
sysmatch = 1
DO WHILE (sysmatch.LE.maxSys(it))
IF (sysA(it)%sys(sysmatch)%grp.EQ.igrp) EXIT
sysmatch = sysmatch+1
END DO
IF (sysmatch.LE.maxSys(it)) THEN
! Match found: fill the output matrix with this data
leng = sysA(it)%sys(sysmatch)%nPoints
DO ind = 1, leng
dum(sysA(it)%sys(sysmatch)%i(ind), &
sysA(it)%sys(sysmatch)%j(ind)) = &
sysA(it)%sys(sysmatch)%hs(ind)
END DO
ELSE
leng = 0
END IF
IF(outputType == 1) THEN
WRITE(22,'(I6,69x,A)') igrp,'System number'
WRITE(22,'(I6,69x,A)') leng,'Number of points in system'
DO J = maxJ,1,-1
DO i = 1,maxI
WRITE(22,'(F8.2)',ADVANCE='NO') dum(i,j)
END DO
WRITE(22,'(A)',ADVANCE='YES') ''
END DO
ELSE
!/TRKNC DO J = maxJ,1,-1
!/TRKNC DO i = 1,maxI
!/TRKNC dum2nc(i,j,igrp,it)=dum(i,j)
!/TRKNC END DO
!/TRKNC END DO
ENDIF
END DO
END DO
!/TRKNC IF(outputType == 3 .OR. outputType == 4 ) THEN
!/TRKNC call t2netcdf(longitude_nc,latitude_nc,dum2nc,maxI,maxJ,&
!/TRKNC maxGroup,date1,date2,dt,NTIME_NC,1,outputType)
!/TRKNC ENDIF
IF(outputType.EQ.1) CLOSE(22)
!-----Final SYSTEM output: tp
IF(outputType == 1) THEN
OPEN(unit=23,file='sys_tp.ww3',status='unknown')
WRITE(23,'(I6,69X,A)') maxJ,'Number of rows'
WRITE(23,'(I6,69X,A)') maxI,'Number of cols'
ENDIF
DO it = 1,SIZE(sysA)
! Loop through identified groups, limiting the output in file to ulimGroup
IF(outputType == 1) THEN
WRITE(23,'(F15.6,60x,A)') wsdat(it)%date,'Time'
WRITE(23,'(I6,69X,A)') MIN(ulimGroup,maxGroup), &
'Tot number of systems'
ENDIF
DO igrp = 1,MIN(ulimGroup,maxGroup)
dum(1:maxI,1:maxJ) = 9999.00
! Find system with this group tag
sysmatch = 1
DO WHILE (sysmatch.LE.maxSys(it))
IF (sysA(it)%sys(sysmatch)%grp.EQ.igrp) EXIT
sysmatch = sysmatch+1
END DO
IF (sysmatch.LE.maxSys(it)) THEN
! Match found: fill the output matrix with this data
leng = sysA(it)%sys(sysmatch)%nPoints
DO ind = 1, leng
dum(sysA(it)%sys(sysmatch)%i(ind), &
sysA(it)%sys(sysmatch)%j(ind)) = &
sysA(it)%sys(sysmatch)%tp(ind)
END DO
ELSE
leng = 0
END IF
IF(outputType == 1) THEN
WRITE(23,'(I6,69X,A)') igrp,'System number'
WRITE(23,'(I6,69X,A)') leng,'Number of points in system'
DO J = maxJ,1,-1
DO i = 1,maxI
WRITE(23,'(F8.2)',ADVANCE='NO') dum(i,j)
END DO
WRITE(23,'(A)',ADVANCE='YES') ''
END DO
ELSE
!/TRKNC DO J = maxJ,1,-1
!/TRKNC DO i = 1,maxI
!/TRKNC dum2nc(i,j,igrp,it)=dum(i,j)
!/TRKNC END DO
!/TRKNC END DO
ENDIF
END DO
END DO
!/TRKNC IF(outputType.EQ.3 .OR. outputType.EQ. 4 ) THEN
!/TRKNC call t2netcdf(longitude_nc,latitude_nc,dum2nc,maxI,maxJ,&
!/TRKNC maxGroup,date1,date2,dt,NTIME_NC,2,outputType)
!/TRKNC ENDIF
IF(outputType.EQ.1) CLOSE(23)
!-----Final SYSTEM output: dir
IF(outputType == 1) THEN
OPEN(unit=24,file='sys_dir.ww3',status='unknown')
WRITE(24,'(I6,69X,A)') maxJ,'Number of rows'
WRITE(24,'(I6,69X,A)') maxI,'Number of cols'
ENDIF
DO it = 1,SIZE(sysA)
! Loop through identified groups, limiting the output in file to
! ulimGroup
IF(outputType == 1) THEN
WRITE(24,'(F15.6,60x,A)') wsdat(it)%date,'Time'
WRITE(24,'(I6,69X,A)') MIN(ulimGroup,maxGroup), &
'Tot number of systems'
ENDIF
DO igrp = 1,MIN(ulimGroup,maxGroup)
dum(1:maxI,1:maxJ) = 9999.00
! Find system with this group tag
sysmatch = 1
DO WHILE (sysmatch.LE.maxSys(it))
IF (sysA(it)%sys(sysmatch)%grp.EQ.igrp) EXIT
sysmatch = sysmatch+1
END DO
IF (sysmatch.LE.maxSys(it)) THEN
! Match found: fill the output matrix with this data
leng = sysA(it)%sys(sysmatch)%nPoints
DO ind = 1, leng
dum(sysA(it)%sys(sysmatch)%i(ind), &
sysA(it)%sys(sysmatch)%j(ind)) = &
sysA(it)%sys(sysmatch)%dir(ind)
END DO
ELSE
leng = 0
END IF
IF(outputType == 1) THEN
WRITE(24,'(I6,69X,A)') igrp,'System number'
WRITE(24,'(I6,69X,A)') leng,'Number of points in system'
DO J = maxJ,1,-1
DO i = 1,maxI
WRITE(24,'(F8.2)',ADVANCE='NO') dum(i,j)
END DO
WRITE(24,'(A)',ADVANCE='YES') ''
END DO
ELSE
!/TRKNC DO J = maxJ,1,-1
!/TRKNC DO i = 1,maxI
!/TRKNC dum2nc(i,j,igrp,it)=dum(i,j)
!/TRKNC END DO
!/TRKNC END DO
END IF
END DO
END DO
!/TRKNC IF(outputType.EQ.3 .OR. outputType.EQ.4 ) THEN
!/TRKNC call t2netcdf(longitude_nc,latitude_nc,dum2nc,maxI,maxJ,&
!/TRKNC maxGroup,date1,date2,dt,NTIME_NC,3,outputType)
!/TRKNC ENDIF
IF(outputType.EQ.1) CLOSE(24)
!-----Final SYSTEM output: dspr
IF(outputType == 1) THEN
OPEN(unit=25,file='sys_dspr.ww3',status='unknown')
WRITE(25,'(I6,69X,A)') maxJ,'Number of rows'
WRITE(25,'(I6,69X,A)') maxI,'Number of cols'
ENDIF
DO it = 1,SIZE(sysA)
! Loop through identified groups, limiting the output in file to ulimGroup
IF(outputType == 1) THEN
WRITE(25,'(F15.6,60x,A)') wsdat(it)%date,'Time'
WRITE(25,'(I6,69X,A)') MIN(ulimGroup,maxGroup), &
'Tot number of systems'
ENDIF
DO igrp = 1,MIN(ulimGroup,maxGroup)
dum(1:maxI,1:maxJ) = 9999.00
! Find system with this group tag
sysmatch = 1
DO WHILE (sysmatch.LE.maxSys(it))
IF (sysA(it)%sys(sysmatch)%grp.EQ.igrp) EXIT
sysmatch = sysmatch+1
END DO
IF (sysmatch.LE.maxSys(it)) THEN
! Match found: fill the output matrix with this data
leng = sysA(it)%sys(sysmatch)%nPoints
DO ind = 1, leng
dum(sysA(it)%sys(sysmatch)%i(ind), &
sysA(it)%sys(sysmatch)%j(ind)) = &
sysA(it)%sys(sysmatch)%dspr(ind)
END DO
ELSE
leng = 0
END IF
IF(outputType == 1) THEN
WRITE(25,'(I6,69X,A)') igrp,'System number'
WRITE(25,'(I6,69X,A)') leng,'Number of points in system'
DO J = maxJ,1,-1
DO i = 1,maxI
WRITE(25,'(F8.2)',ADVANCE='NO') dum(i,j)
END DO
WRITE(25,'(A)',ADVANCE='YES') ''
END DO
ELSE
!/TRKNC DO J = maxJ,1,-1
!/TRKNC DO i = 1,maxI
!/TRKNC dum2nc(i,j,igrp,it)=dum(i,j)
!/TRKNC END DO
!/TRKNC END DO
ENDIF
END DO
END DO
!/TRKNC IF(outputType.EQ.3 .OR. outputType.EQ.4 ) THEN
!/TRKNC call t2netcdf(longitude_nc,latitude_nc,dum2nc,maxI,maxJ,&
!/TRKNC maxGroup,date1,date2,dt,NTIME_NC,4,outputType)
!/TRKNC ENDIF
IF(outputType.EQ.1) CLOSE(25)
IF (ALLOCATED(DUM)) DEALLOCATE(dum)
!/TRKNC IF (ALLOCATED(dum2nc)) DEALLOCATE(dum2nc)
!/TRKNC IF(outputType.EQ.3.OR.outputType.EQ.4) THEN
!/TRKNC ALLOCATE( hsprt_nc(10,noutp,NTIME_NC) )
!/TRKNC ALLOCATE( tpprt_nc(10,noutp,NTIME_NC) )
!/TRKNC ALLOCATE( dirprt_nc(10,noutp,NTIME_NC) )
!/TRKNC ALLOCATE( lonprt_nc(noutp) )
!/TRKNC ALLOCATE( latprt_nc(noutp) )
!/TRKNC ENDIF
!-----Final SYSTEM output: point output
IF(outputType == 1) THEN
OPEN(unit=26,file='sys_pnt.ww3',status='unknown')
WRITE(26,'(A)') '%'
WRITE(26,'(A)') '%'
WRITE(26,'(A)') '% WW3 Wave tracking point output'
WRITE(26,'(A)') '%'
WRITE(26,'(10A)') '% Xp Yp ', &
'HsSY01 HsSY02 HsSY03 HsSY04 ', &
'HsSY05 HsSY06 HsSY07 HsSY08 ', &
'HsSY09 HsSY10 ', &
'TpSY01 TpSY02 TpSY03 TpSY04 ', &
'TpSY05 TpSY06 TpSY07 TpSY08 ', &
'TpSY09 TpSY10 ', &
'DrSY01 DrSY02 DrSY03 DrSY04 ', &
'DrSY05 DrSY06 DrSY07 DrSY08 ', &
'DrSY09 DrSY10'
WRITE(26,'(10A)') '% [degr] [degr] ', &
'[m] [m] [m] [m] ', &
'[m] [m] [m] [m] ', &
'[m] [m] ', &
'[sec] [sec] [sec] [sec] ', &
'[sec] [sec] [sec] [sec] ', &
'[sec] [sec] ', &
'[degr] [degr] [degr] [degr] ', &
'[degr] [degr] [degr] [degr] ', &
'[degr] [degr]'
WRITE(26,'(A)') '%'
ENDIF
DO it = 1,SIZE(sysA)
IF(outputType == 1) THEN
WRITE(26,'(A,F15.6)') 'Time : ',wsdat(it)%date
ENDIF
DO ip = 1,noutp
hsprt(1:10) = 999.9999
tpprt(1:10) = 999.9999
dirprt(1:10) = 999.9999
lonprt = 999.9999
latprt = 999.9999
BL_hsprt(1:10) = 999.9999
BL_tpprt(1:10) = 999.9999
BL_dirprt(1:10) = 999.9999
BR_hsprt(1:10) = 999.9999
BR_tpprt(1:10) = 999.9999
BR_dirprt(1:10) = 999.9999
TL_hsprt(1:10) = 999.9999
TL_tpprt(1:10) = 999.9999
TL_dirprt(1:10) = 999.9999
TR_hsprt(1:10) = 999.9999
TR_tpprt(1:10) = 999.9999
TR_dirprt(1:10) = 999.9999
BL_lonprt = 999.9999
BL_latprt = 999.9999
BR_lonprt = 999.9999
BR_latprt = 999.9999
TL_lonprt = 999.9999
TL_latprt = 999.9999
TR_lonprt = 999.9999
TR_latprt = 999.9999
BL_W = 999
BR_W = 999
TR_W = 999
TL_W = 999
DO j = 1, (maxJ-1)
DO i = 1, (maxI-1)
IF ( ( ((lonout(ip).GE. &
wsdat(1)%lon(i,j)).AND. &
(lonout(ip).LT. &
wsdat(1)%lon(i+1,j))).OR. &
((lonout(ip).GT. &
wsdat(1)%lon(i,j)).AND. &
(lonout(ip).LE. &
wsdat(1)%lon(i+1,j))) ).AND. &
( ((latout(ip).GE. &
wsdat(1)%lat(i,j)).AND. &
(latout(ip).LT. &
wsdat(1)%lat(i,j+1))).OR. &
((latout(ip).GT. &
wsdat(1)%lat(i,j)).AND. &
(latout(ip).LE. &
wsdat(1)%lat(i,j+1))) ) ) &
THEN
BL_lonprt = wsdat(1)%lon(i,j)
BL_latprt = wsdat(1)%lat(i,j)
BR_lonprt = wsdat(1)%lon(i+1,j)
BR_latprt = wsdat(1)%lat(i+1,j)
TL_lonprt = wsdat(1)%lon(i,j+1)
TL_latprt = wsdat(1)%lat(i,j+1)
TR_lonprt = wsdat(1)%lon(i+1,j+1)
TR_latprt = wsdat(1)%lat(i+1,j+1)
! Compute weights for this point
t = (lonout(ip)-BL_lonprt)/(BR_lonprt-BL_lonprt)
u = (latout(ip)-BL_latprt)/(TL_latprt-BL_latprt)
BL_W = (1-t)*(1-u)
BR_W = t*(1-u)
TR_W = t*u
TL_W = (1-t)*u
! Compute output values using weights
lonprt = BL_W*BL_lonprt + BR_W*BR_lonprt + &
TL_W*TL_lonprt + TR_W*TR_lonprt
latprt = BL_W*BL_latprt + BR_W*BR_latprt + &
TL_W*TL_latprt + TR_W*TR_latprt
END IF
END DO
END DO
! Loop through identified groups, limiting the output in file to 10
DO igrp = 1,MIN(10,maxGroup)
! Find system with this group tag
sysmatch = 1
DO WHILE (sysmatch.LE.maxSys(it))
IF (sysA(it)%sys(sysmatch)%grp.EQ.igrp) EXIT
sysmatch = sysmatch+1
END DO
IF (sysmatch.LE.maxSys(it)) THEN
! Match found: fill the output matrix with this data
leng = sysA(it)%sys(sysmatch)%nPoints
DO ind = 1, leng
! Write output point data with bilinear interpolation
IF ( (sysA(it)%sys(sysmatch)%lon(ind).EQ.&
BL_lonprt).AND.&
(sysA(it)%sys(sysmatch)%lat(ind).EQ.&
BL_latprt) ) THEN
BL_hsprt(igrp) = sysA(it)%sys(sysmatch)%hs(ind)
BL_tpprt(igrp) = sysA(it)%sys(sysmatch)%tp(ind)
BL_dirprt(igrp) = sysA(it)%sys(sysmatch)%dir(ind)
ELSE IF ( (sysA(it)%sys(sysmatch)%lon(ind).EQ.&
BR_lonprt).AND.&
(sysA(it)%sys(sysmatch)%lat(ind).EQ.&
BR_latprt)) THEN
BR_hsprt(igrp) = sysA(it)%sys(sysmatch)%hs(ind)
BR_tpprt(igrp) = sysA(it)%sys(sysmatch)%tp(ind)
BR_dirprt(igrp) = sysA(it)%sys(sysmatch)%dir(ind)
ELSE IF ( (sysA(it)%sys(sysmatch)%lon(ind).EQ.&
TL_lonprt).AND.&
(sysA(it)%sys(sysmatch)%lat(ind).EQ.&
TL_latprt)) THEN
TL_hsprt(igrp) = sysA(it)%sys(sysmatch)%hs(ind)
TL_tpprt(igrp) = sysA(it)%sys(sysmatch)%tp(ind)
TL_dirprt(igrp) = sysA(it)%sys(sysmatch)%dir(ind)
ELSE IF ( (sysA(it)%sys(sysmatch)%lon(ind).EQ.&
TR_lonprt).AND.&
(sysA(it)%sys(sysmatch)%lat(ind).EQ.&
TR_latprt)) THEN
TR_hsprt(igrp) = sysA(it)%sys(sysmatch)%hs(ind)
TR_tpprt(igrp) = sysA(it)%sys(sysmatch)%tp(ind)
TR_dirprt(igrp) = sysA(it)%sys(sysmatch)%dir(ind)
END IF
END DO
! Compute output value using weights
! (only if output point is surrounded by valid points)
IF ( (BL_hsprt(igrp).NE.999.9999).AND. &
(BR_hsprt(igrp).NE.999.9999).AND. &
(TL_hsprt(igrp).NE.999.9999).AND. &
(TR_hsprt(igrp).NE.999.9999) ) THEN
hsprt(igrp) = BL_W * BL_hsprt(igrp) + &
BR_W * BR_hsprt(igrp) + &
TL_W * TL_hsprt(igrp) + &
TR_W * TR_hsprt(igrp)
tpprt(igrp) = BL_W * BL_tpprt(igrp) + &
BR_W * BR_tpprt(igrp) + &
TL_W * TL_tpprt(igrp) + &
TR_W * TR_tpprt(igrp)
BL_dirx = COS((270-BL_dirprt(igrp))*PI/180.)
BR_dirx = COS((270-BR_dirprt(igrp))*PI/180.)
TR_dirx = COS((270-TR_dirprt(igrp))*PI/180.)
TL_dirx = COS((270-TL_dirprt(igrp))*PI/180.)
BL_diry = SIN((270-BL_dirprt(igrp))*PI/180.)
BR_diry = SIN((270-BR_dirprt(igrp))*PI/180.)
TR_diry = SIN((270-TR_dirprt(igrp))*PI/180.)
TL_diry = SIN((270-TL_dirprt(igrp))*PI/180.)
dirprt(igrp)=270 - 180./PI* &
ATAN2(BL_W*BL_diry+BR_W*BR_diry+ &
TL_W*TL_diry+TR_W*TR_diry, &
BL_W*BL_dirx+BR_W*BR_dirx+ &
TL_W*TL_dirx+TR_W*TR_dirx)
IF (dirprt(igrp).GT.360.) THEN
dirprt(igrp) = dirprt(igrp) - 360.
END IF
ELSE
hsprt(igrp) = 999.9999
tpprt(igrp) = 999.9999
dirprt(igrp) = 999.9999
END IF
END IF
END DO
IF(outputType == 1) THEN
WRITE(26,'(32F14.4)') lonprt,latprt, &
hsprt(1:10),tpprt(1:10),dirprt(1:10)
ENDIF
!/TRKNC IF(outputType.EQ.3.OR.outputType.EQ.4) THEN
!/TRKNC lonprt_nc(ip)=lonprt
!/TRKNC latprt_nc(ip)=latprt
!/TRKNC do igrp=1,10
!/TRKNC hsprt_nc(igrp,ip,it)=hsprt(igrp)
!/TRKNC tpprt_nc(igrp,ip,it)=tpprt(igrp)
!/TRKNC dirprt_nc(igrp,ip,it)=dirprt(igrp)
!/TRKNC enddo
!/TRKNC ENDIF
END DO
END DO
!/TRKNC IF(outputType.EQ.3.OR.outputType.EQ.4) THEN
!/TRKNC call pt2netcdf(lonprt_nc,latprt_nc,hsprt_nc,tpprt_nc, &
!/TRKNC dirprt_nc,noutp,date1,date2,dt,NTIME_NC,outputType)
!/TRKNC ENDIF
IF(outputType.EQ.1) CLOSE(26)
!-----Final SYSTEM output: point output (Nearest neighbor, as a double check)
IF (testout) THEN
OPEN(unit=28,file='sys_pnt_nn.ww3',status='unknown')
WRITE(28,'(A)') '%'
WRITE(28,'(A)') '%'
WRITE(28,'(A)') '% WW3 Wave tracking point output'
WRITE(28,'(A)') '%'
WRITE(28,'(10A)') '% Xp Yp ', &
'HsSY01 HsSY02 HsSY03 HsSY04 ', &
'HsSY05 HsSY06 HsSY07 HsSY08 ', &
'HsSY09 HsSY10 ', &
'TpSY01 TpSY02 TpSY03 TpSY04 ', &
'TpSY05 TpSY06 TpSY07 TpSY08 ', &
'TpSY09 TpSY10 ', &
'DrSY01 DrSY02 DrSY03 DrSY04 ', &
'DrSY05 DrSY06 DrSY07 DrSY08 ', &
'DrSY09 DrSY10'
WRITE(28,'(10A)') '% [degr] [degr] ', &
'[m] [m] [m] [m] ', &
'[m] [m] [m] [m] ', &
'[m] [m] ', &
'[sec] [sec] [sec] [sec] ', &
'[sec] [sec] [sec] [sec] ', &
'[sec] [sec] ', &
'[degr] [degr] [degr] [degr] ', &
'[degr] [degr] [degr] [degr] ', &
'[degr] [degr]'
WRITE(28,'(A)') '%'
DO it = 1,SIZE(sysA)
WRITE(28,'(A,F15.6)') 'Time : ',wsdat(it)%date
DO ip = 1,noutp
hsprt(1:10) = 999.9999
tpprt(1:10) = 999.9999
dirprt(1:10) = 999.9999
lonprt = 999.9999
latprt = 999.9999
DO j = 1, maxJ
DO i = 1, maxI
! Write nearest nearbor output (no bilinear interpolation)
IF ( (lonout(ip).GE. &
(wsdat(1)%lon(i,j)-dlon/2)).AND. &
(lonout(ip).LT. &
(wsdat(1)%lon(i,j)+dlon/2)).AND. &
(latout(ip).GE. &
(wsdat(1)%lat(i,j)-dlat/2)).AND. &
(latout(ip).LT. &
(wsdat(1)%lat(i,j)+dlat/2)) ) &
THEN
lonprt = wsdat(1)%lon(i,j)
latprt = wsdat(1)%lat(i,j)
END IF
END DO
END DO
! Loop through identified groups, limiting the output in file to 10
DO igrp = 1,MIN(10,maxGroup)
! Find system with this group tag
sysmatch = 1
DO WHILE (sysmatch.LE.maxSys(it))
IF (sysA(it)%sys(sysmatch)%grp.EQ.igrp) EXIT
sysmatch = sysmatch+1
END DO
IF (sysmatch.LE.maxSys(it)) THEN
! Match found: fill the output matrix with this data
leng = sysA(it)%sys(sysmatch)%nPoints
DO ind = 1, leng
! Write nearest nearbor output (no bilinear interpolation)
IF ( (lonout(ip).GE. &
(sysA(it)%sys(sysmatch)%lon(ind)-dlon/2)).AND. &
(lonout(ip).LT. &
(sysA(it)%sys(sysmatch)%lon(ind)+dlon/2)).AND. &
(latout(ip).GE. &
(sysA(it)%sys(sysmatch)%lat(ind)-dlat/2)).AND. &
(latout(ip).LT. &
(sysA(it)%sys(sysmatch)%lat(ind)+dlat/2)) ) &
THEN
hsprt(igrp) = sysA(it)%sys(sysmatch)%hs(ind)
tpprt(igrp) = sysA(it)%sys(sysmatch)%tp(ind)
dirprt(igrp) = sysA(it)%sys(sysmatch)%dir(ind)
END IF
END DO
END IF
END DO
WRITE(28,'(32F14.4)') lonprt,latprt, &
hsprt(1:10),tpprt(1:10),dirprt(1:10)
END DO
END DO
CLOSE(28)
END IF
!-------------------------------------------------------------------------
WRITE(20,*) 'In ww3_systrk: Deallocating wsdat ...'
DO it=1,size(wsdat)
IF (ASSOCIATED(wsdat(it)%lat)) DEALLOCATE(wsdat(it)%lat)
IF (ASSOCIATED(wsdat(it)%lon)) DEALLOCATE(wsdat(it)%lon)
IF (ASSOCIATED(wsdat(it)%par)) DEALLOCATE(wsdat(it)%par)
IF (ASSOCIATED(wsdat(it)%wnd)) DEALLOCATE(wsdat(it)%wnd)
END DO
IF (ASSOCIATED(wsdat)) DEALLOCATE(wsdat)
WRITE(20,*) ' Deallocating sysA ...'
DO i=1,size(sysA)
DO j=1,size(sysA(i)%sys)
IF (ASSOCIATED(sysA(i)%sys(j)%i)) DEALLOCATE(sysA(i)%sys(j)%i)
IF (ASSOCIATED(sysA(i)%sys(j)%j)) DEALLOCATE(sysA(i)%sys(j)%j)
IF (ASSOCIATED(sysA(i)%sys(j)%lon)) &
DEALLOCATE(sysA(i)%sys(j)%lon)
IF (ASSOCIATED(sysA(i)%sys(j)%lat)) &
DEALLOCATE(sysA(i)%sys(j)%lat)
IF (ASSOCIATED(sysA(i)%sys(j)%hs)) &
DEALLOCATE(sysA(i)%sys(j)%hs)
IF (ASSOCIATED(sysA(i)%sys(j)%tp)) &
DEALLOCATE(sysA(i)%sys(j)%tp)
IF (ASSOCIATED(sysA(i)%sys(j)%dir)) &
DEALLOCATE(sysA(i)%sys(j)%dir)
IF (ASSOCIATED(sysA(i)%sys(j)%dspr)) &
DEALLOCATE(sysA(i)%sys(j)%dspr)
END DO
END DO
IF (ASSOCIATED(sysA)) DEALLOCATE(sysA)
WRITE(20,*) ' Deallocating maxSys ...'
IF (ASSOCIATED(maxSys)) DEALLOCATE(maxSys)
CLOSE(20)
WRITE(6,*) '... ww3_systrk completed successfully.'
WRITE(6,999)
!/MPI END IF !/IF (rank.EQ.0)
!/MPI CALL MPI_FINALIZE(IERR)
!/MPI/! End of parallel region
!/F90 998 FORMAT ( ' ... finished. Elapsed time : ',F10.2,' s')
993 FORMAT (/' *** WAVEWATCH III ERROR IN WW3_SYSTRK : '/ &
' OutputType=3 needs TRKNC switch ')
994 FORMAT (/' *** WAVEWATCH III ERROR IN WW3_SYSTRK : '/ &
' OutputType=4 needs TRKNC and NC4 switch ')
995 FORMAT (/' *** WAVEWATCH III ERROR IN WW3_SYSTRK : '/ &
' OutputType,',I3,'not valid. Options: 1,3,4')
999 FORMAT (/15X,'End of program '/ &
15X,'==============================================='/ &
15X,' *** WAVEWATCH III Wave system tracking *** ')
2000 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/ &
' ERROR IN OPENING INPUT FILE')
2200 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/ &
' ERROR IN OPENING PARTITION FILE : ',A)
END PROGRAM WW3_SYSTRK
!
!/TRKNC subroutine t2netcdf(lons,lats,data_in,nlons,nlats,nsys,date1,date2,&
!/TRKNC dt,ntime,ivar, outputType)
!/TRKNC USE W3TIMEMD
!/TRKNC use netcdf
!/TRKNC implicit none
!/TRKNC character (len = 15) :: file_name
!/TRKNC integer, parameter :: ndims = 4
!/TRKNC integer, parameter :: deflate = 1
!/TRKNC integer :: outputType, ncid, oldMode
!/TRKNC integer :: nlons,nlats,nsys,rec,ntime,ivar
!/TRKNC double precision :: date1,date2,timenc
!/TRKNC real :: data_in(nlons, nlats, nsys,ntime)
!/TRKNC real :: lats(nlats), lons(nlons),dt
!/TRKNC double precision :: times(ntime)
!/TRKNC integer :: iyc,imc,idc,ihc,iminc,isc,Jday,Jday0
!/TRKNC integer :: iret
!
!/TRKNC integer :: lon_varid, lat_varid, rec_varid
!/TRKNC character (len = *), parameter :: lsys_name = "system_index"
!/TRKNC character (len = *), parameter :: lat_name = "latitude"
!/TRKNC character (len = *), parameter :: lon_name = "longitude"
!/TRKNC character (len = *), parameter :: time_name = "time"
!/TRKNC integer :: sys_dimid, lon_dimid, lat_dimid, rec_dimid
!/TRKNC integer :: start(ndims), count(ndims)
!
!/TRKNC character (len = *), parameter :: var1_name="hs"
!/TRKNC character (len = *), parameter :: var2_name="tp"
!/TRKNC character (len = *), parameter :: var3_name="dir"
!/TRKNC character (len = *), parameter :: var4_name="dspr"
!/TRKNC integer :: var1_varid, var2_varid, var3_varid,var4_varid
!/TRKNC integer :: dimids(ndims)
!
!/TRKNC character (len = *), parameter :: units = "units"
!/TRKNC character (len = *), parameter :: var1_units = "m"
!/TRKNC character (len = *), parameter :: var2_units = "s"
!/TRKNC character (len = *), parameter :: var3_units = "degrees"
!/TRKNC character (len = *), parameter :: var4_units = "degrees"
!/TRKNC character (len = *), parameter :: lat_units = "degrees_north"
!/TRKNC character (len = *), parameter :: lon_units = "degrees_east"
!/TRKNC iyc=date1/10000
!/TRKNC imc=(date1-iyc*10000)/100
!/TRKNC idc=int(date1-DBLE(iyc*10000)-DBLE(imc*100))
!/TRKNC ihc=date2/10000
!/TRKNC iminc=(date2-ihc*10000)/100
!/TRKNC isc=date2-ihc*10000-100*iminc
!/TRKNC timenc=DBLE(julday(idc,imc,iyc))+(DBLE(ihc)+(DBLE(iminc)+ &
!/TRKNC (DBLE(isc)/60.0D0))/60.0D0)/24.0D0
!/TRKNC Jday0=julday(1,1,1990)
!/TRKNC timenc=timenc-Jday0
!/TRKNC do rec=1,ntime
!/TRKNC times(rec)=timenc+DBLE( (rec-1)*dt)/3600.0D0/24.0D0
!/TRKNC enddo
!/TRKNC if( ivar == 1) then
!/TRKNC file_name = "sys_hs.ww3.nc"
!/TRKNC else if( ivar == 2) then
!/TRKNC file_name = "sys_tp.ww3.nc"
!/TRKNC else if( ivar == 3) then
!/TRKNC file_name = "sys_dir.ww3.nc"
!/TRKNC else
!/TRKNC file_name = "sys_dspr.ww3.nc"
!/TRKNC endif
!
!/TRKNC! create the netcdf file.
!/TRKNC if (outputType.EQ.3) then
!/TRKNC call check( nf90_create(file_name, NF90_CLOBBER, ncid) )
!/TRKNC endif
!/NC4 if(outputType.EQ.4) call check( nf90_create(file_name,NF90_NETCDF4,ncid))
!/TRKNC call check ( nf90_set_fill(ncid,nf90_nofill,oldMode) )
!/TRKNC call check( nf90_def_dim(ncid, lsys_name, nsys, sys_dimid) )
!/TRKNC call check( nf90_def_dim(ncid, lat_name, nlats, lat_dimid) )
!/TRKNC call check( nf90_def_dim(ncid, lon_name, nlons, lon_dimid) )
!/TRKNC call check( nf90_def_dim(ncid, time_name, ntime, rec_dimid) )
!/TRKNC call check( nf90_def_var(ncid, lat_name, NF90_REAL, lat_dimid,lat_varid))
!/NC4 call check( nf90_def_var_deflate(ncid,lat_varid,1,1,deflate) )
!/TRKNC call check( nf90_def_var(ncid, lon_name, NF90_REAL, lon_dimid,lon_varid))
!/NC4 call check( nf90_def_var_deflate(ncid,lon_varid,1,1,deflate) )
!/TRKNC call check( nf90_def_var(ncid,time_name,NF90_DOUBLE,rec_dimid,rec_varid))
!/NC4 call check( nf90_def_var_deflate(ncid,rec_varid,1,1,deflate) )
!
!/TRKNC call check( nf90_put_att(ncid, lat_varid, units, lat_units) )
!/TRKNC call check( nf90_put_att(ncid, lat_varid, 'long_name', 'latitude') )
!/TRKNC call check( nf90_put_att(ncid, lat_varid, 'standard_name', 'latitude') )
!/TRKNC call check( nf90_put_att(ncid, lat_varid, 'axis','Y'))
!/TRKNC call check( nf90_put_att(ncid, lon_varid, units, lon_units) )
!/TRKNC call check( nf90_put_att(ncid, lon_varid, 'long_name', 'longitude') )
!/TRKNC call check( nf90_put_att(ncid, lon_varid, 'standard_name', 'longitude') )
!/TRKNC call check( nf90_put_att(ncid, lon_varid, 'axis','X'))
!/TRKNC call check(nf90_put_att(ncid,rec_varid,units,&
!/TRKNC 'days since 1990-01-01 00:00:00'))
!/TRKNC call check(nf90_put_att(ncid,rec_varid,'long_name','julian day (UT)'))
!/TRKNC call check( nf90_put_att(ncid, rec_varid,'standard_name', 'time') )
!/TRKNC call check( nf90_put_att(ncid, rec_varid, 'conventions',&
!/TRKNC 'relative julian day with decimal part (as part of the day)' ) )
!/TRKNC call check( nf90_put_att(ncid, rec_varid, 'axis','T'))
!
!/TRKNC dimids = (/ lon_dimid, lat_dimid, sys_dimid, rec_dimid /)
!/TRKNC if( ivar == 1) then
!/TRKNC call check( nf90_def_var(ncid, var1_name, NF90_REAL, dimids,var1_varid) )
!/NC4 call check( nf90_def_var_deflate(ncid,var1_varid,1,1,deflate) )
!/TRKNC call check( nf90_put_att(ncid, var1_varid, units, var1_units) )
!/TRKNC call check( nf90_put_att(ncid, var1_varid,'long_name','significant_wave_height') )
!/TRKNC call check( nf90_put_att(ncid, var1_varid,'missing_value','9999.00'))
!/TRKNC else if( ivar == 2) then
!/TRKNC call check( nf90_def_var(ncid, var2_name, NF90_REAL, dimids, var2_varid) )
!/NC4 call check( nf90_def_var_deflate(ncid,var2_varid,1,1,deflate) )
!/TRKNC call check( nf90_put_att(ncid, var2_varid, units, var2_units) )
!/TRKNC call check( nf90_put_att(ncid, var2_varid,'long_name','peak_period') )
!/TRKNC call check( nf90_put_att(ncid, var2_varid,'missing_value','9999.00') )
!/TRKNC else if ( ivar ==3 ) then
!/TRKNC call check( nf90_def_var(ncid, var3_name, NF90_REAL, dimids, var3_varid) )
!/NC4 call check( nf90_def_var_deflate(ncid,var3_varid,1,1,deflate) )
!/TRKNC call check( nf90_put_att(ncid, var3_varid, units, var3_units) )
!/TRKNC call check( nf90_put_att(ncid, var3_varid,'long_name','peak_direction') )
!/TRKNC call check( nf90_put_att(ncid, var3_varid,'missing_value','9999.00') )
!/TRKNC else
!/TRKNC call check( nf90_def_var(ncid, var4_name, NF90_REAL, dimids, var4_varid) )
!/NC4 call check( nf90_def_var_deflate(ncid,var4_varid,1,1,deflate) )
!/TRKNC call check( nf90_put_att(ncid, var4_varid, units, var4_units) )
!/TRKNC call check( nf90_put_att(ncid,var4_varid,'long_name','directional_spread') )
!/TRKNC call check( nf90_put_att(ncid, var4_varid,'missing_value','9999.00') )
!/TRKNC endif
!/TRKNC call check( nf90_enddef(ncid) )
!
!/TRKNC call check( nf90_put_var(ncid, lat_varid, lats) )
!/TRKNC call check( nf90_put_var(ncid, lon_varid, lons) )
!/TRKNC call check( nf90_put_var(ncid, rec_varid, times) )
!
!/TRKNC count = (/ nlons, nlats, nsys, ntime /)
!/TRKNC start = (/ 1, 1, 1, 1 /)
!/TRKNC if( ivar == 1) then
!/TRKNC call check( nf90_put_var(ncid, var1_varid, data_in, start = start, &
!/TRKNC count = count) )
!/TRKNC else if( ivar == 2) then
!/TRKNC call check( nf90_put_var(ncid, var2_varid, data_in, start = start, &
!/TRKNC count = count) )
!/TRKNC else if( ivar == 3) then
!/TRKNC call check( nf90_put_var(ncid, var3_varid, data_in, start = start, &
!/TRKNC count = count) )
!/TRKNC else
!/TRKNC call check( nf90_put_var(ncid, var4_varid, data_in, start = start, &
!/TRKNC count = count) )
!/TRKNC endif
!/TRKNC call check( nf90_close(ncid) )
!/TRKNC end subroutine t2netcdf
!
!/TRKNC subroutine check(status)
!/TRKNC use netcdf
!/TRKNC integer, intent ( in) :: status
!/TRKNC if(status /= nf90_noerr) then
!/TRKNC write(6,996)
!/TRKNC 996 FORMAT (/' *** WAVEWATCH III ERROR IN WW3_SYSTRK:'/ &
!/TRKNC 'netCDF error:')
!/TRKNC print *, trim(nf90_strerror(status))
!/TRKNC stop "Stopped in netcdf output part"
!/TRKNC endif
!/TRKNC end subroutine check
!
!/TRKNC subroutine pt2netcdf(longitude,latitude,hs,tp,&
!/TRKNC dir,npoints,date1,date2,dt,ntime,outputType)
!/TRKNC USE W3TIMEMD
!/TRKNC use netcdf
!/TRKNC implicit none
!/TRKNC integer :: ntime,npoints,outputType
!/TRKNC integer, parameter :: deflate = 1
!/TRKNC integer :: iret, oldMode
!/TRKNC integer :: ncid
!/TRKNC integer :: system_index_dim
!/TRKNC integer :: point_dim,rec_dim
!/TRKNC integer :: nsys
!/TRKNC integer :: start(3), count(3)
!/TRKNC parameter (nsys = 10)
!/TRKNC integer :: latitude_id
!/TRKNC integer :: longitude_id
!/TRKNC integer :: time_id
!/TRKNC integer :: hs_id
!/TRKNC integer :: tp_id
!/TRKNC integer :: dir_id
!/TRKNC integer :: time_rank
!/TRKNC integer :: hs_rank
!/TRKNC integer :: tp_rank
!/TRKNC integer :: dir_rank
!/TRKNC parameter (time_rank = 1)
!/TRKNC parameter (hs_rank = 3)
!/TRKNC parameter (tp_rank = 3)
!/TRKNC parameter (dir_rank = 3)
!
!/TRKNC integer :: hs_dims(hs_rank)
!/TRKNC integer :: tp_dims(tp_rank)
!/TRKNC integer :: dir_dims(dir_rank)
!/TRKNC real :: latitude(npoints),dt
!/TRKNC real :: longitude(npoints)
!/TRKNC real :: hs(nsys, npoints, ntime)
!/TRKNC real :: tp(nsys, npoints, ntime)
!/TRKNC real :: dir(nsys, npoints, ntime)
!/TRKNC integer :: iyc,imc,idc,ihc,iminc,isc,Jday,Jday0,rec
!/TRKNC double precision date1,date2,timenc
!/TRKNC double precision times(ntime)
!
!/TRKNC iyc=date1/10000
!/TRKNC imc=(date1-iyc*10000)/100
!/TRKNC idc=int(date1-DBLE(iyc*10000)-DBLE(imc*100))
!/TRKNC ihc=date2/10000
!/TRKNC iminc=(date2-ihc*10000)/100
!/TRKNC isc=date2-ihc*10000-100*iminc
!/TRKNC timenc=DBLE(julday(idc,imc,iyc))+(DBLE(ihc)+(DBLE(iminc)+ &
!/TRKNC (DBLE(isc)/60.0D0))/60.0D0)/24.0D0
!/TRKNC Jday0=julday(1,1,1990)
!/TRKNC timenc=timenc-Jday0
!/TRKNC do rec=1,ntime
!/TRKNC times(rec)=timenc+DBLE( (rec-1)*dt)/3600.0D0/24.0D0
!/TRKNC enddo
!
!/TRKNC if(outputType.EQ.3) then
!/TRKNC iret = nf90_create('sys_pnt.ww3.nc', NF90_CLOBBER, ncid)
!/TRKNC endif
!/TRKNC if (outputType.EQ.4) iret = nf90_create('sys_pnt.ww3.nc',NF90_NETCDF4, ncid)
!/TRKNC call check(iret)
!/TRKNC iret = nf90_set_fill(ncid,nf90_nofill,oldMode)
!/TRKNC call check(iret)
!/TRKNC! define dimensions
!/TRKNC iret = nf90_def_dim(ncid, 'system_index', nsys, system_index_dim)
!/TRKNC call check(iret)
!/TRKNC iret = nf90_def_dim(ncid, 'point', npoints, point_dim)
!/TRKNC call check(iret)
!/TRKNC iret = nf90_def_dim(ncid, 'time', ntime, rec_dim)
!/TRKNC call check(iret)
!/TRKNC! define variables
!/TRKNC iret = nf90_def_var(ncid, 'latitude', NF90_REAL, point_dim, &
!/TRKNC latitude_id)
!/TRKNC call check(iret)
!/NC4 call check( nf90_def_var_deflate(ncid,latitude_id,1,1,deflate))
!/TRKNC iret = nf90_def_var(ncid, 'longitude', NF90_REAL, point_dim, &
!/TRKNC longitude_id)
!/TRKNC call check(iret)
!/NC4 call check( nf90_def_var_deflate(ncid,longitude_id,1,1,deflate))
!/TRKNC iret = nf90_def_var(ncid, 'time', NF90_DOUBLE, rec_dim, &
!/TRKNC time_id)
!/TRKNC call check(iret)
!/NC4 call check( nf90_def_var_deflate(ncid,time_id,1,1,deflate) )
!/TRKNC hs_dims(3) = rec_dim
!/TRKNC hs_dims(2) = point_dim
!/TRKNC hs_dims(1) = system_index_dim
!/TRKNC iret = nf90_def_var(ncid, 'hs', NF90_REAL, &
!/TRKNC hs_dims, hs_id)
!/TRKNC call check(iret)
!/NC4 call check( nf90_def_var_deflate(ncid,hs_id,1,1,deflate))
!/TRKNC tp_dims(3) = rec_dim
!/TRKNC tp_dims(2) = point_dim
!/TRKNC tp_dims(1) = system_index_dim
!/TRKNC iret = nf90_def_var(ncid, 'tp', NF90_REAL, &
!/TRKNC tp_dims, tp_id)
!/TRKNC call check(iret)
!/NC4 call check( nf90_def_var_deflate(ncid,tp_id,1,1,deflate))
!/TRKNC dir_dims(3) = rec_dim
!/TRKNC dir_dims(2) = point_dim
!/TRKNC dir_dims(1) = system_index_dim
!/TRKNC iret = nf90_def_var(ncid, 'dir', NF90_REAL, &
!/TRKNC dir_dims, dir_id)
!/TRKNC call check(iret)
!/NC4 call check( nf90_def_var_deflate(ncid,dir_id,1,1,deflate))
!/TRKNC! assign attributes
!/TRKNC iret = nf90_put_att(ncid, latitude_id, 'units', 'degrees_north')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, latitude_id, 'long_name', 'latitude')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, latitude_id, 'standard_name', 'latitude')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, latitude_id, 'axis', 'Y')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, longitude_id, 'units', 'degrees_east')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, longitude_id,'long_name','longitude')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, longitude_id,'standard_name','longitude')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, longitude_id, 'axis', 'X')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, time_id, 'units', &
!/TRKNC 'days since 1990-01-01 00:00:00')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, time_id, 'long_name','julian day(UT)')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, time_id, 'standard_name','time')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, time_id, 'conventions', &
!/TRKNC 'relative julian day with decimal part (as part of the day)')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, time_id, 'axis', 'T')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, hs_id, 'units', 'm')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, hs_id,'long_name','significant_wave_height')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, hs_id, 'missing_value', &
!/TRKNC '999.9999')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, tp_id, 'units', 's')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, tp_id,'long_name','peak_period')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, tp_id, 'missing_value', &
!/TRKNC '999.9999')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, dir_id, 'units', 'degrees')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, dir_id,'long_name','peak_direction')
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_att(ncid, dir_id, 'missing_value',&
!/TRKNC '999.9999')
!/TRKNC call check(iret)
!/TRKNC! leave define mode
!/TRKNC iret = nf90_enddef(ncid)
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_var(ncid, latitude_id, latitude)
!/TRKNC call check(iret)
!
!/TRKNC iret = nf90_put_var(ncid, longitude_id, longitude)
!/TRKNC call check(iret)
!
!/TRKNC iret = nf90_put_var(ncid, time_id, times)
!/TRKNC call check(iret)
!
!/TRKNC start = (/ 1, 1, 1 /)
!/TRKNC count = (/ nsys,npoints,ntime /)
!
!/TRKNC iret = nf90_put_var(ncid, hs_id, hs,&
!/TRKNC start = start, count = count )
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_var(ncid, tp_id, tp, &
!/TRKNC start = start, count = count )
!/TRKNC call check(iret)
!/TRKNC iret = nf90_put_var(ncid, dir_id, dir,&
!/TRKNC start = start, count = count )
!/TRKNC call check(iret)
!/TRKNC iret = nf90_close(ncid)
!/TRKNC call check(iret)
!/TRKNC return
!/TRKNC end subroutine pt2netcdf