WW3/w3sbt4md_8F90_source.html

#include "w3macros.h"

!/ ------------------------------------------------------------------- /

MODULE w3sbt4md

  !/

  !/                  +-----------------------------------+

  !/                  | WAVEWATCH III           NOAA/NCEP |

  !/                  |    F. Ardhuin and J. Lepesqueur   |

  !/                  |                        FORTRAN 90 |

  !/                  | Last update :         14-Mar-2012 |

  !/                  +-----------------------------------+

  !/

  !/    20-Dec-2004 : Origination.                        ( version 3.06 )

  !/    23-Jun-2006 : Formatted for submitting code for   ( version 3.09 )

  !/                  inclusion in WAVEWATCH III.

  !/    29-May-2009 : Preparing distribution version.     ( version 3.14 )

  !/    14-Mar-2012 : Preparing distribution version.     ( version 4.05 )

  !/

  !/    Copyright 2009 National Weather Service (NWS),

  !/       National Oceanic and Atmospheric Administration.  All rights

  !/       reserved.  WAVEWATCH III is a trademark of the NWS.

  !/       No unauthorized use without permission.

  !/

  !  1. Purpose :

  !

  !     SHOWEX bottom friction source term (Ardhuin et al. 2003),

  !     using a subgrid depth parameterization based on Tolman (CE 1995).

  !

  !  2. Variables and types :

  !

  !      Name      Type  Scope    Description

  !     ----------------------------------------------------------------

  !     ----------------------------------------------------------------

  !

  !  3. Subroutines and functions :

  !

  !      Name      Type  Scope    Description

  !     ----------------------------------------------------------------

  !      W3SBT4    Subr. Public   SHOWEX bottom friction (movable bed)

  !      INSBT4    Subr. Public   Corresponding initialization routine.

  !      TABU_ERF  Subr. Public   Tabulation of ERF function

  !     ----------------------------------------------------------------

  !

  !  4. Subroutines and functions used :

  !

  !      Name      Type  Module   Description

  !     ----------------------------------------------------------------

  !      STRACE    Subr. W3SERVMD Subroutine tracing.

  !     ----------------------------------------------------------------

  !

  !  5. Remarks :

  !

  !     WAVEWATCH III is designed as a highly plug-compatible code.

  !     Source term modules can be included as self-contained modules,

  !     with limited changes needed to the interface of routine calls

  !     in W3SRCE, and in the point postprocessing programs only.

  !     Codes submitted for inclusion in WAVEWATCH III should be

  !     self-contained in the way described below, and might be

  !     provided with distributions fully integrated in the data

  !     structure, or as an optional version of this module to be

  !     included by the user.

  !

  !     Rules for preparing a module to be included in or distributed

  !     with WAVEWATCH III :

  !

  !      - Fully document the code following the outline given in this

  !        file, and according to all other WAVEWATCH III routines.

  !      - Provide a file with necessary modifications to W3SRCE and

  !        all other routines that require modification.

  !      - Provide a test case with expected results.

  !      - It is strongly recommended that the programming style used

  !        in WAVEWATCH III is followed, in particular

  !          a) for readability, write as if in fixed FORTRAN format

  !             regarding column use, even though all files are F90

  !             free format.

  !          b) I prefer upper case programming for permanent code,

  !             as I use lower case in debugging and temporary code.

  !

  !     This module needs to be self-contained in the following way.

  !

  !      a) All saved variables connected with this source term need

  !         to be declared in the module header. Upon acceptance as

  !         permanent code, they will be converted to the WAVEWATCH III

  !         dynamic data structure.

  !      b) Provide a separate computation and initialization routine.

  !         In the submission, the initialization should be called

  !         from the computation routine upon the first call to the

  !         routine. Upon acceptance as permanent code, the

  !         initialization routine will be moved to a more appropriate

  !         location in the code (i.e., being absorbed in ww3_grid or

  !         being moved to W3IOGR).

  !

  !     See notes in the file below where to add these elements.

  !

  !  6. Switches :

  !

  !     !/S  Enable subroutine tracing.

  !

  !  7. Source code :

  !/

  !/ ------------------------------------------------------------------- /

  !/

  !


  PUBLIC

  !

  ! Parameters for ERF function

  !

  INTEGER, PARAMETER      :: sizeerftable=300

  REAL                    :: erftable(0:sizeerftable)

  REAL                    :: delxerf

  REAL,    PARAMETER      :: xerfmax =  4. ! number of stdev

  !/

