#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 (XP1XP) or (XP1>XP and XP2YP ! 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