C/ ------------------------------------------------------------------- /
PROGRAM W3SPLT
C/
C/ +-----------------------------------+
C/ | WAVEWATCH III NOAA/NCEP |
C/ | H. L. Tolman |
C/ | FORTRAN 77 |
C/ | Last update : 19-May-1999 |
C/ +-----------------------------------+
C/
C 1. Purpose :
C
C Read a WAVEWATCH III version 1.17 point output data file and
C produces a table of mean parameters for all individual wave
C systems.
C
C 2. Method :
C
C Partitioning as devised by Gerling and similar to used in
C WAM model. From standard input the following data is read :
C
C name of output point to be considered.
C name of run of model.
C name of input file.
C logical to indentify FORMATTED input file.
C name of output file.
C
C All strings are read as characters, therefore no quotes are
C needed.
C
C 3. Parameters :
C
C Parameter statements
C ----------------------------------------------------------------
C NFR Int. Number of frequencies in spectrum.
C NTH Int. Number of directions in spectrum.
C NPMAX Int. Maximum number of peaks to be looked for.
C NPTAB Int. Number of columns in table.
C HSMIN Real Minimum wave height for includion in table.
C HSDROP Real Minimum wave high for system to be considered
C as separate wave system.
C DHSMAX Real Max. change in Hs for system to be considered
C related to previous time.
C DTPAX Real Id. Tp.
C DDMMAX Real Id. Dm.
C DDWMAX Real Maximum differences in wind and wave direction
C for marking of system as under the influence
C of the local wind,
C AGEMIN Real Id. wave age.
C ----------------------------------------------------------------
C
C 4. Subroutines used :
C
C MPARS Calculate mean parameters
C FNDPRT Get mask for partition.
C WAVNU2 Solve dispersion relation.
C
C 5. Called by :
C
C None, main program,
C
C 6. Error messages :
C
C See "Error escape locations" at end of code.
C
C 7. Remarks :
C
C 8. Structure :
C
C See source code.
C
C 10. Source code :
C
C/ ------------------------------------------------------------------- /
C/ Parameter statements
C/
INTEGER NFR, NTH, NPMAX, NPTAB
REAL XFR, HSMIN, HSDROP,
& DHSMAX, DTPAX, DDMMAX, DDWMAX, AGEMIN
*
PARAMETER ( NFR = 25 )
PARAMETER ( NTH = 24 )
PARAMETER ( NPMAX = 20 )
PARAMETER ( NPTAB = 6 )
PARAMETER ( XFR = 1.1 )
PARAMETER ( HSMIN = 0.15 )
PARAMETER ( HSDROP = 0.05 )
PARAMETER ( DHSMAX = 1.5 )
PARAMETER ( DTPMAX = 1.5 )
PARAMETER ( DDMMAX = 15. )
PARAMETER ( DDWMAX = 30. )
PARAMETER ( AGEMIN = 0.8 )
*
INTEGER NDSI, NDSO, IERR, NFRF, NTHF, IFR, ITH
INTEGER NPNTS, IP, TIME(2), IREAD
INTEGER IFL, IFH, ITHL, ITHH
INTEGER NPEAK, NFRP(NPMAX), NTHP(NPMAX), IPNOW, IOUT, ITAB
INTEGER NZERO
REAL PI, RADE, FR(NFR), TH(NTH), DTH, DFR(NFR),
& TAILF, SPEC(NFR,NTH), W1(NFR,NTH), W2(NFR,NTH),
& SP1D(NFR), HS, TP, DM, XFRV, STMAX, HMAX
REAL HSTOT, HSP(NPMAX), TPP(NPMAX), DMP(NPMAX)
REAL HST(NPMAX,2), TPT(NPMAX,2), DMT(NPMAX,2)
REAL D, UA, UD, UABS, UDIR, DELDW, DELHS, DELTP, DELDM
REAL AFR, WN, CG, AGE, Y, X
LOGICAL FORMI, DATA, FLAG(NPMAX), HEADER
CHARACTER POINT*10, PID*10, FNAMEI*40, FNAMEO*40, IDSTR*21
CHARACTER MODEL*40, RUN*20, GNAME*30, IDLAT*1, IDLON*1
CHARACTER*129 BLANK, TAIL, STRING
CHARACTER*15 PART
C/
C/ ------------------------------------------------------------------- /
C/
*
* 1. Initializations ------------------------------------------------ *
* 1.a Constants etc.
*
NDSI = 10
NDSO = 50
*
PI = 4. * ATAN(1.)
RADE = 180. / PI
XFRV = XFR
*
IREAD = 0
*
TAIL ( 1: 40) = '+-------+-----------+-----------------+-'
TAIL ( 41: 80) = '----------------+-----------------+-----'
TAIL ( 81:120) = '------------+-----------------+---------'
TAIL (120:129) = '---------+'
BLANK( 1: 40) = '| nn nn | nn | | '
BLANK( 41: 80) = ' | | '
BLANK( 81:120) = ' | | '
BLANK(120:129) = ' |'
STRING = BLANK
*
DO 100, IP=1, NPTAB
HST(IP,1) = -1.
TPT(IP,1) = -1.
DMT(IP,1) = -1.
100 CONTINUE
*
* 1.b Initial I/O
*
WRITE (*,900)
READ (*,'(A)') POINT
READ (*,'(A)') RUN
WRITE (*,901) POINT, RUN
READ (*,'(A)') FNAMEI
READ (*,*) FORMI
WRITE (*,902) FNAMEI, FORMI
READ (*,'(A)') FNAMEO
WRITE (*,903) FNAMEO
*
* 1.c Open input file and process header
*
IF ( FORMI ) THEN
OPEN (NDSI,FILE=FNAMEI,STATUS='OLD',ERR=800,IOSTAT=IERR)
ELSE
OPEN (NDSI,FILE=FNAMEI,STATUS='OLD',
& FORM='UNFORMATTED',ERR=800,IOSTAT=IERR)
ENDIF
WRITE (*,904)
*
IF ( FORMI ) THEN
READ (NDSI,*,END=801,ERR=802,IOSTAT=IERR)
& IDSTR, NFRF, NTHF, NPNTS, GNAME
ELSE
READ (NDSI , END=801,ERR=802,IOSTAT=IERR)
& IDSTR, NFRF, NTHF, NPNTS, GNAME
ENDIF
*
IF ( IDSTR .NE. 'WAVEWATCH III SPECTRA' ) GOTO 803
IF ( NFR.NE.NFRF .OR. NTH.NE.NTHF ) GOTO 804
WRITE (*,905)
*
IF ( FORMI ) THEN
READ (NDSI,*,END=801,ERR=802,IOSTAT=IERR) FR
ELSE
READ (NDSI , END=801,ERR=802,IOSTAT=IERR) FR
ENDIF
WRITE (*,906)
*
IF ( FORMI ) THEN
READ (NDSI,*,END=801,ERR=802,IOSTAT=IERR) TH
ELSE
READ (NDSI , END=801,ERR=802,IOSTAT=IERR) TH
ENDIF
WRITE (*,907)
*
DTH = 2. * PI / REAL(NTH)
TAILF = 0.25 * FR(NFR)
*
DO 110, IFR=1, NFR
DFR(IFR) = 0.5 * (XFR-1./XFR) * FR(IFR)
110 CONTINUE
DFR(NFR) = 0.5 * DFR(NFR)
WRITE (*,908)
*
* 1.d Open output file and process header
*
OPEN (NDSO,FILE=FNAMEO,ERR=810,IOSTAT=IERR)
WRITE (*,909)
HEADER = .TRUE.
*
* === LOOP OVER DATA ================================================= *
*
200 CONTINUE
*
* 2. Read next data ------------------------------------------------- *
*
IF ( FORMI ) THEN
READ (NDSI,*,END=888,ERR=802,IOSTAT=IERR) TIME
ELSE
READ (NDSI , END=888,ERR=802,IOSTAT=IERR) TIME
ENDIF
*
DATA = .FALSE.
DO 210, IP=1, NPNTS
IF ( FORMI ) THEN
READ (NDSI,*,END=801,ERR=802,IOSTAT=IERR)
& PID, Y, X, D, UA, UD
ELSE
READ (NDSI , END=801,ERR=802,IOSTAT=IERR)
& PID, Y, X, D, UA, UD
ENDIF
IF ( PID .EQ. POINT ) THEN
IF ( FORMI ) THEN
READ (NDSI,*,END=801,ERR=802,IOSTAT=IERR) SPEC
ELSE
READ (NDSI , END=801,ERR=802,IOSTAT=IERR) SPEC
ENDIF
DATA = .TRUE.
UABS = UA
UDIR = MOD( UD+180., 360. )
ELSE
IF ( FORMI ) THEN
READ (NDSI,*,END=801,ERR=802,IOSTAT=IERR) W1
ELSE
READ (NDSI , END=801,ERR=802,IOSTAT=IERR) W1
ENDIF
ENDIF
210 CONTINUE
*
IF ( .NOT. DATA ) GOTO 805
IREAD = IREAD + 1
*
IF ( HEADER ) THEN
X = MOD ( X+720. , 360. )
IF ( X .LE. 180. ) THEN
IDLON = 'E'
ELSE
X = 360. - X
IDLON = 'W'
ENDIF
IF ( ABS(Y) .LE. 0.0049 ) THEN
IDLAT = ' '
ELSE IF ( Y .GT. 0. ) THEN
IDLAT = 'N'
ELSE
IDLAT = 'S'
X = -X
ENDIF
WRITE (NDSO,950) POINT, Y, IDLAT, X, IDLON, GNAME, RUN
WRITE (NDSO,951) TAIL
WRITE (NDSO,952)
WRITE (NDSO,951) TAIL
HEADER = .FALSE.
ENDIF
*
* 3. Get overall wave height ---------------------------------------- *
*
DO 310, ITH=1, NTH
DO 300, IFR=1, NFR
W1(IFR,ITH) = 1.
300 CONTINUE
310 CONTINUE
*
CALL MPARS ( NFR, NTH, SPEC, W1, FR, TH, DFR, DTH, TAILF, XFRV,
& SP1D, HS, TP, DM )
*
HSTOT = HS
*
* 4. Determine maxima ----------------------------------------------- *
*
NPEAK = 0
*
DO 410, IFR=NFR, 1, -1
IFL = MAX ( 1 , IFR-1 )
IFH = MIN ( NFR , IFR+1 )
DO 400, ITH=1, NTH
ITHL = 1 + MOD(NTH+ITH-2,NTH)
ITHH = 1 + MOD(ITH,NTH)
IF ( SPEC(IFR,ITH) .GE. SPEC(IFL,ITH ) .AND.
& SPEC(IFR,ITH) .GE. SPEC(IFH,ITH ) .AND.
& SPEC(IFR,ITH) .GE. SPEC(IFL,ITHL) .AND.
& SPEC(IFR,ITH) .GE. SPEC(IFR,ITHL) .AND.
& SPEC(IFR,ITH) .GE. SPEC(IFH,ITHL) .AND.
& SPEC(IFR,ITH) .GE. SPEC(IFL,ITHH) .AND.
& SPEC(IFR,ITH) .GE. SPEC(IFR,ITHH) .AND.
& SPEC(IFR,ITH) .GE. SPEC(IFH,ITHH) .AND.
& SPEC(IFR,ITH) .GT. 0. ) THEN
NPEAK = NPEAK + 1
NFRP(NPEAK) = IFR
NTHP(NPEAK) = ITH
ENDIF
W1(IFR,ITH) = 0.
400 CONTINUE
410 CONTINUE
*
IF ( NPEAK .LE. NPMAX ) THEN
WRITE (*,910) TIME, NPEAK
ELSE
NPEAK = NPMAX
WRITE (*,911) TIME, NPEAK
ENDIF
*
* 5. Process all partial fields ------------------------------------- *
*
NZERO = 0
*
DO 500, IP=1, NPEAK
*
CALL FNDPRT ( NFR, NTH, NFRP(IP), NTHP(IP), SPEC, W1, W2 )
*
CALL MPARS ( NFR, NTH, SPEC, W2, FR, TH, DFR, DTH, TAILF, XFRV,
& SP1D, HSP(IP), TPP(IP), DMP(IP) )
IF ( HSP(IP) .LE. HSDROP ) NZERO = NZERO + 1
*
500 CONTINUE
*
DO 510, IP=NPEAK+1, NPMAX
HSP(IP) = 0.00
TPP(IP) = -999.99
DMP(IP) = -999.99
510 CONTINUE
*
DO 520, IP=1, NPTAB
HST(IP,2) = HST(IP,1)
TPT(IP,2) = TPT(IP,1)
DMT(IP,2) = DMT(IP,1)
HST(IP,1) = -1.
TPT(IP,1) = -1.
DMT(IP,1) = -1.
520 CONTINUE
*
* 6. Generate output table ------------------------------------------ *
* 6.a Time and overall wave height to string
*
STRING = BLANK
*
WRITE (STRING(3:4),'(I2)') MOD(TIME(1),100)
WRITE (STRING(6:7),'(I2)') TIME(2)/10000
IF ( HSTOT .GT. 0. ) WRITE (STRING(11:14),'(F4.1)') HSTOT
WRITE (STRING(16:17),'(I2)') NPEAK - NZERO
IF ( NPEAK.EQ.0 .OR. HSTOT.LT.0.1 ) GOTO 699
*
* 6.b Switch off peak with too low wave height
*
DO 600, IP=1, NPEAK
FLAG(IP) = HSP(IP) .GT. HSMIN
600 CONTINUE
*
* 6.c Find next highest wave height
*
IOUT = 0
601 CONTINUE
*
HMAX = 0.
IPNOW = 0
DO 610, IP=1, NPEAK
IF ( HSP(IP).GT.HMAX .AND. FLAG(IP) ) THEN
IPNOW = IP
HMAX = HSP(IP)
ENDIF
610 CONTINUE
*
* 6.d No more peaks, skip to output
*
IF ( IPNOW .EQ. 0 ) GOTO 699
*
* 6.e Find matching field
*
ITAB = 0
*
DO 620, IP=1, NPTAB
IF ( TPT(IP,2) .GT. 0. ) THEN
DELHS = ABS ( HST(IP,2) - HSP(IPNOW) )
DELTP = ABS ( TPT(IP,2) - TPP(IPNOW) )
DELDM = ABS ( DMT(IP,2) - DMP(IPNOW) )
IF ( DELDM .GT. 180. ) DELDM = 360. - DELDM
IF ( DELHS.LT.DHSMAX .AND.
& DELTP.LT.DTPMAX .AND.
& DELDM.LT.DDMMAX ) ITAB = IP
ENDIF
620 CONTINUE
*
* 6.f No matching field, find empty fields
*
IF ( ITAB .EQ. 0 ) THEN
DO 630, IP=NPTAB, 1, -1
IF ( TPT(IP,1).LT.0. .AND. TPT(IP,2).LT.0. ) ITAB = IP
630 CONTINUE
ENDIF
*
* 6.g Slot in table found, write
*
IF ( ITAB .NE. 0 ) THEN
*
WRITE (PART,'(2X,F4.1,F5.1,I4)')
& HSP(IPNOW), TPP(IPNOW), NINT(DMP(IPNOW))
DELDW = MOD ( ABS ( UDIR - DMP(IPNOW) ) , 360. )
IF ( DELDW .GT. 180. ) DELDW = 360. - DELDW
AFR = 2.*PI/TPP(IPNOW)
CALL WAVNU2 ( AFR, D, WN, CG, 1.E-5, 15, ICON )
AGE = UABS * WN / AFR
IF ( DELDW.LT.DDWMAX .AND. AGE.GT.AGEMIN ) PART(1:1) = '*'
*
STRING(5+ITAB*18:19+ITAB*18) = PART
*
HST(ITAB,1) = HSP(IPNOW)
TPT(ITAB,1) = TPP(IPNOW)
DMT(ITAB,1) = DMP(IPNOW)
*
* 6.g No slot in table found, write
*
ELSE
*
IOUT = IOUT + 1
WRITE (STRING(19:19),'(I1)') IOUT
*
ENDIF
*
FLAG(IPNOW) = .FALSE.
GOTO 601
*
* 6.h End of processing, write line in table
*
699 CONTINUE
*
WRITE (NDSO,951) STRING
*
* === BRANCH BACK IN DATA LOOP ======================================= *
*
GOTO 200
*
* Error escape locations ----------------------------------------- *
*
c GOTO 888
*
800 CONTINUE
WRITE (*,1000) FNAMEI, IERR
STOP
*
801 CONTINUE
WRITE (*,1001) FNAMEI
STOP
*
802 CONTINUE
WRITE (*,1002) FNAMEI, IERR
STOP
*
803 CONTINUE
WRITE (*,1003) IDSTR, 'WAVEWATCH III SPECTRA'
STOP
*
804 CONTINUE
WRITE (*,1004) NFR, NTH, NFRF, NTHF
STOP
*
805 CONTINUE
WRITE (*,1005) POINT
STOP
*
810 CONTINUE
WRITE (*,1010) FNAMEO, IERR
STOP
*
888 CONTINUE
WRITE (NDSO,951) TAIL
WRITE (NDSO,953) HSDROP, HSMIN
WRITE (*,999)
*
* Format statements ---------------------------------------------- *
*
900 FORMAT (/' Splitting WAVEWATCH III spectra : '/
& ' ----------------------------------')
901 FORMAT ( ' Location : ',A/
& ' Run ID : ',A)
902 FORMAT ( ' Input file : ',A,' FORMATTED : ',L1)
903 FORMAT ( ' Output file : ',A)
*
904 FORMAT (/' Input file opened ...')
905 FORMAT ( ' Header read and processed.')
906 FORMAT ( ' Frequency info read.')
907 FORMAT ( ' Direction info read.')
908 FORMAT ( ' Bin sizes preprocessed.')
*
909 FORMAT (/' Output file opened ...')
910 FORMAT ( ' Processing ',I8.8,' ',I6.6,
& ' Number of peaks :',I3)
911 FORMAT ( ' Processing ',I8.8,' ',I6.6,
& ' Number of peaks :',I3,' (TRUNCATED)')
*
950 FORMAT ( ' Location : ',A,' (',F5.2,A,1X,F6.2,A,')'/
& ' Model : ',A/
& ' Cycle : ',A)
951 FORMAT (1X,A)
952 FORMAT (' | day & | Hst n x | Hs Tp dir |',
& ' Hs Tp dir |',
& ' Hs Tp dir |',
& ' Hs Tp dir |',
& ' Hs Tp dir |',
& ' Hs Tp dir |'/
& ' | hour | (m) - - | (m) (s) (d) |',
& ' (m) (s) (d) |',
& ' (m) (s) (d) |',
& ' (m) (s) (d) |',
& ' (m) (s) (d) |',
& ' (m) (s) (d) |')
953 FORMAT (
& 75X,'Hst : Total sigificant wave height.'/
& 75X,'n : Number of fields with Hs > ',f6.2,
& ' in 2-D spectrum.'/
& 75X,'x : Number of fields with Hs > ',f6.2,
& ' not in table.'/
& 75X,'Hs : Significant wave height of separate wave field.'/
& 75X,'Tp : Peak period of separate wave field.'/
& 75X,'dir : Mean direction of separate wave field.'/
& 75X,'* : Wave generation due to local wind probable.')
*
999 FORMAT (/' End of program '/)
*
1000 FORMAT (/' *** ERROR IN OPENING INPUT FILE ',A/
& ' IOSTAT = ',I8/)
1001 FORMAT (/' *** PREMATURE END OF INPUT FILE ',A/)
1002 FORMAT (/' *** ERROR IN READING FROM INPUT FILE ',A/
& ' IOSTAT = ',I8/)
1003 FORMAT (/' *** UNEXPECTED HEADER IN FILE : '/
& ' READ : ',A/
& ' EXPECTED : ',A)
1004 FORMAT (/' *** UNEXPECTED SPECTRAL DIMENSIONS : '/
& ' IN PROGRAM :',2I6/
& ' IN FILE :',2I6/)
1005 FORMAT (/' *** NO DATA FOUND FOR POINT ',A/)
*
1010 FORMAT (/' *** ERROR IN OPENING OUTPUT FILE ',A/
& ' IOSTAT = ',I8/)
C/
C/ End of W3SPLT ----------------------------------------------------- /
C/
END
C/ ------------------------------------------------------------------- /
SUBROUTINE MPARS ( NFR, NTH, SPEC, MASK, FREQ, DIR, DFR, DTH,
& TAILF, XFR, SP1D, HS, TP, DM )
C/
C/ +-----------------------------------+
C/ | WAVEWATCH III NOAA/NCEP |
C/ | H. L. Tolman |
C/ | FORTRAN 77 |
C/ | Last update : 15-Oct-1998 |
C/ +-----------------------------------+
C/
C 1. Purpose :
C
C Calculate mean parameters of a spectrum using a mask.
C
C 2. Method :
C
C Wave height HS : Integral over spectrum using par. tail.
C Peak period TP : Parabolic fit to 1-D spectrum.
C Mean direction DM : From first Fourier components.
C
C 3. Parameters :
C
C Parameter list
C ----------------------------------------------------------------
C NFR Int. I Number of frequencies.
C NTH Int. I Number of directions.
C SPEC R.A. I Spectrum.
C MASK R.A. I Mask for spectrumi (real weight function).
C FREQ R.A. I Frequencies.
C DIR R.A. I Directions (radians).
C DFR R.A. I Bind width info for all frequencies.
C DTH R.A. I Directionsl increment.
C TAILF Real I Tail factor.
C XFR Real I Frequency increment factor.
C SP1D R.A. O 1-D spectrum.
C HS Real O Wave height -999.99 if undef.
C TP Real O Peak wave period -999.99 if undef.
C DM Real O Mean wave direction -999.99 if undef.
C ----------------------------------------------------------------
C
C 4. Subroutines used :
C
C None.
C
C 5. Called by :
C
C Any.
C
C 10. Source code :
C
C/ ------------------------------------------------------------------- /
C/ Parameter list
C/
INTEGER NFR, NTH
REAL SPEC(NFR,NTH), MASK(NFR,NTH), FREQ(NFR), DIR(NTH),
& DFR(NFR), DTH, TAILF, XFR, SP1D(NFR), HS, TP, DM
C/
C/ ------------------------------------------------------------------- /
C/ Local parameters
C/
INTEGER IFR, ITH, IMAX, ILOW, IHGH
REAL COSMOM, SINMOM, TOTMOM, COS0, SIN0, PI, RADE, EMAX
REAL XL, XH, XL2, XH2, EL, EH, DENOM, FP
C/
C/ ------------------------------------------------------------------- /
*
PI = 4. * ATAN(1.)
RADE = 180./PI
*
XL = 1./XFR - 1.
XH = XFR - 1.
XL2 = XL**2
XH2 = XH**2
*
COSMOM = 0.
SINMOM = 0.
TOTMOM = 0.
*
DO 110, IFR=1, NFR
*
SP1D(IFR) = 0.
COS0 = 0.
SIN0 = 0.
*
DO 100, ITH=1, NTH
SP1D(IFR) = SP1D(IFR) + MASK(IFR,ITH)*SPEC(IFR,ITH)
COS0 = COS0 + MASK(IFR,ITH)*SPEC(IFR,ITH) * COS(DIR(ITH))
SIN0 = SIN0 + MASK(IFR,ITH)*SPEC(IFR,ITH) * SIN(DIR(ITH))
100 CONTINUE
*
SP1D(IFR) = SP1D(IFR) * DTH
COSMOM = COSMOM + COS0 * DTH * DFR(IFR)
SINMOM = SINMOM + SIN0 * DTH * DFR(IFR)
TOTMOM = TOTMOM + SP1D(IFR)*DFR(IFR)
*
110 CONTINUE
*
TOTMOM = TOTMOM + SP1D(NFR)*TAILF
TOTMOM = MAX ( 0. , TOTMOM )
*
* 2. Mean parameters ------------------------------------------------ *
* 2.a Set to undefined
*
HS = 0.00
TP = -999.99
DM = -999.99
*
IF ( TOTMOM .GT. 0. ) THEN
*
* 2.b Wave height
*
HS = 4. * SQRT(TOTMOM)
*
* 2.c Mean direction
*
DM = RADE * ATAN2(SINMOM,COSMOM)
IF ( DM .LT. 0. ) DM = DM + 360.
*
* 2.d Peak period
*
EMAX = 0.
*
DO 200, IFR=1, NFR
IF ( SP1D(IFR) .GT. EMAX ) THEN
EMAX = SP1D(IFR)
IMAX = IFR
ENDIF
200 CONTINUE
*
ILOW = MAX ( 1 , IMAX-1 )
IHGH = MIN ( NFR , IMAX+1 )
EL = SP1D(ILOW) - SP1D(IMAX)
EH = SP1D(IHGH) - SP1D(IMAX)
DENOM = XL*EH - XH*EL
FP = FREQ(IMAX) * ( 1. + 0.5 * ( XL2*EH - XH2*EL )
& / SIGN ( MAX(ABS(DENOM),1.E-15) , DENOM ) )
*
TP = 1. / FP
*
ENDIF
*
RETURN
C/
C/ End of MPARS ------------------------------------------------------ /
C/
END
C/ ------------------------------------------------------------------- /
SUBROUTINE FNDPRT ( NFR, NTH, IFRC, ITHC, SPEC, W1, W2 )
C/
C/ +-----------------------------------+
C/ | WAVEWATCH III NOAA/NCEP |
C/ | H. L. Tolman |
C/ | FORTRAN 77 |
C/ | Last update : 09-Sep-1998 |
C/ +-----------------------------------+
C/
C 1. Purpose :
C
C Find partition starting at given peak and given mask.
C
C 2. Method :
C
C 3. Parameters :
C
C Parameter list
C ----------------------------------------------------------------
C NFR Int. I Number of frequencies.
C NTH Int. I Number of directions.
C IFRC Int. I Peak discrete frequency.
C ITHC Int. I Peak discrete direction.
C SPEC R.A. I Spectrum.
C W1 R.A. I/O Map of bins used so far.
C W2 R.A. O Map of bins used now.
C ----------------------------------------------------------------
C
C 4. Subroutines used :
C
C None.
C
C 5. Called by :
C
C Any.
C
C 10. Source code :
C
C/ ------------------------------------------------------------------- /
C/ Parameter list
C/
INTEGER NFR, NTH, IFRC, ITHC
REAL SPEC(NFR,NTH), W1(NFR,NTH), W2(NFR,NTH)
C/
C/ ------------------------------------------------------------------- /
C/ Local parameters
C/
INTEGER IFR, ITH, IFRL, ITHL
LOGICAL CHANGE, ADD
C/
C/ ------------------------------------------------------------------- /
*
* 1. Set up the W2 map ---------------------------------------------- *
*
DO 110, IFR=1, NFR
DO 100, ITH=1, NTH
W2(IFR,ITH) = 0.
100 CONTINUE
110 CONTINUE
*
IFRL = MAX ( 1 , IFRC-1 )
IFRH = MIN ( NFR , IFRC+1 )
ITHL = 1 + MOD(NTH+ITHC-2,NTH)
ITHH = 1 + MOD(ITHC,NTH)
*
DO 130, ITH=ITHC-1, ITHC+1
ITHL = 1 + MOD(NTH+ITH-1,NTH)
DO 120, IFR=IFRC-1, IFRC+1
IFRL = MAX(1,MIN(NFR,IFR))
IF ( W1(IFRL,ITHL) .LE. 0.5 ) W2(IFRL,ITHL) = 0.5
120 CONTINUE
130 CONTINUE
*
IF ( W1(IFRC,ITHC) .LT. 0.25 ) W2(IFRC,ITHC) = 1.0
*
* 2. Itterate search ------------------------------------------------ *
*
NITT = 0
*
* 2.a Branch point
*
200 CONTINUE
NITT = NITT + 1
CHANGE = .FALSE.
*
* 2.b Determine central points
*
DO 240, IFR=1, NFR
DO 230, ITH=1, NTH
*
IF ( W2(IFR,ITH).EQ.0.5 .AND. W1(IFR,ITH).LT.0.5 ) THEN
ADD = .TRUE.
DO 220, ITHR=ITH-1, ITH+1
ITHL = 1 + MOD(NTH+ITHR-1,NTH)
DO 210, IFRL=MAX(1,IFR-1), MIN(NFR,IFR+1)
IF ( W2 (IFRL,ITHL).EQ.0. .AND.
& SPEC(IFRL,ITHL).GT.SPEC(IFR,ITH) ) ADD = .FALSE.
210 CONTINUE
220 CONTINUE
IF ( ADD ) W2(IFR,ITH) = 1.
CHANGE = CHANGE .OR. ADD
ENDIF
*
230 CONTINUE
240 CONTINUE
*
* 2.c Determine central points
*
DO 280, IFR=1, NFR
DO 270, ITH=1, NTH
*
IF ( W2(IFR,ITH).EQ.0. ) THEN
ADD = .FALSE.
DO 260, ITHR=ITH-1, ITH+1
ITHL = 1 + MOD(NTH+ITHR-1,NTH)
DO 250, IFRL=MAX(1,IFR-1), MIN(NFR,IFR+1)
IF ( W2 (IFRL,ITHL).EQ.1. ) ADD = .TRUE.
250 CONTINUE
260 CONTINUE
IF ( ADD ) W2(IFR,ITH) = 0.5
CHANGE = CHANGE .OR. ADD
ENDIF
*
270 CONTINUE
280 CONTINUE
*
* 2.d Branch back ?
*
IF ( CHANGE .AND. NITT.LT.25 ) GOTO 200
*
* 3 Update the overall map ----------------------------------------- *
*
DO 310, IFR=1, NFR
DO 300, ITH=1, NTH
W1(IFR,ITH) = W1(IFR,ITH) + W2(IFR,ITH)
300 CONTINUE
310 CONTINUE
*
RETURN
C/
C/ End of FNDPRT ----------------------------------------------------- /
C/
END
C/ ------------------------------------------------------------------- /
SUBROUTINE WAVNU2 (W,H,K,CG,EPS,NMAX,ICON)
C/
C/ +-----------------------------------+
C/ | WAVEWATCH III NOAA/NCEP |
C/ | H. L. Tolman |
C/ | FORTRAN 77 |
C/ | Last update : 17-Jul-1990 |
C/ +-----------------------------------+
C/
C 1. Purpose :
C
C Calculation of wavenumber K from a given angular
C frequency W and waterdepth H.
C
C 2. Method :
C
C Used equation :
C 2
C W = G*K*TANH(K*H)
C
C Because of the nature of the equation, K is calculated
C with an itterative procedure.
C
C 3. Parameters :
C
C Parameter list
C ----------------------------------------------------------------
C W Real I Angular frequncy
C H Real I Waterdepth
C K Real O Wavenumber ( same sign as W )
C CG Real O Group velocity (same sign as W)
C EPS Real I Wanted max. difference between K and Kold
C NMAX Int. I Max number of repetitions in calculation
C ICON Int. O Contol counter ( See error messages )
C ----------------------------------------------------------------
C
C 9. Switches :
C
C C/S Enable subroutine tracing.
C
C 10. Source code :
C/
C/ ------------------------------------------------------------------- /
C/ Parameter list
C/
INTEGER ICON, NMAX
REAL EPS, CG, K, H, W
C/
C/ ------------------------------------------------------------------- /
C/ Local parameters
C/
INTEGER I
C/S INTEGER IENT
REAL G, F, W0, FD, DIF, RDIF, KOLD
C/
C/ ------------------------------------------------------------------- /
C/
DATA G / 9.81 /
*
* Initialisations :
*
CG = 0
KOLD = 0
ICON = 0
W0 = ABS(W)
*
* 1st approach :
*
IF (W0.LT.SQRT(G/H)) THEN
K = W0/SQRT(G*H)
ELSE
K = W0*W0/G
END IF
*
* Refinement :
*
DO 7, I=1, NMAX
DIF = ABS(K-KOLD)
IF (K.NE.0) THEN
RDIF = DIF/K
ELSE
RDIF = 0
END IF
IF (DIF .LT. EPS .AND. RDIF .LT. EPS) THEN
ICON = 1
GOTO 8
ELSE
KOLD = K
F = G*KOLD*TANH(KOLD*H)-W0**2
IF (KOLD*H.GT.25) THEN
FD = G*TANH(KOLD*H)
ELSE
FD = G*TANH(KOLD*H) + G*KOLD*H/((COSH(KOLD*H))**2)
END IF
K = KOLD - F/FD
END IF
7 CONTINUE
DIF = ABS(K-KOLD)
RDIF = DIF/K
IF (DIF .LT. EPS .AND. RDIF .LT. EPS) ICON = 1
8 CONTINUE
IF (2*K*H.GT.25) THEN
CG = W0/K * 0.5
ELSE
CG = W0/K * 0.5*(1+(2*K*H/SINH(2*K*H)))
END IF
IF (W.LT.0.0) THEN
K = (-1)*K
CG = CG*(-1)
END IF
*
RETURN
C/
C/ End of WAVNU2 ----------------------------------------------------- /
C/
END