CONTAINS


  !/ ------------------------------------------------------------------- /

  SUBROUTINE insbt4

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |                         SHOM      |

    !/                  |            F. Ardhuin             |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         14-Mar-2012 |

    !/                  +-----------------------------------+

    !/

    !/    14-Mar-2012 : Origination.                        ( version 4.05 )

    !

    !  1. Purpose :

    !

    !     Initialization for bottom friction source term routine.

    !

    !  2. Method :

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !     ----------------------------------------------------------------

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      STRACE    Subr. W3SERVMD Subroutine tracing.

    !     ----------------------------------------------------------------

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      W3SBT4    Subr. W3SRC3MD Corresponding source term.

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !       None.

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

    !

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    !      NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    !/

    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_S

    CALL strace (ient, 'INSIN3')

#endif

    !

    ! 1.  .... ----------------------------------------------------------- *

    !

    CALL tabu_erf   !tabulates ERF function

    !/

    !/ End of INSBT4 ----------------------------------------------------- /

    !/

  END SUBROUTINE insbt4

  ! ----------------------------------------------------------------------


  SUBROUTINE tabu_erf

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |        J. Lepesqueur              |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         14-Mar-2012 |

    !/                  +-----------------------------------+

    !/

    !/    14-Mar-2012 : Origination.                        ( version 3.13 )

    !/

    !  1. Purpose :

    !     Tabulation of ERF function, which is used in bottom friction subgrid modelling

    !

    !     Initialization for source term routine.

    !

    !  2. Method :

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !     ----------------------------------------------------------------

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      STRACE    Subr. W3SERVMD Subroutine tracing.

    !     ----------------------------------------------------------------

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      W3SIN3    Subr. W3SRC3MD Corresponding source term.

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !       None.

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

    IMPLICIT NONE

    INTEGER :: I

    REAL :: x,y


    delxerf   = (2*xerfmax)/real(sizeerftable)

    DO i=0,sizeerftable

      x=-1.*xerfmax+i*delxerf

      if(x.lt.0.)then

        y=2**(1/2)*(1-abs(erf(x)))/2

      else

        y=2**(1/2)*(1+erf(x))/2

      end if

      erftable(i)=y

    END DO

    RETURN

    !/ ------------------------------------------------------------------- /

  END SUBROUTINE tabu_erf

  !/ ------------------------------------------------------------------- /


  !/ ------------------------------------------------------------------- /

  SUBROUTINE w3sbt4 (A, CG, WN, DEPTH, D50, PSIC, TAUBBL, BEDFORM, S, D, IX, IY )

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |            F. Ardhuin             |

    !/                  !            J. Lepesqueur          !

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         15-Mar-2012 |

    !/                  +-----------------------------------+

    !/

    !/    23-Jun-2011 : Origination.                        ( version 4.04 )

    !/    04-Jul-2011 : Adding momentum flux TAUBBL         ( version 4.05 )

    !/    15-Mar-2012 : Adding subgrid treatment for depth  ( version 4.05 )

    !/

    !  1. Purpose :

    !

    !     Computes the SHOWEX bottom friction with movable bed effects

    !

    !  2. Method :

    !     Uses a Gaussian distribution for friction factors, and estimates

    !     the contribution of rippled and non-rippled fractions based on

    !     the bayesian approach of Tolman (1995).

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !       A       R.A.  I   Action density spectrum.

    !       CG      R.A.  I   Group velocities.

    !       WN      R.A.  I   Wavenumbers.

    !       DEPTH   Real  I   Water depth.

    !       D50     Real  I   Median grain size.

    !       PSIC    Real  I   Critical Shields parameter

    !       BEFORMS Real I/O  Ripple parameters (roughness and wavelength).

    !       TAUBBL  Real  O   Components of stress leaking to the bottom.

    !       S       R.A.  O   Source term (1-D version).

    !       D       R.A.  O   Diagonal term of derivative.             *)

    !       IX,IY   Int. I   Spatial grid indices

    !     ----------------------------------------------------------------

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      STRACE    Subr. W3SERVMD Subroutine tracing.

    !     ----------------------------------------------------------------

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      W3SRCE    Subr. W3SRCEMD Source term integration.

    !      W3EXPO    Subr.   N/A    Point output post-processor.

    !      GXEXPO    Subr.   N/A    GrADS point output post-processor.

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !       None.

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /


    USE constants

    USE w3odatmd, ONLY: ndse

    USE w3servmd, ONLY: extcde

    USE w3gdatmd, ONLY: nk, nth, nspec, sig, dden,  &

         sbtcx, ecos, esin, dth


