FVCOMin/doxygen/swanser_8f90_source.html

 !
 !     SWAN - SERVICE ROUTINES
 !
 !  Contents of this file
 !
 !     READXY
 !     REFIXY
 !     INFRAM
 !     DISTR
 !     KSCIP1
 !     AC2TST
 !     CVCHEK                                                              30.60
 !     CVMESH                                                              30.60
 !     NEWTON                                                              30.60
 !     EVALF                                                               30.60
 !     SWOBST                                                              30.60
 !     TCROSS                                                              40.04
 !     SWTRCF
 !     SSHAPE                                                              40.00
 !     SINTRP                                                              40.00
 !     HSOBND                                                              32.01
 !     CHGBAS                                                              40.00
 !     GAMMA                                                               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
 !     SWI2B                                                               40.30
 !     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
    USE swcomm2
 !
 !
 !   --|-----------------------------------------------------------|--
 !     | Delft University of Technology                            |
 !     | Faculty of Civil Engineering                              |
 !     | Environmental Fluid Mechanics Section                     |
 !     | P.O. Box 5048, 2600 GA  Delft, The Netherlands            |
 !     |                                                           |
 !     | Programmers: R.C. Ris, N. Booij,                          |
 !     |              IJ.G. Haagsma, A.T.M.M. Kieftenburg,         |
 !     |              M. Zijlema, E.E. Kriezi,                     |
 !     |              R. Padilla-Hernandez, L.H. Holthuijsen       |
 !   --|-----------------------------------------------------------|--
 !
 !
 !     SWAN (Simulating WAves Nearshore); a third generation wave model
 !     Copyright (C) 2004-2005  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')
 !
 !JQI   CALL INDBLE (NAMX, XTMP, KONT, DBLE(XSTA)+DBLE(XOFFS))
    IF(chgval)THEN
 !       a valid value was given for XX
 !JQI     CALL INDBLE (NAMY, YTMP, 'REQ', DBLE(YSTA)+DBLE(YOFFS))
    ELSE
 !JQI     CALL INDBLE (NAMY, YTMP, KONT, DBLE(YSTA)+DBLE(YOFFS))
    ENDIF
    IF(.NOT.lxoffs)THEN
      xoffs = real(xtmp)
      yoffs = real(ytmp)
      lxoffs = .true.
    ENDIF
 !JQI   IF(.NOT.EQREAL(XOFFS,REAL(XTMP)))THEN
 !JQI     XX = REAL(XTMP-DBLE(XOFFS))
 !JQI   ELSE IF(OPTG == 3)THEN
 !JQI     XX = 1.E-5
 !JQI   ELSE
 !JQI     XX = 0.
 !JQI   END IF
 !JQI   IF(.NOT.EQREAL(YOFFS,REAL(YTMP)))THEN
 !JQI     YY = REAL(YTMP-DBLE(YOFFS))
 !JQI   ELSE IF(OPTG == 3)THEN
 !JQI     YY = 1.E-5
 !JQI   ELSE
 !JQI     YY = 0.
 !JQI   END IF
 !
    RETURN
    END SUBROUTINE readxy


 !****************************************************************
 !
    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: M. Zijlema                                    |
 !   --|-----------------------------------------------------------|--
 !
 !
 !     SWAN (Simulating WAves Nearshore); a third generation wave model
 !     Copyright (C) 2004-2005  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('  ')
 !UNIX      ITABVL = ICHAR(' ')
    lentxt = len(text)
    IF = 1
    found = .false.
 100 IF(IF <= lentxt .AND. .NOT. found)THEN
      IF(.NOT. (text(if:if) ==  ' ' .OR.                   &
        ichar(text(if:if)) == itabvl))THEN
        found = .true.
      ELSE
        IF = IF + 1
      ENDIF
      GOTO 100
    ENDIF
    il = lentxt + 1
    found = .false.
 200 IF(il > 1 .AND. .NOT. found)THEN
     il = il - 1
      IF(.NOT. (text(il:il) ==  ' ' .OR.                   &
        ichar(text(il:il)) == itabvl))THEN
        found = .true.
      ENDIF
      GOTO 200
    ENDIF

    RETURN
    END SUBROUTINE txpbla

 !****************************************************************
 !
     CHARACTER*20 FUNCTION numstr ( IVAL, RVAL, FORM )
 !
 !****************************************************************
 !
     USE ocpcomm4
 !
     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: M. Zijlema                                    |
 !   --|-----------------------------------------------------------|--
 !
 !
 !     SWAN (Simulating WAves Nearshore); a third generation wave model
 !     Copyright (C) 2004-2005  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 /= inan)THEN
       WRITE(numstr,form) ival
     ELSE IF( rval /= rnan)THEN
       WRITE (numstr,form) rval
     ELSE
       numstr = ''
     END IF

     RETURN
     END FUNCTION numstr
 !****************************************************************
 !
    SUBROUTINE swcopi ( IARR1, IARR2, LENGTH )
 !
 !****************************************************************
 !
    USE ocpcomm4
 !
    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: M. Zijlema                                    |
 !   --|-----------------------------------------------------------|--
 !
 !
 !     SWAN (Simulating WAves Nearshore); a third generation wave model
 !     Copyright (C) 2004-2005  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 <= 0)THEN
      CALL msgerr( 3, 'Array length should be positive' )
    END IF

 !  --- copy elements of array IARR1 to IARR2

    DO i = 1, length
      iarr2(i) = iarr1(i)
    END DO

    RETURN
    END SUBROUTINE swcopi


 !***********************************************************************
 !                                                                      *
     SUBROUTINE kscip1(MMT,SIG,D,K,CG,N,ND)
 !                                                                      *
 !***********************************************************************
 !
     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: R.C. Ris, N. Booij,                          |
 !     |              IJ.G. Haagsma, A.T.M.M. Kieftenburg,         |
 !     |              M. Zijlema, E.E. Kriezi,                     |
 !     |              R. Padilla-Hernandez, L.H. Holthuijsen       |
 !   --|-----------------------------------------------------------|--
 !
 !  0. Authors
 !
 !  1. Updates
 !
 !  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_w)
     wgd    = rootdg*grav_w
     DO is = 1, mmt
 !     SND is dimensionless frequency
       snd = sig(is) * rootdg
       IF(snd >= 2.5)THEN
 !     ******* deep water *******
         k(is)  = sig(is) * sig(is) / grav_w
         cg(is) = 0.5 * grav_w / sig(is)
         n(is)  = 0.5
         nd(is) = 0.
       ELSE IF(snd < 1.e-6)THEN
 !     *** very shallow water ***
 !        print*,'IN VERY SHALLOW WATER'
 !   stop
         k(is)  = snd/d
         cg(is) = wgd
         n(is)  = 1.
         nd(is) = 0.
       ELSE
         snd2  = snd*snd
         c     = sqrt(grav_w*d/(snd2+1./(1.+0.666*snd2+0.445*snd2**2      &
             -0.105*snd2**3+0.272*snd2**4)))
         k(is) = sig(is)/c
         knd   = k(is)*d
         fac1  = 2.*knd/sinh(2.*knd)
         n(is) = 0.5*(1.+fac1)
         cg(is)= n(is)*c
         fac2  = snd2/knd
         fac3  = 2.*fac2/(1.+fac2*fac2)
         nd(is)= fac1*(0.5/d - k(is)/fac3)
       ENDIF
     END DO

     RETURN
     END SUBROUTINE kscip1
 !***********************************************************************
 !                                                                      *
     SUBROUTINE kscip2(S,D,CG)
 !                                                                      *
 !***********************************************************************
 !
     USE ocpcomm4
     USE swcomm3
 !
     IMPLICIT NONE
       REAL      CG(MSC), D, K,S(MSC)
       INTEGER   IENT, IS
       REAL      KND, ROOTDG, SND, WGD, SND2, C, FAC1,N
     rootdg = sqrt(d/grav_w)
     wgd    = rootdg*grav_w
 !     SND is dimensionless frequency
      DO is=1,msc
       snd = s(is) * rootdg
       IF(snd >= 2.5)THEN
 !     ******* deep water *******
         k = s(is) * s(is) / grav_w
         cg(is) = 0.5 * grav_w / s(is)
       ELSE IF(snd < 1.e-6)THEN
 !     *** very shallow water ***
 !        print*,'IN VERY SHALLOW WATER'
 !       stop
         k  = snd/d
         cg(is) = wgd
       ELSE
         snd2  = snd*snd
         c     = sqrt(grav_w*d/(snd2+1./(1.+0.666*snd2+0.445*snd2**2      &
                 -0.105*snd2**3+0.272*snd2**4)))
         k = s(is)/c
         knd   = k*d
         fac1  = 2.*knd/sinh(2.*knd)
         n = 0.5*(1.+fac1)
         cg(is)= n*c
       ENDIF
     END DO
     RETURN
     END SUBROUTINE kscip2


 !****************************************************************

 !****************************************************************
 !
       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: M. Zijlema                                    |
 !   --|-----------------------------------------------------------|--
 !
 !
 !     SWAN (Simulating WAves Nearshore); a third generation wave model
 !     Copyright (C) 2004-2005  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 FUNCTION intstr

 !
 !********************************************************************
 !                                                                   *
       REAL function gamma(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
       SAVE ient, abig
       DATA ient /0/, abig /30./
       CALL strace (ient, 'GAMMA')
       yy = gammln(xx)
       IF (yy.GT.abig) yy = abig
       IF (yy.LT.-abig) yy = -abig
       gamma = exp(yy)
       RETURN
       END FUNCTION gamma

 !********************************************************************
 !                                                                   *
       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 FUNCTION gammln

 !
 !***********************************************************************
 !                                                                      *
       SUBROUTINE sshape (ACLOC, SPCSIG, SPCDIR, FSHAPL, DSHAPL)
 !                                                                      *
 !***********************************************************************
 !
       USE ocpcomm4
       USE swcomm3
 !
 !  0. Authors
 !
 !  1. Updates
 !
 !  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)
       REAL    SPCSIG(MSC)
 !
 ! 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
 !
 !  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
 !
 !     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 >= 80) WRITE (prtest, 8) fshapl, dshapl,          &
                                          (spparm(jj), jj = 1,4)
    8  FORMAT (' entry SSHAPE ', 2i3, 4e12.4)
       IF(fshapl < 0)THEN
         lshape = - fshapl
         logpm  = .false.
       ELSE
         lshape = fshapl
         logpm  = .true.
       ENDIF

 !
       IF(spparm(1) <= 0.)                                         &
          CALL msgerr(1,' sign. wave height at boundary is not positive')
 !
       pkper = spparm(2)
       itper = 0
       IF(lshape == 3)THEN
 !       select bin closest to given period
         difper = 1.e10
         DO is = 1, msc
           IF(abs(pkper - pi2_w/spcsig(is)) < difper)THEN
             isp = is
             difper = abs(pkper - pi2_w/spcsig(is))
           ENDIF
         ENDDO
       ENDIF
 !
 !     compute spectral shape using peak period PKPER
 !
       fac  = 1.
  100  fpk  = (1./pkper)
       fpk4 = fpk**4
       IF(lshape == 1)THEN
         salpha = ((spparm(1) ** 2) * (fpk4)) * 5. / 16.
       ELSE IF(lshape == 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 == 4)THEN
         aux1 = spparm(1)**2 / ( 16.* sqrt(pi2_w) * pshape(2))
         aux3 = 2. * pshape(2)**2
       ENDIF
 !
       ctott = 0.
       DO is = 1, msc
 !
         IF(lshape == 1)THEN
 !         *** LSHAPE = 1 : Pierson and Moskowitz ***
           sf = spcsig(is) / pi2_w
           sf4 = sf**4
           sf5 = sf**5
           ra = (salpha/sf5)*exp(-(5.*fpk4)/(4.*sf4))/(pi2_w*spcsig(is))
           acloc(mdc,is) = ra
         ELSE IF(lshape == 2)THEN
 !         *** LSHAPE = 2 : JONSWAP ***
           sf = spcsig(is)/(pi2_w)
           sf4 = sf**4
           sf5 = sf**5
           cpshap = 1.25 * fpk4 / sf4
           IF(cpshap > 10.)THEN
             ra = 0.
           ELSE
             ra = (salpha/sf5) * exp(-cpshap)
           ENDIF
           IF(sf < fpk)THEN
             coeff = 0.07
           ELSE
             coeff = 0.09
           ENDIF
           apshap =  0.5 * ((sf-fpk) / (coeff*fpk)) **2
           IF(apshap > 10.)THEN
             syf = 1.
           ELSE
             ppshap = exp(-apshap)
             syf = pshape(1)**ppshap
           ENDIF
           ra = syf*ra/(spcsig(is)*pi2_w)
           acloc(mdc,is) = ra
           IF(itest >= 120) WRITE (prtest, 112)                  &
                 sf, salpha, cpshap, apshap, syf, ra
  112      FORMAT (' SSHAPE freq. ', 8e12.4)
         ELSE IF(lshape == 3)THEN
 !
 !         *** all energy concentrated in one BIN ***
 !
           IF(is == isp)THEN
             acloc(mdc,is) = ( spparm(1)**2 ) /                   &
                         ( 16. * spcsig(is)**2 * frintf )
           ELSE
             acloc(mdc,is) = 0.
           ENDIF
         ELSE IF(lshape == 4)THEN
 !
 !         *** energy Gaussian distributed (wave-current tests) ***
 !
           aux2 = ( spcsig(is) - ( pi2_w / pkper ) )**2
           ra = aux1 * exp( -1. * aux2 / aux3 ) / spcsig(is)
           acloc(mdc,is) = ra
         ELSE
           IF(is == 1)THEN
             CALL msgerr (2,'Wrong type for frequency shape')
             WRITE (printf, *) ' -> ', fshapl, lshape
           ENDIF
         ENDIF
         IF (itest >= 10)                                         &
             ctott = ctott + frintf * acloc(mdc,is) * spcsig(is)**2
       END DO   !END DO IS
       IF(itest >= 10)THEN
         IF(spparm(1) > 0.01)THEN
           hstmp = 4. * sqrt(ctott)
           IF(abs(hstmp-spparm(1)) > 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 < 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.
         aptail = 1. / (pptail * (1. + pptail * (frinth-1.)))
         am0 = am0 * frintf + aptail * as2
         pptail = pwtail(1) - 2.
         eptail = 1. / (pptail * (1. + pptail * (frinth-1.)))
         am1 = am1 * frintf + eptail * as3
 !       Mean period:
         IF( am1 /= 0.)THEN
            mper = pi2_w * am0 / am1
         ELSE
            CALL msgerr(3, ' first moment is zero in calculating the')
            CALL msgerr(3, ' spectrum at boundary using param. bc.')
         END IF
         IF(itest >= 80) WRITE (prtest, 72) itper, spparm(2), mper,     &
                                         pkper
   72    FORMAT (' SSHAPE iter=', i2, '  period values:', 3f7.2)
         IF(abs(mper-spparm(2)) > 0.01*spparm(2))THEN
 !         modification suggested by Mauro Sclavo
           pkper = (spparm(2) / mper) * pkper
           GOTO 100
         ENDIF
       ENDIF
 !
       IF (itper >= 10)THEN
         CALL msgerr(3, 'No convergence calculating the spectrum')
         CALL msgerr(3, 'at the boundary using parametric bound. cond.')
       ENDIF
 !
 !     now introduce distribution over directions
 !
       adir = pi_w * degcnv(spparm(3)) / 180.
       IF(dshapl == 1)THEN
         dspr = pi_w * spparm(4) / 180.
         ms = max(dspr**(-2) - 2., 1.)
       ELSE
         ms = spparm(4)
       ENDIF
       IF(ms < 12.)THEN
         ctot = (2.**ms) * (gamma(0.5*ms+1.))**2 / (pi_w * gamma(ms+1.))
       ELSE
         ctot =  sqrt(0.5*ms/pi_w) / (1. - 0.25/ms)
       ENDIF
       IF(itest >= 100)THEN
         esom = 0.
         DO is = 1, msc
           esom = esom + frintf * spcsig(is)**2 * acloc(mdc,is)
         ENDDO
         WRITE (prtest, *) ' SSHAPE dir ', 4.*sqrt(abs(esom)),         &
            spparm(1), ctot, ms, gamma(0.5*ms+1.), gamma(ms+1.),   &
            ctot
       ENDIF
       dverif = .false.
       ctott = 0.
       DO id = 1, mdc
         acos = cos(spcdir(id,1) - adir)
         IF(acos > 0.)THEN
           cdir = ctot * max(acos**ms, 1.e-10)
           IF(.NOT.fulcir)THEN
             IF(acos >= cos(ddir)) dverif = .true.
           ENDIF
         ELSE
           cdir = 0.
         ENDIF
         IF(itest >= 10) ctott = ctott + cdir * ddir
         IF(itest >= 100) WRITE (prtest, 360) id,spcdir(id,1),cdir
  360    FORMAT (' ID Spcdir Cdir: ',i3,3(1x,e11.4))
         DO is = 1, msc
           acloc(id,is) = cdir * acloc(mdc,is)
         ENDDO
       ENDDO
       IF(itest >= 10)THEN
         IF(abs(ctott-1.) > 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 SUBROUTINE sshape

 !***********************************************************************
 !                                                                      *
       REAL FUNCTION DEGCNV (DEGREE)
 !                                                                      *
 !***********************************************************************
 !
       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: R.C. Ris, N. Booij,                          |
 !     |              IJ.G. Haagsma, A.T.M.M. Kieftenburg,         |
 !     |              M. Zijlema, E.E. Kriezi,                     |
 !     |              R. Padilla-Hernandez, L.H. Holthuijsen       |
 !   --|-----------------------------------------------------------|--
 !
 !
 !     SWAN (Simulating WAves Nearshore); a third generation wave model
 !     Copyright (C) 2004-2005  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
 !
 !***********************************************************************
 !
       INTEGER :: ient
       REAL    :: degree

       SAVE ient
       DATA ient /0/
       CALL strace(ient,'DEGCNV')
 !
       IF ( bnaut ) THEN
         degcnv = 180. + dnorth - degree
       ELSE
         degcnv = degree
       ENDIF
 !
       IF (degcnv >= 360.) THEN
         degcnv = mod(degcnv, 360.)
       ELSE IF (degcnv < 0.) THEN
         degcnv = mod(degcnv, 360.) + 360.
       ELSE
 !       DEGCNV between 0 and 360; do nothing
       ENDIF
 !
       RETURN
       END FUNCTION degcnv

 !
 !***********************************************************************
 !                                                                      *
       SUBROUTINE swanout(DEP2)    !AC2   ,SPCSIG,HSIBC,KGRPNT)
       RETURN
       END SUBROUTINE swanout
mod_time::abs
Definition: mod_time.f90:64

strace
subroutine strace(IENT, SUBNAM)
Definition: ocpmix.f90:468

swcomm3::bnaut
logical bnaut
Definition: swmod1.f90:2144

ocpcomm4::ltrace
logical ltrace
Definition: swmod1.f90:537

readxy
subroutine readxy(NAMX, NAMY, XX, YY, KONT, XSTA, YSTA)
Definition: swanser.f90:63

intstr
character *20 function intstr(IVAL)
Definition: swanser.f90:696

ocpcomm4::prtest
integer prtest
Definition: swmod1.f90:517

swcomm3::grav_w
real grav_w
Definition: swmod1.f90:1721

swcomm3::pshape
real, dimension(mshape) pshape
Definition: swmod1.f90:2134

ocpcomm4::rnan
real rnan
Definition: swmod1.f90:485

swcomm3::spparm
real, dimension(msppar) spparm
Definition: swmod1.f90:2138

swcomm3::pi2_w
real pi2_w
Definition: swmod1.f90:1722

numstr
character *20 function numstr(IVAL, RVAL, FORM)
Definition: swanser.f90:305

swcomm3::dnorth
real dnorth
Definition: swmod1.f90:1721

mod_time::mod
Definition: mod_time.f90:69

swcomm3::fulcir
logical fulcir
Definition: swmod1.f90:1683

swcomm3::pi_w
real pi_w
Definition: swmod1.f90:1721

ocpcomm4::printf
integer printf
Definition: swmod1.f90:517

swcomm3::frintf
real frintf
Definition: swmod1.f90:1678

swcomm2::yoffs
real yoffs
Definition: swmod1.f90:1364

ocpcomm4::itest
integer itest
Definition: swmod1.f90:536

sshape
subroutine sshape(ACLOC, SPCSIG, SPCDIR, FSHAPL, DSHAPL)
Definition: swanser.f90:840

gamma
real function gamma(XX)
Definition: swanser.f90:788

ocpcomm1
Definition: swmod1.f90:26

swcomm2::lxoffs
logical lxoffs
Definition: swmod1.f90:1365

degcnv
real function degcnv(DEGREE)
Definition: swanser.f90:1225

swcomm3::frinth
real frinth
Definition: swmod1.f90:1678

kscip1
subroutine kscip1(MMT, SIG, D, K, CG, N, ND)
Definition: swanser.f90:504

msgerr
subroutine msgerr(LEV, STRING)
Definition: ocpmix.f90:582

gammln
function gammln(XX)
Definition: swanser.f90:814

swcomm3::pwtail
real, dimension(10) pwtail
Definition: swmod1.f90:1722

ocpcomm4::inan
integer inan
Definition: swmod1.f90:484

ocpcomm1::chgval
logical chgval
Definition: swmod1.f90:154

swcomm2
Definition: swmod1.f90:1216

swanout
subroutine swanout(DEP2)
Definition: swanser.f90:1356

txpbla
subroutine txpbla(TEXT, IF, IL)
Definition: swanser.f90:202

swcomm2::xoffs
real xoffs
Definition: swmod1.f90:1363

swcopi
subroutine swcopi(IARR1, IARR2, LENGTH)
Definition: swanser.f90:387

swcomm3
Definition: swmod1.f90:1417

kscip2
subroutine kscip2(S, D, CG)
Definition: swanser.f90:650

ocpcomm4
Definition: swmod1.f90:418

swcomm3::ddir
real ddir
Definition: swmod1.f90:1676