WW3/w3swldmd_8F90_source.html

#include "w3macros.h"

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

MODULE w3swldmd

  !/

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

  !/                  | WAVEWATCH III           NOAA/NCEP |

  !/                  |           H. L. Tolman            |

  !/                  |                        FORTRAN 90 |

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

  !/

  !/    21-Nov-2011 : Origination.                        ( version 4.07 )

  !/

  !/    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 :

  !

  !     Source term module for swell dissipation based on different

  !     physics that can be independently selected form the input

  !     and whitecapping dissipation terms in the model setup.

  !

  !  2. Variables and types :

  !

  !      Name      Type  Scope    Description

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

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

  !

  !  3. Subroutines and functions :

  !

  !      Name      Type  Scope    Description

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

  !      W3SWL4    Subr. Public   Ardhuin et al (2010+) swell dissipation

  !      W3SWL6    Subr. Public   Babanin (2011) swell dissipation

  !

  !      IRANGE    Func. Private  Generate a sequence of integer values

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

  !

  !  4. Subroutines and functions used :

  !

  !      Name      Type  Module   Description

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

  !      STRACE    Subr. W3SERVMD Subroutine tracing.

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

  !

  !  5. Remarks :

  !

  !

  !  6. Switches :

  !

  !     !/S  Enable subroutine tracing.

  !

  !  7. Source code :

  !/

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

  !/

  PUBLIC  :: w3swl4, w3swl6

  PRIVATE :: irange

  !/

CONTAINS

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

  SUBROUTINE w3swl4 (A, CG, WN, DAIR, S, D)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           H. L. Tolman            |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         13-Aug-2021 |

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

    !/

    !/    29-May-2009 : Origination (w3srcxmd.ftn)          ( version 3.14 )

    !/    06-Jan-2012 : Implementation                         (S. Zieger)

    !/    13-Aug-2021 : Consider DAIR a variable           ( version x.xx )

    !/

    !  1. Purpose :

    !

    !     FIXME

    !

    !  2. Method :

    !

    !  3. Parameters :

    !

    !      Parameter list

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

    !      A¹      R.A. I  Action density spectrum

    !      CG      R.A. I  Group velocities

    !      WN      R.A. I  Wavenumbers

    !      DAIR    R.A. I   Air density

    !      S¹      R.A. O  Source term

    !      D¹      R.A. O  Diagonal term of derivative

    !      ¹ Stored as 1-D array with dimension NTH*NK (column by column).

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

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !      IRANGE    Func. W3SWLDMD

    !      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 comments in source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

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

    USE constants, ONLY: grav, dwat

    USE w3gdatmd,  ONLY: nk, nth, nspec, sig2, dden, fte, swl6b1

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    IMPLICIT NONE

    !/

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

    !/ Parameter list

    REAL, INTENT(IN)  :: a(nspec), cg(nk), wn(nk), dair

    REAL, INTENT(OUT) :: s(nspec), d(nspec)

    !/

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

    !/ Local parameters

#ifdef W3_S

    INTEGER, SAVE     :: ient = 0

#endif

    INTEGER           :: ikn(nk), ith

    REAL, PARAMETER   :: va = 1.4e-5 ! Air kinematic viscosity (used in WAM).

    REAL              :: eb(nk), wn2(nspec), emean

    REAL              :: fe, aorb, re, recrit, uosig, cdsv

    !/

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

    !/

#ifdef W3_S

    CALL strace (ient, 'W3SWL4')

#endif

    !

    ikn = irange(1,nspec,nth)

    d   = 0.

    wn2 = 0.

    !

    DO ith = 1, nth

      wn2(ikn+(ith-1)) = wn        ! Wavenumbers to all directions.

    END DO

    !

    eb     = sum(reshape(a,(/ nth,nk /)),1) * dden(1:nk) / cg

    emean  = sum(eb) + (eb(nk) / dden(nk)) * fte

    !

    aorb   = 2.0*sqrt(emean)

    !

    eb     = sum(reshape(a*sig2**2,(/ nth,nk /)),1) * dden(1:nk) / cg

    uosig  = 2.0*sqrt(sum(eb))


    fe     = swl6b1    ! (from NAMELIST)

    !     FE     = 0.001     ! (from NAMELIST)

    !/             0.001 - 0.019 with median value 0.007 (Ardhuin et al 2009, Babanin 2011)

    cdsv   = 1.2

    !

    recrit = 1.0e5

    re     = 4.0 * uosig * aorb / va

    !

    IF (re .GT. recrit) THEN

      d = -(16.0/grav) * (dair/dwat) * fe * (sig2**2) *uosig

    ELSE

      d = -2.0 * (dair/dwat) * cdsv * wn2 * sqrt(2.0 * va * sig2)

    END IF

    !

    s = d * a

    !

    !  WRITE(*,*) ' FE       =',FE

    !  WRITE(*,*) ' HS       =',4.*SQRT(EMEAN)

    !  WRITE(*,*) ' UOSIG    =',UOSIG

    !  WRITE(*,*) ' AORB     =',AORB

    !  WRITE(*,*) ' RE/RECRIT=',RE/RECRIT

    !  WRITE(*,*) ' SWL4_tot =',SUM(SUM(RESHAPE(S,(/ NTH,NK /)),1)*DDEN/CG)

    !/

    !/ End of W3SWL4 ----------------------------------------------------- /

    !/

  END SUBROUTINE w3swl4

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

  SUBROUTINE w3swl6 (A, CG, WN, S, D)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           H. L. Tolman            |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         16-Feb-2012 |

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

    !/

    !/    29-May-2009 : Origination (w3srcxmd.ftn)          ( version 3.14 )

    !/    16-Feb-2012 : Implementation                      ( version 4.07 )

    !/                                                         (S. Zieger)

    !/

    !  1. Purpose :

    !

    !     Turbulent dissipation of narrow-banded swell as described in

    !     Babanin (2011, Section 7.5).

    !

    !     Babanin 2011: Cambridge Press, 295-321, 463pp.

    !

    !  2. Method :

    !

    !     S = D * A

    !

    !  3. Parameters :

    !

    !      Parameter list

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

    !      A¹      R.A. I  Action density spectrum

    !      CG      R.A. I  Group velocities

    !      WN      R.A. I  Wavenumbers

    !      S¹      R.A. O  Source term

    !      D¹      R.A. O  Diagonal term of derivative

    !      ¹ Stored as 1-D array with dimension NTH*NK (column by column).

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

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !      IRANGE    Func. W3SWLDMD

    !      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 comments in source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

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

    USE constants, ONLY: grav

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

    USE w3gdatmd,  ONLY: swl6cstb1, swl6b1, fte, ftwn