#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    !/

    LOGICAL, SAVE           :: FIRST = .true.

    REAL, INTENT(IN)        :: CG(NK), WN(NK), DEPTH, A(NSPEC), D50

    REAL, INTENT(IN)        :: PSIC

    INTEGER, INTENT(IN)     :: IX, IY

    REAL, INTENT(OUT)       :: S(NSPEC), D(NSPEC), TAUBBL(2)

    REAL, INTENT(INOUT)     :: BEDFORM(3)

    REAL                    :: CBETA(NK)

    REAL :: UORB2,UORB,AORB, EBX, EBY, AX, AY, LX, LY

    REAL :: CONST2, TEMP2

    REAL :: FW, KSUBN, KSUBS, KSUBR, MINADIM

    REAL :: SHIELDS(3), PSI, DELI1, DELI2, EB, XI, VARU, DD50

    INTEGER :: IK, ITH, IS, IND, INDE, ISUB


    REAL :: KRR, DSUB

    REAL DSUM(NK)

    ! These are the 3-point Gauss-Hermitte quadrature coefficients

    REAL, PARAMETER :: WSUB(3) = (/ 0.1666667,   0.1666666  , 0.6666667/)

    REAL, PARAMETER :: XSUB(3) = (/ -0.001,  0.001 , 0. /)


    REAL :: PROBA1, PROBA2, PSIX, PSIXT, PSIN2, DPSI , FACTOR

    REAL :: BACKGROUND


    !/

    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_S

    CALL strace (ient, 'W3SBT4')

