#include "swancpp.h"
!
! SWAN - SERVICE ROUTINES
!
! Contents of this file
!
! READXY
! REFIXY
! INFRAM
! DISTR
! KSCIP1
! KSCIP2 40.59
! NG 40.59
! AC2TST
! CVCHEK 30.60
! CVMESH 30.60
! NEWTON 30.60
! EVALF 30.60
! SWOBST 30.60
! SWOBSTO 40.86
! TCROSS 40.04
! SWTRCF
! SSHAPE 40.00
! SINTRP 40.00
! HSOBND 32.01
! CHGBAS 40.00
! GAMMAF 40.00
! WRSPEC 40.00
!TIMG! SWTSTA 40.23
!TIMG! SWTSTO 40.23
!TIMG! SWPRTI 40.23
! TXPBLA 40.23
! INTSTR 40.23
! NUMSTR 40.23
! SWCOPI 40.23
! SWCOPR 40.23
!MatL4! SWI2B 40.30
!MatL4! SWR2B 40.30
!
! functions:
! ----------
! DEGCNV (converts from cartesian convention to nautical and 32.01
! vice versa) 32.01
! ANGRAD (converts radians to degrees) 32.01
! ANGDEG (converts degrees to radians) 32.01
!
! subroutines:
! ------------
! HSOBND (Hs is calculated after a SWAN computation at all sides. 32.01
! The calculated wave height from SWAN is then compared with 32.01
! the wave heigth as provided by the user 32.01
!
!***********************************************************************
! *
SUBROUTINE READXY (NAMX, NAMY, XX, YY, KONT, XSTA, YSTA)
! *
!***********************************************************************
!
USE OCPCOMM1 40.41
USE SWCOMM2 40.41
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
! 40.22: John Cazes and Tim Campbell
! 40.13: Nico Booij
! 40.51: Marcel Zijlema
!
! 1. UPDATE
!
! Nov. 1996 offset values are added to standard values
! because they will be subtracted later
! 40.13, Nov. 01: a valid value for YY is required if a valid value
! for XX has been given; ocpcomm1.inc reactivated
! 40.51, Feb. 05: correction to location points equal to offset values
!
! 2. PURPOSE
!
! Read x and y, initialize offset values XOFFS and YOFFS
!
! 3. METHOD
!
! ---
!
! 4. PARAMETERLIST
!
! NAMX, NAMY inp char names of the two coordinates as given in
! the user manual
! XX, YY out real values of x and y taking into account offset
! KONT inp char what to be done if values are missing
! see doc. of INDBLE (Ocean Pack doc.)
! XSTA, YSTA inp real standard values of x and y
!
! 5. SUBROUTINES CALLING
!
!
!
! 6. SUBROUTINES USED
!
! INDBLE (Ocean Pack)
LOGICAL EQREAL
!
! 7. ERROR MESSAGES
!
! ---
!
! 8. REMARKS
!
! ---
!
! 9. STRUCTURE
!
! ----------------------------------------------------------------
! Read x and y in double prec.
! If this is first couple of values
! Then assign values to XOFFS and YOFFS
! make LXOFFS True
! ---------------------------------------------------------------
! make XX and YY equal to x and y taking into account offset
! ----------------------------------------------------------------
!
! 10. SOURCE TEXT
!
DOUBLE PRECISION XTMP, YTMP
CHARACTER NAMX *(*), NAMY *(*), KONT *(*)
SAVE IENT
DATA IENT /0/
CALL STRACE (IENT,'READXY')
!
CALL INDBLE (NAMX, XTMP, KONT, DBLE(XSTA)+DBLE(XOFFS))
IF (CHGVAL) THEN 40.13
! a valid value was given for XX 40.13
CALL INDBLE (NAMY, YTMP, 'REQ', DBLE(YSTA)+DBLE(YOFFS)) 40.13
ELSE 40.13
CALL INDBLE (NAMY, YTMP, KONT, DBLE(YSTA)+DBLE(YOFFS))
ENDIF 40.13
IF (.NOT.LXOFFS) THEN
XOFFS = REAL(XTMP)
YOFFS = REAL(YTMP)
LXOFFS = .TRUE.
ENDIF
IF (.NOT.EQREAL(XOFFS,REAL(XTMP))) THEN 40.51
XX = REAL(XTMP-DBLE(XOFFS))
ELSE IF (OPTG.EQ.3) THEN 40.51
XX = 1.E-5 40.51
ELSE 40.51
XX = 0. 40.51
END IF 40.51
IF (.NOT.EQREAL(YOFFS,REAL(YTMP))) THEN 40.51
YY = REAL(YTMP-DBLE(YOFFS))
ELSE IF (OPTG.EQ.3) THEN 40.51
YY = 1.E-5 40.51
ELSE 40.51
YY = 0. 40.51
END IF 40.51
!
RETURN
! * end of subroutine READXY *
END
!***********************************************************************
! *
SUBROUTINE REFIXY (NDS, XX, YY, IERR)
! *
!***********************************************************************
!
USE SWCOMM2 40.41
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
! 40.22: John Cazes and Tim Campbell
! 40.51: M. Zijlema
!
! 1. UPDATE
!
! first version: 10.18 (Sept 1994)
!
! 2. PURPOSE
!
! initialize offset values XOFFS and YOFFS, and shift XX and YY
!
! 3. METHOD
!
! ---
!
! 4. PARAMETERLIST
!
! NDS in int file reference number
! XX, YY out real values of x and y taking into account offset
! IERR out int error indicator: IERR=0: no error, =-1: end-
! of-file, =-2: read error
!
! 5. SUBROUTINES CALLING
!
!
!
! 6. SUBROUTINES USED
!
LOGICAL EQREAL
!
! 7. ERROR MESSAGES
!
! ---
!
! 8. REMARKS
!
! ---
!
! 9. STRUCTURE
!
! ----------------------------------------------------------------
! If this is first couple of values
! Then assign values to XOFFS and YOFFS
! make LXOFFS True
! ---------------------------------------------------------------
! make XX and YY equal to x and y taking into account offset
! ----------------------------------------------------------------
!
! 10. SOURCE TEXT
!
DOUBLE PRECISION XTMP, YTMP
REAL XX, YY
SAVE IENT
DATA IENT /0/
CALL STRACE (IENT,'REFIXY')
!
READ (NDS, *, END=10, ERR=20) XTMP, YTMP
IF (.NOT.LXOFFS) THEN
XOFFS = REAL(XTMP)
YOFFS = REAL(YTMP)
LXOFFS = .TRUE.
ENDIF
IF (.NOT.EQREAL(XOFFS,REAL(XTMP))) THEN 40.51
XX = REAL(XTMP-DBLE(XOFFS))
ELSE IF (OPTG.EQ.3) THEN 40.51
XX = 1.E-5 40.51
ELSE 40.51
XX = 0. 40.51
END IF 40.51
IF (.NOT.EQREAL(YOFFS,REAL(YTMP))) THEN 40.51
YY = REAL(YTMP-DBLE(YOFFS))
ELSE IF (OPTG.EQ.3) THEN 40.51
YY = 1.E-5 40.51
ELSE 40.51
YY = 0. 40.51
END IF 40.51
!
IERR = 0
RETURN
! end of file
10 IERR = -1
RETURN
! read error
20 IERR = -2
RETURN
! * end of subroutine REFIXY *
END
!***********************************************************************
! *
LOGICAL FUNCTION INFRAM (XQQ, YQQ)
! *
!***********************************************************************
!
USE OCPCOMM4 40.41
USE SWCOMM1 40.41
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
! 40.41: Marcel Zijlema
!
! 1. UPDATE
!
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. PURPOSE
!
! Checking whether a point given in frame coordinates is located
! in the plotting frame (INFRAM = .TRUE.) or not (INFRAM = .FALSE.)
!
! 3. METHOD
!
! ---
!
! 4. PARAMETERLIST
!
! XQQ REAL input X-coordinate (output grid) of the point
! YQQ REAL input Y-coordinate (output grid) of the point
!
! 5. SUBROUTINES CALLING
!
! SPLSIT, PLNAME
!
! 6. SUBROUTINES USED
!
! none
!
! 7. ERROR MESSAGES
!
! ---
!
! 8. REMARKS
!
! ---
!
! 9. STRUCTURE
!
! ----------------------------------------------------------------
! Give INFRAM initial value true
! IF XQQ < 0, XQQ > XQLEN, YQQ < 0 OR YQQ > YQLEN, THEN
! INFRAM = false
! ----------------------------------------------------------------
!
! 10. SOURCE TEXT
!
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'INFRAM')
!
INFRAM = .TRUE.
IF (XQQ .LT. 0.) INFRAM = .FALSE.
IF (XQQ .GT. XQLEN) INFRAM = .FALSE.
IF (YQQ .LT. 0.) INFRAM = .FALSE.
IF (YQQ .GT. YQLEN) INFRAM = .FALSE.
!
RETURN
! * end of function INFRAM *
END
!***********************************************************************
! *
SUBROUTINE DISTR (CDIR, DIR, COEF, SPCDIR) 20.43
! *
!***********************************************************************
!
USE OCPCOMM4 40.41
USE SWCOMM3 40.41
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
!
! 0. Authors
!
! 30.82: IJsbrand Haagsma
! 40.22: John Cazes and Tim Campbell
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! 0.1 , Jul. 87: Standard heading added
! 0.2 , Dec. 89: Value for energy outside of the sector changed
! from 0. to 1.E-6
! Oct. 90: Value for energy outside the sector changed to 1.E-10
! logical BDIR introduced to take care for case where
! none of the values is positive
! 30.82, Oct. 98: Updated description of several variables
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. PURPOSE
!
! Computation of the distribution of the wave energy over the
! sectors, according to the given directional spread.
!
! 3. METHOD
!
! ---
!
! 4. Argument variables
!
! i SPCDIR: (*,1); spectral directions (radians) 30.82
! (*,2); cosine of spectral directions 30.82
! (*,3); sine of spectral directions 30.82
! (*,4); cosine^2 of spectral directions 30.82
! (*,5); cosine*sine of spectral directions 30.82
! (*,6); sine^2 of spectral directions 30.82
!
REAL SPCDIR(MDC,6) 30.82
!
! CDIR REAL output array containing the coefficients of
! (energy) distribution
! DIR REAL input main wave direction, in radians
! COEF REAL input coefficient of the directional distri-
! bution (cos**COEF)
!
! 5. SUBROUTINES CALLING
!
! REINVA (HISWA/SWREAD), STARTB and SWIND (both HISWA/COMPU)
!
! 6. SUBROUTINES USED
!
! none
!
! 7. Common blocks used
!
!
! 8. REMARKS
!
! ---
!
! 9. STRUCTURE
!
! -----------------------------------------------------------------
! For every direction of the grid do
! If the direction deviates less than PI/2 from the main wave
! direction, then
! Compute the coefficient cos**n
! Else
! Coefficient is 1.E-10
! -----------------------------------------------------------------
! If any of the directions deviated less than PI/2
! Then Compute the total of the coefficients
! For every direction of the grid do
! Divide the fraction of the distribution by the total
! -----------------------------------------------------------------
!
! 13. Source text
!
LOGICAL BDIR
REAL CDIR(*)
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT, 'DISTR')
!
SOMC = 0.
BDIR = .FALSE.
DO 10 ID0 = 1, MDC
TETA = SPCDIR(ID0,1)
ACOS = COS(TETA-DIR)
IF (ACOS .GT. 0.) THEN
BDIR = .TRUE.
CDIR(ID0) = MAX (ACOS**COEF, 1.E-10)
SOMC = SOMC + CDIR(ID0)
ELSE
CDIR(ID0) = 1.E-10
ENDIF
10 CONTINUE
IF (BDIR) THEN
CNORM = 1./(SOMC*DDIR)
DO 20 ID0 = 1, MDC
CDIR(ID0) = CDIR(ID0) * CNORM
20 CONTINUE
ENDIF
!
! ***** test *****
! IF(TESTFL .AND. ITEST .GE. 200)
IF( ITEST .GE. 200)
&WRITE(PRINTF,6010) CNORM,(CDIR(JJ) , JJ=1,MDC)
6010 FORMAT (' Test DISTR',F10.3/(10E12.4))
!
RETURN
! * end of subroutine DISTR *
END
!***********************************************************************
! *
SUBROUTINE KSCIP1 (MMT, SIG, D, K, CG, N, ND)
! *
!***********************************************************************
!
USE OCPCOMM4 40.41
USE SWCOMM3 40.41
!
IMPLICIT NONE
!
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 30.81: Annette Kieftenburg
! 40.22: John Cazes and Tim Campbell
! 40.41: Marcel Zijlema
! 41.16: Greg Wilson
!
! 1. Updates
!
! Aug. 94, ver. 10.10: arguments N and ND added
! Dec. 98, ND corrected, argument list adjusted and IMPLICIT NONE added
! 40.41, Aug. 04: tables replaced by Pade and other formulas
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
! 41.16, Mar. 11: correction: add dk/dh to dn/dh
!
! 2. Purpose
!
! Calculation of the wave number, group velocity, group number N
! and the derivative of N w.r.t. depth (=ND)
!
! 3. Method
!
! --
!
! 4. Argument variables
!
! MMT input number of frequency points
!
INTEGER MMT
!
! CG output group velocity
! D input local depth
! K output wave number
! N output ratio of group and phase velocity
! ND output derivative of N with respect to D
! SIG input rel. frequency for which wave parameters
! must be determined
!
REAL CG(MMT), D,
& K(MMT), N(MMT), ND(MMT), SIG(MMT)
!
! 6. Local variables
!
! C phase velocity
! FAC1 auxiliary factor
! FAC2 auxiliary factor
! FAC3 auxiliary factor
! IENT number of entries
! IS counter in frequency (sigma-space)
! KND dimensionless wave number
! ROOTDG square root of D/GRAV
! WGD square root of GRAV*D
! SND dimensionless frequency
! SND2 = SND*SND
!
INTEGER IENT, IS
REAL KND, ROOTDG, SND, WGD, SND2, C, FAC1, FAC2, FAC3
!
! 8. Subroutines used
!
! --
!
! 9. Subroutines calling
!
! SWOEXA, SWOEXF (Swan/Output)
!
! 10. Error messages
!
! --
!
! 11. Remarks
!
! --
!
! 12. Structure
!
! -----------------------------------------------------------------
! Compute non-dimensional frequency SND
! IF SND >= 2.5, then
! Compute wave number K, group velocity CGO, ratio of group
! and phase velocity N and its derivative ND according to
! deep water theory
! ELSE IF SND =< 1.e-6
! Compute wave number K, group velocity CGO, ratio of group
! and phase velocity N and its derivative ND
! according to extremely shallow water
! ELSE
! Compute wave number K, group velocity CGO and the ratio of
! group and phase velocity N by Pade and other simple formulas.
! Compute the derivative of N w.r.t. D = ND.
! -----------------------------------------------------------------
!
! 13. Source text
!
SAVE IENT
DATA IENT /0/
IF (LTRACE) CALL STRACE (IENT, 'KSCIP1')
!
ROOTDG = SQRT(D/GRAV) 30.81
WGD = ROOTDG*GRAV 30.81
DO 200 IS = 1, MMT
! SND is dimensionless frequency
SND = SIG(IS) * ROOTDG
IF (SND .GE. 2.5) THEN
! ******* deep water *******
K(IS) = SIG(IS) * SIG(IS) / GRAV 30.81
CG(IS) = 0.5 * GRAV / SIG(IS) 30.81
N(IS) = 0.5
ND(IS) = 0.
ELSE IF (SND.LT.1.E-6) THEN
! *** very shallow water *** 30.81
K(IS) = SND/D 30.81
CG(IS) = WGD
N(IS) = 1.
ND(IS) = 0.
ELSE
SND2 = SND*SND 40.41
C = SQRT(GRAV*D/(SND2+1./(1.+0.666*SND2+0.445*SND2**2 40.41
& -0.105*SND2**3+0.272*SND2**4))) 40.41
K(IS) = SIG(IS)/C 40.41
KND = K(IS)*D 40.41
FAC1 = 2.*KND/SINH(2.*KND) 40.41
N(IS) = 0.5*(1.+FAC1) 40.41
CG(IS)= N(IS)*C 40.41
FAC2 = SND2/KND 40.41
FAC3 = 2.*FAC2/(1.+FAC2*FAC2) 40.41
FAC2 = -K(IS)*(2.*N(IS)-1.)/(2.*D*N(IS)) 41.16
ND(IS)= FAC1*(0.5/D - K(IS)/FAC3 + FAC2*(0.5/K(IS) - D/FAC3)) 41.16 40.41
ENDIF
200 CONTINUE
!
RETURN
! end of subroutine KSCIP1 *
END
!***********************************************************************
! *
SUBROUTINE KSCIP2 (MMT, SIG, D, K, CG, N, ND, DMW, DM)
! *
!***********************************************************************
!
USE OCPCOMM4
USE SWCOMM3
!
IMPLICIT NONE
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 40.59: Erick Rogers
!
! 1. Updates
!
! 40.59, Aug. 07: copy from KSCIP1 and adaption to muddy bottom
!
! 2. Purpose
!
! Calculation of the wave number, group velocity, group number N
! and the derivative of N w.r.t. depth (=ND) in case of muddy bottom
!
! 3. Method
!
! Wave number based on Ng (2000) and other quantities derived thereof
!
! 4. Argument variables
!
! MMT input number of frequency points
!
INTEGER MMT
!
! CG output group velocity
! D input local depth
! DM input local depth of mud layer
! DMW output mud dissipation rate
! (is the imaginary part of the wave number)
! K in/out (the real part of) the wave number
! N output ratio of group and phase velocity
! ND output derivative of N with respect to D
! SIG input rel. frequency for which wave parameters
! must be determined
!
REAL CG(MMT), D, DM,
& K(MMT), N(MMT), ND(MMT), SIG(MMT), DMW(MMT)
!
! 6. Local variables
!
! C phase velocity
! DTILDE normalized mud depth = mud depth / delta_m,
! delta is the sblt= sqrt(2*visc/sigma)
! FAC1 auxiliary factor
! FAC2 auxiliary factor
! FAC3 auxiliary factor
! GAMMA this is the gamma used in Ng pg. 238. this is
! density(water)/density(mud)
! IENT number of entries
! IS counter in frequency space
! KINVISM kinematic viscosity of mud
! KINVISW kinematic viscosity of water
! KMUD (local) muddy wave number
! KND dimensionless wave number
! RHOM density of mud
! RHOW density of water
! SBLTM a function of viscosity and frequency
! ZETA this is zeta as used in Ng pg. 238. it is the
! ratio of Stokes boundary layer thicknesses,
! or delta_m/delta_w
!
INTEGER IENT, IS
REAL C, DTILDE, FAC1, FAC2, FAC3, GAMMA, KINVISM,
& KINVISW, KMUD, KND, RHOM, RHOW, SBLTM, ZETA
!
! 8. Subroutines used
!
! --
!
! 9. Subroutines calling
!
! --
!
! 10. Error messages
!
! --
!
! 11. Remarks
!
! --
!
! 12. Structure
!
! --
!
! 13. Source text
!
SAVE IENT
DATA IENT /0/
IF (LTRACE) CALL STRACE (IENT, 'KSCIP2')
!
RHOM = PMUD(2)
KINVISM = PMUD(3)
RHOW = PMUD(4)
KINVISW = PMUD(5)
!
ZETA = SQRT(KINVISM/KINVISW)
GAMMA = RHOW/RHOM
!
DO IS = 1, MMT
KND = K(IS)*D
IF ( KND.LT.10. .AND. DM.GT.1.E-5 ) THEN
SBLTM = SQRT(2.*KINVISM/SIG(IS))
DTILDE = DM/SBLTM
! calculate muddy wave number and dissipation rate using Ng (2000)
CALL NG(SIG(IS),D,DTILDE,ZETA,SBLTM,GAMMA,K(IS),KMUD,DMW(IS))
K(IS) = KMUD
! calculate ratio N, CG and derivative of N w.r.t. D
KND = K(IS)*D
IF ( KND.LT.35. ) THEN
FAC1 = 2.*KND/SINH(2.*KND)
ELSE
FAC1 = 2.E-30*KND
ENDIF
N(IS) = 0.5*(1.+FAC1)
C = SIG(IS)/K(IS)
CG(IS) = N(IS)*C
FAC2 = SIG(IS) * C / GRAV
FAC3 = 2.*FAC2/(1.+FAC2*FAC2)
FAC2 = -K(IS)*(2.*N(IS)-1.)/(2.*D*N(IS))
ND(IS) = FAC1*(0.5/D - K(IS)/FAC3 + FAC2*(0.5/K(IS) - D/FAC3))
ELSE
DMW(IS)= 0.
ENDIF
ENDDO
!
RETURN
! end of subroutine KSCIP2 *
END
!************************************************************************
! *
SUBROUTINE NG(SIGMA,H_WDEPTH,DTILDE,ZETA,SBLTM,GAMMA,WK,WKDR,DISS)
! *
!************************************************************************
!
USE OCPCOMM4
!
IMPLICIT NONE
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 1.00: Jim Kaihatu
! 40.??: Erick Rogers
! 40.59: Erick Rogers
!
! 1. Updates
!
! 1.00, Dec. 04: original
! 40.??, Mar. 05: adapted for SWAN
! 40.59, Aug. 07: included patch provided by Jim Kaihatu
! and finalized
!
! 2. Purpose
!
! Computes the muddy wavenumber using Ng (2000)
!
! 3. Method
!
! Reference: Ng, C.O., 2000. Water waves over a muddy bed:
! a two-layer Stokes boundary layer model
! Coastal Engineering 40(3), 221-242
!
! 4. Argument variables
!
REAL, INTENT(IN) :: SIGMA ! radian frequency (rad)
REAL, INTENT(IN) :: H_WDEPTH! water depth, denoted "h" in Ng (m)
REAL, INTENT(IN) :: DTILDE ! normalized mud depth = mud depth / SBLTM,
! delta is the sblt= sqrt(2*visc/sigma)
REAL, INTENT(IN) :: ZETA ! this is zeta as used in Ng pg. 238. it is
! the ratio of Stokes boundary layer
! thicknesses, or SBLTM/delta_w
REAL, INTENT(IN) :: GAMMA ! this is the gamma used in Ng pg. 238.
! this is density(water)/density(mud)
REAL, INTENT(IN) :: SBLTM ! SBLTM is what you get if you calculate
! sblt using the viscosity of the mud,
! SBLTM=sqrt(2*visc_m/sigma)
! .....also delta_m
REAL, INTENT(IN) :: WK ! unmuddy wavenumber
REAL, INTENT(OUT) :: WKDR ! muddy wavenumber
REAL, INTENT(OUT) :: DISS ! dissipation rate
!
! 5. Local variables
!
INTEGER :: IENT
REAL :: B1 ! an Ng coefficient
REAL :: B2 ! an Ng coefficient
REAL :: B3 ! an Ng coefficient
REAL :: BR ! an Ng coefficient
REAL :: BI ! an Ng coefficient
REAL :: BRP ! an Ng coefficient
REAL :: BIP ! an Ng coefficient
REAL :: DM ! mud depth, added June 2 2006
!
! 11. Remarks
!
! Calculations for the "B coefficients" came from a code by Jim Kaihatu
!
! 13. Source text
!
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'NG')
!
DM=DTILDE*SBLTM ! DTILDE=DM/SBLTM
! now calculate Ng's B coefficients : see Ng pg 238
B1=GAMMA*(-2.*GAMMA**2+2.*GAMMA-1.-ZETA**2)*SINH(DTILDE)*
& COSH(DTILDE)-GAMMA**2*ZETA*((COSH(DTILDE))**2+
& (SINH(DTILDE))**2)-(GAMMA-1.)**2*ZETA*((COSH(DTILDE))**2
& *(COS(DTILDE))**2+(SINH(DTILDE))**2*(SIN(DTILDE))**2)-2.
& *GAMMA*(1.-GAMMA)*(ZETA*COSH(DTILDE)+GAMMA*SINH(DTILDE))
& *COS(DTILDE)
B2=GAMMA*(-2.*GAMMA**2+2.*GAMMA-1.+ZETA**2)*SIN(DTILDE)*
& COS(DTILDE) -2.*GAMMA*(1.-GAMMA)*(ZETA*SINH(DTILDE)+GAMMA
& *COSH(DTILDE))*SIN(DTILDE)
B3=(ZETA*COSH(DTILDE)+GAMMA*SINH(DTILDE))**2*(COS(DTILDE))**2
& +(ZETA*SINH(DTILDE)+GAMMA*COSH(DTILDE))**2*(SIN(DTILDE))**2
BR=WK*SBLTM*(B1-B2)/(2.*B3)+GAMMA*WK*DM
BI=WK*SBLTM*(B1+B2)/(2.*B3)
BRP=B1/B3 ! "B_R PRIME"
BIP=B2/B3 ! "B_I PRIME"
! now calculate dissipation rate and wavenumber
DISS=-SBLTM*(BRP+BIP)*WK**2/(SINH(2.*WK*H_WDEPTH)+2.*WK*H_WDEPTH)
WKDR=WK-BR*WK/(SINH(WK*H_WDEPTH)*COSH(WK*H_WDEPTH)+WK*H_WDEPTH)
RETURN
END SUBROUTINE NG
!***********************************************************************
! *
SUBROUTINE AC2TST (XYTST, AC2,KGRPNT)
! *
!***********************************************************************
!
USE OCPCOMM4 40.41
USE SWCOMM2 41.13
USE SWCOMM3 40.41
USE SWCOMM4 40.41
USE M_PARALL 40.31
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
!
!
!
!
INTEGER XYTST(*) ,KGRPNT(MXC,MYC)
REAL AC2(MDC,MSC,MCGRD) 30.21
!.................................................................
IF ( ITEST .GE. 100 .AND. TESTFL) THEN
DO II = 1, NPTST
IF (OPTG.NE.5) THEN 40.80
IX = XYTST(2*II-1)
IY = XYTST(2*II)
INDEX = KGRPNT(IX,IY)
WRITE (PRINTF, 618) IX+MXF-2, IY+MYF-2, KGRPNT(IX,IY) 40.30
ELSE 40.80
INDEX = XYTST(II) 40.80
WRITE (PRINTF, 619) INDEX 40.80
ENDIF 40.80
DO ID = 1, MDC
WRITE (PRINTF, 620) (AC2(ID,IS,INDEX), IS=1,MIN(10,MSC)) 30.21
ENDDO
ENDDO
ENDIF
618 FORMAT(/,'Spectrum for test point(index):', 2I5,2X,'(',I5,')')
619 FORMAT(/,'Spectrum for test point: (',I5,')') 40.80
620 FORMAT (10(1X,E12.4))
RETURN
END
!****************************************************************
!
SUBROUTINE CVCHEK (KGRPNT, XCGRID, YCGRID) 30.72
!
!****************************************************************
!
USE OCPCOMM4 40.41
USE SWCOMM2 40.41
USE SWCOMM3 40.41
USE SWCOMM4 40.41
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
! 0. Authors
!
! 30.72: IJsbrand Haagsma
! 40.13: Nico Booij
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! May 96: New subroutine
! 30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN
! 40.13, Mar. 01: messages corrected and extended
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. Purpose
!
! Checks whether the given curvilinear grid is correct
! also set the value of CVLEFT.
!
! 3. Method
!
! Going around a mesh in the same direction the interior
! of the mesh must be always in the same side if the
! coordinates are correct
!
! 4. Argument variables
!
! KGRPNT: input Array of indirect addressing
!
INTEGER KGRPNT(MXC,MYC) 30.72
!
! XCGRID: input Coordinates of computational grid in x-direction 30.72
! YCGRID: input Coordinates of computational grid in y-direction 30.72
!
REAL XCGRID(MXC,MYC), YCGRID(MXC,MYC) 30.72
!
!
! 5. SUBROUTINES CALLING
!
! SWRBC
!
! 6. SUBROUTINES USED
!
! ---
!
! 7. ERROR MESSAGES
!
! ---
!
! 8. REMARKS
!
!
! 9. STRUCTURE
!
! ------------------------------------------------------------
! FIRST = True
! For ix=1 to MXC-1 do
! For iy=1 to MYC-1 do
! For iside=1 to 4 do
! Case iside=
! 1: K1 = KGRPNT(ix,iy), K2 = KGRPNT(ix+1,iy),
! K3 = KGRPNT(ix+1,iy+1)
! 2: K1 = KGRPNT(ix+1,iy), K2 = KGRPNT(ix+1,iy+1),
! K3 = KGRPNT(ix,iy+1)
! 3: K1 = KGRPNT(ix+1,iy+1), K2 = KGRPNT(ix,iy+1),
! K3 = KGRPNT(ix,iy)
! 4: K1 = KGRPNT(ix,iy+1), K2 = KGRPNT(ix,iy),
! K3 = KGRPNT(ix+1,iy)
! ---------------------------------------------------
! If K1>1 and K2>1 and K3>1
! Then Det = (xpg(K3)-xpg(K1))*(ypg(K2)-ypg(K1)) -
! (ypg(K3)-ypg(K1))*(xpg(K2)-xpg(K1))
! If FIRST
! Then Make FIRST = False
! If Det>0
! Then Make CVleft = False
! Else Make CVleft = True
! ----------------------------------------------
! If ((CVleft and Det<0) or (not CVleft and Det>0))
! Then Write error message with IX, IY, ISIDE
! ------------------------------------------------------------
!
! 10. SOURCE
!
!****************************************************************
!
!
LOGICAL FIRST
!
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'CVCHEK')
!
! test output
!
IF (ITEST .GE. 150 .OR. INTES .GE. 30) THEN
WRITE(PRINTF,186)
186 FORMAT(/,' ... Subroutine CVCHEK...',
& /,2X,'POINT( IX, IY), INDEX, COORDX, COORDY')
ICON = 0
DO 5 IIY = 1, MYC
DO 6 IIX = 1, MXC
ICON = ICON + 1
WRITE(PRINTF,7)IIX-1,IIY-1,KGRPNT(IIX,IIY),
& XCGRID(IIX,IIY)+XOFFS, YCGRID(IIX,IIY)+YOFFS 30.72 40.13
6 CONTINUE
5 CONTINUE
ENDIF
7 FORMAT(4X,I5,1X,I5,3X,I4,5X,F10.2,4X,F10.2)
!
FIRST = .TRUE.
!
DO 10 IX = 1,MXC-1
DO 15 IY = 1,MYC-1
DO 20 ISIDE = 1,4
IF (ISIDE .EQ. 1) THEN
IX1 = IX 40.13
IY1 = IY 40.13
IX2 = IX+1 40.13
IY2 = IY 40.13
IX3 = IX+1 40.13
IY3 = IY+1 40.13
ELSE IF (ISIDE .EQ. 2) THEN
IX1 = IX+1 40.13
IY1 = IY 40.13
IX2 = IX+1 40.13
IY2 = IY+1 40.13
IX3 = IX 40.13
IY3 = IY+1 40.13
ELSE IF (ISIDE .EQ. 3) THEN
IX1 = IX+1 40.13
IY1 = IY+1 40.13
IX2 = IX 40.13
IY2 = IY+1 40.13
IX3 = IX 40.13
IY3 = IY 40.13
ELSE IF (ISIDE .EQ. 4) THEN
IX1 = IX 40.13
IY1 = IY+1 40.13
IX2 = IX 40.13
IY2 = IY 40.13
IX3 = IX+1 40.13
IY3 = IY 40.13
ENDIF
K1 = KGRPNT(IX1,IY1) 40.13
XC1 = XCGRID(IX1,IY1) 40.13 30.72
YC1 = YCGRID(IX1,IY1) 40.13 30.72
K2 = KGRPNT(IX2,IY2) 40.13
XC2 = XCGRID(IX2,IY2) 40.13 30.72
YC2 = YCGRID(IX2,IY2) 40.13 30.72
K3 = KGRPNT(IX3,IY3) 40.13
XC3 = XCGRID(IX3,IY3) 30.72
YC3 = YCGRID(IX3,IY3) 30.72
DET = 0.
IF (K1 .GE. 2 .AND. K2 .GE. 2 .AND. K3 .GE. 2) THEN
DET = ((XC3 - XC1) * (YC2 - YC1)) -
& ((YC3 - YC1) * (XC2 - XC1))
IF (DET .EQ. 0.) THEN
! three grid points on one line 40.13
CALL MSGERR (2,'3 comp. grid points on one line') 40.13
WRITE (PRINTF, 112)
& IX1-1, IY1-1, XC1+XOFFS, YC1+YOFFS, 40.13
& IX2-1, IY2-1, XC2+XOFFS, YC2+YOFFS, 40.13
& IX3-1, IY3-1, XC3+XOFFS, YC3+YOFFS 40.13
112 FORMAT (3(1X, 2I3, 2(1X, F14.4))) 40.13
ENDIF
!
IF (FIRST) THEN
FIRST = .FALSE.
IF (DET .GT. 0.) THEN
CVLEFT = .FALSE.
ELSE
CVLEFT = .TRUE.
ENDIF
ENDIF
IF ( ( CVLEFT .AND. DET .GT. 0.)
& .OR. (.NOT. CVLEFT .AND. DET .LT. 0.)) THEN
! crossing grid lines in a mesh 40.13
CALL MSGERR (2,'Grid angle <0 or >180 degrees') 40.13
WRITE (PRINTF, 112)
& IX1-1, IY1-1, XC1+XOFFS, YC1+YOFFS, 40.13
& IX2-1, IY2-1, XC2+XOFFS, YC2+YOFFS, 40.13
& IX3-1, IY3-1, XC3+XOFFS, YC3+YOFFS 40.13
ENDIF
ENDIF
20 CONTINUE
15 CONTINUE
10 CONTINUE
RETURN
! *** end of subroutine CVCHEK ***
END
!***********************************************************************
! *
SUBROUTINE CVMESH (XP, YP, XC, YC, KGRPNT, XCGRID ,YCGRID, KGRBND)
! *
!***********************************************************************
!
USE OCPCOMM4 40.41
USE SWCOMM2 40.41
USE SWCOMM3 40.41
USE SWCOMM4 40.41
#ifdef SWAN_MODEL
USE M_CVMESH
USE M_COUPLING
USE M_PARALL
#endif
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 30.72: IJsbrand Haagsma
! 40.00, 40.13: Nico Booij
! 40.02: IJsbrand Haagsma
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! 30.21, Jun. 96: New for curvilinear version
! 30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN
! 40.00, May 98: procedure for points outside grid accelerated
! 40.00, Feb 99: procedure extended for 1D case
! XOFFS and YOFFS added in write statements
! 40.02, Mar. 00: Fixed bug that placed dry testpoints outside computational grid
! 40.13, Mar. 01: message "CVMESH 2nd attempt .." suppressed
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
! 40.41, Nov. 04: search for boundary points improved
!
! 2. Purpose
!
! procedure to find location in curvilinear grid for a point
! given in problem coordinates
!
! 3. Method
!
! First attempt: use Newton-Raphson method to find XC and YC
! (Note: in the program XC and YC indicate the mesh and position in
! the mesh) in a few steps; this may be most efficient if a series of
! points is processed, because the previous point provides a good
! first estimate.
! This procedure may fail if the number of iterations is larger than
! a previously set limit (default=5).
!
! If the first attempt fails then determine whether the points (XP,YP)
! is inside the mesh. If so, then the Newton-Raphson procedure is used
! again with the pivoting point like first guess. Otherwise, scan the
! boundaries whether the point is on the boundaries. If this fails, it
! may be concluded that the point (XP,YP) is outside the grid.
!
! 4. Argument variables
!
! XCGRID input Coordinates of computational grid in x-direction 30.72
! YCGRID input Coordinates of computational grid in y-direction 30.72
! XP, YP input a point given in problem coordinates
! XC, YC outp same point in computational grid coordinates
!
REAL XCGRID(MXC,MYC), YCGRID(MXC,MYC) 30.72
REAL XP, YP, XC, YC
!
! KGRPNT input array(MXC,MYC) grid numbers
! if KGRPNT <= 1, point is not in comp. grid.
! KGRBND input lists all boundary grid points consecutively
!
INTEGER KGRPNT(MXC,MYC), KGRBND(*) 40.00
!
! Local variables
!
! MXITNR number of iterations in Newton-Raphson procedure
! IX, IY counter of computational grid point
! K1 address of grid point
! IXMIN counter of grid point closest to (XP,YP)
! IYMIN counter of grid point closest to (XP,YP)
! IBND counter of boundary grid points
!
INTEGER :: IX, IY, K1, IXMIN, IYMIN, IBND
#if !defined SWAN_MODEL
INTEGER, SAVE :: MXITNR = 0
#endif
INTEGER, SAVE :: IENT = 0
!
! INMESH if True, point (XP,YP) is inside the computational grid
! FINDXY if True, Newton-Raphson procedure succeeded
! ONBND if True, given point is on boundary 40.41
!
LOGICAL INMESH ,FINDXY, ONBND
!
! DISMIN minimal distance found 40.41
! XPC1 user coordinate of a computational grid point
! YPC1 user coordinate of a computational grid point
! XC0 grid coordinate of grid point closest to (XP,YP)
! YC0 grid coordinate of grid point closest to (XP,YP)
!
REAL :: DISMIN 40.41
REAL :: XPC1, YPC1, XC0, YC0
#ifdef SWAN_MODEL
MXITNR=>MXITNR_G(Gridnum)
#endif
!
! 5. SUBROUTINES CALLING
!
! SINCMP
!
! 6. SUBROUTINES USED
!
! NEWTON
!
! 7. ERROR MESSAGES
!
! ---
!
! 8. REMARKS
!
! ---
!
! 9. STRUCTURE
!
! --------------------------------------------------------------
! Determine XC and YC from XP and YP using Newton-Raphson iteration
! process
! If (XC and YC were found) then
! Procedure is ready; Return values of XC and YC
! return
! else
! ---------------------------------------------------------------------
! For ix=1 to MXC-1 do
! For iy=1 to MYC-1 do
! Inmesh = True
! For iside=1 to 4 do
! Case iside=
! 1: K1 = KGRPNT(ix,iy), K2 = KGRPNT(ix+1,iy)
! 2: K1 = KGRPNT(ix+1,iy), K2 = KGRPNT(ix+1,iy+1)
! 3: K1 = KGRPNT(ix+1,iy+1), K2 = KGRPNT(ix,iy+1)
! 4: K1 = KGRPNT(ix,iy+1), K2 = KGRPNT(ix,iy)
! ----------------------------------------------------------
! If K1>0 and K2>0
! Then Det = (xp-xpg(K1))*(ypg(K2)-ypg(K1)) -
! (yp-ypg(K1))*(xpg(K2)-xpg(K1))
! If ((CVleft and Det>0) or (not CVleft and Det<0))
! Then Make Inmesh = False
! Else Inmesh = true and XC = IX and YC = IY
! Else Make Inmesh = False
! --------------------------------------------------------
! If Inmesh
! Then Determine XC and YC using Newton-Raphson iteration
! process
! Procedure is ready; Return values of XC and YC
! ---------------------------------------------------------------------
! No mesh is found: Make XC and YC = exception value
! Return values of XC and YC
! ---------------------------------------------------------------------
!
!****************************************************************
!
!
IF (LTRACE) CALL STRACE (IENT,'CVMESH')
!
IF (ONED) THEN
CALL NEWT1D (XP, YP, XCGRID, YCGRID, KGRPNT, 40.00
& XC ,YC ,FINDXY)
IF (.NOT.FINDXY) THEN
XC = -99.
YC = -99.
IF (ITEST .GE. 150 .OR. INTES .GE. 20) THEN
WRITE(PRINTF, 85) XP+XOFFS, YP+YOFFS 40.00
ENDIF
ENDIF
GOTO 99
ELSE
! two-dimensional computation
XC = 1.
YC = 1.
! --- First attempt, to find XC,YC with Newton-Raphson method
MXITNR = 5
CALL NEWTON (XP, YP, XCGRID, YCGRID, 40.00
& MXITNR ,ITER, XC ,YC ,FINDXY) 40.41 40.02
IF ((ITEST .GE. 150 .OR. INTES .GE. 20) .AND. FINDXY) THEN 40.02
WRITE(PRINTF,25) XP+XOFFS ,YP+YOFFS ,XC ,YC 40.03
ENDIF
25 FORMAT (' CVMESH: (XP,YP)=','(',F12.4,',',F12.4,
& '), (XC,YC)=','(',F9.2,',',F9.2,')')
IF (FINDXY) GOTO 80
!
IF (INMESH (XP, YP, XCGRID ,YCGRID, KGRBND)) THEN
! --- select grid point closest to (XP,YP)
DISMIN = 1.E20
DO 50 IX = 1,MXC
DO 40 IY = 1,MYC
K1 = KGRPNT(IX,IY)
IF (K1.GT.1) THEN
XPC1 = XCGRID(IX,IY)
YPC1 = YCGRID(IX,IY)
DISXY = SQRT ((XP-XPC1)**2 + (YP-YPC1)**2)
IF (DISXY .LT. DISMIN) THEN
IXMIN = IX
IYMIN = IY
DISMIN = DISXY
ENDIF
ENDIF
40 CONTINUE
50 CONTINUE
! second attempt using closest grid point as first guess
MXITNR = 20
XC0 = REAL(IXMIN)
YC0 = REAL(IYMIN)
! ITEST condition changed from 20 to 120 40.13
IF (ITEST.GE.120) WRITE (PRTEST, 55) XP+XOFFS ,YP+YOFFS , 40.13
& XC0-1. ,YC0-1.
55 FORMAT (' CVMESH 2nd attempt, (XP,YP)=','(',F12.4,',',F12.4,
& '), (XC,YC)=','(',F9.2,',',F9.2,')')
DO KORNER = 1, 4
IF (KORNER.EQ.1) THEN
XC = XC0 + 0.2
YC = YC0 + 0.2
ELSE IF (KORNER.EQ.2) THEN
XC = XC0 - 0.2
YC = YC0 + 0.2
ELSE IF (KORNER.EQ.3) THEN
XC = XC0 - 0.2
YC = YC0 - 0.2
ELSE
XC = XC0 + 0.2
YC = YC0 - 0.2
ENDIF
CALL NEWTON (XP, YP, XCGRID, YCGRID, 40.00
& MXITNR ,ITER, XC ,YC ,FINDXY) 40.41 40.02
IF (FINDXY) THEN
IF (ITEST .GE. 150 .OR. INTES .GE. 20) THEN
WRITE(PRINTF,25) XP+XOFFS ,YP+YOFFS ,XC ,YC 40.00
ENDIF
GOTO 80
ENDIF
ENDDO
IF (ITER.GE.MXITNR) THEN 40.41
WRITE (PRINTF, 75) XP+XOFFS, YP+YOFFS, MXITNR 40.00
75 FORMAT (' search for point with location ', 2F12.4, 40.41
& ' fails in', I3, ' iterations') 40.41
END IF 40.41
ELSE
! scan boundary to see whether the point is close to the boundary
DISMIN=99999. 40.41
ONBND =.FALSE. 40.41
IX1 = 0 40.41
IY1 = 0 40.51
IX2 = 0
DO IBND = 1, NGRBND
IF (IX2.NE.0) THEN 40.41
IX1 = IX2
IY1 = IY2
XP1 = XP2
YP1 = YP2
END IF
IX2 = KGRBND(2*IBND-1)
IY2 = KGRBND(2*IBND)
IF (IX2.NE.0 .AND. (ABS(IX2-IX1).GT.1 .OR. 40.51
& ABS(IY2-IY1).GT.1)) IX1 = 0 40.51
IF (IX2.GT.0) THEN
XP2 = XCGRID(IX2,IY2)
YP2 = YCGRID(IX2,IY2)
IF (IBND.GT.1 .AND. IX1.GT.0) THEN 40.51
! --- determine relative distance from boundary segment
! with respect to the length of that segment
SLEN2 = (XP2-XP1)**2 + (YP2-YP1)**2
RELDIS = ABS((XP-XP1)*(YP2-YP1)-(YP-YP1)*(XP2-XP1)) /
& SLEN2
IF (RELDIS.LT.0.01) THEN 40.41
! --- determine location on the boundary section
IF (RELDIS-DISMIN.LE.0.01) THEN 40.41
DISMIN = RELDIS 40.41
RELLOC = ((XP-XP1)*(XP2-XP1)+(YP-YP1)*(YP2-YP1)) /
& SLEN2
IF (RELLOC.GE.-0.001 .AND. RELLOC.LE.1.001) THEN 40.41
RELLCM = RELLOC 40.41
IF (RELLCM.LT.0.01) RELLCM=0. 40.41
IF (RELLCM.GT.0.99) RELLCM=1. 40.41
IX1M = IX1 40.41
IX2M = IX2 40.41
IY1M = IY1 40.41
IY2M = IY2 40.41
ONBND = .TRUE. 40.41
ENDIF 40.41
ENDIF
ENDIF
ENDIF
ENDIF
ENDDO
IF (ONBND) THEN 40.41
XC = FLOAT(IX1M) + RELLCM * FLOAT(IX2M-IX1M) - 1. 40.41
YC = FLOAT(IY1M) + RELLCM * FLOAT(IY2M-IY1M) - 1. 40.41
IF (ITEST .GE. 150 .OR. INTES .GE. 20) THEN
WRITE(PRINTF, 65) XP+XOFFS, YP+YOFFS, XC, YC 40.00
65 FORMAT (' CVMESH: (XP,YP)=','(',F12.4,',',F12.4,
& ') is on the boundary, (XC,YC)=(',
& F9.2,',',F9.2,')')
ENDIF
GOTO 80
ENDIF
XC = -99.
YC = -99.
IF (ITEST .GE. 150 .OR. INTES .GE. 20) THEN
WRITE(PRINTF, 85) XP+XOFFS, YP+YOFFS 40.00
85 FORMAT (' CVMESH: (XP,YP)=','(',F12.4,',',F12.4,
& ') is outside grid')
ENDIF
GOTO 99
ENDIF
ENDIF 40.00
80 IF (KGRPNT(INT(XC+3.001)-2,INT(YC+3.001)-2).LE.1) THEN
WRITE (PRINTF, 90) XP+XOFFS, YP+YOFFS 40.41
90 FORMAT (' point with location ',2F12.4,' is not active') 40.41
XC = -99.
YC = -99.
ENDIF
99 RETURN
END
!***********************************************************************
! *
LOGICAL FUNCTION INMESH (XP, YP, XCGRID ,YCGRID, KGRBND)
! *
!***********************************************************************
!
USE OCPCOMM4 40.41
USE SWCOMM2 40.41
USE SWCOMM3 40.41
USE M_PARALL 40.41
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! Nico Booij
! 40.41: Marcel Zijlema
! 40.51: Marcel Zijlema
!
! 1. Updates
!
! New function for curvilinear version (ver. 40.00). May '98
! 40.03, Dec 99: test output added; commons swcomm2 and ocpcomm4 added
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
! 40.41, Nov. 04: search for points restricted to subdomain
! 40.51, Feb. 05: determining number of crossing points improved
!
! 2. Purpose
!
! procedure to find whether a given location is
! in the (curvilinear) computational grid
!
! 3. Method suggested by Gerbrant van Vledder
!
! draw a line from the point (XP,YP) in vertical direction
! determine the number of crossings with the boundary of the
! grid; if this number is even the point is outside
!
! 4. Argument variables
!
!
! KGRBND int input array containing boundary grid points
!
INTEGER KGRBND(*)
!
! XP, YP real, input a point given in problem coordinates
! XCGRID real, input array(IX,IY) x-coordinate of a grid point
! YCGRID real, input array(IX,IY) y-coordinate of a grid point
!
REAL XCGRID(MXC,MYC) ,YCGRID(MXC,MYC),
& XP, YP
!
! 5. Parameter variables
!
! 6. Local variables
!
! NUMCRS number of crossings with boundary outline
!
INTEGER NUMCRS, IX1, IY1, IX2, IY2
REAL XP1, XP2, YP1, YP2, YPS, RELDIS, RELDO
!
! 8. Subroutines used
!
! 9. Subroutines calling
!
! CVMESH
!
! 10. Error messages
!
! 11. Remarks
!
! 12. Structure
!
! --------------------------------------------------------------
! numcros = 0
! For all sections of the boundary do
! determine coordinates of end points (XP1,YP1) and (XP2,YP2)
! If (XP1<XP and XP2>XP) or (XP1>XP and XP2<XP)
! then If not (YP1<YP and YP2<YP)
! if YPS>YP
! then numcros = numcros + 1
! ---------------------------------------------------------------
! If numcros is even
! Then Inmesh = False
! Else Inmesh = True
! ---------------------------------------------------------------
!
! 13. Source text
!
SAVE IENT
DATA IENT/0/
CALL STRACE (IENT,'INMESH')
!
IF (XP.LT.XCLMIN .OR. XP.GT.XCLMAX .OR. 40.41
& YP.LT.YCLMIN .OR. YP.GT.YCLMAX) THEN 40.41
IF (ITEST.GE.70) WRITE (PRTEST, 22) XP+XOFFS, YP+YOFFS, 40.03
& XCLMIN+XOFFS, XCLMAX+XOFFS, YCLMIN+YOFFS, YCLMAX+YOFFS 40.41
22 FORMAT (1X, 2F12.4, ' is outside region ', 4F12.4)
INMESH = .FALSE.
GOTO 90
ENDIF
!
IF (NGRBND.LE.0) THEN
CALL MSGERR (3, 'grid outline not yet determined')
RETURN
ENDIF
!
NUMCRS = 0
IX1 = 0
IY1 = 0
IX2 = 0
RELDIS = -1.
! loop over the boundary of the computational grid
DO IBND = 1, NGRBND
IF (IX2.NE.0) THEN
IX1 = IX2
IY1 = IY2
XP1 = XP2
YP1 = YP2
END IF
IX2 = KGRBND(2*IBND-1)
IY2 = KGRBND(2*IBND)
IF (IX2.NE.0 .AND. (ABS(IX2-IX1).GT.1 .OR. 40.51
& ABS(IY2-IY1).GT.1)) IX1 = 0 40.51
IF (IX2.GT.0) THEN
XP2 = XCGRID(IX2,IY2)
YP2 = YCGRID(IX2,IY2)
IF (ITEST.GE.180) WRITE (PRTEST, 28) XP2+XOFFS, 40.03
& YP2+YOFFS
28 FORMAT (' boundary point ', 2F12.4)
IF (IBND.GT.1 .AND. IX1.GT.0) THEN 40.51
IF (((XP1.GT.XP).AND.(XP2.LE.XP)).OR.
& ((XP1.LE.XP).AND.(XP2.GT.XP))) THEN
IF (YP1.GT.YP .OR. YP2.GT.YP) THEN
! determine y-coordinate of crossing point
YPS = YP1 + (XP-XP1) * (YP2-YP1) / (XP2-XP1)
! determine relative distance from boundary segment 40.51
! with respect to the length of that segment 40.51
RELDO = RELDIS 40.51
RELDIS = ABS(YP-YPS) / SQRT((XP2-XP1)**2 + (YP2-YP1)**2) 40.51
IF (YPS.GT.YP.AND.ABS(RELDIS-RELDO).GT.0.1) THEN 40.51
NUMCRS = NUMCRS + 1
IF (ITEST.GE.70) WRITE (PRTEST, 32) NUMCRS, 40.03
& XP+XOFFS, YP+YOFFS, YPS+YOFFS 40.03
32 FORMAT (' crossing ', I1, ' point ', 3F12.4)
ENDIF
ENDIF
ENDIF
ENDIF
ENDIF
ENDDO
! point is inside the grid is number of crossings is odd
IF (MOD(NUMCRS,2) .EQ. 1) THEN
INMESH = .TRUE.
ELSE
INMESH = .FALSE.
ENDIF
90 RETURN
END
!***********************************************************************
! *
SUBROUTINE NEWTON (XP, YP, XCGRID, YCGRID, 40.00
& MXITNR, ITER, XC, YC, FIND) 40.41 40.02
! *
!***********************************************************************
!
USE OCPCOMM4 40.41
USE SWCOMM3 40.41
USE SWCOMM4 40.41
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
! 0. Authors
!
! 30.72: IJsbrand Haagsma
! 30.80: Nico Booij
! 30.82: IJsbrand Haagsma
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! 30.21, Jun. 96: New for curvilinear version
! 30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN
! 30.82, Oct. 98: Updated description of several variables
! 30.80, Oct. 98: computation of update of XC,YC modified to avoid
! division by 0
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. Purpose
!
! Solve eqs. and find a point (XC,YC) in a curvilinear grid (compt.
! grid) for a given point (XP ,YP) in a cartesian grid (problem coord).
!
! 3. Method
!
! In this subroutine the next equations are solved :
!
! @XP @XP
! XP(xc,yc) + --- * @XC + --- * @YC - XP(x,y) = 0
! @XC @YC
!
! @YP @YP
! YP(xc,yc) + --- * @XC + --- * @YC - YP(x,y) = 0
! @XC @YC
!
! In the subroutine, next notation is used for the previous eqs.
! XVC + DXDXC * DXC + DXDYC * DYC - XP = 0.
! YVC + DYDXC * DXC + DYDYC * DYC - YP = 0.
!
!
! 4. Argument variables
!
! i MXITNR: Maximum number of iterations 30.82
!
INTEGER MXITNR, ITER 40.41 40.02
!
! o XC : X-coordinate in computational coordinates 30.82
! i XCGRID: Coordinates of computational grid in x-direction 30.72
! i XP : X-coordinate in problem coordinates 30.82
! o YC : Y-coordinate in computational coordinates 30.82
! i YCGRID: Coordinates of computational grid in y-direction 30.72
! i YP : Y-coordinate in problem coordinates 30.82
!
REAL XC, XCGRID(MXC,MYC), XP 30.82
REAL YC, YCGRID(MXC,MYC), YP 30.82
!
! o FIND : Whether XC and YC are found 30.82
!
LOGICAL FIND 30.82
!
! 6. SUBROUTINES USED
!
! STRACE
!
!
! 7. ERROR MESSAGES
!
! ---
!
! 8. REMARKS
!
! ---
!
! 9. STRUCTURE
!
! -----------------------------------------------------------------
! -----------------------------------------------------------------
!
! 13. Source text
!
INTEGER, SAVE :: IENT = 0
IF (LTRACE) CALL STRACE (IENT,'NEWTON')
!
DXC = 1000.
DYC = 1000.
TOLDC = 0.001
FIND = .FALSE.
!
IF (ITEST .GE. 200) THEN
WRITE(PRINTF,*) ' Coordinates in subroutine NEWTON '
DO J = 1, MYC
DO I = 1, MXC
WRITE(PRINTF,30) I ,J ,XCGRID(I,J) ,YCGRID(I,J) 30.72
ENDDO
ENDDO
ENDIF
30 FORMAT(2(2X,I5),2(2X,E12.4))
!
DO 14 K = 1 ,MXITNR
ITER = K 40.41
I1 = INT(XC) 40.00
J1 = INT(YC)
IF (I1 .EQ. MXC) I1 = I1 - 1
IF (J1 .EQ. MYC) J1 = J1 - 1
I2 = I1 + 1
J2 = J1 + 1
FJ1 = FLOAT(J1)
FI1 = FLOAT(I1)
FJ2 = FLOAT(J2)
FI2 = FLOAT(I2)
!
XVC = (YC-FJ1)*((XC-FI1)*XCGRID(I2,J2) +
& (FI2-XC)*XCGRID(I1,J2)) +
& (FJ2-YC)*((XC-FI1)*XCGRID(I2,J1) +
& (FI2-XC)*XCGRID(I1,J1))
YVC = (YC-FJ1)*((XC-FI1)*YCGRID(I2,J2) +
& (FI2-XC)*YCGRID(I1,J2)) +
& (FJ2-YC)*((XC-FI1)*YCGRID(I2,J1) +
& (FI2-XC)*YCGRID(I1,J1))
DXDXC = (YC -FJ1)*(XCGRID(I2,J2) - XCGRID(I1,J2)) +
& (FJ2-YC )*(XCGRID(I2,J1) - XCGRID(I1,J1))
DXDYC = (XC -FI1)*(XCGRID(I2,J2) - XCGRID(I2,J1)) +
& (FI2-XC )*(XCGRID(I1,J2) - XCGRID(I1,J1))
DYDXC = (YC -FJ1)*(YCGRID(I2,J2) - YCGRID(I1,J2)) +
& (FJ2-YC )*(YCGRID(I2,J1) - YCGRID(I1,J1))
DYDYC = (XC -FI1)*(YCGRID(I2,J2) - YCGRID(I2,J1)) +
& (FI2-XC )*(YCGRID(I1,J2) - YCGRID(I1,J1))
!
IF (ITEST .GE. 150)
& WRITE(PRINTF,35) K, XC-1., YC-1., XP, YP, XVC, YVC 40.00
35 FORMAT(' NEWTON iter=', I2, ' (XC,YC)=', 2(1X,F10.2),/, 40.00
& ' (XP,YP)=', 2(1X,F10.2),
& ' X,Y(XC,YC) = ', 2(1X,F10.2))
IF (ITEST .GE. 180) WRITE(PRINTF,36)
& XCGRID(I1,J1), XCGRID(I1,J2), XCGRID(I2,J1), XCGRID(I2,J2),
& YCGRID(I1,J1), YCGRID(I1,J2), YCGRID(I2,J1), YCGRID(I2,J2),
& DXDXC, DXDYC, DYDXC, DYDYC 40.00
36 FORMAT(' NEWTON grid coord:', 8(1x, F10.0), /
& ' deriv=', 4(1X,F10.2)) 40.00
!
! *** the derivated terms of the eqs. are evaluated and ***
! *** the eqs. are solved ***
DDEN = DXDXC*DYDYC - DYDXC*DXDYC 30.80
DXP = XP - XVC 30.80
DYP = YP - YVC 30.80
IF ( DDEN.NE.0. ) THEN
DXC = ( DYDYC*DXP - DXDYC*DYP) / DDEN 30.80
DYC = (-DYDXC*DXP + DXDXC*DYP) / DDEN 30.80
ENDIF
!
XC = XC + DXC
YC = YC + DYC
!
! *** If the guess point (XC,YC) is outside of compt. ***
! *** grid, put that point in the closest boundary ***
IF (XC .LT. 1. ) XC = 1.
IF (YC .LT. 1. ) YC = 1.
IF (XC .GT. MXC) XC = FLOAT(MXC)
IF (YC .GT. MYC) YC = FLOAT(MYC)
!
IF (ITEST .GE. 120 .OR. INTES .GE. 50 .OR. IOUTES .GE. 50)
& WRITE(PRINTF,42) DXC, DYC, XC-1., YC-1. 40.00
42 FORMAT(' (DXC,DYC)=', 2(1X,F10.2), ' (XC,YC)=', 2(1X,F10.2)) 40.00
!
! *** If the accuracy is reached stop the iteration, ***
IF (ABS(DXC) .LE. TOLDC .AND. ABS(DYC) .LE. TOLDC) THEN
!
FIND = .TRUE.
XC = XC -1.
YC = YC -1.
RETURN
ENDIF
!
14 CONTINUE
RETURN
! *** end of subroutine NEWTON ***
END
!***********************************************************************
! *
SUBROUTINE NEWT1D (XP, YP, XCGRID, YCGRID, KGRPNT, 40.00
& XC, YC, FIND)
! *
!***********************************************************************
!
USE OCPCOMM4 40.41
USE SWCOMM2 40.41
USE SWCOMM3 40.41
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
! 0. Authors
!
! 40.00, 40.13: Nico Booij
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! 40.00, Feb. 99: New (adaptation from subr NEWTON for 1D case)
! 40.13, Feb. 01: DX and DY renamed to DELX and DELY (DX and DY are
! common var.); error in expression for RS corrected
! PRINTF replaced by PRTEST in test output
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. Purpose
!
! Finds broken coordinate XC for a given point XP in a rectilinear grid
!
! 3. Method
!
! In this subroutine the step on the computational grid is selected
! for which
!
! (X-X1).(X2-X1)
! 0 <= --------------- <= 1
! (X2-X1).(X2-X1)
!
! where X, X1 and X2 are vectors; X corresponds to (Xp,Yp)
! X1 and X2 are two neighbouring grid points
!
! 4. Argument variables
!
! i KGRPNT: Grid adresses 40.00
!
INTEGER KGRPNT(MXC,MYC) 40.00
!
! o XC : X-coordinate in computational coordinates 30.82
! i XCGRID: Coordinates of computational grid in x-direction 30.72
! i XP : X-coordinate in problem coordinates 30.82
! o YC : Y-coordinate in computational coordinates 30.82
! i YCGRID: Coordinates of computational grid in y-direction 30.72
! i YP : Y-coordinate in problem coordinates 30.82
!
REAL XC, XCGRID(MXC,MYC), XP 30.82
REAL YC, YCGRID(MXC,MYC), YP 30.82
!
! o FIND : Whether XC and YC are found 30.82
!
LOGICAL FIND 30.82
!
! Local variables:
REAL :: DELX, DELY ! grid line 40.13
!
! 6. SUBROUTINES USED
!
! ---
!
! 7. ERROR MESSAGES
!
! ---
!
! 8. REMARKS
!
! ---
!
! 9. STRUCTURE
!
! -----------------------------------------------------------------
! -----------------------------------------------------------------
!
! 13. Source text
!
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'NEWT1D')
!
IF (ITEST .GE. 120) THEN 40.13
WRITE(PRTEST,*) ' Coordinates in subroutine NEWT1D ' 40.13
DO I = 1, MXC
WRITE(PRTEST,30) I, XCGRID(I,1)+XOFFS ,YCGRID(I,1)+YOFFS 40.13
ENDDO
ENDIF
30 FORMAT(2X,I5,2(2X,E12.4))
!
FIND = .FALSE.
DO 40 IX = 2 ,MXC
IF (KGRPNT(IX-1,1).GT.1) THEN
X1 = XCGRID(IX-1,1)
Y1 = YCGRID(IX-1,1)
ELSE
GOTO 40
ENDIF
IF (KGRPNT(IX,1).GT.1) THEN
X2 = XCGRID(IX,1)
Y2 = YCGRID(IX,1)
ELSE
GOTO 40
ENDIF
! both ends of the step are valid grid points
! now verify whether projection of (Xp,Yp) is within the step
DELX = X2 - X1 40.13
DELY = Y2 - Y1 40.13
RS = ((XP - X1) * DELX + (YP - Y1) * DELY) / 40.13
& (DELX * DELX + DELY * DELY) 40.13
IF (RS.GE.0. .AND. RS.LE.1.) THEN
FIND = .TRUE.
XC = REAL(IX-2) + RS 40.00
YC = 0.
GOTO 50
ENDIF
40 CONTINUE
50 RETURN
! *** end of subroutine NEWT1D ***
END
!***********************************************************************
! *
SUBROUTINE EVALF (XC ,YC ,XVC ,YVC ,XCGRID ,YCGRID) 30.72
! *
!***********************************************************************
!
USE OCPCOMM4 40.41
USE SWCOMM3 40.41
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 30.72: IJsbrand Haagsma
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! 30.21, Jun. 96: New for curvilinear version
! 30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. Purpose
!
! Evaluate the coordinates (in problem coordinates) of point (XC,YC)
! given in computational coordinates
!
! 3. Method
!
! Bilinear interpolation
!
! 4. Argument variables
!
! XCGRID: input Coordinates of computational grid in x-direction 30.72
! YCGRID: input Coordinates of computational grid in y-direction 30.72
!
REAL XCGRID(MXC,MYC), YCGRID(MXC,MYC) 30.72
!
! XC, YC real, outp point in computational grid coordinates
! XVC, YCV real, OUTP same point but in problem coordinates
!
! 6. SUBROUTINES USED
!
! none
!
! 7. ERROR MESSAGES
!
! ---
!
! 8. REMARKS
!
!
! 9. STRUCTURE
!
! -----------------------------------------------------------------
! -----------------------------------------------------------------
!
! 10. SOURCE TEXT
!
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'EVALF')
!
I = INT(XC)
J = INT(YC)
!
! *** If the guess point (XC,YC) is in the boundary ***
! *** where I = MXC or/and J = MYC the interpolation ***
! *** is done in the mesh with pivoting point ***
! *** (MXC-1, J) or/and (I,MYC-1) ***
!
IF (I .EQ. MXC) I = I - 1
IF (J .EQ. MYC) J = J - 1
T = XC - FLOAT(I)
U = YC - FLOAT(J)
! *** For x-coord. ***
P1 = XCGRID(I,J) 30.72
P2 = XCGRID(I+1,J) 30.72
P3 = XCGRID(I+1,J+1) 30.72
P4 = XCGRID(I,J+1) 30.72
XVC = (1.-T)*(1.-U)*P1+T*(1.-U)*P2+T*U*P3+(1.-T)*U*P4
! *** For y-coord. ***
P1 = YCGRID(I,J) 30.72
P2 = YCGRID(I+1,J) 30.72
P3 = YCGRID(I+1,J+1) 30.72
P4 = YCGRID(I,J+1) 30.72
YVC = (1.-T)*(1.-U)*P1+T*(1.-U)*P2+T*U*P3+(1.-T)*U*P4
RETURN
! *** end of subroutine EVALF ***
END
!
!***********************************************************************
! *
SUBROUTINE SWOBST (XCGRID, YCGRID, KGRPNT, CROSS) 40.31 30.70
! *
!***********************************************************************
USE OCPCOMM4 40.41
USE SWCOMM3 40.41
USE M_OBSTA 40.31
IMPLICIT NONE 40.04
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 30.70
! 30.72 IJsbrand Haagsma
! 30.74 IJsbrand Haagsma
! 40.04 Annette Kieftenburg
! 40.28 Annette Kieftenburg
! 40.31 Marcel Zijlema
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! 30.70, Feb. 98: check if neighbouring point is a true grid point
! loop over grid points moved from calling routine into this
! argument list changed
! 30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN
! 40.04, Nov. 99: IMPLICIT NONE added, header updated
! : Removed include files that are not used
! 40.28, Feb. 02: Adjustments for extended REFLECTION option
! 40.31, Oct. 03: changes w.r.t. obstacles
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. Purpose
!
! Obtains all the data required to find obstacles and
! use subroutine TCROSS to find them
!
! 3. Method
!
! 4. Argument variables
!
! CROSS output Array which contains 0's if there is no
! obstacle crossing
! if an obstacle is crossing between the
! central point and its neighbour CROSS is equal
! to the number of the obstacle
! KGRPNT input Indirect addressing for computational grid points
! XCGRID input Coordinates of computational grid in x-direction 30.72
! YCGRID input Coordinates of computational grid in y-direction 30.72
!
INTEGER KGRPNT(MXC,MYC), CROSS(2,MCGRD)
REAL XCGRID(MXC,MYC), YCGRID(MXC,MYC) 30.72
!
! 5. Parameter variables
!
! 6. Local variables
!
! ICC index
! ICGRD index
! IENT number of entries of this subroutine
! ILINK indicates which link is analyzed: 1 -> neighbour in x
! 2 -> neighbour in y
! IX counter of gridpoints in x-direction
! IY counter of gridpoints in y-direction
! JJ counter for number of obstacles
! JP counter for number of corner points of obstacles
! NUMCOR number of corner points of obstacle
! X1, Y1 user coordinates of one end of grid link
! X2, Y2 user coordinates of other end of grid link
! X3, Y3 user coordinates of one end of obstacle side
! X4, Y4 user coordinates of other end of obstacle side
!
INTEGER ICC, ICGRD, IENT, ILINK, IX, IY, JJ, JP
INTEGER NUMCOR
REAL X1, X2, X3, X4, Y1, Y2, Y3, Y4
LOGICAL XONOBST 40.04
TYPE(OBSTDAT), POINTER :: COBST 40.31
!
! 8. Subroutines used
!
! TCROSS 40.04
! STRACE
!
LOGICAL TCROSS 40.04
!
! 9. Subroutines calling
!
! SWPREP
!
! 10. Error messages
!
! 11. Remarks
!
! 12. Structure
! ----------------------------------------------------------------
! Read number of obstacles from array OBSTA
! For every obstacle do
! Read number of corners of the obstacle
! For every corner of the obstacle do
! For every grid point do
! call function TCROSS to search if there is crossing 40.04
! point 40.04
! between the line of two points of the stencil and the
! line of the corners of the obstacle.
! If there is crossing point then
! then CROSS(link,kcgrd) = number of the crossing obstacle
! else CROSS(link,kcgrd) = 0
! ----------------------------------------------------------------
!
! 13. Source text
! ======================================================================
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'SWOBST')
!
IF (NUMOBS .GT. 0) THEN
! NUMOBS is the number of obstacles ***
COBST => FOBSTAC 40.31
DO 120 JJ = 1, NUMOBS
! number of corner points of the obstacle
NUMCOR = COBST%NCRPTS
IF (ITEST.GE. 120) THEN
WRITE(PRINTF,50) JJ, NUMCOR
50 FORMAT( ' Obstacle number : ', I4,' has ',I4,' corners')
ENDIF
! *** X1 X2 X3 ETC. are the coordinates of point according ***
! *** with the scheme in the subroutine TCROSS header *** 40.04
X3 = COBST%XCRP(1) 40.31
Y3 = COBST%YCRP(1) 40.31
IF (ITEST.GE. 120) WRITE(PRINTF,30) 1,X3,Y3
DO 110 JP = 2, NUMCOR
X4 = COBST%XCRP(JP) 40.31
Y4 = COBST%YCRP(JP) 40.31
IF (ITEST.GE. 120) WRITE(PRINTF,30) JP,X4,Y4
30 FORMAT(' Corner number:', I4,' XP: ',E10.4,' YP: ',E11.4)
DO 100 IX = 1, MXC
DO 90 IY = 1, MYC
ICC = KGRPNT(IX,IY)
IF (ICC .GT. 1) THEN
X1 = XCGRID(IX,IY) 30.72
Y1 = YCGRID(IX,IY) 30.72
!
! *** "ILINK" indicates which link is analyzed. Initial ***
! *** neighbour in x , second link with neighbouring in y***
DO 80 ILINK = 1, 2
IF (ILINK.EQ.1 .AND. IX.GT.1) THEN
X2 = XCGRID(IX-1,IY) 30.72
Y2 = YCGRID(IX-1,IY) 30.72
ICGRD = KGRPNT(IX-1,IY) 30.70
ELSE IF (ILINK.EQ.2 .AND. IY.GT.1) THEN
X2 = XCGRID(IX,IY-1) 30.72
Y2 = YCGRID(IX,IY-1) 30.72
ICGRD = KGRPNT(IX,IY-1) 30.70
ELSE
ICGRD = 0
ENDIF
IF (ICGRD.GT.1) THEN 30.70
!
! *** All links are analyzed in each point otherwise the ***
! *** boundaries can be excluded ***
!
IF (TCROSS(X1, X2, X3, X4, Y1, Y2, Y3, Y4, 40.04
& XONOBST)) THEN 40.04
CROSS(ILINK,ICC) = JJ
ENDIF
ENDIF
80 CONTINUE
ENDIF
90 CONTINUE
100 CONTINUE
X3 = X4
Y3 = Y4
110 CONTINUE
IF (.NOT.ASSOCIATED(COBST%NEXTOBST)) EXIT
COBST => COBST%NEXTOBST
120 CONTINUE
ENDIF
!
RETURN
! * end of subroutine SWOBST *
END
!
!***********************************************************************
! *
SUBROUTINE SWOBSTO (XCGRID, YCGRID, XP, YP, XC, YC, KGRPNT, CROSS,
& MIP)
! *
!***********************************************************************
USE OCPCOMM4
USE SWCOMM3
USE SWCOMM4
USE M_OBSTA
IMPLICIT NONE
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering and Geosciences |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmer: Nico Booij |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 40.86: Nico Booij
!
! 1. Updates
!
! 40.86, Feb. 08: new subroutine, based on SWOBST
!
! 2. Purpose
!
! Find out whether a line segment defined in user coordinates
! crosses one or more obstacles
! use subroutine TCROSS to find crossings
!
! 3. Method
!
! 4. Argument variables
!
! CROSS output false if no obstacle crossing
! true if an obstacle is crossing between the
! output point and neighbouring grid point
! KGRPNT input Indirect addressing for computational grid points
! MIP input number of output points
! XC, YC input array containing indices of output points
! XP, YP input user coordinates of output point
! XCGRID input Coordinates of computational grid in x-direction
! YCGRID input Coordinates of computational grid in y-direction
!
INTEGER MIP
INTEGER KGRPNT(MXC,MYC)
LOGICAL CROSS(4,MIP)
REAL XC(MIP), YC(MIP), XP(MIP), YP(MIP)
REAL XCGRID(MXC,MYC), YCGRID(MXC,MYC)
!
! 5. Parameter variables
!
! 6. Local variables
!
INTEGER :: IP ! counter for number of output points
INTEGER :: JJ ! counter for number of obstacles
INTEGER :: JP ! counter for number of corner points of obstacles
INTEGER :: NUMCOR ! number of corner points of obstacle
REAL :: X1, Y1 ! user coordinates of one end of line segment
REAL :: X2, Y2 ! user coordinates of other end of line segment
REAL :: X3, Y3 ! user coordinates of one end of obstacle side
REAL :: X4, Y4 ! user coordinates of other end of obstacle side
INTEGER :: JX(1:4), JY(1:4) ! grid counters for the 4 corners
INTEGER :: JC ! corner counter
INTEGER :: JX1, JY1, JX2, JY2
LOGICAL :: OUTSID
!
TYPE(OBSTDAT), POINTER :: COBST
!
! 8. Subroutines used
!
! STRACE
!
LOGICAL :: TCROSS ! determines whether two line segments cross
!
! 9. Subroutines calling
!
! SWOUTP
!
! 10. Error messages
!
! 11. Remarks
!
! 12. Structure
!
! 13. Source text
! ======================================================================
LOGICAL :: XONOBST ! not used
INTEGER, SAVE :: IENT=0 ! number of entries of this subroutine
IF (LTRACE) CALL STRACE (IENT,'SWOBSTO')
!
CROSS = .FALSE.
!
IF (NUMOBS .GT. 0) THEN
! NUMOBS is the number of obstacles ***
COBST => FOBSTAC
DO JJ = 1, NUMOBS
! number of corner points of the obstacle
NUMCOR = COBST%NCRPTS
IF (ITEST.GE. 120) THEN
WRITE(PRINTF,50) JJ, NUMCOR
50 FORMAT( ' Obstacle number : ', I4,' has ',I4,' corners')
ENDIF
! *** X1 X2 X3 ETC. are the coordinates of point according ***
! *** with the scheme in the subroutine TCROSS header ***
X3 = COBST%XCRP(1)
Y3 = COBST%YCRP(1)
IF (ITEST.GE. 120) WRITE(PRINTF,30) 1,X3,Y3
DO JP = 2, NUMCOR
X4 = COBST%XCRP(JP)
Y4 = COBST%YCRP(JP)
IF (ITEST.GE. 120) WRITE(PRINTF,30) JP,X4,Y4
30 FORMAT(' Corner number:', I4,' XP: ',E10.4,
& ' YP: ',E11.4)
DO IP = 1, MIP
X1 = XP(IP)
Y1 = YP(IP)
!
IF (XC(IP) .LE. -0.5 .OR. YC(IP) .LE. -0.5) CYCLE
OUTSID = .FALSE.
JX1 = INT(XC(IP)+3.001) - 2
JX2 = JX1 + 1
IF (JX1.LT.0) OUTSID = .TRUE.
IF (KREPTX .EQ. 0) THEN
IF (JX1.GT.MXC) OUTSID = .TRUE.
IF (JX1.EQ.MXC) JX2 = MXC
IF (JX1.EQ.0) JX1 = 1
ELSE
JX1 = 1 + MODULO (JX1-1,MXC)
JX2 = 1 + MODULO (JX2-1,MXC)
ENDIF
IF (ONED) THEN
JY1 = 1
JY2 = 1
ELSE
JY1 = INT(YC(IP)+3.001) - 2
JY2 = JY1 + 1
IF (JY1.LT.0) OUTSID = .TRUE.
IF (JY1.GT.MYC) OUTSID = .TRUE.
IF (JY1.EQ.MYC) JY2 = MYC
IF (JY1.EQ.0) JY1 = 1
ENDIF
IF (.NOT.OUTSID) THEN
JX(1) = JX1
JY(1) = JY1
JX(2) = JX2
JY(2) = JY1
JX(3) = JX1
JY(3) = JY2
JX(4) = JX2
JY(4) = JY2
!
DO JC = 1, 4
IF (KGRPNT(JX(JC),JY(JC)).GT.1) THEN
X2 = XCGRID(JX(JC),JY(JC))
Y2 = YCGRID(JX(JC),JY(JC))
IF (TCROSS(X1, X2, X3, X4, Y1, Y2, Y3, Y4,
& XONOBST)) CROSS(JC,IP) = .TRUE.
ENDIF
ENDDO
ENDIF
ENDDO
!
X3 = X4
Y3 = Y4
ENDDO
IF (.NOT.ASSOCIATED(COBST%NEXTOBST)) EXIT
COBST => COBST%NEXTOBST
ENDDO
ENDIF
!
RETURN
! end of subroutine SWOBSTO
END
!***********************************************************************
! *
LOGICAL FUNCTION TCROSS (X1, X2, X3, X4, Y1, Y2, Y3, Y4, X1ONOBST)
! *
!***********************************************************************
!
USE OCPCOMM4 40.41
!
IMPLICIT NONE 40.04
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 40.00 Gerbrant van Vledder
! 40.04 Annette Kieftenburg
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! 30.70, Feb 98: argument list simplified
! subroutine changed into logical function
! 40.00, Aug 98: division by zero prevented
! 40.04, Aug 99: method corrected, IMPLICIT NONE added, XCONOBST added,
! introduced TINY and EPSILON (instead of comparing to 0)
! replaced 0 < LMBD,MIU by 0 <= LMBD,MIU
! XCONOBST added to argument list
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. Purpose
!
! Find if there is an obstacle crossing the stencil in used
!
! 3. Method
!
! For the next situation (A, B and C are the points in the stencil,
! D and E are corners of the obstacle
!
!
! obstacle --> D(X3,Y3)
! *
! *
! *
! (X2,Y2) * (XC,YC)
! B-----------@--------------------------A (X1,Y1)
! ^* /
! _____| * /
! | * /
! | * /
! crossing point * /
! * /
! E /
! (X4,Y4) /
! C
!
!
! The crossing point (@) should be found solving the next eqs.
! for LMBD and MIU.
!
! | XC | | X1 | | X2 - X1 |
! | | = | | + LMBD * | |
! | YC | | Y1 | | Y2 - Y1 |
!
!
! | XC | | X3 | | X4 - X3 | 40.04
! | | = | | + MIU * | |
! | YC | | Y3 | | Y4 - Y3 | 40.04
!
!
! If solution exist and (0 <= LMBD <= 1 and 0 <= MIU <= 1) 40.04
! there is an obstacle crossing the stencil
!
! 4. Argument variables
!
! X1, Y1 inp user coordinates of one end of grid link
! X2, Y2 inp user coordinates of other end of grid link
! X3, Y3 inp user coordinates of one end of obstacle side
! X4, Y4 inp user coordinates of other end of obstacle side
! X1ONOBST outp boolean which tells whether (X1,Y1) is on obstacle
!
REAL EPS, X1, X2, X3, X4, Y1, Y2, Y3, Y4
LOGICAL X1ONOBST
!
! 5. Parameter variables
!
! 6. Local variables
!
! A,B,C,D dummy variables
! DIV1 denominator of value of LMBD (or MIU)
! E,F dummy variables
! IENT number of entries
! LMBD coefficient in vector equation for stencil points (or obstacle)
! MIU coefficient in vector equation for obstacle (or stencil points)
!
INTEGER IENT
REAL A, B, C, D, DIV1, E, F, LMBD, MIU
!
! 8. Subroutines used
!
! 9. Subroutines calling
!
! SWOBST
! SWTRCF
! OBSTMOVE
!
! 10. Error messages
!
! 11. Remarks
!
! 12. Structure
!
! Calculate MIU and LMBD
! If 0 <= MIU, LMBD <= 1 40.04
! Then TCROSS is .True.
! Else TCROSS is .False.
!
! 13. Source text
! ======================================================================
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'TCROSS')
!
EPS = EPSILON(X1)*SQRT((X2-X1)*(X2-X1)+(Y2-Y1)*(Y2-Y1)) 40.04
IF (EPS ==0.) EPS = TINY(X1) 40.04
A = X2 - X1
! A not equal to zero
IF (ABS(A) .GT. TINY(X1)) THEN 40.04
B = X4 - X3
C = X3 - X1
D = Y2 - Y1
E = Y4 - Y3
F = Y3 - Y1
ELSE
! exchange MIU and LMBD 40.04
A = X4 - X3
B = X2 - X1
C = X1 - X3
D = Y4 - Y3
E = Y2 - Y1
F = Y1 - Y3
ENDIF
DIV1 = ((A*E) - (D*B)) 40.00
!
! DIV1 = 0 means that obstacle is parallel to line through 40.04
! stencil points, or (X3,Y3) = (X4,Y4); 40.04
! A = 0 means trivial set of equations X4= X3 and X2 =X1 40.04
!
IF ((ABS(DIV1).LE.TINY(X1)) .OR. 40.04
& (ABS(A).LE.TINY(X1))) THEN 40.04
MIU = -1. 40.00
LMBD = -1. 40.04
ELSE 40.00
MIU = ((D*C) - (A*F)) / DIV1
LMBD = (C + (B*MIU)) / A
END IF 40.00
!
IF (MIU .GE. 0. .AND. MIU .LE. 1. .AND. 40.04
& LMBD .GE. 0. .AND. LMBD .LE. 1.) THEN 40.04
!
! Only (X1,Y1) is checked, because of otherwise possible double 40.04
! counting 40.04
IF ((LMBD.LE.EPS .AND. ABS(X2-X1).GT.EPS).OR. 40.04
& (MIU .LE.EPS .AND. ABS(X2-X1).LE.EPS))THEN 40.04
X1ONOBST = .TRUE. 40.04
ELSE 40.04
X1ONOBST = .FALSE. 40.04
ENDIF 40.04
!
! *** test output ***
IF (ITEST .GE. 120) THEN
WRITE(PRINTF,70)X1,Y1,X2,Y2,X3,Y3,X4,Y4
70 FORMAT(' Obstacle crossing :',/,
& ' Coordinates of comp grid points and corners of obstacle:',/,
& ' P1(',E10.4,',',E10.4,')',' P2(',E10.4,',',E10.4,')',/,
& ' P3(',E10.4,',',E10.4,')',' P4(',E10.4,',',E10.4,')')
ENDIF
!
TCROSS = .TRUE.
ELSE
TCROSS = .FALSE.
ENDIF
!
! End of subroutine TCROSS
RETURN
END
!
!***********************************************************************
! *
RECURSIVE SUBROUTINE OBSTMOVE (XCGRID, YCGRID, KGRPNT) 40.31
! *
!***********************************************************************
!
USE OCPCOMM4 40.41
USE SWCOMM2 40.81
USE SWCOMM3 40.41
USE M_OBSTA 40.31
IMPLICIT NONE
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 40.09 Annette Kieftenburg
! 40.28 Annette Kieftenburg
! 40.31 Marcel Zijlema
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! 40.09, July 00: new subroutine
! 40.28, Feb. 02: Adjustments for extended REFLECTION option
! 40.31, Oct. 03: changes w.r.t. obstacles
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. Purpose
!
! Move OBSTACLE points (X3,Y3) and (X4,Y4) a bit if computational gridcell
! (X1,Y1) is on the OBSTACLE line piece.
!
! 3. Method
!
! Add EPS*(dY,-dX) to OBSTACLE line piece coordinates so that movement of
! these OBSTACLE points is perpendicular to the direction
!
! 4. Argument variables
!
! KGRPNT input Indirect addressing for computational grid points
! XCGRID input Coordinates of computational grid in x-direction
! YCGRID input Coordinates of computational grid in y-direction
!
INTEGER KGRPNT(MXC,MYC)
REAL XCGRID(MXC,MYC), YCGRID(MXC,MYC)
!
! 5. Parameter variables
!
! 6. Local variables
!
! DISTA distance between (X1,Y1) and (X2A,Y2A)
! DISTB distance between (X1,Y1) and (X2 ,Y2 )
! DXA, DYA difference X1 - X2A respectively Y1 - Y2A
! DXB, DYB difference X2 - X1 respectively Y2 - Y1
! DXO, DYO difference X4 - X3 respectively Y4 - Y3
! DXYO distance between (X3,Y3) and (X4,Y4)
! EPS multiplication factor
! ICC index
! ICGRD index
! IENT number of entries of this subroutine
! ILINK indicates which link is analyzed: 1 -> neighbour in x
! 2 -> neighbour in y
! IX counter of gridpoints in x-direction
! IY counter of gridpoints in y-direction
! JJ counter for number of obstacles
! JP counter for number of corner points of obstacles
! MOVED boolean which tells whether OBSTACLE has been moved
! NUMCOR number of corner points of obstacle
! X1, Y1 computational grid coordinates of one end of grid link
! X2, Y2 computational grid coordinates of other end of grid link
! i.e. neighbouring point of X1,Y1 associated with linknumber
! X2A,Y2A other neighbouring point of X1,Y1 associated with other
! linknumber, or if invalid: = X2,Y2
! X3, Y3 user coordinates of one end of obstacle side
! X4, Y4 user coordinates of other end of obstacle side
! XCONOBST boolean variable which tells whether XC in on OBSTACLE
! line piece
! XYEPS displacement factor relative to local computational grid
!
INTEGER ICC, ICGRD, IENT, ILINK, IX, IY, JJ, JP
INTEGER NUMCOR
REAL DISTA, DISTB, DXA, DXB, DYA, DYB,
& DXO, DXYO, DYO, EPS, X1, X2, X2A, X3, X4,
& XYEPS, Y1, Y2, Y2A, Y3, Y4
LOGICAL MOVED, XCONOBST
TYPE(OBSTDAT), POINTER :: COBST 40.31
!
! 8. Subroutines used
!
! STRACE
! TCROSS
! MSGERR
!
LOGICAL TCROSS
!
! 9. Subroutines calling
!
! SWPREP
!
! 10. Error messages
!
! 11. Remarks
!
! 12. Structure
! ----------------------------------------------------------------
! Read number of obstacles from array OBSTA
! For every obstacle do
! Read number of corners of the obstacle
! For every corner of the obstacle do
! For every grid point do
! call function TCROSS to search if there is crossing point
! between the line of two points of the stencil and the
! line of the corners of the obstacle.
! If there is crossing point then
! If there computational gridpoint is on the obstacle
! then move corner points of obstacle perpendicular
! to obstacle with factor
! (XYEPS*DYO/DXYO,-XYEPS*DXO/DXYO)
! Moved = .True.
! If Moved then call OBSTMOVE again to check whether
! there are still computational grid points on OBSTACLE
! ----------------------------------------------------------------
!
! 13. Source text
! ======================================================================
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'OBSTMOVE')
MOVED = .FALSE.
EPS =1.E-2
!
IF (NUMOBS .GT. 0) THEN
COBST => FOBSTAC 40.31
DO JJ = 1, NUMOBS
! number of corner points of the obstacle
NUMCOR = COBST%NCRPTS
IF (ITEST.GE. 120) THEN
WRITE(PRINTF,50) JJ, NUMCOR
50 FORMAT( ' Obstacle number : ', I4,' has ',I4,' corners')
ENDIF
! *** X1 X2 X3 ETC. are the coordinates of point according ***
! *** with the scheme in the subroutine TCROSS header ***
X3 = COBST%XCRP(1) 40.31
Y3 = COBST%YCRP(1) 40.31
DO JP = 2, NUMCOR
X4 = COBST%XCRP(JP) 40.31
Y4 = COBST%YCRP(JP) 40.31
IF (ITEST.GE. 120) WRITE(PRINTF,30) JP,X4,Y4
30 FORMAT(' Corner number:', I4,' XP: ',E10.4,' YP: ',E11.4)
DO IX = 1, MXC
DO IY = 1, MYC
ICC = KGRPNT(IX,IY)
IF (ICC .GT. 1) THEN
X1 = XCGRID(IX,IY)
Y1 = YCGRID(IX,IY)
!
! *** "ILINK" indicates which link is analyzed. Initial ***
! *** neighbour in x , second link with neighbouring in y***
DO ILINK = 1, 2
IF (ILINK.EQ.1 .AND. IX.GT.1) THEN
X2 = XCGRID(IX-1,IY)
Y2 = YCGRID(IX-1,IY)
IF (IY.GT.1) THEN
X2A = XCGRID(IX,IY-1)
Y2A = YCGRID(IX,IY-1)
ELSE
X2A = X2
Y2A = Y2
ENDIF
ICGRD = KGRPNT(IX-1,IY)
ELSE IF (ILINK.EQ.2 .AND. IY.GT.1) THEN
X2 = XCGRID(IX,IY-1)
Y2 = YCGRID(IX,IY-1)
IF (IX.GT.1) THEN
X2A = XCGRID(IX-1,IY)
Y2A = YCGRID(IX-1,IY)
ELSE
X2A = X2
Y2A = Y2
ENDIF
ICGRD = KGRPNT(IX,IY-1)
ELSE
ICGRD = 0
ENDIF
IF (ICGRD.GT.1) THEN
!
! *** All links are analyzed in each point otherwise the ***
! *** boundaries can be excluded ***
!
IF (TCROSS(X1,X2,X3,X4,Y1,Y2,Y3,Y4,XCONOBST))THEN
IF (XCONOBST) THEN
DXA=(X1-X2A)
DYA=(Y1-Y2A)
DXB=(X2-X1)
DYB=(Y2-Y1)
DISTA=SQRT(DXA*DXA+DYA*DYA)
DISTB=SQRT(DXB*DXB+DYB*DYB)
XYEPS = EPS*MIN(DISTA,DISTB)
DXO = X4-X3
DYO = Y4-Y3
DXYO = SQRT(DXO*DXO+DYO*DYO)
! -DXO/DXYO and DYO/DXYO are used (instead of -DXO
! and DYO) because otherwise displacement is dependent
! on length of OBSTACLE line piece
COBST%XCRP(JP-1) = X3 + XYEPS*DYO/DXYO 40.31
COBST%YCRP(JP-1) = Y3 - XYEPS*DXO/DXYO 40.31
COBST%XCRP(JP ) = X4 + XYEPS*DYO/DXYO 40.31
COBST%YCRP(JP ) = Y4 - XYEPS*DXO/DXYO 40.31
CALL MSGERR (1, 'Obstacle points moved')
WRITE(PRINTF, 17) X3+XOFFS, Y3+YOFFS, 40.81
& X4+XOFFS, Y4+YOFFS, 40.81
& X3+XYEPS*DYO/DXYO+XOFFS, 40.81
& Y3-XYEPS*DXO/DXYO+YOFFS, 40.81
& X4+XYEPS*DYO/DXYO+XOFFS, 40.81
& Y4-XYEPS*DXO/DXYO+YOFFS, 40.81
& X1+XOFFS,Y1+YOFFS 40.81
17 FORMAT ('OBSTACLE POINTS (', F11.2, ',', F11.2,
& '), and (', F11.2,',', F11.2,'),',
& 'moved to: (', F11.2,',',
& F11.2,'), and (', F11.2,',', F11.2,
& '), because OBSTACLE line piece ',
& 'was on computational grid point (',
& F11.2,',', F11.2,').')
X3 = X3 + XYEPS * DYO/DXYO
Y3 = Y3 - XYEPS * DXO/DXYO
X4 = X4 + XYEPS * DYO/DXYO
Y4 = Y4 - XYEPS * DXO/DXYO
MOVED = .TRUE.
ENDIF
ENDIF
ENDIF
END DO
ENDIF
END DO
END DO
X3 = X4
Y3 = Y4
END DO
IF (.NOT.ASSOCIATED(COBST%NEXTOBST)) EXIT
COBST => COBST%NEXTOBST
END DO
ENDIF
IF (MOVED) CALL OBSTMOVE(XCGRID, YCGRID, KGRPNT) 40.31
RETURN
! * end of subroutine OBSTMOVE *
END
!
!***********************************************************************
! *
SUBROUTINE SWTRCF (WLEV2 , CHS , 40.31 40.00
& LINK , OBREDF, 40.03
& AC2 , REFLSO, KGRPNT, XCGRID, 40.41 40.09
& YCGRID, CAX, CAY, RDX, RDY, ANYBIN, 40.09
& SPCSIG, SPCDIR) 40.13 40.28
! *
!***********************************************************************
USE OCPCOMM4 40.41
USE SWCOMM2 40.66
USE SWCOMM3 40.41
USE SWCOMM4 40.41
USE M_OBSTA 40.31
USE M_PARALL 40.31
USE SwanGriddata 40.80
USE SwanCompdata 40.80
IMPLICIT NONE 40.09
!
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 30.70
! 40.03 Nico Booij
! 40.08 Erick Rogers
! 40.09 Annette Kieftenburg
! 40.13 Nico Booij
! 40.14 Annette Kieftenburg
! 40.18 Annette Kieftenburg
! 40.28 Annette Kieftenburg
! 40.30 Marcel Zijlema
! 40.31 Marcel Zijlema
! 40.41: Marcel Zijlema
! 40.66: Marcel Zijlema
! 40.80: Marcel Zijlema
! 41.65: Marcel Zijlema
!
! 1. Updates
!
! 30.70, Feb. 98: water level (WLEV2) replaced depth
! incident wave height introduced using argument
! CHS (sign. wave height in whole comput. grid)
! 40.03, Jul. 00: LINK1 and LINK2 in argumentlist replaced by LINK
! 40.09, Nov. 99: IMPLICIT NONE added, Method corrected
! Reflection option for obstacle added
! 40.14, Dec. 00: Reflection call corrected: reduced to neighbouring
! linepiece of obstacle (bug fix 40.11D)
! Jan. 01: Constant waterlevel taken into account as well (bug fix 40.11E)
! 40.18, Apr. 01: Scattered reflection against obstacles added 40.18
! 40.28, Dec. 01: Frequency dependent reflection added
! 40.13, Aug. 02: subroutine restructured:
! loop in reflection procedure changed to avoid double
! reflection
! argument list of subr REFLECT revised
! argument SPCDIR added
! 40.30, Mar. 03: correcting indices of test point with offsets MXF, MYF
! 40.08, Mar. 03: Dimensioning of RDX, RDX changed to be consistent
! with other subroutines
! 40.31, Oct. 03: changes w.r.t. obstacles
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
! 40.66, Mar. 07: extension with d'Angremond and Van der Meer transmission
! 40.80, Mar. 08: extension to unstructured grids
! 41.65, Jun. 16: extension frequency and direction dependent tranmission coefficients
!
! 2. Purpose
!
! take the value of transmission coefficient given
! by the user in case obstacle TRANSMISSION
!
! or
!
! compute the transmision coeficient in case obstacle DAM
! based on Goda (1967) [from Seelig (1979)]
! or d'Angremond and Van der Meer formula's (1996) 40.66
!
! if reflections are switched on, calculate sourceterm in 40.09
! subroutine REFLECT 40.09
!
! 3. Method
!
! Calculate transmission coefficient based on Goda (1967) 40.09
! from Seelig (1979) 40.09
! Kt = 0.5*(1-sin {pi/(2*alpha)*(WATHIG/Hi +beta)})
! where
! Kt transmission coefficient
!
! alpha,beta coefficients dependent on structure of obstacle
! and waves
! WATHIG = F = h-d is the freeboard of the dam, where h is the 40.09
! crest level of the dam above the reference level and d 40.09
! is the mean water level (relative to reference level) 40.09
! Hi incident (significant) wave height 40.09
! 40.09
! If reflection are switched on and obstacle is not exactly on line 40.09
! of two neighbouring gridpoints, calculate reflections 40.09
!
! 4. Argument variables
!
! AC2 input Action density array 40.09
! ANYBIN input Set a particular bin TRUE or FALSE depending on 40.09
! SECTOR 40.09
! CAX input Propagation velocity 40.09
! CAY input Propagation velocity 40.09
! CHS input Hs in all computational grid points
! KCGRD input Grid address of points of computational stencil
! LINK input indicates whether link in stencil 40.03
! crosses an obstacle 40.03
! OBREDF output Array of action density reduction coefficients
! (reduction at the obstacle)
! REFLSO inp/outp contribution to the source term of action 40.41
! balance equation due to reflection 40.41
! RDX,RDY input Array containing spatial derivative coefficients 40.09
! WLEV2 input Water level in grid points
!
INTEGER KGRPNT(MXC,MYC)
INTEGER LINK(2) 40.03
REAL CHS(MCGRD), OBREDF(MDC,MSC,2), WLEV2(MCGRD) 30.70
REAL :: AC2(MDC,MSC,MCGRD) 40.09 40.22
! Changed ICMAX to MICMAX, since MICMAX doesn't vary over gridpoint 40.22
REAL :: CAX(MDC,MSC,MICMAX), CAY(MDC,MSC,MICMAX) 40.09 40.22
REAL :: REFLSO(MDC,MSC), RDX(MICMAX), RDY(MICMAX) 40.41 40.09 40.22 40.08
REAL :: SPCSIG(MSC), SPCDIR(MDC,6) 40.13 40.28
LOGICAL :: ANYBIN(MDC,MSC) 40.09
!
! 5. Parameter variables
!
! 6. Local variables
!
! ALOW Lower limit for FVH
! BK crest width 40.66
! BUPL Upper limit for FVH
! BVH Bk/Hsin 40.66
! FD1 Coeff. for freq. dep. reflection: vertical displacement 40.28
! FD2 Coeff. for freq. dep. reflection: shape parameter 40.28
! FD3 Coeff. for freq. dep. reflection: directional coefficient 40.28
! FD4 Coeff. for freq. dep. reflection: bending point of freq. 40.28
! FVH WATHIG/Hsin
! HGT elevation of top of obstacle above reference level
! HSIN incoming significant wave height
! ID counter in directional space
! IENT number of entries of this subroutine
! ILINK indicates which link is analyzed: 1 -> neighbour in x
! 2 -> neighbour in y
! IS counter in frequency space
! ITRAS indicates kind of obstacle: 0 -> constant transm
! 1 -> dam, Goda
! 2 -> dam, d'Angremond and 40.66
! Van der Meer 40.66
! JP counter for number of corner points of obstacles
! L0P wave length in deep water 40.66
! LREFDIFF indicates whether reflected energy should be 40.18
! scattered (1) or not (0) 40.18
! LREFL if LREFL=0, no reflection; if LREFL=1, constant 40.13
! reflection coeff. 40.13
! LRFRD Indicates whether frequency dependent reflection is 40.28
! active (#0.) or not (=0.) 40.28
! NMPO link number
! NUMCOR number of corner points of obstacle
! OBET user defined coefficient (beta) in formulation of
! Goda/Seelig (1967/1979)
! OBHKT transmission coefficient in terms of wave height
! OGAM user defined coefficient (alpha) in formulation of
! Goda/Seelig (1967/1979)
! POWN user defined power of redistribution function 40.28
! REFLCOEF reflection coefficient in terms of action density
! REFLTST used to test Refl^2+Transm^2 <=1 40.13
! SLOPE slope of obstacle 40.66
! SQRTREF dummy variable
! TRCF transmission coefficient in terms of action density
! (user defined or calculated (in terms of waveheight))
! X1, Y1 user coordinates of one end of grid link
! X2, Y2 user coordinates of other end of grid link
! X3, Y3 user coordinates of one end of obstacle side
! X4, Y4 user coordinates of other end of obstacle side
! XCGRID Coordinates of computational grid in x-direction
! XI0P breaker parameter 40.66
! XOBS x-coordinate of obstacle point 40.80
! XONOBST Indicates whether computational point (X1,Y1) is on 40.14
! obstacle 40.14
! XV x-coordinate of vertex of face 40.80
! YCGRID Coordinates of computational grid in y-direction
! YOBS y-coordinate of obstacle point 40.80
! YV y-coordinate of vertex of face 40.80
! WATHIG freeboard of the dam (= HGT-waterlevel)
!
INTEGER ID, IENT, ILINK, ITRAS, IS, JP, ICGRD, LREFL,
& NUMCOR, NMPO, ISIGM
INTEGER LREFDIFF, LRFRD 40.31
REAL ALOW, BUPL, FVH, HGT, HSIN, OBET, OBHKT,
& SLOPE, BK, L0P, XI0P, BVH, EMAX, ETD, TP, 40.66
& FAC1, FAC2,
& POWN, OGAM, REFLCOEF, 40.18 40.09
& X1, X2, X3, X4, Y1, Y2, Y3, Y4, WATHIG
REAL TRCF(MSC,MDC) 41.65
REAL FD1, FD2, FD3, FD4 40.31 40.28
REAL SQRTREF 40.09
LOGICAL XONOBST 40.14
LOGICAL :: REFLTST 40.13
REAL XCGRID(MXC,MYC), YCGRID(MXC,MYC)
INTEGER ICC, JJ
REAL :: XOBS(2), XV(2), YOBS(2), YV(2) 40.80
LOGICAL :: SwanCrossObstacle 40.80
TYPE(OBSTDAT), POINTER :: COBST 40.31
!
! 8. Subroutines used
!
! MSGERR Writes error message
! REFLECT Computes effect of reflection
! TCROSS Searches for crossing point if exist 40.14
!
LOGICAL TCROSS 40.14
!
! 9. Subroutines calling
!
! SWOMPU 30.70
!
! 10. Error messages
!
! 11. Remarks
!
! Here the formulation of the transmission coefficients concerns the 40.09
! ratio of action densities! 40.09
!
! 12. Structure
!
! ------------------------------------------------------------------
! For both links from grid point (X1,Y1) do 40.13
! calculate transmission coefficients 40.09
! assign values to OBREDF 40.13
! If there is reflection 40.13
! Then select obstacle side which crosses the grid link 40.13
! calculate reflection source terms 40.09
! ------------------------------------------------------------------ 40.13
!
! 13. Source text
! ======================================================================
!
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'SWTRCF')
!
REFLTST = .TRUE. 40.13
DO 90 ILINK = 1 ,2
! default transmission coefficient
TRCF = 1.
NMPO = LINK(ILINK)
IF (NMPO .EQ. 0) THEN
OBREDF(1:MDC,1:MSC,ILINK) = 1. 40.13
GOTO 90
END IF
! incoming wave height
HSIN = CHS(KCGRD(ILINK+1)) 40.03
IF (HSIN.LT.0.1E-4) HSIN = 0.1E-4
COBST => FOBSTAC 40.31
DO JJ = 1, NMPO-1 40.31
IF (.NOT.ASSOCIATED(COBST%NEXTOBST)) EXIT 40.31
COBST => COBST%NEXTOBST 40.31
END DO 40.31
ITRAS = COBST%TRTYPE 40.31
IF (ITRAS .EQ. 0) THEN
! constant transmission coefficient
! User defined transmission coefficient concerns ratio of 40.09
! wave heights, so 40.09
OBHKT = COBST%TRCOEF(1) 40.31
TRCF(1,1) = OBHKT * OBHKT 40.09
ELSE IF (ITRAS .EQ. 1) THEN
! transmission coefficient according to Goda and Seelig
HGT = COBST%TRCOEF(1) 40.31
OGAM = COBST%TRCOEF(2) 40.31
OBET = COBST%TRCOEF(3) 40.31
! level of dam above the water surface (freeboard):
WATHIG = HGT - WLEV2(KCGRD(1)) - WLEV 40.14 30.70
!
! *** Here the transmission coeff. is that of Goda and Seelig ***
FVH = WATHIG/HSIN
ALOW = -OBET-OGAM 40.09
BUPL = OGAM-OBET
!
IF (FVH.LT.ALOW) FVH = ALOW
IF (FVH.GT.BUPL) FVH = BUPL
OBHKT = 0.5*(1.0-SIN(PI*(FVH+OBET)/(2.0*OGAM)))
IF (TESTFL) WRITE (PRTEST, 20) IXCGRD(1)+MXF-2, 40.30
& IYCGRD(1)+MYF-2, 40.30
& ILINK, HGT, WATHIG, HSIN, OBHKT 40.01
20 FORMAT (' test SWTRCF ', 2X, 3I5, ' dam level=', F6.2,
& ' depth=', F6.2, ' Hs=', F6.2, ' transm=', F6.3) 40.01
IF (TESTFL .AND. ITEST.GE.140) WRITE (PRTEST, 22)
& OGAM, OBET, ALOW, BUPL, FVH 40.01
22 FORMAT (8X, 5E12.4) 40.01
!
! Formulation of Goda/Seelig concerns ratio of waveheights. 40.09
! Here we use action density, so 40.09
TRCF(1,1) = OBHKT * OBHKT
ELSE IF (ITRAS.EQ.2) THEN 40.66
! d'Angremond and Van der Meer formulae (1996) 40.66
HGT = COBST%TRCOEF(1)
SLOPE = COBST%TRCOEF(2)
BK = COBST%TRCOEF(3)
!
! level of dam above the water surface:
WATHIG = HGT - WLEV2(KCGRD(1)) - WLEV
!
! compute peak frequency of incoming wave
EMAX = 0.
ISIGM = -1
DO IS = 1, MSC
ETD = 0.
DO ID = 1, MDC
ETD = ETD + SPCSIG(IS)*AC2(ID,IS,KCGRD(ILINK+1))*DDIR
END DO
IF (ETD.GT.EMAX) THEN
EMAX = ETD
ISIGM = IS
END IF
END DO
IF (ISIGM.LE.0) ISIGM=MSC
TP=2.*PI/SPCSIG(ISIGM)
!
! compute breaker parameter
L0P = MAX(1.E-8,1.5613*TP*TP)
XI0P = TAN(SLOPE*PI/180.)/SQRT(HSIN/L0P)
!
! compute transmission coefficient
FVH = WATHIG/HSIN
BVH = BK/HSIN
IF (BVH.EQ.0.) THEN
OBHKT= -0.40*FVH
IF (OBHKT.LT.0.075) OBHKT = 0.075
IF (OBHKT.GT.0.900) OBHKT = 0.9
ELSE IF (BVH.LT.8.) THEN
OBHKT= -0.40*FVH + 0.64*(BVH**(-0.31))*(1.-EXP(-0.50*XI0P))
IF (OBHKT.LT.0.075) OBHKT = 0.075
IF (OBHKT.GT.0.900) OBHKT = 0.9
ELSE IF (BVH.GT.12.) THEN
OBHKT= -0.35*FVH + 0.51*(BVH**(-0.65))*(1.-EXP(-0.41*XI0P))
IF (OBHKT.GT.0.93-0.006*BVH) OBHKT = 0.93-0.006*BVH
IF (OBHKT.LT.0.05 ) OBHKT = 0.05
ELSE
! linear interpolation
FAC1 = -0.40*FVH + 0.64*( 8.**(-0.31))*(1.-EXP(-0.50*XI0P))
IF (FAC1.LT.0.075) FAC1 = 0.075
IF (FAC1.GT.0.900) FAC1 = 0.9
FAC2 = -0.35*FVH + 0.51*(12.**(-0.65))*(1.-EXP(-0.41*XI0P))
IF (FAC2.LT.0.050) FAC2 = 0.050
IF (FAC2.GT.0.858) FAC2 = 0.858
OBHKT = 3.*FAC1 - 2.*FAC2 + BVH*(FAC2 - FAC1)/4.
END IF
IF(OPTG.NE.5) THEN
X1 = XCGRID(IXCGRD(1),IYCGRD(1))+XOFFS
Y1 = YCGRID(IXCGRD(1),IYCGRD(1))+YOFFS
ELSE
X1 = xcugrd(vs(1))+XOFFS
Y1 = ycugrd(vs(1))+YOFFS
ENDIF
WRITE (PRINTF, 23) X1, Y1,
& ILINK, HGT, WATHIG, HSIN, SQRT(L0P/1.5613), XI0P, OBHKT 40.66
23 FORMAT (' Transmission: ', 2X, 2F12.4, I5, ' dam level=',F6.2,
& ' board=', F6.2, ' Hs=', F6.2, ' Tp=', F6.2,
& ' Xi0p=', F6.3, ' Kt=', F6.3) 40.66
IF (TESTFL .AND. ITEST.GE.140) WRITE (PRTEST, 22)
& TAN(SLOPE*PI/180.), FVH, BVH, L0P, XI0P
!
! Formulation of d'Angremond concerns ratio of waveheights
! Here we use action density, so
TRCF(1,1) = OBHKT * OBHKT
ELSE IF (ITRAS .EQ. 11) THEN 41.65
! frequency dependent transmission coefficients 41.65
! user defined values concerns ratio of wave heights, so 41.65
DO IS = 1, MSC
OBHKT = COBST%TRCF1D(IS) 41.65
TRCF(IS,1) = OBHKT * OBHKT 41.65
ENDDO
ELSE IF (ITRAS .EQ. 12) THEN 41.65
! frequency and direction dependent transmission coefficients 41.65
! user defined values concerns ratio of wave heights, so 41.65
DO ID = 1, MDC
DO IS = 1, MSC
OBHKT = COBST%TRCF2D(ID,IS) 41.65
TRCF(IS,ID) = OBHKT * OBHKT 41.65
ENDDO
ENDDO
ENDIF
!
! assign values to array OBREDF
IF (ITRAS.EQ.11) THEN 41.65
DO IS = 1, MSC
DO ID = 1, MDC
OBREDF(ID,IS,ILINK) = TRCF(IS,1)
ENDDO
ENDDO
ELSE IF (ITRAS.EQ.12) THEN
DO IS = 1, MSC
DO ID = 1, MDC
OBREDF(ID,IS,ILINK) = TRCF(IS,ID)
ENDDO
ENDDO
ELSE
DO IS = 1, MSC
DO ID = 1, MDC
OBREDF(ID,IS,ILINK) = TRCF(1,1)
ENDDO
ENDDO
ENDIF
!
! *** REFLECTION ****
! *** X1 X2 X3 ETC. are the coordinates of point according ***
! *** with the scheme in the function TCROSS header *** 40.04
LREFL = COBST%RFTYP1 40.13
IF ( LREFL.GT.0. ) THEN 40.31
! Reflections are activated 40.09
SQRTREF = COBST%RFCOEF(1) 40.31
REFLCOEF = SQRTREF * SQRTREF 40.09
LREFDIFF = COBST%RFTYP2 40.31
POWN = COBST%RFCOEF(2) 40.31
FD1 = COBST%RFCOEF(3) 40.31
FD2 = COBST%RFCOEF(4) 40.31
FD3 = COBST%RFCOEF(5) 40.31
FD4 = COBST%RFCOEF(6) 40.31
LRFRD = COBST%RFTYP3 40.31
!
IF (OPTG.NE.5) THEN 40.80
! determine (X1,Y1) and (X2,Y2)
ICC = KCGRD(1) 40.09
ICGRD = 0 40.09
IF ( ICC.GT.1 ) THEN 40.09
X1 = XCGRID(IXCGRD(1),IYCGRD(1)) 40.09
Y1 = YCGRID(IXCGRD(1),IYCGRD(1)) 40.09
IF (KGRPNT(IXCGRD(ILINK+1),IYCGRD(ILINK+1)).GT.1) THEN 40.09
X2 = XCGRID(IXCGRD(ILINK+1),IYCGRD(ILINK+1)) 40.09
Y2 = YCGRID(IXCGRD(ILINK+1),IYCGRD(ILINK+1)) 40.09
ICGRD = KCGRD(ILINK+1) 40.09
ENDIF
ENDIF
IF (ICGRD.EQ.0) GOTO 90 40.13
! select obstacle side crossing the grid link
X3 = COBST%XCRP(1) 40.31
Y3 = COBST%YCRP(1) 40.31
NUMCOR = COBST%NCRPTS 40.31
DO JP = 2, NUMCOR 40.09
X4 = COBST%XCRP(JP) 40.31
Y4 = COBST%YCRP(JP) 40.31
IF (TCROSS(X1,X2,X3,X4,Y1,Y2,Y3,Y4,XONOBST)) GOTO 70 40.14
X3 = X4 40.09
Y3 = Y4 40.09
ENDDO
ELSE 40.80
! determine begin and end points of link 40.80
X1 = xcugrd(vs(1)) 40.80
Y1 = ycugrd(vs(1)) 40.80
X2 = xcugrd(vs(ILINK+1)) 40.80
Y2 = ycugrd(vs(ILINK+1)) 40.80
XV(1) = X1 40.80
YV(1) = Y1 40.80
XV(2) = X2 40.80
YV(2) = Y2 40.80
! select obstacle side crossing the grid link 40.80
X3 = COBST%XCRP(1) 40.80
Y3 = COBST%YCRP(1) 40.80
XOBS(1) = X3 40.80
YOBS(1) = Y3 40.80
DO JP = 2, COBST%NCRPTS 40.80
X4 = COBST%XCRP(JP) 40.80
Y4 = COBST%YCRP(JP) 40.80
XOBS(2) = X4 40.80
YOBS(2) = Y4 40.80
IF ( SwanCrossObstacle( XV, YV, XOBS, YOBS ) ) GOTO 70 40.80
X3 = X4 40.80
Y3 = Y4 40.80
XOBS(1) = X3 40.80
YOBS(1) = Y3 40.80
ENDDO 40.80
ENDIF 40.80
! no crossing found, skip procedure
GOTO 90
70 CALL REFLECT(AC2, REFLSO, X1, Y1, X2, Y2, 40.13
& X3, Y3, X4, Y4, CAX, 40.13
& CAY, RDX, RDY, ILINK, 40.13
& REFLCOEF, LREFDIFF, POWN, ANYBIN, 40.13
& LRFRD, SPCSIG, SPCDIR, FD1, FD2, FD3, FD4, 40.13
& OBREDF, REFLTST) 40.13
END IF
!
IF (ITEST .GE. 120) WRITE (PRTEST,10)
& IXCGRD(1)-1, IYCGRD(1)-1, NMPO, TRCF(1,1)
10 FORMAT(' SWTRCF: Point=', 2I5, ' NMPO = ', I5, ' transm ', 40.03
& F8.3) 40.03
90 CONTINUE
IF (.NOT.REFLTST) THEN 40.13
CALL MSGERR(3,'Kt^2 + Kr^2 > 1 ') 40.13
IF (ITEST.LT.50) THEN 40.13
WRITE (PRTEST, 74) IXCGRD(1)-1, IYCGRD(1)-1 40.13
74 FORMAT (' Kt^2 + Kr^2 > 1 in grid point:', 2I4) 40.13
ENDIF 40.13
ENDIF 40.13
RETURN
! * end of SUBROUTINE SWTRCF
END
!
!************************************************************************
! *
SUBROUTINE REFLECT (AC2, REFLSO, X1, Y1, X2, Y2, X3, Y3, 40.13
& X4, Y4, CAX, CAY, RDX, RDY, 40.13
& ILINK, REF0, LREFDIFF, POWN, ANYBIN, 40.13
& LRFRD, SPCSIG, SPCDIR, FD1, FD2, FD3, FD4, 40.13
& OBREDF, REFLTST) 40.13
! *
!************************************************************************
! 40.09
USE OCPCOMM4 40.41
USE SWCOMM3 40.41
USE SWCOMM4 40.41
! 40.09
IMPLICIT NONE 40.09
!
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors 40.09
! 40.09
! 40.09 Annette Kieftenburg 40.09
! 40.13 Nico Booij
! 40.18 Annette Kieftenburg 40.18
! 40.28 Annette Kieftenburg 40.28
! 40.38 Annette Kieftenburg 40.38
! 40.41: Marcel Zijlema
! 40.09
! 1. Updates 40.09
! 40.09
! 40.09, Nov. 99: Subroutine created 40.09
! 40.18, Apr. 01: Scattered reflection against obstacles added 40.18
! 40.28, Dec. 01: Frequency dependent reflection added 40.28
! 40.38, Feb. 02: Diffuse reflection against obstacles added 40.38
! 40.13, Sep. 02: Subroutine restructured
! assumptions changed: reflected energy can come
! from Th_norm-PI/2 to Th_norm+PI/2
! 40.08, Mar. 03: Dimensioning of RDX, RDX changed to be consistent 40.08
! with other subroutines 40.08
! 40.13, Nov. 03: test on refl + transm added
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
! 40.09
! 2. Purpose 40.09
! 40.09
! Computation of REFLECTIONS near obstacles 40.09
! 40.09
! 3. Method 40.09
! 40.09
! Determine the angle of the obstacle, 40.09
! Determine the angles between which reflections should be taken 40.09
! into account 40.09
! Determine redistribution function 40.18
! determine expression of reflection coefficient for frequency
! dependency, if appropriate 40.28
! Determine reflected action density (corrected for angle obstacle 40.09
! and if option is on: redistribute energy) 40.18
! Add reflected spectrum to contribution for the right hand side 40.41 40.09
! of matrix equation 40.09
! 40.09
! 4. Modules used 40.18
! 40.18
! -- 40.18
! 40.18
! 5. Argument variables 40.09
! 40.09
! AC2 inp action density
! ANYBIN inp Determines whether a bin fall within a sweep
! CAX inp Propagation velocity in x-direction 40.09
! CAY inp Propagation velocity in y-direction 40.09
! FD1 inp Coeff. freq. dep. reflection: vertical displacement 40.28
! FD2 inp Coeff. freq. dep. reflection: shape parameter 40.28
! FD3 inp Coeff. freq. dep. reflection: directional coefficient 40.28
! FD4 inp Coeff. freq. dep. reflection: bending point of freq. 40.28
! ILINK inp Indicates which link is analyzed: 1 -> neighbour in x 40.09
! 2 -> neighbour in y 40.09
! LREFDIFF inp Indicates whether reflected energy should be 40.18
! scattered (1) or not (0) 40.18
! LRFRD inp Indicates whether frequency dependent reflection is 40.28
! active (#0.) or not (=0.) 40.28
! OBREDF inp transmission coefficients
! POWN inp User defined power of redistribution function 40.18
! REF0 inp reflection coefficient in terms of action density
! REFLSO i/o contribution to the source term due to reflection 40.41
! REFLTST i/o used to test Refl^2+Transm^2 <=1 40.13
! RDX,RDY inp Array containing spatial derivative coefficients 40.09
! SPCDIR(*,1) spectral directions (radians)
! SPCSIG inp Relative frequency (= 2*PI*Freq.) 40.28
! X1, Y1 inp Coordinates of computational grid point under 40.09
! consideration 40.09
! X2, Y2 inp Coordinates of computational grid point neighbour 40.09
! X3, Y3 inp User coordinates of one end of obstacle side 40.09
! X4, Y4 inp User coordinates of other end of obstacle side 40.09
! 40.09
REAL :: AC2(MDC,MSC,MCGRD) 40.18
! Changed ICMAX to MICMAX, since MICMAX doesn't vary over gridpoint 40.22
REAL :: CAX(MDC,MSC,MICMAX), CAY(MDC,MSC,MICMAX) 40.18 40.22
REAL :: REFLSO(MDC,MSC), OBREDF(MDC,MSC,2) 40.41 40.18
REAL :: RDX(MICMAX), RDY(MICMAX) 40.18 40.08
REAL :: FD1, FD2, FD3, FD4, SPCSIG(MSC), SPCDIR(MDC,6) 40.41 40.28
REAL :: REF0 40.18
REAL :: X1, X2, X3, X4, Y1, Y2, Y3, Y4 40.18
LOGICAL :: ANYBIN(MDC,MSC) 40.18
INTEGER :: ILINK 40.18
REAL :: POWN 40.41 40.18
INTEGER :: LREFDIFF, LRFRD 40.41
LOGICAL :: REFLTST 40.13
! 40.18
INTENT (IN) AC2, CAX, CAY, OBREDF, 40.18
& FD1, FD2, FD3, FD4, LRFRD, 40.28
& RDX, RDY, SPCSIG, SPCDIR, X1, X2, 40.18
& X3, X4, Y1, Y2, Y3, Y4, ANYBIN, ILINK, 40.18
& POWN, LREFDIFF 40.18
INTENT (IN OUT) REF0, REFLSO 40.41 40.18
! 40.09
! 6. Parameter variables 40.09
! 40.18
! 7. Local variables 40.09
! 40.09
REAL :: AC2REF ! reflected action density of one spectral bin
REAL :: BETA ! local angle of obstacle 40.09
REAL :: EPS 40.09
REAL, ALLOCATABLE :: PRDIF(:) ! scattering filter 40.13
REAL :: TH_INC ! direction of incident wave
REAL :: TH_NORM ! direction of normal to obstacle
REAL :: TH_OUT ! direction of outgoing wave
REAL :: IANG ! angle divided by DTheta (DDIR)
REAL :: SUMRD ! sum of PRDIF array
REAL :: W1, W2 ! interpolation coefficients
INTEGER :: ID ! counter of directions
INTEGER :: IS ! counter of frequencies
INTEGER :: MAXIDR ! width of scattering filter
INTEGER :: IDR ! relative directional counter
INTEGER :: ID_I1, ID_I2 ! counters of incoming directions
INTEGER :: IDA, IDB ! counters of incoming directions
INTEGER :: IENT ! number of entries
! 40.18
! 8. Subroutines used 40.09
! 40.09
! 9. Subroutines calling 40.09
! 40.09
! SWTRCF 40.09
! 40.09
! 10. Error messages 40.09
! 40.09
! if obstacle linepiece is of length < EPS 40.09
! 40.09
! 11. Remarks 40.09
! 40.09
! -In case the obstacle cuts exactly through computational grid point, 40.09
! the obstacle should be moved a bit with subroutine OBSTMOVE. 40.09
! -The length of the obstacle linepiece is assumed to be 40.09
! 'long enough' compared to grid resolution (> 0.5*sqrt(dx^2+dy^2)) 40.09
! (if this restriction is violated, the reflections due to an obsta- 40.09
! cle of one straight line can be very different from a similar line 40.09
! consisting of several pieces (because only the directions of the 40.09
! spectrum that are directed towards the obstacle linepiece are 40.09
! reflected). 40.09
! -There should be only one intersection per computational gridcell. 40.09
! Therefore it is better to avoid sharp edges in obstacles. 40.18
! 40.09
! 12. Structure 40.09
! 40.09
! -----------------------------------------------------------------
! Determine angle of obstacle, Beta 40.09
! Determine angle of normal from (X1,Y1) to obstacle 40.13
! If there is constant diffuse reflection
! Then determine scattering distribution
! -----------------------------------------------------------------
! For all frequencies do
! If amount of reflection varies with frequency
! Then determine reflection coefficient
! -------------------------------------------------------------
! If there is diffuse reflection
! Then if scattering varies with frequency
! Then determine power of cos
! --------------------------------------------------------
! Determine distribution
! -------------------------------------------------------------
! For all active directions do
! Determine specular incoming direction
! For directions of scattering filter do
! multiply incoming action with scattering coefficient
! add this to array AC2REF
! -----------------------------------------------------------------
! 40.09
! Add reflected spectrum to right hand side of matrix equation 40.09
! 40.09
! 13. Source text 40.09
! 40.18
SAVE IENT 40.09
DATA IENT /0/ 40.09
CALL STRACE (IENT, 'REFLECT') 40.09
! 40.09
EPS = EPSILON(X1)*SQRT((X2-X1)*(X2-X1)+(Y2-Y1)*(Y2-Y1)) 40.09
IF (EPS ==0) EPS = TINY(X1) 40.09
IF (LREFDIFF .EQ. 0) THEN
ALLOCATE (PRDIF(0:0))
MAXIDR = 0
PRDIF(0) = 1.
ELSE
MAXIDR = MDC/2
ALLOCATE (PRDIF(0:MDC/2))
ENDIF
! 40.09
! Determine angle of obstacle BETA, and related 40.09
! 40.09
IF (.NOT. ((ABS(Y4-Y3).LE.EPS) .AND. (ABS(X4-X3).LE.EPS)) ) THEN 40.09
BETA = ATAN2((Y4-Y3),(X4-X3)) 40.09
ELSE 40.09
CALL MSGERR (3, 'Obstacle contains line piece of length 0!') 40.09
END IF 40.09
! determine direction of normal (4)
! this is the normal from (X1,Y1) |
! towards the obstacle |
! see sketch to establish sign (2)----+------(1)
! |
! (3)
IF ((X1-X3)*(Y4-Y3)-(Y1-Y3)*(X4-X3) .GT. 0.) THEN 40.13
TH_NORM = BETA + 0.5*PI
ELSE
TH_NORM = BETA - 0.5*PI
ENDIF 40.13
! prepare directional filter in case of diffuse reflection
IF (LREFDIFF.EQ.1) THEN
PRDIF(1:MAXIDR) = 0.
PRDIF(0) = 1.
SUMRD = 1.
DO ID = 1, MAXIDR
PRDIF(ID) = (COS(ID*DDIR))**POWN
IF (PRDIF(ID) .GT. 0.01) THEN
SUMRD = SUMRD + 2.*PRDIF(ID)
ELSE
MAXIDR = ID-1
GOTO 23
ENDIF
ENDDO
23 DO ID = 0, MAXIDR
PRDIF(ID) = PRDIF(ID) / SUMRD
ENDDO
IF (TESTFL .AND. ITEST.GE.50) THEN
WRITE (PRTEST, 27) POWN, MAXIDR
27 FORMAT (' power scattering filter:', F4.1, I3)
IF (ITEST.GE.130) WRITE (PRTEST, 28)
& (PRDIF(IDR), IDR=0, MAXIDR)
28 FORMAT (10 F7.3)
ENDIF
ENDIF
DO IS = 1, MSC
IF (LRFRD.EQ.1) THEN
! amount of reflection varies with wave frequency
REF0 = FD1 + 40.13
& FD2/PI * ATAN2(PI*FD3*(SPCSIG(IS)-FD4),FD2) 40.13
IF (REF0 > 1.) REF0 = 1. 40.28
IF (REF0 < 0.) REF0 = 0. 40.28
! >>> should REF0 not be squared? <<<
ENDIF
! check whether reflection + transmission <= 1 40.13
DO ID = 1, MDC 40.13
IF ((REF0 + OBREDF(ID,IS,ILINK)) .GT. 1.) THEN 40.13
REFLTST = .FALSE. 40.13
IF (ITEST.GE.50) THEN 40.13
WRITE (PRTEST, 72) IXCGRD(1)-1, IYCGRD(1)-1, ILINK, 40.13
& IS, ID, REF0, OBREDF(ID,IS,ILINK) 40.13
72 FORMAT (' Refl+Transm>1 in ', 2I4, 2X, 3I3, 2X, 2F6.2) 40.13
ENDIF 40.13
ENDIF 40.13
ENDDO 40.13
IF (LREFDIFF.EQ.2) THEN
! spreading varies with frequency; not yet implemented
ENDIF
DO ID = 1, MDC
IF (ANYBIN(ID,IS)) THEN
AC2REF = 0.
TH_OUT = SPCDIR(ID,1)
! corresponding incident direction (assuming specular reflection)
TH_INC = 2.*BETA-TH_OUT
! determine counter for which direction is TH_INC:
IANG = MOD (TH_INC-SPCDIR(1,1), 2.*PI) / DDIR
! incident angle is between ID_I1 and ID_I2
ID_I1 = 1 + INT (IANG)
ID_I2 = ID_I1+1
! W1 and W2 are weighting coefficients for the above directions
! by linear interpolation
W2 = IANG + 1. - REAL(ID_I1)
W1 = 1. - W2
DO IDR = -MAXIDR, MAXIDR
IDA = ID_I1 + IDR
IDB = ID_I2 + IDR
IF (FULCIR) THEN
IDA = 1+MOD(2*MDC+IDA-1,MDC)
! only outgoing reflected waves, i.e. not towards obstacle
IF (COS(TH_NORM-SPCDIR(IDA,1)) .GT. 0.)
& AC2REF = AC2REF +
& REF0 * W1 * PRDIF(ABS(IDR)) * AC2(IDA,IS,KCGRD(1))
IDB = 1+MOD(2*MDC+IDB-1,MDC)
IF (COS(TH_NORM-SPCDIR(IDB,1)) .GT. 0.)
& AC2REF = AC2REF +
& REF0 * W2 * PRDIF(ABS(IDR)) * AC2(IDB,IS,KCGRD(1))
ELSE
IF (IDA.GE.1 .AND. IDA.LE.MDC) THEN
IF (COS(TH_NORM-SPCDIR(IDA,1)) .GT. 0.)
& AC2REF = AC2REF +
& REF0 * W1 * PRDIF(ABS(IDR)) * AC2(IDA,IS,KCGRD(1))
ENDIF
IF (IDB.GE.1 .AND. IDB.LE.MDC) THEN
IF (COS(TH_NORM-SPCDIR(IDB,1)) .GT. 0.)
& AC2REF = AC2REF +
& REF0 * W2 * PRDIF(ABS(IDR)) * AC2(IDB,IS,KCGRD(1))
ENDIF
ENDIF
ENDDO
! add reflected energy to right hand side of matrix
REFLSO(ID,IS) = REFLSO(ID,IS) + AC2REF *
& (RDX(ILINK)*CAX(ID,IS,1) + RDY(ILINK)*CAY(ID,IS,1)) 40.13
END IF 40.09
END DO 40.09
END DO 40.09
! 40.18
IF (LREFDIFF .GT. 0) DEALLOCATE (PRDIF)
RETURN 40.09
! End of subroutine REFLECT 40.09
END 40.09
!
!***********************************************************************
! *
SUBROUTINE SSHAPE (ACLOC, SPCSIG, SPCDIR, FSHAPL, DSHAPL)
! *
!***********************************************************************
!
USE OCPCOMM4 40.41
USE SWCOMM3 40.41
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! Roeland Ris
! Roberto Padilla
! 30.73: Nico Booij
! 30.80: Nico Booij
! 30.82: IJsbrand Haagsma
! 40.02: IJsbrand Haagsma
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! Dec. 92: new for SWAN
! Dec. 96: option MEAN freq. introduced see LOGPM
! 30.73, Nov. 97: revised in view of new boundary treatment
! 30.82, Sep. 98: Added error message in case of non-convergence
! 30.80, Oct. 98: correction suggested by Mauro Sclavo, and renames
! computation of tail added to improve accuracy
! 30.82, Oct. 98: Updated description of several variables
! 40.02, Oct. 00: Modified test write statement to avoid division by MS=0
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. Purpose
!
! Calculating of energy density at boundary point (x,y,sigma,theta)
!
! 3. Method (updated...)
!
! see: M. Yamaguchi: Approximate expressions for integral properties
! of the JONSWAP spectrum; Proc. JSCE, No. 345/II-1, pp. 149-152,
! 1984.
!
! computation of mean period: see Swan system documentation
!
! 4. Argument variables
!
! o ACLOC : Energy density at a point in space
! i SPCDIR: (*,1); spectral directions (radians) 30.82
! (*,2); cosine of spectral directions 30.82
! (*,3); sine of spectral directions 30.82
! (*,4); cosine^2 of spectral directions 30.82
! (*,5); cosine*sine of spectral directions 30.82
! (*,6); sine^2 of spectral directions 30.82
! i SPCSIG: Relative frequencies in computational domain in sigma-space 30.82
!
REAL ACLOC(MDC,MSC)
REAL SPCDIR(MDC,6) 30.82
REAL SPCSIG(MSC) 30.82
!
! i DSHAPL: Directional distribution
! i FSHAPL: Shape of spectrum:
! =1; Pierson-Moskowitz spectrum
! =2; Jonswap spectrum
! =3; bin
! =4; Gauss curve
! (if >0: period is interpreted as peak per.
! if <0: period is interpreted as mean per.)
!
INTEGER FSHAPL, DSHAPL 40.00
!
! 5. Parameter variables
!
! 6. Local variables
!
! ID counter of directions
! IS counter of frequencies
! LSHAPE absolute value of FSHAPL
!
INTEGER ID, IS, LSHAPE
!
! PKPER peak period 30.80
! APSHAP aux. var. used in computation of spectrum
! AUX1 auxiliary variable
! AUX2 auxiliary variable
! AUX3 auxiliary variable
! COEFF coefficient for behaviour around the peak (Jonswap)
! CPSHAP aux. var. used in computation of spectrum
! CTOT total energy
! CTOTT total energy (used for comparison)
! DIFPER auxiliary variable used to select bin closest
! to given frequency
! MPER
! MS power in directional distribution
! RA action density
! SALPHA
! SF frequency (Hz)
! SF4 SF**4
! SF5 SF**5
! FPK frequency corresponding to peak period (1/PKPER) 30.80
! FPK4 FPK**4
! SYF peakedness parameter
!
REAL APSHAP, AUX1, AUX2, AUX3
REAL COEFF ,SYF ,MPER ,CTOT ,CTOTT,PKPER ,DIFPER
REAL MS
REAL RA ,SALPHA,SF ,SF4 ,SF5 ,FPK ,FPK4, FAC
!
! LOGPM indicates whether peak or mean frequency is used
! DVERIF logical used in verification of incident direction
!
LOGICAL LOGPM, DVERIF 40.00
!
! PSHAPE coefficients of spectral distribution (see remarks)
! SPPARM array containing integral wave parameters (see remarks)
!
! 8. Subroutines used
!
! ---
!
! 9. Subroutines calling
!
! 10. Error messages
!
! 11. Remarks
!
! PSHAPE(1): SY0, peak enhancement factor (gamma) in Jonswap spectrum
! PSHAPE(2): spectral width in case of Gauss spectrum in rad/s
!
! SPPARM real input incident wave parameters (Hs, Period,
! direction, Ms (dir. spread))
! SPPARM(1): Hs, sign. wave height
! SPPARM(2): Wave period given by the user (either peak or mean) 30.80
! SPPARM(3): average direction
! SPPARM(4): directional spread
!
! ---------------------------------------------------------------------
!
! In the case of a JONSWAP spectrum the initial conditions are given by
! _ _ _ _ _
! | _ _ -4 | | | S - S |
! 2 | | S | | | | p |
! a g | | _ | | exp|-1/2 * |________ |* 2/pi COS(T-T )
! E(S,D )= ___ exp|-5/4 *| S | | G | | e * S | wi
! wa 5 | | p | | |_ |_ p _|
! S | |_ _| |
! |_ _|
!
! where
! S : rel. frequency
!
! D : Dir. of wave component
! wa
!
! a : equili. range const. (Phillips' constant)
! g : gravity acceleration
!
! S : Peak frequency
! p
!
! G : Peak enhancement factor
! e : Peak width
!
! T : local wind direction
! wi
!
! 12. Structure
!
! ----------------------------------------------------------------
! case shape
! =1: calculate value of Pierson-Moskowitz spectrum
! =2: calculate value of Jonswap spectrum
! =3: calculate value of bin spectrum
! =4: calculate value of Gauss spectrum
! else: Give error message because of wrong shape
! ----------------------------------------------------------------
! if LOGPM is True
! then calculate average period
! if it differs from given average period
! then recalculate peak period
! restart procedure to compute spectral shape
! ----------------------------------------------------------------
! for all spectral bins do
! multiply all action densities by directional distribution
! ----------------------------------------------------------------
!
! 13. Source text
!
SAVE IENT
DATA IENT/0/
CALL STRACE(IENT,'SSHAPE')
!
IF (ITEST.GE.80) WRITE (PRTEST, 8) FSHAPL, DSHAPL,
& (SPPARM(JJ), JJ = 1,4)
8 FORMAT (' entry SSHAPE ', 2I3, 4E12.4)
IF (FSHAPL.LT.0) THEN
LSHAPE = - FSHAPL
LOGPM = .FALSE.
ELSE
LSHAPE = FSHAPL
LOGPM = .TRUE.
ENDIF
!
IF (SPPARM(1).LE.0.) 40.31
& CALL MSGERR(1,'sign. wave height at boundary is not positive') 40.31
!
PKPER = SPPARM(2)
ITPER = 0
IF (LSHAPE.EQ.3) THEN
! select bin closest to given period
DIFPER = 1.E10
DO IS = 1, MSC
IF (ABS(PKPER - PI2/SPCSIG(IS)) .LT. DIFPER) THEN
ISP = IS
DIFPER = ABS(PKPER - PI2/SPCSIG(IS))
ENDIF
ENDDO
ENDIF
!
! compute spectral shape using peak period PKPER 30.80
!
FAC = 1.
100 FPK = (1./PKPER) 30.80
FPK4 = FPK**4
IF (LSHAPE.EQ.1) THEN
SALPHA = ((SPPARM(1) ** 2) * (FPK4)) * 5. / 16.
ELSE IF (LSHAPE.EQ.2) THEN
! *** SALPHA = alpha*(grav**2)/(2.*pi)**4)
SALPHA = (SPPARM(1)**2 * FPK4) /
& ((0.06533*(PSHAPE(1)**0.8015)+0.13467)*16.)
ELSE IF (LSHAPE.EQ.4) THEN
AUX1 = SPPARM(1)**2 / ( 16.* SQRT (PI2) * PSHAPE(2))
AUX3 = 2. * PSHAPE(2)**2
ENDIF
!
CTOTT = 0.
DO 300 IS = 1, MSC 30.80
!
IF (LSHAPE.EQ.1) THEN
! *** LSHAPE = 1 : Pierson and Moskowitz ***
SF = SPCSIG(IS) / PI2
SF4 = SF**4
SF5 = SF**5
RA = (SALPHA/SF5)*EXP(-(5.*FPK4)/(4.*SF4))/(PI2*SPCSIG(IS))
ACLOC(MDC,IS) = RA
ELSE IF (LSHAPE.EQ.2) THEN
! *** LSHAPE = 2 : JONSWAP ***
SF = SPCSIG(IS)/(PI2)
SF4 = SF**4
SF5 = SF**5
CPSHAP = 1.25 * FPK4 / SF4
IF (CPSHAP.GT.10.) THEN 30.50
RA = 0.
ELSE
RA = (SALPHA/SF5) * EXP(-CPSHAP)
ENDIF
IF (SF .LT. FPK) THEN
COEFF = 0.07
ELSE
COEFF = 0.09
ENDIF
APSHAP = 0.5 * ((SF-FPK) / (COEFF*FPK)) **2
IF (APSHAP.GT.10.) THEN 30.50
SYF = 1.
ELSE
PPSHAP = EXP(-APSHAP)
SYF = PSHAPE(1)**PPSHAP
ENDIF
RA = SYF*RA/(SPCSIG(IS)*PI2)
ACLOC(MDC,IS) = RA
IF (ITEST.GE.120) WRITE (PRTEST, 112)
& SF, SALPHA, CPSHAP, APSHAP, SYF, RA
112 FORMAT (' SSHAPE freq. ', 8E12.4)
ELSE IF (LSHAPE.EQ.3) THEN
!
! *** all energy concentrated in one BIN ***
!
IF (IS.EQ.ISP) THEN
ACLOC(MDC,IS) = ( SPPARM(1)**2 ) /
& ( 16. * SPCSIG(IS)**2 * FRINTF )
ELSE
ACLOC(MDC,IS) = 0.
ENDIF
ELSE IF (LSHAPE.EQ.4) THEN
!
! *** energy Gaussian distributed (wave-current tests) ***
!
AUX2 = ( SPCSIG(IS) - ( PI2 / PKPER ) )**2
RA = AUX1 * EXP ( -1. * AUX2 / AUX3 ) / SPCSIG(IS)
ACLOC(MDC,IS) = RA
ELSE
IF (IS.EQ.1) THEN
CALL MSGERR (2,'Wrong type for frequency shape')
WRITE (PRINTF, *) ' -> ', FSHAPL, LSHAPE
ENDIF
ENDIF
IF (ITEST.GE.10)
& CTOTT = CTOTT + FRINTF * ACLOC(MDC,IS) * SPCSIG(IS)**2
300 CONTINUE
IF (ITEST.GE.10) THEN
IF (SPPARM(1).GT.0.01) THEN
HSTMP = 4. * SQRT(CTOTT)
IF (ABS(HSTMP-SPPARM(1)) .GT. 0.1*SPPARM(1))
& WRITE (PRINTF, 303) SPPARM(1), HSTMP
303 FORMAT (' SSHAPE, deviation in Hs, should be ', F8.3,
& ', calculated ', F8.3)
ENDIF
ENDIF
!
! if mean frequency was given recalculate PKPER and restart
!
IF (.NOT.LOGPM .AND. ITPER.LT.10) THEN
ITPER = ITPER + 1
! calculate average frequency
AM0 = 0.
AM1 = 0.
DO IS = 1, MSC
AS2 = ACLOC(MDC,IS) * (SPCSIG(IS))**2
AS3 = AS2 * SPCSIG(IS)
AM0 = AM0 + AS2
AM1 = AM1 + AS3
ENDDO
! contribution of tail to total energy density
PPTAIL = PWTAIL(1) - 1. 30.80
APTAIL = 1. / (PPTAIL * (1. + PPTAIL * (FRINTH-1.))) 30.80
AM0 = AM0 * FRINTF + APTAIL * AS2 30.80
PPTAIL = PWTAIL(1) - 2. 30.80
EPTAIL = 1. / (PPTAIL * (1. + PPTAIL * (FRINTH-1.))) 30.80
AM1 = AM1 * FRINTF + EPTAIL * AS3 30.80
! Mean period:
IF ( AM1.NE.0. ) THEN 40.31
MPER = PI2 * AM0 / AM1
ELSE 40.31
CALL MSGERR(3, ' first moment is zero in calculating the') 40.31
CALL MSGERR(3, ' spectrum at boundary using param. bc.') 40.31
END IF 40.31
IF (ITEST.GE.80) WRITE (PRTEST, 72) ITPER, SPPARM(2), MPER,
& PKPER
72 FORMAT (' SSHAPE iter=', I2, ' period values:', 3F7.2)
IF (ABS(MPER-SPPARM(2)) .GT. 0.01*SPPARM(2)) THEN
! modification suggested by Mauro Sclavo
PKPER = (SPPARM(2) / MPER) * PKPER 30.80
GOTO 100
ENDIF
ENDIF
!
IF (ITPER.GE.10) THEN
CALL MSGERR(3, 'No convergence calculating the spectrum') 30.82
CALL MSGERR(3, 'at the boundary using parametric bound. cond.') 30.82
ENDIF
!
! now introduce distribution over directions
!
ADIR = PI * DEGCNV(SPPARM(3)) / 180. 40.00
IF (DSHAPL.EQ.1) THEN
DSPR = PI * SPPARM(4) / 180.
MS = MAX (DSPR**(-2) - 2., 1.)
ELSE
MS = SPPARM(4)
ENDIF
IF (MS.LT.12.) THEN
CTOT = (2.**MS) * (GAMMAF(0.5*MS+1.))**2 / (PI * GAMMAF(MS+1.))
ELSE
CTOT = SQRT (0.5*MS/PI) / (1. - 0.25/MS)
ENDIF
IF (ITEST.GE.100) THEN
ESOM = 0.
DO IS = 1, MSC
ESOM = ESOM + FRINTF * SPCSIG(IS)**2 * ACLOC(MDC,IS)
ENDDO
GAM1 = GAMMAF(0.5*MS+1.) 40.84
GAM2 = GAMMAF(MS+1.) 40.84
WRITE (PRTEST, *) ' SSHAPE dir ', 4.*SQRT(ABS(ESOM)),
& SPPARM(1), CTOT, MS, GAM1, GAM2, CTOT 40.84 40.02
ENDIF
DVERIF = .FALSE.
CTOTT = 0.
DO ID = 1, MDC
ACOS = COS(SPCDIR(ID,1) - ADIR)
IF (ACOS .GT. 0.) THEN
CDIR = CTOT * MAX (ACOS**MS, 1.E-10)
IF (.NOT.FULCIR) THEN
IF (ACOS .GE. COS(DDIR)) DVERIF = .TRUE.
ENDIF
ELSE
CDIR = 0.
ENDIF
IF (ITEST.GE.10) CTOTT = CTOTT + CDIR * DDIR
IF (ITEST.GE.100) WRITE (PRTEST, 360) ID,SPCDIR(ID,1),CDIR
360 FORMAT (' ID Spcdir Cdir: ',I3,3(1X,E10.4))
DO IS = 1, MSC
ACLOC(ID,IS) = CDIR * ACLOC(MDC,IS)
ENDDO
ENDDO
IF (ITEST.GE.10) THEN
IF (ABS(CTOTT-1.) .GT. 0.1) WRITE (PRINTF, 363) CTOTT
363 FORMAT (' SSHAPE, integral of Cdir is not 1, but:', F6.3)
ENDIF
IF (.NOT.FULCIR .AND. .NOT.DVERIF)
& CALL MSGERR (1, 'incident direction is outside sector')
!
RETURN
!
! End of subroutine SSHAPE
END
!*******************************************************************
! *
SUBROUTINE SINTRP (W1, W2, FL1, FL2, FL, SPCDIR, SPCSIG)
! *
!*******************************************************************
!
USE OCPCOMM4 40.41
USE SWCOMM3 40.41
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering, Fluid Mechanics Group |
! | P.O. Box 5048, 2600 GA Delft, the Netherlands |
! | |
! | Authors : Weimin Luo, Roeland Ris, Nico Booij |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 30.73: Nico Booij
! 30.82: IJsbrand Haagsma
! 40.00: Nico Booij
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! 30.01, Jan. 96: New subroutine for SWAN Ver. 30.01
! 30.73, Nov. 97: revised
! 40.00, Apr. 98: procedure to maintain peakedness introduced
! 30.82, Oct. 98: Update description of several variables
! 30.82, Oct. 98: Made arguments in ATAN2 double precision to prevent
! underflows
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. Purpose
!
! interpolation of spectra
!
! 3. Method (updated...)
!
! linear interpolation with peakedness maintained
! interpolated average direction and frequency are determined
! average direction and frequency of interpolated spectrum are determ.
! shifts in frequency and direction are determined from spectrum 1 and
! 2 to the interpolated spectrum
! bilinear interpolation in spectral space is used to calculate
! contributions from spectrum 1 and 2.
! in full circle cases interpolation crosses the boundary 0-360 degr.
!
! 4. Argument variables
!
! o FL : Interpolated spectrum.
! i FL1 : Input spectrum 1.
! i FL2 : Input spectrum 2.
! i SPCDIR: (*,1); spectral directions (radians) 30.82
! (*,2); cosine of spectral directions 30.82
! (*,3); sine of spectral directions 30.82
! (*,4); cosine^2 of spectral directions 30.82
! (*,5); cosine*sine of spectral directions 30.82
! (*,6); sine^2 of spectral directions 30.82
! i SPCSIG: Relative frequencies in computational domain in sigma-space 30.82
! i W1 : Weighting coefficient for spectrum 1.
! i W2 : Weighting coefficient for spectrum 2.
!
REAL FL1(MDC,MSC), FL2(MDC,MSC), FL(MDC, MSC)
REAL SPCDIR(MDC,6) 30.82
REAL SPCSIG(MSC) 30.82
REAL W1, W2
!
! 5. Parameter variables
!
! 6. Local variables
!
! ID counter of directions
! IS counter of frequencies
!
INTEGER ID, IS
!
! DOADD indicates whether or not values have to be added
!
LOGICAL DOADD
!
! ATOT1 integral over spectrum 1
! ATOT2 integral over spectrum 2
! AXTOT1 integral over x-component of spectrum 1
! AXTOT2 integral over x-component of spectrum 2
! AYTOT1 integral over y-component of spectrum 1
! AYTOT2 integral over y-component of spectrum 2
! ASTOT1 integral over Sigma * spectrum 1
! ASTOT2 integral over Sigma * spectrum 2
! ASIG1 average Sigma of spectrum 1
! ASIG2 average Sigma of spectrum 2
! DELD1 difference in direction between spectrum 1 and
! the interpolated spectrum in number of directional steps
! DELD2 same for spectrum 2
! DELSG1 shift in frequency between spectrum 1 and interpolated
! spectrum in number of frequency steps
! DELSG2 same for spectrum 2
!
REAL ATOT1, ATOT2, AXTOT1, AXTOT2, AYTOT1, AYTOT2,
& ASTOT1, ASTOT2
REAL ASIG1, ASIG2
REAL DELD1, DELD2, DELSG1, DELSG2
!
! 8. Subroutines used
!
! 9. Subroutines calling
!
! SNEXTI, RBFILE
!
! 10. Error messages
!
! 11. Remarks
!
! 12. Structure
!
! -----------------------------------------------------------------
! If W1 close to 1
! Then copy FL from FL1
! Else If W2 close to 1
! Then copy FL from FL2
! Else determine total energy in FL1 and FL2
! If energy of FL1 = 0
! Then make FL = W2 * FL2
! Else If energy of FL2 = 0
! Then make FL = W1 * FL1
! Else determine average direction of FL1 and FL2
! make ADIR = W1 * ADIR1 + W2 * ADIR2
! determine average frequency of FL1 and FL2
! make ASIG = W1 * ASIG1 + W2 * ASIG2
! determine directional shift from FL1
! determine directional shift from FL2
! determine frequency shift from FL1
! determine frequency shift from FL2
! For all spectral components do
! compose FL from components of FL1 and FL2
! -----------------------------------------------------------------
!
! 13. Source text
!
SAVE IENT
DATA IENT/0/
CALL STRACE(IENT,'SINTRP')
!
! interpolation of spectra
! ------------------------
!
IF (W1.GT.0.99) THEN
DO 101 ID=1,MDC
DO 102 IS=1,MSC
FL(ID,IS) = FL1(ID,IS)
102 CONTINUE
101 CONTINUE
ELSE IF (W1.LT.0.01) THEN
DO 201 ID=1,MDC
DO 202 IS=1,MSC
FL(ID,IS) = FL2(ID,IS)
202 CONTINUE
201 CONTINUE
ELSE
ATOT1 = 0.
ATOT2 = 0.
AXTOT1 = 0.
AXTOT2 = 0.
AYTOT1 = 0.
AYTOT2 = 0.
ASTOT1 = 0.
ASTOT2 = 0.
DO 301 ID=1,MDC
DO 302 IS=1,MSC
ATOT1 = ATOT1 + FL1(ID,IS)
AXTOT1 = AXTOT1 + FL1(ID,IS) * SPCDIR(ID,2)
AYTOT1 = AYTOT1 + FL1(ID,IS) * SPCDIR(ID,3)
ASTOT1 = ASTOT1 + FL1(ID,IS) * SPCSIG(IS)
ATOT2 = ATOT2 + FL2(ID,IS)
AXTOT2 = AXTOT2 + FL2(ID,IS) * SPCDIR(ID,2)
AYTOT2 = AYTOT2 + FL2(ID,IS) * SPCDIR(ID,3)
ASTOT2 = ASTOT2 + FL2(ID,IS) * SPCSIG(IS)
302 CONTINUE
301 CONTINUE
IF (ATOT1.LT.1.E-9) THEN
DO 401 ID=1,MDC
DO 402 IS=1,MSC
FL(ID,IS) = W2*FL2(ID,IS)
402 CONTINUE
401 CONTINUE
ELSE IF (ATOT2.LT.1.E-9) THEN
DO 501 ID=1,MDC
DO 502 IS=1,MSC
FL(ID,IS) = W1*FL1(ID,IS)
502 CONTINUE
501 CONTINUE
ELSE
! determine interpolation factors in Theta space
AXTOT = W1 * AXTOT1 + W2 * AXTOT2
AYTOT = W1 * AYTOT1 + W2 * AYTOT2
IF (ITEST.GE.80) THEN
WRITE (PRTEST, 509) ATOT1, ATOT2, 40.02
& AXTOT, AXTOT1, AXTOT2, AYTOT, AYTOT1, AYTOT2
509 FORMAT (' SINTRP factors ', 8E11.4, /, 15X, 4F7.3)
ENDIF
! DELD1 is the difference in direction between spectrum 1 and
! the interpolated spectrum in number of directional steps
DELD1 = REAL(ATAN2(DBLE(AXTOT*AYTOT1 - AYTOT*AXTOT1), 30.82
& DBLE(AXTOT*AXTOT1 + AYTOT*AYTOT1))) / DDIR 30.82
! DELD2 is the difference between spectrum 2 and
! the interpolated spectrum
DELD2 = REAL(ATAN2(DBLE(AXTOT*AYTOT2 - AYTOT*AXTOT2), 30.82
& DBLE(AXTOT*AXTOT2 + AYTOT*AYTOT2))) / DDIR 30.82
IDD1A = NINT(DELD1)
RDD1B = DELD1 - REAL(IDD1A)
IF (RDD1B .LT. 0.) THEN
IDD1A = IDD1A - 1
RDD1B = RDD1B + 1.
ENDIF
IDD1B = IDD1A + 1
RDD1B = W1 * RDD1B
RDD1A = W1 - RDD1B
IDD2A = NINT(DELD2)
RDD2B = DELD2 - REAL(IDD2A)
IF (RDD2B .LT. 0.) THEN
IDD2A = IDD2A - 1
RDD2B = RDD2B + 1.
ENDIF
IDD2B = IDD2A + 1
RDD2B = W2 * RDD2B
RDD2A = W2 - RDD2B
!
! determine interpolation factors in Sigma space
ASIG1 = ASTOT1 / ATOT1
ASIG2 = ASTOT2 / ATOT2
ATOT = W1 * ATOT1 + W2 * ATOT2
ASTOT = W1 * ASTOT1 + W2 * ASTOT2
ASIG = ASTOT / ATOT
!
! DELSG1 is shift in frequency between spectrum 1 and interpolated
! spectrum in number of frequency steps
DELSG1 = ALOG (ASIG1 / ASIG) / FRINTF
IDS1A = NINT(DELSG1)
RDS1B = DELSG1 - REAL(IDS1A)
IF (RDS1B .LT. 0.) THEN
IDS1A = IDS1A - 1
RDS1B = RDS1B + 1.
ENDIF
IDS1B = IDS1A + 1
RDS1A = 1. - RDS1B
!
! DELSG2 is shift in frequency between spectrum 2 and interpolated
! spectrum in number of frequency steps
DELSG2 = ALOG (ASIG2 / ASIG) / FRINTF
IDS2A = NINT(DELSG2)
RDS2B = DELSG2 - REAL(IDS2A)
IF (RDS2B .LT. 0.) THEN
IDS2A = IDS2A - 1
RDS2B = RDS2B + 1.
ENDIF
IDS2B = IDS2A + 1
RDS2A = 1. - RDS2B
! test output
IF (ITEST.GE.80) THEN
WRITE (PRTEST, 510) ATOT, ATOT1, ATOT2,
& AXTOT, AXTOT1, AXTOT2, AYTOT, AYTOT1, AYTOT2,
& DELD1, DELD2, DELSG1, DELSG2
510 FORMAT (' SINTRP factors ', 9E11.4, /, 15X, 4F7.3)
WRITE (PRTEST, 512) IDS1A, RDS1A, IDS1B, RDS1B,
& IDS2A, RDS2A, IDS2B, RDS2B,
& IDD1A, RDD1A, IDD1B, RDD1B,
& IDD2A, RDD2A, IDD2B, RDD2B
512 FORMAT (' SINTRP ', 8(I2, F7.3))
ENDIF
!
DO 601 ID=1,MDC
DO 602 IS=1,MSC
FL(ID,IS) = 0.
602 CONTINUE
601 CONTINUE
DO 611 ID=1,MDC
DOADD = .TRUE.
ID1A = ID + IDD1A
IF (FULCIR) THEN
IF (ID1A.LT.1) ID1A = ID1A + MDC
IF (ID1A.GT.MDC) ID1A = ID1A - MDC
ELSE
IF (ID1A.LT.1) DOADD = .FALSE.
IF (ID1A.GT.MDC) DOADD = .FALSE.
ENDIF
IF (DOADD) THEN
DO 612 IS = MAX(1,1-IDS1A), MIN(MSC,MSC-IDS1A)
FL(ID,IS) = FL(ID,IS) +
& RDD1A * RDS1A * FL1(ID1A,IS+IDS1A)
612 CONTINUE
DO 613 IS = MAX(1,1-IDS1B), MIN(MSC,MSC-IDS1B)
FL(ID,IS) = FL(ID,IS) +
& RDD1A * RDS1B * FL1(ID1A,IS+IDS1B)
613 CONTINUE
ENDIF
611 CONTINUE
DO 621 ID=1,MDC
DOADD = .TRUE.
ID1B = ID + IDD1B
IF (FULCIR) THEN
IF (ID1B.LT.1) ID1B = ID1B + MDC
IF (ID1B.GT.MDC) ID1B = ID1B - MDC
ELSE
IF (ID1B.LT.1) DOADD = .FALSE.
IF (ID1B.GT.MDC) DOADD = .FALSE.
ENDIF
IF (DOADD) THEN
DO 622 IS = MAX(1,1-IDS1A), MIN(MSC,MSC-IDS1A)
FL(ID,IS) = FL(ID,IS) +
& RDD1B * RDS1A * FL1(ID1B,IS+IDS1A)
622 CONTINUE
DO 623 IS = MAX(1,1-IDS1B), MIN(MSC,MSC-IDS1B)
FL(ID,IS) = FL(ID,IS) +
& RDD1B * RDS1B * FL1(ID1B,IS+IDS1B)
623 CONTINUE
ENDIF
621 CONTINUE
DO 631 ID=1,MDC
DOADD = .TRUE.
ID2A = ID + IDD2A
IF (FULCIR) THEN
IF (ID2A.LT.1) ID2A = ID2A + MDC
IF (ID2A.GT.MDC) ID2A = ID2A - MDC
ELSE
IF (ID2A.LT.1) DOADD = .FALSE.
IF (ID2A.GT.MDC) DOADD = .FALSE.
ENDIF
IF (DOADD) THEN
DO 632 IS = MAX(1,1-IDS2A), MIN(MSC,MSC-IDS2A)
FL(ID,IS) = FL(ID,IS) +
& RDD2A * RDS2A * FL2(ID2A,IS+IDS2A)
632 CONTINUE
DO 633 IS = MAX(1,1-IDS2B), MIN(MSC,MSC-IDS2B)
FL(ID,IS) = FL(ID,IS) +
& RDD2A * RDS2B * FL2(ID2A,IS+IDS2B)
633 CONTINUE
ENDIF
631 CONTINUE
DO 641 ID=1,MDC
DOADD = .TRUE.
ID2B = ID + IDD2B
IF (FULCIR) THEN
IF (ID2B.LT.1) ID2B = ID2B + MDC
IF (ID2B.GT.MDC) ID2B = ID2B - MDC
ELSE
IF (ID2B.LT.1) DOADD = .FALSE.
IF (ID2B.GT.MDC) DOADD = .FALSE.
ENDIF
IF (DOADD) THEN
DO 642 IS = MAX(1,1-IDS2A), MIN(MSC,MSC-IDS2A)
FL(ID,IS) = FL(ID,IS) +
& RDD2B * RDS2A * FL2(ID2B,IS+IDS2A)
642 CONTINUE
DO 643 IS = MAX(1,1-IDS2B), MIN(MSC,MSC-IDS2B)
FL(ID,IS) = FL(ID,IS) +
& RDD2B * RDS2B * FL2(ID2B,IS+IDS2B)
643 CONTINUE
ENDIF
641 CONTINUE
ENDIF
ENDIF
!
! Test output
IF (ITEST.GE.80) THEN
A1 = 0.
A2 = 0.
AA = 0.
DO 801 ID=1,MDC
DO 802 IS=1,MSC
A1 = MAX(A1,FL1(ID,IS))
A2 = MAX(A2,FL2(ID,IS))
AA = MAX(AA,FL(ID,IS))
802 CONTINUE
801 CONTINUE
WRITE (PRTEST, *) ' SINTRP, maxima ', A1, A2, AA
ENDIF
!
RETURN
! end of subroutine of SINTRP
END
!***********************************************************************
! *
REAL FUNCTION DEGCNV (DEGREE)
! *
!***********************************************************************
!
USE SWCOMM3 40.41
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 1. UPDATE
!
! SEP 1997: New for SWAN 32.01
! Cor van der Schelde - Delft Hydraulics
! 30.70, Feb. 98: test output suppressed (causes problem if subr
! is used during output
!
! 2. PURPOSE
!
! Transform degrees from nautical to cartesian or vice versa.
!
! 3. METHOD
!
! DEGCNV = 180 + dnorth - degree
!
! 4. PARAMETERLIST
!
! DEGCNV direction in cartesian or nautical degrees.
! DEGREE direction in nautical or cartesian degrees.
!
! 5. SUBROUTINES CALLING
!
! ---
!
! 6. SUBROUTINES USED
!
! NONE
!
! 7. ERROR MESSAGES
!
! NONE
!
! 8. REMARKS
!
! Nautical convention Cartesian convention
!
! 0 90
! | |
! | |
! | |
! | |
! 270 --------+-------- 90 180 --------+-------- 0
! | |
! | |
! | |
! | |
! 180 270
!
! 9. STRUCTURE
!
! ---------------------------------
! IF (NAUTICAL DEGREES) THEN
! CONVERT DEGREES
! IF (DEGREES > 360 OR < 0) THEN
! CORRECT DEGREES WITHIN 0 - 360
! ---------------------------------
!
! 10. SOURCE TEXT
!
!***********************************************************************
!
SAVE IENT
DATA IENT /0/
CALL STRACE(IENT,'DEGCNV')
!
IF ( BNAUT ) THEN
DEGCNV = 180. + DNORTH - DEGREE
ELSE
DEGCNV = DEGREE
ENDIF
!
IF (DEGCNV .GE. 360.) THEN
DEGCNV = MOD (DEGCNV, 360.)
ELSE IF (DEGCNV .LT. 0.) THEN
DEGCNV = MOD (DEGCNV, 360.) + 360.
ELSE
! DEGCNV between 0 and 360; do nothing
ENDIF
!
!
! *** end of subroutine DEGCNV ***
!
RETURN
END
!
!***********************************************************************
! *
REAL FUNCTION ANGRAD (DEGREE)
! *
!***********************************************************************
!
USE SWCOMM3 40.41
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 1. UPDATE
!
! SEP 1997: New for SWAN 32.01
! Cor van der Schelde - Delft Hydraulics
! 30.70, Feb. 98: test output suppressed (causes problem if subr
! is used during output
!
! 2. PURPOSE
!
! Transform degrees to radians
!
! 3. METHOD
!
! ANGRAD = DEGREE * PI / 180
!
! 4. PARAMETERLIST
!
! ANGRAD radians
! DEGREE degrees
!
! 5. SUBROUTINES CALLING
!
! ---
!
! 6. SUBROUTINES USED
!
! NONE
!
! 7. ERROR MESSAGES
!
! NONE
!
! 8. REMARKS
!
! NONE
!
! 9. STRUCTURE
!
! ---------------------------------
! ANGLE[radian] = ANGLE[degrees} * PI / 180
! ---------------------------------
!
! 10. SOURCE TEXT
!
!***********************************************************************
!
SAVE IENT
DATA IENT /0/
CALL STRACE(IENT,'ANGRAD')
!
ANGRAD = DEGREE * PI / 180.
!
!
! *** end of subroutine ANGRAD ***
!
RETURN
END
!
!***********************************************************************
! *
REAL FUNCTION ANGDEG (RADIAN)
! *
!***********************************************************************
!
USE SWCOMM3 40.41
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 1. UPDATE
!
! SEP 1997: New for SWAN 32.01
! Cor van der Schelde - Delft Hydraulics
! 30.70, Feb. 98: test output suppressed (causes problem if subr
! is used during output
!
! 2. PURPOSE
!
! Transform radians to degrees
!
! 3. METHOD
!
! ANGDEG = RADIAN * 180 / PI
!
! 4. PARAMETERLIST
!
! RADIAN radians
! ANGDEG degrees
!
! 5. SUBROUTINES CALLING
!
! ---
!
! 6. SUBROUTINES USED
!
! NONE
!
! 7. ERROR MESSAGES
!
! NONE
!
! 8. REMARKS
!
! NONE
!
! 9. STRUCTURE
!
! ---------------------------------
! ANGLE[degrees] = ANGLE[radians} * 180 / PI
! ---------------------------------
!
! 10. SOURCE TEXT
!
!***********************************************************************
!
SAVE IENT
DATA IENT /0/
CALL STRACE(IENT,'ANGDEG')
!
ANGDEG = RADIAN * 180. / PI
!
!
! *** end of subroutine ANGDEG ***
!
RETURN
END
!
!***********************************************************************
! *
SUBROUTINE HSOBND (AC2 ,SPCSIG,HSIBC ,KGRPNT) 40.00
! *
!***********************************************************************
!
USE OCPCOMM4 40.41
USE SWCOMM3 40.41
USE M_PARALL 40.31
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 30.72: IJsbrand Haagsma
! 32.01: Roeland Ris
! 30.70: Nico Booij
! 40.00: Nico Booij
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! 32.01, Sep. 97: new for SWAN
! 30.72, Jan. 98: Changed number of elements for HSI to MCGRD
! 30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN
! 30.70, Feb. 98: structure scheme corrected
! 40.00, Mar. 98: integration method changed (as in SNEXTI)
! structure corrected
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. Purpose
!
! Compare computed significant wave height with the value of
! the significant wave height as predescribed by the user. If
! the values differ more than e.g. 10 % give an error message
! and the gridpoints where the error has been located
!
! 3. Method
!
! ---
!
! 4. Argument variables
!
! SPCSIG: input Relative frequencies in computational domain in 30.72
! sigma-space 30.72
!
REAL SPCSIG(MSC) 30.72
!
! REALS:
! ------
! AC2 action density
! HSI significant wave height at boundary (using SWAN
! resolution (has thus not to be equal to the WAVEC
! significant wave height )
! ETOT total energy in a gridpoint
! DS increment in frequency space
! DDIR increment in directional space
! HSC computed wave height after SWAN computation
! EFTAIL contribution of tail to spectrum
!
! INTEGERS:
! ---------
! KGRPNT values of grid indices
!
! 5. SUBROUTINES CALLING
!
! ---
!
! 6. SUBROUTINES USED
!
! TRACE
!
! 7. ERROR MESSAGES
!
! NONE
!
! 8. REMARKS
!
! NONE
!
! 9. STRUCTURE
!
! ------------------------------------------------------------------
! for all computational grid points do 30.70
! if HSI is non-zero
! then compute Hs from action density array
! if relative difference is large than HSRERR
! then write error message
! -------------------------------------------------------------------
!
! 10. SOURCE TEXT
!
!***********************************************************************
!
REAL AC2(MDC,MSC,MCGRD) ,HSIBC(MCGRD) 30.72
!
REAL ETOT ,HSC 40.00
!
INTEGER ID ,IS ,IX ,IY ,INDX
!
LOGICAL HSRR
!
INTEGER KGRPNT(MXC,MYC)
!
SAVE IENT, HSRR
DATA IENT/0/, HSRR/.TRUE./
CALL STRACE (IENT, 'HSOBND')
!
! *** initializing ***
!
HSRR = .TRUE. 40.03
!
DO IY = MYC, 1, -1
DO IX = 1, MXC
INDX = KGRPNT(IX,IY)
IF ( HSIBC(INDX) .GT. 1.E-25 ) THEN
! *** compute Hs for boundary point (without tail) ***
ETOT = 0.
DO ID = 1, MDC
DO IS = 1, MSC 40.00
ETOT = ETOT + SPCSIG(IS)**2 * AC2(ID,IS,INDX) 40.00
ENDDO
ENDDO
IF (ETOT .GT. 1.E-8) THEN
HSC = 4. * SQRT(ETOT*FRINTF*DDIR)
ELSE
HSC = 0.
ENDIF
HSREL = ABS(HSIBC(INDX) - HSC) / HSIBC(INDX) 40.51 40.00
IF (HSREL .GT. HSRERR) THEN
IF ( HSRR ) THEN
WRITE (PRINTF,*) ' ** WARNING : ',
& 'Differences in wave height at the boundary'
WRITE (PRINTF,802) HSRERR
802 FORMAT (' Relative difference between input and ',
& 'computation >= ', F6.2)
WRITE (PRINTF,*) ' Hs[m]',
& ' Hs[m] Hs[-]'
WRITE (PRINTF,*) ' ix iy index (input)',
& ' (computed) (relative)'
WRITE (PRINTF,*) ' ----------------------------',
& '----------------------'
HSRR = .FALSE.
ENDIF
WRITE (PRINTF,'(2(1x,I5),I7,3(1x,F10.2))')
& IX+MXF-1, IY+MYF-1, INDX, HSIBC(INDX), HSC, HSREL
ENDIF
ENDIF
ENDDO
ENDDO
WRITE(PRINTF,*)
!
IF ( ITEST .GE. 150 ) THEN
WRITE(PRINTF,*) 'Values of wave height at boundary (HSOBND)'
WRITE(PRINTF,*) '------------------------------------------'
DO IY = MYC, 1, -1
WRITE (PRINTF,'(13F8.3)') ( HSIBC(KGRPNT(IX,IY)), IX=1 , MXC)
ENDDO
ENDIF
!
! *** end of subroutine HSOBND ***
!
RETURN
END
!
!*****************************************************************
! *
SUBROUTINE CHGBAS (X1, X2, PERIOD, Y1, Y2, N1, N2,
& ITEST, PRTEST)
! *
!*****************************************************************
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering, Fluid Mechanics Group |
! | P.O. Box 5048, 2600 GA Delft, the Netherlands |
! | |
! | Authors : G. van Vledder, N. Booij |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Update history
!
! ver 20.48: also accomodates periodic variables such as directions
!
! 1. Purpose
!
! change x-basis of a discretized y-function
!
! 2. Method
!
! A piecewise constant representation of the functions is assumed
!
! first boundaries of a cell in X1 are determined
! then it is determined whether there are overlaps with cells
! in X2. if so Y1*common length is added to Y2
! Finally Y2 values are divided by cell lengths
!
! 3. Parameter list
!
! Name I/O Type Description
!
! X1 i ra x-coordinates of input grid
! X2 i ra x-coordinates of output grid
! PERIOD i r period, i.e. x-axis is periodic if period>0
! e.g. spectral directions
! Y1 i ra function values of input grid
! Y2 o ra function values of output grid
! N1 i i number of x-values of input grid
! N2 i i number of x-values of output grid
!
! 4. Subroutines used
!
! ---
!
! 5. Error messages
!
! 6. Remarks
!
! Cell boundaries in X1 are: X1A and X1B
! X2 is assumed to be monotonically increasing; this is checked
! X1 is assumed to be monotonous but not necessarily increasing
!
! 7. Structure
!
! ------------------------------------------------------------------
! Make all values of Y2 = 0
! For each cell in X1 do
! determine boundaries of cell in X1
! --------------------------------------------------------------
! For each cell in X2 do
! determine overlap with cell in X1; limits: RLOW and RUPP
! add to Y2: Y1 * length of overlapping interval
! ------------------------------------------------------------------
! For each cell in X2 do
! divide Y2 value by cell length
! ------------------------------------------------------------------
!
! 8. Source text
!
INTEGER I1, I2, N1, N2, ITEST, PRTEST
REAL X1(N1), Y1(N1), X2(N2), Y2(N2), PERIOD
REAL X1A, X1B, X2A, X2B, RLOW, RUPP
LOGICAL TWICE
SAVE IENT
DATA IENT /0/
CALL STRACE (IENT, 'CHGBAS')
!
! initialize output data
!
DO I2 = 1, N2
Y2(I2) = 0.
ENDDO
DO I2 = 2, N2
IF (X2(I2).LE.X2(I2-1))
& CALL MSGERR (2, 'subr. CHGBAS: values of X2 not increasing')
ENDDO
! boundaries of the range in X2
X2LO = 1.5 * X2(1) - 0.5 * X2(2)
X2HI = 1.5 * X2(N2) - 0.5 * X2(N2-1)
TWICE = .FALSE.
!
! loop over cells in X1
!
DO 300 I1 = 1, N1
IF (ABS(Y1(I1)) .LT. 1.E-20) GOTO 300
!
! determine cell boundaries in X1
!
IF (I1.EQ.1) THEN
X1A = 1.5 * X1(1) - 0.5 * X1(2)
ELSE
X1A = 0.5 * (X1(I1) + X1(I1-1))
ENDIF
IF (I1.EQ.N1) THEN
X1B = 1.5 * X1(N1) - 0.5 * X1(N1-1)
ELSE
X1B = 0.5 * (X1(I1) + X1(I1+1))
ENDIF
!
! swap X1A and X1B if X1A > X1B
!
IF (X1A.GT.X1B) THEN
RR = X1A
X1A = X1B
X1B = RR
ENDIF
IF (PERIOD.LE.0.) THEN
IF (X1A.GT.X2HI) GOTO 300
IF (X1B.LT.X2LO) GOTO 300
ELSE
! X is periodic; move interval in X1 if necessary
TWICE = .FALSE.
IADD = 0
60 IF (X1B.GT.X2HI) THEN
X1A = X1A - PERIOD
X1B = X1B - PERIOD
IADD = IADD + 1
IF (IADD.GT.99)
& CALL MSGERR (2, 'endless loop in CHGBAS')
GOTO 60
ENDIF
70 IF (X1A.LT.X2LO) THEN
X1A = X1A + PERIOD
X1B = X1B + PERIOD
IADD = IADD + 1
IF (IADD.GT.99)
& CALL MSGERR (2, 'endless loop in CHGBAS')
GOTO 70
ENDIF
IF (X1A.GT.X2HI) GOTO 300
IF (X1B.LT.X2LO) GOTO 300
IF (X1A.LT.X2LO .AND. X1A+PERIOD.LT.X2HI) TWICE = .TRUE.
IF (X1B.GT.X2HI .AND. X1B-PERIOD.GT.X2LO) TWICE = .TRUE.
ENDIF
!
! loop over cells in X2
!
100 DO 200 I2 = 1, N2
IF (I2.EQ.1) THEN
X2A = X2LO
ELSE
X2A = 0.5 * (X2(I2) + X2(I2-1))
ENDIF
IF (I2.EQ.N2) THEN
X2B = X2HI
ELSE
X2B = 0.5 * (X2(I2) + X2(I2+1))
ENDIF
!
! (RLOW,RUPP) is overlapping interval of (X1A,X1B) and (X2A,X2B)
!
IF (X1A.LT.X2B) THEN
RLOW = MAX (X1A, X2A)
ELSE
GOTO 200
ENDIF
IF (X1B.GT.X2A) THEN
RUPP = MIN (X1B, X2B)
ELSE
GOTO 200
ENDIF
IF (RUPP.LT.RLOW) THEN
CALL MSGERR (3, 'interpolation error')
WRITE (PRTEST, 140) I1, X1A, X1B, I2, X2A, X2B
140 FORMAT (' I, XA, XB ', 2(I3, 2(1X,E12.4)))
ELSE
Y2(I2) = Y2(I2) + Y1(I1) * (RUPP-RLOW)
ENDIF
200 CONTINUE
!
! Cell in X1 covers both ends of sector boundary
IF (TWICE) THEN
IF (X1A.LT.X2LO) THEN
X1A = X1A + PERIOD
X1B = X1B + PERIOD
ENDIF
IF (X1B.GT.X2HI) THEN
X1A = X1A - PERIOD
X1B = X1B - PERIOD
ENDIF
TWICE = .FALSE.
GOTO 100
ENDIF
300 CONTINUE
!
DO I2 = 1, N2
IF (I2.EQ.1) THEN
CELLEN = X2(2) - X2(1)
ELSE IF (I2.EQ.N2) THEN
CELLEN = X2(N2) - X2(N2-1)
ELSE
CELLEN = 0.5 * (X2(I2+1) - X2(I2-1))
ENDIF
! divide Y2 by cell length
Y2(I2) = Y2(I2) / CELLEN
ENDDO
IF (ITEST.GE.160) THEN
WRITE (PRTEST, 84) N1, N2
84 FORMAT (' test CHGBAS ', 2I5)
WRITE (PRTEST, 85) (X1(II), II = 1, N1)
WRITE (PRTEST, 85) (Y1(II), II = 1, N1)
WRITE (PRTEST, 85) (X2(II), II = 1, N2)
WRITE (PRTEST, 85) (Y2(II), II = 1, N2)
85 FORMAT (10 (1X,E10.3))
ENDIF
!
RETURN
END
!
!********************************************************************
! *
REAL FUNCTION GAMMAF(XX)
! *
!********************************************************************
!
! Updates
! ver 30.70, Oct 1997 by N.Booij: new subroutine
!
! Purpose
! Compute the transcendental function Gamma
!
! Subroutines used
! GAMMLN (Numerical Recipes)
!
REAL XX, YY, ABIG 40.00
SAVE IENT, ABIG
DATA IENT /0/, ABIG /30./
CALL STRACE (IENT, 'GAMMAF')
YY = GAMMLN(XX)
IF (YY.GT.ABIG) YY = ABIG
IF (YY.LT.-ABIG) YY = -ABIG
GAMMAF = EXP(YY)
RETURN
END
!********************************************************************
! *
FUNCTION GAMMLN(XX)
! *
!********************************************************************
!
! Method:
! function is copied from: Press et al., "Numerical Recipes"
!
DOUBLE PRECISION COF(6),STP,HALF,ONE,FPF,X,TMP,SER
DATA COF,STP/76.18009173D0,-86.50532033D0,24.01409822D0,
& -1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/
DATA HALF,ONE,FPF/0.5D0,1.0D0,5.5D0/
X=XX-ONE
TMP=X+FPF
TMP=(X+HALF)*LOG(TMP)-TMP
SER=ONE
DO 11 J=1,6
X=X+ONE
SER=SER+COF(J)/X
11 CONTINUE
GAMMLN=TMP+LOG(STP*SER)
RETURN
END
!***********************************************************************
! *
SUBROUTINE WRSPEC (NREF, ACLOC)
! *
!***********************************************************************
USE OCPCOMM4 40.41
USE SWCOMM3 40.41
USE OUTP_DATA 40.13
IMPLICIT NONE 40.13
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 40.00, 40.13: Nico Booij
! 40.41: Marcel Zijlema
!
! 1. UPDATE
!
! new subroutine, update 40.00
! 40.03, Mar. 00: precision increased; 2 decimals more in output table
! 40.13, July 01: variable format using module OUTP_DATA
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. PURPOSE
!
! Writing of action density spectrum in Swan standard format
!
! 3. METHOD
!
!
! 4. Argument variables
!
! NREF int input unit ref. number of output file
! ACLOC real local 2-D spectrum or source term at one
! output location
INTEGER, INTENT(IN) :: NREF
REAL, INTENT(IN) :: ACLOC(1:MDC,1:MSC)
! 5. Parameter variables
!
! 6. Local variables
!
! ID counter of spectral directions
! IS counter of spectral frequencies
!
INTEGER :: ID, IS
!
! EFAC multiplication factor written to file
!
REAL :: EFAC
!
! 8. Subroutines used
!
! 9. Subroutines calling
!
! SWOUTP (SWAN/OUTP)
!
! 10. Error messages
!
! 11. Remarks
!
! 12. Structure
!
! ----------------------------------------------------------------
! determine maximum value of ACLOC
! if maximum = 0
! then write 'ZERO' to file
! else write 'FACTOR'
! determine multiplication factor, write this to file
! write values of ACLOC/factor to file
! ----------------------------------------------------------------
!
! 13. Source text
!
INTEGER, SAVE :: IENT = 0 40.13
IF (LTRACE) CALL STRACE (IENT, 'WRSPEC')
!
! first determine maximum energy density
EFAC = 0.
DO ID = 1, MDC
DO IS = 1, MSC
IF (ACLOC(ID,IS).GE.0.) THEN
EFAC = MAX (EFAC, ACLOC(ID,IS))
ELSE
EFAC = MAX (EFAC, 10.*ABS(ACLOC(ID,IS)))
ENDIF
ENDDO
ENDDO
IF (EFAC .LE. 1.E-10) THEN
WRITE (NREF, 12) 'ZERO' 40.00
12 FORMAT (A4)
ELSE
EFAC = 1.01 * EFAC * 10.**(-DEC_SPEC) 40.13
! factor PI/180 introduced to account for change from rad to degr
! factor 2*PI to account for transition from rad/s to Hz
WRITE (NREF, 95) EFAC * 2. * PI**2 / 180.
95 FORMAT ('FACTOR', /, E18.8) 40.13
DO IS = 1, MSC
! write spectral energy densities to file
WRITE (NREF, FIX_SPEC) (NINT(ACLOC(ID,IS)/EFAC), ID=1,MDC) 40.13
ENDDO
ENDIF
RETURN
! end of subroutine WRSPEC
END
!****************************************************************
!
SUBROUTINE SWACC(AC2, AC2OLD, ACNRMS, ISSTOP, IDCMIN, IDCMAX)
!
!****************************************************************
!
USE SWCOMM3 40.41
USE OCPCOMM4 40.41
!
IMPLICIT NONE
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering and Geosciences |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmer: Marcel Zijlema |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 40.23: Marcel Zijlema
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! 40.23, Sep. 02: New subroutine
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. Purpose
!
! Determine some infinity norms meant for stop criterion
!
! 4. Argument variables
!
! AC2 action density
! AC2OLD action density at previous iteration
! ACNRMS array containing infinity norms
! IDCMIN integer array containing minimum counter of directions
! IDCMAX integer array containing maximum counter of directions
! ISSTOP maximum frequency counter in this sweep
!
INTEGER IDCMIN(MSC), IDCMAX(MSC), ISSTOP
REAL AC2(MDC,MSC,MCGRD), AC2OLD(MDC,MSC), ACNRMS(2)
!
! 6. Local variables
!
! DIFFAC: difference between AC2 and AC2OLD
! ID : counter of direction
! IDDUM : uncorrected counter of direction
! IENT : number of entries
! IS : counter of frequency
!
INTEGER ID, IDDUM, IS, IENT
REAL DIFFAC
!
! 7. Common blocks used
!
!
! 8. Subroutines used
!
! STRACE Tracing routine for debugging
!
! 9. Subroutines calling
!
! SWOMPU (in SWANCOM1)
!
! 12. Structure
!
! determine infinity norms |ac2 - ac2old| and |ac2|
! in selected sweep
!
! 13. Source text
!
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'SWACC')
DO IS = 1, ISSTOP
DO IDDUM = IDCMIN(IS), IDCMAX(IS)
ID = MOD ( IDDUM - 1 + MDC , MDC ) + 1
!
! *** determine infinity norms |ac2 - ac2old| and |ac2|
!
DIFFAC = ABS(AC2(ID,IS,KCGRD(1)) - AC2OLD(ID,IS))
IF (DIFFAC.GT.ACNRMS(1)) ACNRMS(1) = DIFFAC
IF (ABS(AC2(ID,IS,KCGRD(1))).GT.ACNRMS(2))
& ACNRMS(2) = ABS(AC2(ID,IS,KCGRD(1)))
END DO
END DO
RETURN
END
!TIMG!****************************************************************
!TIMG!
!TIMG SUBROUTINE SWTSTA (ITIMER)
!TIMG!
!TIMG!****************************************************************
!TIMG!
!TIMG USE TIMECOMM 40.41
!TIMG USE OCPCOMM4 40.41
!TIMG!
!TIMG IMPLICIT NONE
!TIMG!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering and Geosciences |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmer: Marcel Zijlema |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!TIMG!
!TIMG! 0. Authors
!TIMG!
!TIMG! 40.23: Marcel Zijlema
!TIMG! 40.41: Marcel Zijlema
!TIMG!
!TIMG! 1. Updates
!TIMG!
!TIMG! 40.23, Aug. 02: New subroutine
!TIMG! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!TIMG!
!TIMG! 2. Purpose
!TIMG!
!TIMG! Start timing
!TIMG!
!TIMG! 3. Method
!TIMG!
!TIMG! Get cpu and wall-clock times and store
!TIMG!
!TIMG! 4. Argument variables
!TIMG!
!TIMG! ITIMER number of timer to be used
!TIMG!
!TIMG INTEGER :: ITIMER
!TIMG!
!TIMG! 6. Local variables
!TIMG!
!TIMG! C : clock count of processor
!TIMG! I : index in LISTTM, loop variable
!TIMG! IFOUND: index in LISTTM, location of ITIMER
!TIMG! IFREE : index in LISTTM, first free position
!TIMG! M : maximum clock count
!TIMG! R : number of clock counts per second
!TIMG!F95! TIMER : current real cpu-time
!TIMG! TIMER1: current cpu-time used
!TIMG! TIMER2: current wall-clock time used
!TIMG!
!TIMG INTEGER :: I, IFOUND, IFREE
!TIMG INTEGER :: C, R, M
!TIMG!F95 REAL :: TIMER
!TIMG DOUBLE PRECISION :: TIMER1, TIMER2
!TIMG!
!TIMG! 7. Common blocks used
!TIMG!
!TIMG!
!TIMG! 8. Subroutines used
!TIMG!
!TIMG!F95! CPU_TIME Returns real value from cpu-time clock
!TIMG! SYSTEM_CLOCK Returns integer values from a real-time clock
!TIMG!
!TIMG! 9. Subroutines calling
!TIMG!
!TIMG! SWMAIN, SWCOMP, SWOMPU
!TIMG!
!TIMG! 12. Structure
!TIMG!
!TIMG! Get and store the cpu and wall-clock times
!TIMG!
!TIMG! 13. Source text
!TIMG!
!TIMG
!TIMG!
!TIMG! --- check whether a valid timer number is given
!TIMG!
!TIMG IF (ITIMER.LE.0 .OR. ITIMER.GT.NSECTM) THEN
!TIMG WRITE(PRINTF,*) 'SWTSTA: ITIMER out of range: ',
!TIMG & ITIMER, 1, NSECTM
!TIMG STOP
!TIMG END IF
!TIMG!
!TIMG! --- check whether timing for ITIMER was started already,
!TIMG! also determine first free location in LISTTM
!TIMG!
!TIMG IFOUND=0
!TIMG IFREE =0
!TIMG I =0
!TIMG 100 IF (I.LT.LASTTM .AND. (IFOUND.EQ.0 .OR. IFREE.EQ.0)) THEN
!TIMG I=I+1
!TIMG IF (LISTTM(I).EQ.ITIMER) THEN
!TIMG IFOUND=I
!TIMG END IF
!TIMG IF (IFREE.EQ.0 .AND. LISTTM(I).EQ.-1) THEN
!TIMG IFREE =I
!TIMG END IF
!TIMG GOTO 100
!TIMG END IF
!TIMG
!TIMG IF (IFOUND.EQ.0 .AND. IFREE.EQ.0 .AND. LASTTM.LT.MXTIMR) THEN
!TIMG LASTTM=LASTTM+1
!TIMG IFREE =LASTTM
!TIMG END IF
!TIMG!
!TIMG! --- produce warning if found in the list
!TIMG!
!TIMG IF (IFOUND.GT.0) THEN
!TIMG WRITE(PRINTF,*)
!TIMG & 'SWTSTA: warning: previous timing for section ',
!TIMG & ITIMER,' not closed properly/will be ignored.'
!TIMG END IF
!TIMG!
!TIMG! --- produce error if not found and no free position available
!TIMG!
!TIMG IF (IFOUND.EQ.0 .AND. IFREE.EQ.0) THEN
!TIMG WRITE(PRINTF,*)
!TIMG & 'SWTSTA: maximum number of simultaneous timers',
!TIMG & ' exceeded:',MXTIMR
!TIMG STOP
!TIMG END IF
!TIMG!
!TIMG! --- register ITIMER in appropriate location of LISTTM
!TIMG!
!TIMG IF (IFOUND.EQ.0) THEN
!TIMG IFOUND=IFREE
!TIMG END IF
!TIMG LISTTM(IFOUND)=ITIMER
!TIMG!
!TIMG! --- get current cpu/wall-clock time and store in TIMERS
!TIMG!
!TIMG TIMER1=0D0
!TIMG!F95 CALL CPU_TIME (TIMER)
!TIMG!F95 TIMER1=DBLE(TIMER)
!TIMG CALL SYSTEM_CLOCK (C,R,M)
!TIMG TIMER2=DBLE(C)/DBLE(R)
!TIMG
!TIMG TIMERS(IFOUND,1)=TIMER1
!TIMG TIMERS(IFOUND,2)=TIMER2
!TIMG
!TIMG RETURN
!TIMG END
!TIMG!****************************************************************
!TIMG!
!TIMG SUBROUTINE SWTSTO (ITIMER)
!TIMG!
!TIMG!****************************************************************
!TIMG!
!TIMG USE TIMECOMM 40.41
!TIMG USE OCPCOMM4 40.41
!TIMG!
!TIMG IMPLICIT NONE
!TIMG!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering and Geosciences |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmer: Marcel Zijlema |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!TIMG!
!TIMG! 0. Authors
!TIMG!
!TIMG! 40.23: Marcel Zijlema
!TIMG! 40.41: Marcel Zijlema
!TIMG!
!TIMG! 1. Updates
!TIMG!
!TIMG! 40.23, Aug. 02: New subroutine
!TIMG! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!TIMG!
!TIMG! 2. Purpose
!TIMG!
!TIMG! Stop timing
!TIMG!
!TIMG! 3. Method
!TIMG!
!TIMG! Get cpu and wall-clock times and store
!TIMG!
!TIMG! 4. Argument variables
!TIMG!
!TIMG! ITIMER number of timer to be used
!TIMG!
!TIMG INTEGER :: ITIMER
!TIMG!
!TIMG! 6. Local variables
!TIMG!
!TIMG! C : clock count of processor
!TIMG! I : index in LISTTM, loop variable
!TIMG! IFOUND: index in LISTTM, location of ITIMER
!TIMG! M : maximum clock count
!TIMG! R : number of clock counts per second
!TIMG!F95! TIMER : current real cpu-time
!TIMG! TIMER1: current cpu-time used
!TIMG! TIMER2: current wall-clock time used
!TIMG!
!TIMG INTEGER :: I, IFOUND
!TIMG INTEGER :: C, R, M
!TIMG!F95 REAL :: TIMER
!TIMG DOUBLE PRECISION :: TIMER1, TIMER2
!TIMG!
!TIMG! 7. Common blocks used
!TIMG!
!TIMG!
!TIMG! 8. Subroutines used
!TIMG!
!TIMG!F95! CPU_TIME Returns real value from cpu-time clock
!TIMG! SYSTEM_CLOCK Returns integer values from a real-time clock
!TIMG!
!TIMG! 9. Subroutines calling
!TIMG!
!TIMG! SWMAIN, SWCOMP, SWOMPU
!TIMG!
!TIMG! 12. Structure
!TIMG!
!TIMG! Get and store the cpu and wall-clock times
!TIMG!
!TIMG! 13. Source text
!TIMG!
!TIMG
!TIMG!
!TIMG! --- check whether a valid timer number is given
!TIMG!
!TIMG IF (ITIMER.LE.0 .OR. ITIMER.GT.NSECTM) THEN
!TIMG WRITE(PRINTF,*) 'SWTSTO: ITIMER out of range: ',
!TIMG & ITIMER, 1, NSECTM
!TIMG STOP
!TIMG END IF
!TIMG!
!TIMG! --- check whether timing for ITIMER was started already,
!TIMG! also determine first free location in LISTTM
!TIMG!
!TIMG IFOUND=0
!TIMG I =0
!TIMG 100 IF (I.LT.LASTTM .AND. IFOUND.EQ.0) THEN
!TIMG I=I+1
!TIMG IF (LISTTM(I).EQ.ITIMER) THEN
!TIMG IFOUND=I
!TIMG END IF
!TIMG GOTO 100
!TIMG END IF
!TIMG!
!TIMG! --- produce error if not found
!TIMG!
!TIMG IF (IFOUND.EQ.0) THEN
!TIMG WRITE(PRINTF,*)
!TIMG & 'SWTSTO: section ',ITIMER,' not found',
!TIMG & ' in list of active timings'
!TIMG STOP
!TIMG END IF
!TIMG!
!TIMG! --- get current cpu/wall-clock time
!TIMG!
!TIMG TIMER1=0D0
!TIMG!F95 CALL CPU_TIME (TIMER)
!TIMG!F95 TIMER1=DBLE(TIMER)
!TIMG CALL SYSTEM_CLOCK (C,R,M)
!TIMG TIMER2=DBLE(C)/DBLE(R)
!TIMG!
!TIMG! --- calculate elapsed time since start of timing,
!TIMG! store in appropriate location in DCUMTM,
!TIMG! increment number of timings for current section
!TIMG!
!TIMG DCUMTM(ITIMER,1)=DCUMTM(ITIMER,1)+(TIMER1-TIMERS(IFOUND,1))
!TIMG DCUMTM(ITIMER,2)=DCUMTM(ITIMER,2)+(TIMER2-TIMERS(IFOUND,2))
!TIMG NCUMTM(ITIMER) =NCUMTM(ITIMER)+1
!TIMG!
!TIMG! --- free appropriate location of LISTTM,
!TIMG! adjust last occupied position of LISTTM
!TIMG!
!TIMG IF (IFOUND.GT.0) THEN
!TIMG LISTTM(IFOUND)=-1
!TIMG END IF
!TIMG 200 IF (LASTTM.GT.1 .AND. LISTTM(LASTTM).EQ.-1) THEN
!TIMG LASTTM=LASTTM-1
!TIMG GOTO 200
!TIMG END IF
!TIMG IF (LISTTM(LASTTM).EQ.-1) LASTTM=0
!TIMG
!TIMG RETURN
!TIMG END
!TIMG!****************************************************************
!TIMG!
!TIMG SUBROUTINE SWPRTI
!TIMG!
!TIMG!****************************************************************
!TIMG!
!TIMG USE OCPCOMM4 40.41
!TIMG USE TIMECOMM 40.41
!TIMG USE M_PARALL 40.31
!TIMG!PUN USE SIZES, ONLY: MNPROC, MYPROC 40.95
!TIMG
!TIMG IMPLICIT NONE
!TIMG!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering and Geosciences |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmer: Marcel Zijlema |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!TIMG!
!TIMG! 0. Authors
!TIMG!
!TIMG! 40.23: Marcel Zijlema
!TIMG! 40.30: Marcel Zijlema
!TIMG! 40.41: Marcel Zijlema
!TIMG!
!TIMG! 1. Updates
!TIMG!
!TIMG! 40.23, Aug. 02: New subroutine
!TIMG! 40.30, Jan. 03: introduction distributed-memory approach using MPI
!TIMG! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!TIMG!
!TIMG! 2. Purpose
!TIMG!
!TIMG! Print timings info
!TIMG!
!TIMG! 6. Local variables
!TIMG!
!TIMG! IDEBUG: level of timing output requested:
!TIMG! 0 - no output for detailed timings
!TIMG! 1 - aggregate output for detailed timings
!TIMG! 2 - complete output for all detailed timings
!TIMG! IENT : number of entries
!TIMG! J : loop counter
!TIMG! K : loop counter
!TIMG! MYPRC : own process number
!TIMG! TABLE : array for computing aggregate cpu- and wallclock-times
!TIMG!
!TIMG INTEGER :: IENT, J, K, IDEBUG, MYPRC
!TIMG DOUBLE PRECISION :: TABLE(30,2)
!TIMG PARAMETER (IDEBUG=0)
!TIMG!
!TIMG! 7. Common blocks used
!TIMG!
!TIMG!
!TIMG! 8. Subroutines used
!TIMG!
!TIMG! STRACE Tracing routine for debugging
!TIMG!
!TIMG! 9. Subroutines calling
!TIMG!
!TIMG! SWMAIN (in SWANMAIN)
!TIMG!
!TIMG! 12. Structure
!TIMG!
!TIMG! Compile table with overview of cpu/wall clock time used in
!TIMG! important parts of SWAN and write to PRINT file
!TIMG!
!TIMG! 13. Source text
!TIMG!
!TIMG SAVE IENT
!TIMG DATA IENT/0/
!TIMG IF (LTRACE) CALL STRACE (IENT,'SWPRTI')
!TIMG!
!TIMG MYPRC = INODE
!TIMG!PUN IF ( MNPROC>1 ) MYPRC = MYPROC
!TIMG!
!TIMG! --- compile table with overview of cpu/wall clock time used in
!TIMG! important parts of SWAN and write to PRINT file
!TIMG!
!TIMG IF ( ITEST.GE.1 .OR. IDEBUG.GE.1 ) THEN
!TIMG!
!TIMG! --- initialise table to zero
!TIMG!
!TIMG DO K = 1, 30
!TIMG DO J = 1, 2
!TIMG TABLE(K,J) = 0D0
!TIMG END DO
!TIMG END DO
!TIMG!
!TIMG! --- compute times for basic blocks
!TIMG!
!TIMG DO J = 1, 2
!TIMG!
!TIMG! --- total run-time
!TIMG!
!TIMG TABLE(1,J) = DCUMTM(1,J)
!TIMG!
!TIMG! --- initialisation, reading, preparation:
!TIMG!
!TIMG DO K = 2, 7
!TIMG TABLE(2,J) = TABLE(2,J) + DCUMTM(K,J)
!TIMG END DO
!TIMG!
!TIMG! --- domain decomposition:
!TIMG!
!TIMG TABLE(2,J) = TABLE(2,J) + DCUMTM(211,J)
!TIMG TABLE(2,J) = TABLE(2,J) + DCUMTM(212,J)
!TIMG!JAC TABLE(2,J) = TABLE(2,J) + DCUMTM(215,J)
!TIMG TABLE(2,J) = TABLE(2,J) + DCUMTM(201,J)
!TIMG!
!TIMG! --- total calculation including communication:
!TIMG!
!TIMG TABLE(3,J) = TABLE(3,J) + DCUMTM(8,J)
!TIMG!
!TIMG! --- output:
!TIMG!
!TIMG TABLE(5,J) = TABLE(5,J) + DCUMTM(9,J)
!TIMG!
!TIMG! --- exchanging data:
!TIMG!
!TIMG TABLE(7,J) = TABLE(7,J) + DCUMTM(213,J)
!TIMG!
!TIMG! --- solving system:
!TIMG!
!TIMG TABLE(9,J) = TABLE(9,J) + DCUMTM(119,J)
!TIMG TABLE(9,J) = TABLE(9,J) + DCUMTM(120,J)
!TIMG!
!TIMG! --- global reductions:
!TIMG!
!TIMG TABLE(10,J) = TABLE(10,J) + DCUMTM(202,J)
!TIMG!
!TIMG! --- collecting data:
!TIMG!
!TIMG TABLE(11,J) = TABLE(11,J) + DCUMTM(214,J)
!TIMG!
!TIMG! --- setup:
!TIMG!
!TIMG TABLE(12,J) = TABLE(12,J) + DCUMTM(106,J)
!TIMG!
!TIMG! --- propagation velocities:
!TIMG!
!TIMG TABLE(14,J) = TABLE(14,J) + DCUMTM(111,J)
!TIMG TABLE(14,J) = TABLE(14,J) + DCUMTM(113,J)
!TIMG TABLE(14,J) = TABLE(14,J) + DCUMTM(114,J)
!TIMG!
!TIMG! --- x-y advection:
!TIMG!
!TIMG TABLE(15,J) = TABLE(15,J) + DCUMTM(140,J)
!TIMG!
!TIMG! --- sigma advection:
!TIMG!
!TIMG TABLE(16,J) = TABLE(16,J) + DCUMTM(141,J)
!TIMG!
!TIMG! --- theta advection:
!TIMG!
!TIMG TABLE(17,J) = TABLE(17,J) + DCUMTM(142,J)
!TIMG!
!TIMG! --- wind:
!TIMG!
!TIMG TABLE(18,J) = TABLE(18,J) + DCUMTM(132,J)
!TIMG!
!TIMG! --- whitecapping:
!TIMG!
!TIMG TABLE(19,J) = TABLE(19,J) + DCUMTM(133,J)
!TIMG!
!TIMG! --- bottom friction:
!TIMG!
!TIMG TABLE(20,J) = TABLE(20,J) + DCUMTM(130,J)
!TIMG!
!TIMG! --- wave breaking:
!TIMG!
!TIMG TABLE(21,J) = TABLE(21,J) + DCUMTM(131,J)
!TIMG!
!TIMG! --- quadruplets:
!TIMG!
!TIMG TABLE(22,J) = TABLE(22,J) + DCUMTM(135,J)
!TIMG!
!TIMG! --- triads:
!TIMG!
!TIMG TABLE(23,J) = TABLE(23,J) + DCUMTM(134,J)
!TIMG!
!TIMG! --- limiter:
!TIMG!
!TIMG TABLE(24,J) = TABLE(24,J) + DCUMTM(122,J)
!TIMG!
!TIMG! --- rescaling:
!TIMG!
!TIMG TABLE(25,J) = TABLE(25,J) + DCUMTM(121,J)
!TIMG!
!TIMG! --- reflections:
!TIMG!
!TIMG TABLE(26,J) = TABLE(26,J) + DCUMTM(136,J)
!TIMG!
!TIMG! --- diffraction:
!TIMG!
!TIMG TABLE(27,J) = TABLE(27,J) + DCUMTM(137,J)
!TIMG!
!TIMG! --- fluid mud:
!TIMG!
!TIMG TABLE(28,J) = TABLE(28,J) + DCUMTM(138,J)
!TIMG!
!TIMG! --- vegetation:
!TIMG!
!TIMG TABLE(29,J) = TABLE(29,J) + DCUMTM(139,J)
!TIMG!
!TIMG! --- turbulence:
!TIMG!
!TIMG TABLE(30,J) = TABLE(30,J) + DCUMTM(143,J)
!TIMG
!TIMG END DO
!TIMG!
!TIMG! --- add up times for some basic blocks
!TIMG!
!TIMG DO J = 1, 2
!TIMG!
!TIMG! --- total calculation:
!TIMG!
!TIMG TABLE(3,J) = TABLE(3,J) - TABLE( 7,J)
!TIMG TABLE(3,J) = TABLE(3,J) - TABLE(10,J)
!TIMG IF ( TABLE(3,J).LT.0D0 ) TABLE(3,J) = 0D0
!TIMG!
!TIMG! --- total communication:
!TIMG! * exchanging data
!TIMG! * global reductions
!TIMG! * collecting data
!TIMG!
!TIMG TABLE(4,J) = TABLE(4,J) + TABLE( 7,J)
!TIMG TABLE(4,J) = TABLE(4,J) + TABLE(10,J)
!TIMG TABLE(4,J) = TABLE(4,J) + TABLE(11,J)
!TIMG!
!TIMG! --- total propagation:
!TIMG! * velocities and derivatives
!TIMG!
!TIMG TABLE(6,J) = TABLE(6,J) + TABLE(14,J)
!TIMG TABLE(6,J) = TABLE(6,J) + TABLE(15,J)
!TIMG TABLE(6,J) = TABLE(6,J) + TABLE(16,J)
!TIMG TABLE(6,J) = TABLE(6,J) + TABLE(17,J)
!TIMG!
!TIMG! --- sources:
!TIMG! * wind, whitecapping, friction, breaking,
!TIMG! * quadruplets, triads, limiter, rescaling,
!TIMG! * reflections
!TIMG!
!TIMG DO K = 18, 26
!TIMG TABLE(8,J) = TABLE(8,J) + TABLE(K,J)
!TIMG END DO
!TIMG!
!TIMG! * diffraction
!TIMG!
!TIMG TABLE(8,J) = TABLE(8,J) + TABLE(27,J)
!TIMG!
!TIMG! * fluid mud
!TIMG!
!TIMG TABLE(8,J) = TABLE(8,J) + TABLE(28,J)
!TIMG!
!TIMG! * vegetation
!TIMG!
!TIMG TABLE(8,J) = TABLE(8,J) + TABLE(29,J)
!TIMG!
!TIMG! * turbulence
!TIMG!
!TIMG TABLE(8,J) = TABLE(8,J) + TABLE(30,J)
!TIMG!
!TIMG! --- other computing:
!TIMG!
!TIMG TABLE(13,J) = TABLE(13,J) + TABLE( 3,J)
!TIMG TABLE(13,J) = TABLE(13,J) - TABLE( 6,J)
!TIMG TABLE(13,J) = TABLE(13,J) - TABLE( 8,J)
!TIMG TABLE(13,J) = TABLE(13,J) - TABLE( 9,J)
!TIMG TABLE(13,J) = TABLE(13,J) - TABLE(12,J)
!TIMG IF ( TABLE(13,J).LT.0D0 ) TABLE(13,J) = 0D0
!TIMG
!TIMG END DO
!TIMG!
!TIMG! --- print CPU-times used in important parts of SWAN
!TIMG!
!TIMG WRITE(PRINTF,'(/)')
!TIMG WRITE(PRINTF,110) MYPRC
!TIMG WRITE(PRINTF,111) MYPRC
!TIMG WRITE(PRINTF,110) MYPRC
!TIMG WRITE(PRINTF,112) MYPRC
!TIMG WRITE(PRINTF,110) MYPRC
!TIMG WRITE(PRINTF,115) MYPRC,'total time:' ,(TABLE(1,J),J=1,2)
!TIMG WRITE(PRINTF,110) MYPRC
!TIMG WRITE(PRINTF,115) MYPRC,'total pre-processing:',
!TIMG & (TABLE(2,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'total calculation:',(TABLE(3,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'total communication:',
!TIMG & (TABLE(4,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'total post-processing:',
!TIMG & (TABLE(5,j),j=1,2)
!TIMG WRITE(PRINTF,110) MYPRC
!TIMG WRITE(PRINTF,113) MYPRC
!TIMG WRITE(PRINTF,110) MYPRC
!TIMG WRITE(PRINTF,115) MYPRC,'calc. propagation:',(TABLE(6,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'exchanging data:' ,(TABLE(7,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'calc. sources:' ,(TABLE(8,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'solving system:' ,(TABLE(9,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'reductions:' ,(TABLE(10,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'collecting data:' ,(TABLE(11,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'calc. setup:' ,(TABLE(12,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'other computing:' ,(TABLE(13,j),j=1,2)
!TIMG WRITE(PRINTF,110) MYPRC
!TIMG WRITE(PRINTF,114) MYPRC
!TIMG WRITE(PRINTF,110) MYPRC
!TIMG WRITE(PRINTF,115) MYPRC,'prop. velocities:',(TABLE(14,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'x-y advection:' ,(TABLE(15,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'sigma advection:' ,(TABLE(16,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'theta advection:' ,(TABLE(17,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'wind:' ,(TABLE(18,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'whitecapping:' ,(TABLE(19,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'bottom friction:' ,(TABLE(20,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'fluid mud:' ,(TABLE(28,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'vegetation:' ,(TABLE(29,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'turbulence:' ,(TABLE(30,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'wave breaking:' ,(TABLE(21,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'quadruplets:' ,(TABLE(22,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'triads:' ,(TABLE(23,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'limiter:' ,(TABLE(24,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'rescaling:' ,(TABLE(25,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'reflections:' ,(TABLE(26,j),j=1,2)
!TIMG WRITE(PRINTF,115) MYPRC,'diffraction:' ,(TABLE(27,j),j=1,2)
!TIMG
!TIMG END IF
!TIMG
!TIMG IF ( IDEBUG.GE.2 ) THEN
!TIMG WRITE(PRINTF,120) MYPRC
!TIMG DO J = 1, NSECTM
!TIMG IF (NCUMTM(J).GT.0)
!TIMG & WRITE(PRINTF,121) MYPRC,J,DCUMTM(J,1),DCUMTM(J,2),
!TIMG & NCUMTM(J)
!TIMG END DO
!TIMG END IF
!TIMG
!TIMG 110 FORMAT(i3,1x,'#')
!TIMG 111 FORMAT(i3,' # Details on timings of the simulation:')
!TIMG 112 FORMAT(i3,1x,'#',26x,'cpu-time',1x,'wall-clock')
!TIMG 113 FORMAT(i3,' # Splitting up calc. + comm. times:')
!TIMG 114 FORMAT(i3,' # Overview source contributions:')
!TIMG 115 FORMAT(i3,1x,'#',1x,a22,2f11.2)
!TIMG
!TIMG 120 FORMAT(/,i3,' # item cpu-time real time count')
!TIMG 121 FORMAT(i3,1x,'#',4x,i4,2f13.4,i10)
!TIMG
!TIMG RETURN
!TIMG END
!****************************************************************
!
SUBROUTINE TXPBLA(TEXT,IF,IL)
!
!****************************************************************
!
IMPLICIT NONE
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering and Geosciences |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmer: Marcel Zijlema |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 40.23: Marcel Zijlema
!
! 1. Updates
!
! 40.23, Feb. 03: New subroutine
!
! 2. Purpose
!
! determines the position of the first and the last non-blank
! (or non-tabulator) character in the text-string
!
! 4. Argument variables
!
! IF position of the first non-blank character in TEXT
! IL position of the last non-blank character in TEXT
! TEXT text string
!
INTEGER IF, IL
CHARACTER*(*) TEXT
!
! 6. Local variables
!
! FOUND : TEXT is found or not
! ITABVL: integer value of tabulator character
! LENTXT: length of TEXT
!
INTEGER LENTXT, ITABVL
LOGICAL FOUND
!
! 12. Structure
!
! Trivial.
!
! 13. Source text
!
!DOS ITABVL = ICHAR(' ')
ITABVL = ICHAR(' ')
LENTXT = LEN (TEXT)
IF = 1
FOUND = .FALSE.
100 IF (IF .LE. LENTXT .AND. .NOT. FOUND) THEN
IF (.NOT. (TEXT(IF:IF) .EQ. ' ' .OR.
& ICHAR(TEXT(IF:IF)) .EQ. ITABVL)) THEN
FOUND = .TRUE.
ELSE
IF = IF + 1
ENDIF
GOTO 100
ENDIF
IL = LENTXT + 1
FOUND = .FALSE.
200 IF (IL .GT. 1 .AND. .NOT. FOUND) THEN
IL = IL - 1
IF (.NOT. (TEXT(IL:IL) .EQ. ' ' .OR.
& ICHAR(TEXT(IL:IL)) .EQ. ITABVL)) THEN
FOUND = .TRUE.
ENDIF
GOTO 200
ENDIF
RETURN
END
!****************************************************************
!
CHARACTER*20 FUNCTION INTSTR ( IVAL )
!
!****************************************************************
!
IMPLICIT NONE
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering and Geosciences |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmer: Marcel Zijlema |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 40.23: Marcel Zijlema
!
! 1. Updates
!
! 40.23, Feb. 03: New subroutine
!
! 2. Purpose
!
! Convert integer to string
!
! 4. Argument variables
!
! IVAL integer to be converted
!
INTEGER IVAL
!
! 6. Local variables
!
! CVAL : character represented an integer of mantisse
! I : counter
! IPOS : position in mantisse
! IQUO : whole quotient
!
INTEGER I, IPOS, IQUO
CHARACTER*1, ALLOCATABLE :: CVAL(:)
!
! 12. Structure
!
! Trivial.
!
! 13. Source text
!
IPOS = 1
100 CONTINUE
IF (IVAL/10**IPOS.GE.1.) THEN
IPOS = IPOS + 1
GO TO 100
END IF
ALLOCATE(CVAL(IPOS))
DO I=IPOS,1,-1
IQUO=IVAL/10**(I-1)
CVAL(IPOS-I+1)=CHAR(INT(IQUO)+48)
IVAL=IVAL-IQUO*10**(I-1)
END DO
WRITE (INTSTR,*) (CVAL(I), I=1,IPOS)
RETURN
END
!****************************************************************
!
CHARACTER*20 FUNCTION NUMSTR ( IVAL, RVAL, FORM )
!
!****************************************************************
!
USE OCPCOMM4 40.41
!
IMPLICIT NONE
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering and Geosciences |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmer: Marcel Zijlema |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 40.23: Marcel Zijlema
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! 40.23, Feb. 03: New subroutine
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. Purpose
!
! Convert integer or real to string with given format
!
! 4. Argument variables
!
! IVAL integer to be converted
! FORM given format
! RVAL real to be converted
!
INTEGER IVAL
REAL RVAL
CHARACTER FORM*20
!
! 6. Local variables
!
! 12. Structure
!
! Trivial.
!
! 13. Source text
!
IF ( IVAL.NE.INAN ) THEN
WRITE (NUMSTR,FORM) IVAL
ELSE IF ( RVAL.NE.RNAN ) THEN
WRITE (NUMSTR,FORM) RVAL
ELSE
NUMSTR = ''
END IF
RETURN
END
!****************************************************************
!
SUBROUTINE SWCOPI ( IARR1, IARR2, LENGTH )
!
!****************************************************************
!
USE OCPCOMM4 40.41
!
IMPLICIT NONE
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering and Geosciences |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmer: Marcel Zijlema |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 40.23: Marcel Zijlema
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! 40.23, Feb. 03: New subroutine
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. Purpose
!
! Copies integer array IARR1 to IARR2
!
! 3. Method
!
! ---
!
! 4. Argument variables
!
! IARR1 source array
! IARR2 target array
! LENGTH array length
!
INTEGER LENGTH
INTEGER IARR1(LENGTH), IARR2(LENGTH)
!
! 6. Local variables
!
! I : loop counter
! IENT : number of entries
!
INTEGER I, IENT
!
! 8. Subroutines used
!
! MSGERR Writes error message
! STRACE Tracing routine for debugging
!
! 9. Subroutines calling
!
! ---
!
! 10. Error messages
!
! ---
!
! 11. Remarks
!
! ---
!
! 12. Structure
!
! Trivial.
!
! 13. Source text
!
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'SWCOPI')
! --- check array length
IF ( LENGTH.LE.0 ) THEN
CALL MSGERR( 3, 'Array length should be positive' )
END IF
! --- copy elements of array IARR1 to IARR2
DO 10 I = 1, LENGTH
IARR2(I) = IARR1(I)
10 CONTINUE
RETURN
END
!****************************************************************
!
SUBROUTINE SWCOPR ( ARR1, ARR2, LENGTH )
!
!****************************************************************
!
USE OCPCOMM4 40.41
!
IMPLICIT NONE
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering and Geosciences |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmer: Marcel Zijlema |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 40.23: Marcel Zijlema
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! 40.23, Feb. 03: New subroutine
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. Purpose
!
! Copies real array ARR1 to ARR2
!
! 3. Method
!
! ---
!
! 4. Argument variables
!
! ARR1 source array
! ARR2 target array
! LENGTH array length
!
INTEGER LENGTH
REAL ARR1(LENGTH), ARR2(LENGTH)
!
! 6. Local variables
!
! I : loop counter
! IENT : number of entries
!
INTEGER I, IENT
!
! 8. Subroutines used
!
! MSGERR Writes error message
! STRACE Tracing routine for debugging
!
! 9. Subroutines calling
!
! ---
!
! 10. Error messages
!
! ---
!
! 11. Remarks
!
! ---
!
! 12. Structure
!
! Trivial.
!
! 13. Source text
!
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'SWCOPR')
! --- check array length
IF ( LENGTH.LE.0 ) THEN
CALL MSGERR( 3, 'Array length should be positive' )
END IF
! --- copy elements of array ARR1 to ARR2
DO 10 I = 1, LENGTH
ARR2(I) = ARR1(I)
10 CONTINUE
RETURN
END
!MatL4!****************************************************************
!MatL4!
!MatL4 SUBROUTINE SWI2B ( IVAL, BVAL )
!MatL4!
!MatL4!****************************************************************
!MatL4!
!MatL4 USE OCPCOMM4 40.41
!MatL4!
!MatL4 IMPLICIT NONE
!MatL4!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering and Geosciences |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmer: Marcel Zijlema |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!MatL4!
!MatL4! 0. Authors
!MatL4!
!MatL4! 40.30: Marcel Zijlema
!MatL4! 40.41: Marcel Zijlema
!MatL4!
!MatL4! 1. Updates
!MatL4!
!MatL4! 40.30, May 03: New subroutine
!MatL4! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!MatL4!
!MatL4! 2. Purpose
!MatL4!
!MatL4! Calculates 32-bit representation of an integer number
!MatL4!
!MatL4! 3. Method
!MatL4!
!MatL4! The representation of an integer number is divided into 4 parts
!MatL4! of 8 bits each, resulting in 32-bit word in memory. Generally,
!MatL4! storage words are represented with bits counted from the right,
!MatL4! making bit 0 the lower-order bit and bit 31 the high-order bit,
!MatL4! which is also the sign bit.
!MatL4!
!MatL4! The integer number is always an exact representation of an
!MatL4! integer of value positive, negative, or zero. Each bit, except
!MatL4! the leftmost bit, corresponds to the actual exponent as power
!MatL4! of two.
!MatL4!
!MatL4! For representing negative numbers, the method called
!MatL4! "excess 2**(m - 1)" is used, which represents an m-bit number by
!MatL4! storing it as the sum of itself and 2**(m - 1). For a 32-bit
!MatL4! machine, m = 32. This results in a positive number, so the
!MatL4! leftmost bit need to be reversed. This method is identical to the
!MatL4! two's complement method.
!MatL4!
!MatL4! An example:
!MatL4!
!MatL4! the 32-bit representation of 5693 is
!MatL4!
!MatL4! decimal : 0 0 22 61
!MatL4! hexidecimal: 0 0 16 3D
!MatL4! binary : 00000000 00000000 00010110 00111101
!MatL4!
!MatL4! since,
!MatL4!
!MatL4! 5693 = 2^12 + 2^10 + 2^9 + 2^5 + 2^4 + 2^3 + 2^2 + 2^0
!MatL4!
!MatL4! 4. Argument variables
!MatL4!
!MatL4! BVAL a byte value as a part of the representation of
!MatL4! integer number
!MatL4! IVAL integer number
!MatL4!
!MatL4 INTEGER BVAL(4), IVAL
!MatL4!
!MatL4! 6. Local variables
!MatL4!
!MatL4! I : loop counter
!MatL4! IENT : number of entries
!MatL4! IQUOT : auxiliary integer with quotient
!MatL4! M : maximal exponent number possible (for 32-bit machine, m=32)
!MatL4!
!MatL4 INTEGER I, IENT, IQUOT, M
!MatL4!
!MatL4! 12. Structure
!MatL4!
!MatL4! initialise 4 parts of the representation
!MatL4! if integer < 0, increased it by 2**(m-1)
!MatL4! compute the actual part of the representation
!MatL4! determine the sign bit
!MatL4!
!MatL4! 13. Source text
!MatL4!
!MatL4 SAVE IENT
!MatL4 DATA IENT/0/, M/32/
!MatL4 IF (LTRACE) CALL STRACE (IENT,'SWI2B')
!MatL4
!MatL4! --- initialise 4 parts of the representation
!MatL4
!MatL4 DO 10 I = 1, 4
!MatL4 BVAL(I) = 0
!MatL4 10 CONTINUE
!MatL4
!MatL4 IQUOT = IVAL
!MatL4
!MatL4! --- if integer < 0, increased it by 2*(m-1)
!MatL4
!MatL4 IF ( IVAL.LT.0 ) IQUOT = IQUOT + 2**(M-1)
!MatL4
!MatL4! --- compute the actual part of the representation
!MatL4
!MatL4 DO 20 I = 4, 1, -1
!MatL4 BVAL(I) = MOD(IQUOT,256)
!MatL4 IQUOT = INT(IQUOT/256)
!MatL4 20 CONTINUE
!MatL4
!MatL4! --- determine the sign bit
!MatL4
!MatL4 IF ( IVAL.LT.0 ) BVAL(1) = BVAL(1) + 128
!MatL4
!MatL4 RETURN
!MatL4 END
!MatL4!****************************************************************
!MatL4!
!MatL4 SUBROUTINE SWR2B ( RVAL, BVAL )
!MatL4!
!MatL4!****************************************************************
!MatL4!
!MatL4 USE OCPCOMM4 40.41
!MatL4!
!MatL4 IMPLICIT NONE
!MatL4!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering and Geosciences |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmer: Marcel Zijlema |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!MatL4!
!MatL4! 0. Authors
!MatL4!
!MatL4! 40.30: Marcel Zijlema
!MatL4! 40.41: Marcel Zijlema
!MatL4!
!MatL4! 1. Updates
!MatL4!
!MatL4! 40.30, May 03: New subroutine
!MatL4! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!MatL4!
!MatL4! 2. Purpose
!MatL4!
!MatL4! Calculates 32-bit representation of a floating-point number
!MatL4!
!MatL4! 3. Method
!MatL4!
!MatL4! The representation of a floating-point number is divided into 4
!MatL4! parts of 8 bits each, resulting in 32-bit word in memory.
!MatL4! Generally, storage words are represented with bits counted from
!MatL4! the right, making bit 0 the lower-order bit and bit 31 the
!MatL4! high-order bit, which is also the sign bit.
!MatL4!
!MatL4! The floating-point number is a processor approximation. Its format
!MatL4! has an 8-bit biased exponent and a 23-bit fraction or mantissa. The
!MatL4! leftmost bit is the sign bit which is zero for plus and 1 for
!MatL4! minus. The biased exponent equals the bias and the actual exponent
!MatL4! (power of two) of the number. For a 32-bit machine, bias=127.
!MatL4!
!MatL4! Furthermore, the floating-point number is usually stored in the
!MatL4! normalized form, i.e. it has a binary point to the left of the
!MatL4! mantissa and an implied leading 1 to the left of the binary point.
!MatL4! Thus, if X is a floating-point number, then it is calculated as
!MatL4! follows:
!MatL4!
!MatL4! X = (-1)**sign bit 1.fraction * 2**(biased exponent-bias)
!MatL4!
!MatL4! There are several exceptions. Let a fraction, biased exponent
!MatL4! and sign bit be denoted as F, E and S, respectively. The following
!MatL4! formats adhere to IEEE standard:
!MatL4!
!MatL4! S = 0, E = 00000000 and F = 00 ... 0 : X = 0
!MatL4! S = 0, E = 00000000 and F <> 00 ... 0 : X = +0.fraction * 2**(1-bias)
!MatL4! S = 1, E = 00000000 and F <> 00 ... 0 : X = -0.fraction * 2**(1-bias)
!MatL4! S = 0, E = 11111111 and F = 00 ... 0 : X = +Inf
!MatL4! S = 1, E = 11111111 and F = 00 ... 0 : X = -Inf
!MatL4! S = 0, E = 11111111 and F <> 00 ... 0 : X = NaN
!MatL4!
!MatL4! A NaN (Not a Number) is a value reserved for signalling an
!MatL4! attempted invalid operation, like 0/0. Its representation
!MatL4! equals the representation of +Inf plus 1, i.e. 2**31 - 2**23 + 1
!MatL4!
!MatL4! An example:
!MatL4!
!MatL4! the 32-bit representation of 23.1 is
!MatL4!
!MatL4! decimal : 65 184 204 205
!MatL4! hexidecimal: 41 B8 CC CD
!MatL4! binary : 01000001 10111000 11001100 11001101
!MatL4!
!MatL4! since,
!MatL4!
!MatL4! 23.1 = 2^4 + 2^2 + 2^1 + 2^0 + 2^-4 + 2^-5 + 2^-8 + 2^-9 +
!MatL4! 2^-12 + 2^-13 + 2^-16 + 2^-17 + 2^-19
!MatL4!
!MatL4! so that the biased exponent = 4 + 127 = 131 = 10000011 = E
!MatL4! and the sign bit = 0 = S. The remaining of the 32-bit word is
!MatL4! the fraction, which is
!MatL4!
!MatL4! 3 2 1 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -10 -11 -12 -13 -14 -15 -16 -17 -18 -19
!MatL4!
!MatL4! F= 0 1 1 1 0 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 1
!MatL4!
!MatL4! 4. Argument variables
!MatL4!
!MatL4! BVAL a byte value as a part of the representation of
!MatL4! floating-point number
!MatL4! RVAL floating-point number
!MatL4!
!MatL4 INTEGER BVAL(4)
!MatL4 REAL RVAL
!MatL4!
!MatL4! 6. Local variables
!MatL4!
!MatL4! ACTEXP: actual exponent in the representation
!MatL4! BEXPO : biased exponent
!MatL4! BIAS : bias (for 32-bit machine, bias=127)
!MatL4! EXPO : calculated exponent of floating-point number
!MatL4! FRAC : fraction of floating-point number
!MatL4! I : loop counter
!MatL4! IENT : number of entries
!MatL4! IPART : i-the part of the representation
!MatL4! IQUOT : auxiliary integer with quotient
!MatL4! LEADNR: leading number of floating-point number
!MatL4! LFRAC : length of fraction in representation
!MatL4! : (for 32-bit machine, lfrac=23)
!MatL4! RFRAC : auxiliary real with fraction
!MatL4!
!MatL4 INTEGER ACTEXP, EXPO, I, IENT, IPART, LFRAC, BIAS, BEXPO
!MatL4 INTEGER*8 LEADNR, IQUOT
!MatL4 REAL FRAC, RFRAC, MAXREAL
!MatL4!
!MatL4! 12. Structure
!MatL4!
!MatL4! initialise 4 parts of the representation and biased exponent
!MatL4! determine leading number and fraction
!MatL4! do while leading number >= 1
!MatL4! calculate positive exponent as power of two
!MatL4! or while fraction > 0
!MatL4! calculate negative exponent as power of two
!MatL4! end do
!MatL4! compute the actual part of the representation
!MatL4!
!MatL4! 13. Source text
!MatL4!
!MatL4 SAVE IENT
!MatL4 DATA IENT/0/, LFRAC/23/, BIAS/127/
!MatL4 IF (LTRACE) CALL STRACE (IENT,'SWR2B')
!MatL4
!MatL4! --- initialise 4 parts of the representation and biased exponent
!MatL4
!MatL4 DO 10 I = 1, 4
!MatL4 BVAL(I) = 0
!MatL4 10 CONTINUE
!MatL4 BEXPO = -1
!MatL4
!MatL4! --- clamp reals to the max integer*8, to prevent overflow
!MatL4
!MatL4 MAXREAL=9.0E+18
!MatL4 IF ( RVAL .GT. MAXREAL ) THEN
!MatL4 RVAL=MAXREAL
!MatL4 ELSEIF ( RVAL .LT. -MAXREAL ) THEN
!MatL4 RVAL=-MAXREAL
!MatL4 END IF
!MatL4
!MatL4! --- determine leading number and fraction
!MatL4
!MatL4 IF ( ABS(RVAL).LT.1.E-7 ) THEN
!MatL4 LEADNR = 0
!MatL4 FRAC = 0.
!MatL4 ELSE
!MatL4 LEADNR = INT(ABS(RVAL),KIND=8)
!MatL4 FRAC = ABS(RVAL) - REAL(LEADNR)
!MatL4 END IF
!MatL4
!MatL4 20 IF ( LEADNR.GE.1 ) THEN
!MatL4
!MatL4! --- calculate positive exponent as power of two
!MatL4
!MatL4 EXPO = 0
!MatL4 IQUOT = LEADNR
!MatL4 30 IF ( IQUOT.GE.2 ) THEN
!MatL4
!MatL4 IQUOT = INT(IQUOT/2,KIND=8)
!MatL4 EXPO = EXPO + 1
!MatL4
!MatL4 GOTO 30
!MatL4 END IF
!MatL4
!MatL4 ELSE IF ( FRAC.GT.0. ) THEN
!MatL4
!MatL4! --- calculate negative exponent as power of two
!MatL4
!MatL4 EXPO = 0
!MatL4 RFRAC = FRAC
!MatL4 40 IF ( RFRAC.LT.1. ) THEN
!MatL4
!MatL4 RFRAC = RFRAC * 2.
!MatL4 EXPO = EXPO - 1
!MatL4
!MatL4 GOTO 40
!MatL4 END IF
!MatL4
!MatL4 ELSE
!MatL4
!MatL4 GOTO 50
!MatL4
!MatL4 END IF
!MatL4
!MatL4! --- compute the actual part of the representation
!MatL4
!MatL4 IF ( BEXPO.EQ.-1 ) THEN
!MatL4
!MatL4! --- determine biased exponent
!MatL4
!MatL4 BEXPO = EXPO + BIAS
!MatL4
!MatL4! --- the first seven bits of biased exponent belong
!MatL4! to first part of the representation
!MatL4
!MatL4 BVAL(1) = INT(BEXPO/2)
!MatL4
!MatL4! --- determine the sign bit
!MatL4
!MatL4 IF ( RVAL.LT.0. ) BVAL(1) = BVAL(1) + 128
!MatL4
!MatL4! --- the eighth bit of biased component is the leftmost
!MatL4! bit of second part of the representation
!MatL4
!MatL4 BVAL(2) = MOD(BEXPO,2)*2**7
!MatL4 IPART = 2
!MatL4
!MatL4 ELSE
!MatL4
!MatL4! --- compute the actual exponent of bit 1 in i-th part of
!MatL4! the representation
!MatL4
!MatL4 ACTEXP = (IPART-2)*8 + 7 - BEXPO + BIAS + EXPO
!MatL4 IF ( ACTEXP.LT.0 ) THEN
!MatL4 ACTEXP = ACTEXP + 8
!MatL4 IPART = IPART + 1
!MatL4 IF ( IPART.GT.4 ) GOTO 50
!MatL4 END IF
!MatL4 BVAL(IPART) = BVAL(IPART) + 2**ACTEXP
!MatL4
!MatL4 END IF
!MatL4
!MatL4 IF ( EXPO.LT.(BEXPO-BIAS-LFRAC) ) GOTO 50
!MatL4 LEADNR = LEADNR - 2.**EXPO
!MatL4 IF ( EXPO.LT.0 ) FRAC = FRAC - 2.**EXPO
!MatL4
!MatL4 GOTO 20
!MatL4 50 CONTINUE
!MatL4
!MatL4 RETURN
!MatL4 END