#ifdef W3_T6

    USE w3odatmd,  ONLY: ndst

#endif

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    IMPLICIT NONE

    !/

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

    !/ Parameter list

    REAL, INTENT(IN)    :: a(nspec), cg(nk), wn(nk)

    REAL, INTENT(OUT)   :: s(nspec), d(nspec)

    !/

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

    !/ Local parameters

#ifdef W3_S

    INTEGER, SAVE     :: ient = 0

#endif

    INTEGER             :: ik, ith, ikn(nk)

    REAL, DIMENSION(NK) :: aband, kmax, anar, bn, aorb, ddis

    REAL                :: k(nth,nk), b1

    !/

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

    !/

#ifdef W3_S

    CALL strace (ient, 'W3SWL6')

#endif

    !

    !/ 0) --- Initialize parameters -------------------------------------- /

    ikn   = irange(1,nspec,nth)            ! Index vector for array access, e.g.

    !                                            ! in form of WN(1:NK) == WN2(IKN).

    aband = sum(reshape(a,(/ nth,nk /)),1) ! action density as function of wavenumber

    ddis  = 0.

    d     = 0.

    b1    = swl6b1                         ! empirical constant from NAMELIST

    !

    !/ 1) --- Choose calculation of steepness a*k ------------------------ /

    !/        Replace the measure of steepness with the spectral

    !         saturation after Banner et al. (2002) ---------------------- /

    k     = reshape(a,(/ nth,nk /))

    kmax  = maxval(k,1)

    DO ik = 1,nk

      IF (kmax(ik).LT.1.0e-34) THEN

        k(1:nth,ik) = 1.

      ELSE

        k(1:nth,ik) = k(1:nth,ik)/kmax(ik)

      END IF

    END DO

    anar  = 1.0/( sum(k,1) * dth )

    bn    = anar * ( aband * sig(1:nk) * dth ) * wn**3

    !

    IF (.NOT.swl6cstb1) THEN

      !

      !/    --- A constant value for B1 attenuates swell too strong in the

      !/        western central Pacific (i.e. cross swell less than 1.0m).

      !/        Workaround is to scale B1 with steepness a*kp, where kp is

      !/        the peak wavenumber. SWL6B1 remains a scaling constant, but

      !/        with different magnitude.  --------------------------------- /

      ik    = maxloc(aband,1)         ! Index for peak

      !         EMEAN = SUM(ABAND * DDEN / CG)  ! Total sea surface variance

      b1    = swl6b1 * ( 2. * sqrt(sum(aband*dden/cg)) * wn(ik) )

      !

    END IF

    !

    !/ 2) --- Calculate the derivative term only (in units of 1/s) ------- /

    DO ik = 1,nk

      IF (aband(ik) .GT. 1.e-30) THEN

        ddis(ik) = -(2./3.) * b1 * sig(ik) * sqrt(bn(ik))

      END IF

    END DO

    !

    !/ 3) --- Apply dissipation term of derivative to all directions ----- /

    DO ith = 1, nth

      d(ikn+(ith-1)) = ddis

    END DO

    !

    s = d * a

    !

    !  WRITE(*,*) ' B1       =',B1

    !  WRITE(*,*) ' DDIS_tot =',SUM(DDIS*ABAND*DDEN/CG)

    !  WRITE(*,*) ' EDENS_tot=',sum(aband*dden/cg)

    !  WRITE(*,*) ' EDENS_tot=',sum(aband*sig*dth*dsii/cg)

    !  WRITE(*,*) ' '

    !  WRITE(*,*) ' SWL6_tot =',sum(SUM(RESHAPE(S,(/ NTH,NK /)),1)*DDEN/CG)

    !

    !/

    !/ End of W3SWL6 ----------------------------------------------------- /

    !/

  END SUBROUTINE w3swl6

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

  !/

  FUNCTION irange(X0,X1,DX) RESULT(IX)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           H. L. Tolman            |

    !/                  |           S. Zieger               |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         15-Feb-2011 |

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

    !/

    !/    15-Feb-2011 : Origination from W3SRC6MD          ( version 4.07 )

    !/                                                        (S. Zieger)

    !/

    !  1. Purpose :

    !         Generate a linear-spaced sequence of integer

    !         numbers. Used for array addressing (indexing).

    !

    !/

    IMPLICIT NONE

    INTEGER, INTENT(IN)  :: x0, x1, dx

    INTEGER, ALLOCATABLE :: ix(:)

    INTEGER              :: n

    INTEGER              :: i

    !

    n = int(real(x1-x0)/real(dx))+1

    ALLOCATE(ix(n))

    DO i = 1, n

      ix(i) = x0+ (i-1)*dx

    END DO

    !/

  END FUNCTION irange

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

  !/

END MODULE w3swldmd