#endif

    !

    ! 0.  Initializations ------------------------------------------------ *

    IF ( first ) THEN

      CALL insbt4

      first  = .false.

    END IF


    !

    ! 1.  Min / Max settings for grain size D50---------------------------- *

    !

    dd50=max(d50,1e-5)

    dd50=min(dd50,1.)

    !

    ! 1.1 Set background roughness when ripples are not active

    !

    background=max(sbtcx(6),sbtcx(7)*dd50)

    !

    ! 2. Subgrid loop

    !

    dsum(:)=0.

    taubbl(:)=0.

    !

    DO isub=1,3

      !

      ! 2.a  Computes bulk parameters : E, Uorb, Aorb------------------------- *

      !

      dsub=depth*(1.+xsub(isub))

      uorb=0.

      aorb=0.

      ax  =0.

      ay  =0.


      DO ik=1, nk

        IF ( wn(ik)*dsub .LT. 6. ) THEN

          eb  = 0.

          ebx = 0.

          eby = 0.

          DO ith=1, nth

            is=ith+(ik-1)*nth

            eb  = eb  + a(is)

            ebx = ebx +a(is)*ecos(ith)

            eby = eby +a(is)*esin(ith)

          END DO

          !

          ! U_bot=sigma * Zeta / sinh(KD)  and CBETA = 0.5*sigma^2 /(g*sinh^(kD))

          ! therefore variance(u_bot)= variance(elevation)*2*CBETA/D

          !

          !            CBETA(IK) = MAX(0., (CG(IK)*WN(IK)/SIG(IK)-0.5) )/DSUB

          cbeta(ik) = 0.5*sig(ik)**2 /(grav*(sinh(wn(ik)*dsub))**2)

          !  N.B.:  could also include shoaling effect on EB ...

          factor= (dden(ik) / cg(ik))*2*cbeta(ik)*grav

          varu= eb * factor

          uorb = uorb + varu

          aorb = aorb + varu/(sig(ik)**2)

          ax   = ax   + (ebx * factor)

          ay   = ay   + (eby * factor)

        ELSE

          cbeta(ik) = 0.

        END IF

      END DO

      !

      ! Computes RMS orbital amplitudes

      !

      uorb2 = 2*uorb

      uorb = sqrt(max(1.0e-7,uorb2))

      aorb = sqrt(max(1.0e-7,2*aorb))

      !

      ! Computes potential ripple wavelength, 1.7 = 2 * sqrt(2) * 0.6

      ! Based on Ardhuin et al. (JGR 2002): lambda = 0.6 * d_1/3

      !

      lx = aorb*1.7*ax/sqrt(ax**2+ay**2+1e-12)

      ly = aorb*1.7*ay/sqrt(ax**2+ay**2+1e-12)

      !

      ! 2.b First use of FWTABLE to get skin roughness and estimate Shields parameter

      !

      xi=max((alog10(max(aorb/dd50,0.3))-abmin)/delab,1.)

      ind  = min(sizefwtable-1, int(xi))

      deli1= min(1. ,xi-float(ind))

      deli2= 1. - deli1

      fw =fwtable(ind)*deli2+fwtable(ind+1)*deli1


      psi=fw*uorb2/(2.*grav*(sed_sg-1)*dd50)

      !

      ! Normalized Shields parameter

      !

      shields(isub)=psi/psic

      !

    END DO ! end of loop on ISUB

    dpsi=(shields(2)-shields(1))/(xsub(2)-xsub(1))*sbtcx(5)

    !

    ! Tests if the variation in psi is large enough to use subgrid

    !

    IF (abs(dpsi).LT.0.0001*shields(3).OR.abs(dpsi).LT.1.e-8) THEN

      !

      ! no subgrid in this case

      !

      IF(shields(3).GT.sbtcx(3)) THEN


        !  ripple roughness, see Ardhuin et al. (2003)

        ksubr=aorb*sbtcx(1)*shields(3)**sbtcx(2)

        !  Sheet flow roughness, see Wilson (1989)

        ksubs=aorb*0.0655*(uorb2/((sed_sg-1)*grav*aorb))**1.4

        ksubn = ksubr + ksubs

        bedform(2)=lx

        bedform(3)=ly

      ELSE

        !  relict roughness, see Ardhuin et al. (2003)

        ksubn=max(background,aorb*sbtcx(4))

        bedform(2)=-lx

        bedform(3)=-ly

      END IF


      bedform(1)=ksubn


    ELSE

      !

      ! subgrid in this case

      !

      psix=(sbtcx(3)-shields(3))/dpsi


      psixt=max((psix + xerfmax)/delxerf,0.)

      inde  = max(min(sizeerftable-1, int(psixt)),0)

      deli1  = min(1. ,psixt-float(inde))

      deli2  = 1. - deli1

      proba2=max(min(erftable(inde)*deli2+erftable(inde+1)*deli1,1.),0.)

      proba1 = 1. - proba2

      ! Mean psi with ripples (Tolman 1995, eq. XX)

      psin2=max(shields(3)+exp(-(0.5*psix**2))/sqrt(tpi)*dpsi/(proba2+0.0001),sbtcx(3))

      ! Sum of relict, ripple and sheet flow roughnesses

      ksubn = proba1*max(background,aorb*sbtcx(4))       &

           +proba2*aorb*(sbtcx(1)*psin2**sbtcx(2)+                   &

           0.0655*(uorb2/((sed_sg-1)*grav*aorb))**1.4)

      !

      IF (proba2.GT.0.5) THEN

        bedform(2)=lx

        bedform(3)=ly

      ELSE

        bedform(2)=-lx

        bedform(3)=-ly

      END IF

      !

    END IF

    bedform(1)=ksubn


    !

    ! 2.c second use of FWTABLE to get FW from the full roughness

    !

    xi=max((alog10(max(aorb/ksubn,0.3))-abmin)/delab,1.)

    ind  = min(sizefwtable-1, int(xi))

    deli1= min(1. ,xi-float(ind))

    deli2= 1. - deli1

    fw =fwtable(ind)*deli2+fwtable(ind+1)*deli1

    !

    ! 5. Fills output arrays and estimates the energy and momentum loss

    !

    DO ik=1, nk

      const2=dden(ik)/cg(ik) &         !Jacobian to get energy in band

           *grav/(sig(ik)/wn(ik))    ! coefficient to get momentum

      dsum(ik)=-fw*uorb*cbeta(ik)      !*WSUB(ISUB)

      DO ith=1,nth

        is=ith+(ik-1)*nth

        d(is)=dsum(ik)

        temp2=const2*d(is)*a(is)

        taubbl(1) = taubbl(1) - temp2*ecos(is)

        taubbl(2) = taubbl(2) - temp2*esin(is)

        s(is)=d(is)*a(is)

      END DO

    END DO

    !

    RETURN

    !/

    !/ End of W3SBT4 ----------------------------------------------------- /

    !/

  END SUBROUTINE w3sbt4


  !/ ------------------------------------------------------------------- /


END MODULE w3sbt4md