WW3/w3uostmd_8F90_source.html

#include "w3macros.h"

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


MODULE w3uostmd

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

  !/                  | WAVEWATCH III           NOAA/NCEP |

  !/                  |             Lorenzo Mentaschi     |

  !/                  |                        FORTRAN 90 |

  !/                  | Last update :         08-Oct-2018 |

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

  !/

  !/    Aug-2018 : Origination.                        ( version 6.07 )

  !/ 18-Sep-2019 : Added UNIT_AB and changed unit      ( version 7.06 )

  !/               number to 110 (C. Bunney, UKMO)

  !/

  !/    Copyright 2010 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 : Parmeterization of the unresoled obstacles

  !

  !  2. Subroutines and functions :

  !

  !      Name               Type  Scope    Description

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

  !    UOST_INITGRID        Subr. Public   allocates UOST variables on each grid,

  !                                        and loads the matrices of alpha and beta coefficient

  !    UOST_SETGRID         Subr. Public   sets the actual computation grid in the source term

  !    UOST_SRCTRMCOMPUTE   Subr. Public   computes the source term for a given input spectrum

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

  !

  !  3. Switches :

  !

  !  4. Source code :

  !/

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

  !/


  USE w3gdatmd, ONLY: grid, sgrd, grids, sgrds

  USE w3odatmd, ONLY: ndso, ndse

  USE w3servmd, ONLY: extcde

#ifdef W3_S

  USE w3servmd, ONLY: strace

#endif


  IMPLICIT NONE


  PUBLIC :: uost_initgrid

  PUBLIC :: uost_setgrid

  PUBLIC :: uost_srctrmcompute


  PRIVATE


  TYPE uost_sourceterm

    REAL, ALLOCATABLE :: COSTH(:), SINTH(:)

    REAL :: GAMMAUP = 10

    REAL :: GAMMADOWN = 20

    ! griddata is a pointer to the grid actually computed

    TYPE(GRID), POINTER :: GRD

    TYPE(SGRD), POINTER :: SGD

  CONTAINS

    !PROCEDURE, PASS, PRIVATE :: COMPUTE_PSI => UOST_SOURCETERM_COMPUTE_PSI


    !compute_ld: estimates the local dissipation (private method)

    PROCEDURE, PASS, PRIVATE :: COMPUTE_LD => uost_sourceterm_compute_ld

    !compute_se: estimates the shadow effect (private method)

    PROCEDURE, PASS, PRIVATE :: COMPUTE_SE => uost_sourceterm_compute_se

    !compute: estimates the whole dissipation

    PROCEDURE, PASS :: COMPUTE => uost_sourceterm_compute

    !setgrid: sets grd pointer and computes some cached structures

    PROCEDURE, PASS :: SETGRID => uost_sourceterm_setgrid

  END TYPE uost_sourceterm


  ! srctrm: global singleton source term

  CLASS(UOST_SOURCETERM), ALLOCATABLE :: SRCTRM


  INTEGER, PARAMETER :: UNIT_AB = 110   ! Unit number for LOAD_ALPHABETA


CONTAINS


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


  SUBROUTINE uost_initgrid(IGRID, FILELOCAL, FILESHADOW, LOCALFACTOR, SHADOWFACTOR)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |          Lorenzo Mentaschi        |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         01-Oct-2018 |

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

    !/

    !/    Aug-2018 : Origination.                        ( version 6.07 )

    !/

    !  1. Purpose : allocate the UOST variables for a given grid, and load

    !               them from file

    !  2. Parameters :

    !

    !     Parameter list

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

    !     IGRID, integer:        id of the grid being initialized

    !     FILELOCAL, string:     file from where the alpha/beta coefficient

    !                            for the local dissipation are loaded

    !     FILESHADOW, string:    file from where the alpha/beta coefficient

    !                            for the shadow effect are loaded

    !     LOCALFACTOR, double:   adjustment parameter for the local

    !                            dissipation alpha and beta

    !     SHADOWFACTOR, double:  adjustment parameter for the shadow

    !                            dissipation alpha and beta

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

    !

    !  3. Called by :

    !

    !      Name      Type  Module     Description

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

    !      W3IOGR    Subr. W3IOGRMD   Initialization of grid objects

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

    !

    !  4. Source code :

    !

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

    IMPLICIT NONE


    INTEGER, INTENT(IN) :: igrid

    CHARACTER(LEN=*), INTENT(IN) :: filelocal, fileshadow

    REAL, INTENT(IN) :: localfactor, shadowfactor

    TYPE(grid), POINTER :: grd

    TYPE(sgrd), POINTER :: sgd

    REAL :: cgmax, minsize

#ifdef W3_S

    INTEGER, SAVE           :: ient   = 0

#endif


#ifdef W3_S

    CALL strace (ient, 'UOST_INITGRID')

#endif


    IF ( (igrid .LE. 0) .OR. (.NOT. ALLOCATED(grids)) ) THEN

      RETURN

    ENDIF


    grd => grids(igrid)

    sgd => sgrds(igrid)

    grd%UOSTFILELOCAL = filelocal

    grd%UOSTFILESHADOW = fileshadow

    grd%UOSTLOCALFACTOR = localfactor

    grd%UOSTSHADOWFACTOR = shadowfactor


    ALLOCATE( grd%UOST_LCL_OBSTRUCTED(grd%NX, grd%NY) )

    grd%UOST_LCL_OBSTRUCTED = .false.

    ALLOCATE( grd%UOST_SHD_OBSTRUCTED(grd%NX, grd%NY) )

    grd%UOST_SHD_OBSTRUCTED = .false.

    ALLOCATE( grd%UOSTCELLSIZE(grd%NX, grd%NY, sgd%NTH) )

    grd%UOSTCELLSIZE = 0

    ALLOCATE( grd%UOSTLOCALALPHA(grd%NX, grd%NY, sgd%NK, sgd%NTH) )

    grd%UOSTLOCALALPHA = 100

    ALLOCATE( grd%UOSTLOCALBETA(grd%NX, grd%NY, sgd%NK, sgd%NTH) )

    grd%UOSTLOCALBETA = 100

    ALLOCATE( grd%UOSTSHADOWALPHA(grd%NX, grd%NY, sgd%NK, sgd%NTH) )

    grd%UOSTSHADOWALPHA = 100

    ALLOCATE( grd%UOSTSHADOWBETA(grd%NX, grd%NY, sgd%NK, sgd%NTH) )

    grd%UOSTSHADOWBETA = 100


    IF ( (igrid .GT. 0) .AND. ( ALLOCATED(srctrm)) ) THEN

      !  loading local/shadow alpha/beta

      CALL load_alphabeta(grd, sgd, unit_ab)


      !  warning the user that for cells too small UOST may be inaccurate

      cgmax = 20 ! simply taking a high value for the max group velocity to give an indication of this threshold

      minsize = cgmax*grd%DTMAX/1000

      WRITE(ndso,*)'*** WAVEWATCH-III WARNING IN W3UOST/UOST_INITGRID'

      WRITE(ndso,*)'UOST: grid ',trim(grd%GNAME),':'

      WRITE(ndso,*)'      global time step == ', grd%DTMAX, ' s'

      WRITE(ndso,*)'      FOR CELLS SMALLER THAN ABOUT ', minsize,    &

           ' KM UOST MAY UNDERESTIMATE THE DISSIPATION'

      WRITE(ndso,*)

    ENDIF

  END SUBROUTINE uost_initgrid


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


  SUBROUTINE uost_setgrid(IGRID)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |         Lorenzo Mentaschi         |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         01-Oct-2018 |

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

    !/

    !/    Aug-2018 : Origination.                        ( version 6.07 )

    !/

    !  1. Purpose : sets the current grid in the sourceterm object

    !  2. Parameters :

    !

    !     Parameter list

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

    !     IGRID, integer:        id of the actual grid

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

    !

    !  3. Called by :

    !

    !      Name      Type  Module     Description

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

    !      W3INIT    Subr. W3INITMD   Initialization of grid objects

    !      W3WAVE    Subr. W3WAVEMD   Initialization of grid objects

    !                                 before computation

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

    !

    !  4. Source code :

    !

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

    IMPLICIT NONE


    INTEGER, INTENT(IN) :: igrid

#ifdef W3_S

    INTEGER, SAVE           :: ient   = 0

#endif


#ifdef W3_S

    CALL strace (ient, 'UOST_SETGRID')

#endif


    IF ( .NOT. ALLOCATED(srctrm) ) THEN

      ALLOCATE(srctrm)

    ENDIF


    CALL srctrm%SETGRID(grids(igrid), sgrds(igrid))

  END SUBROUTINE uost_setgrid

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


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

  SUBROUTINE uost_srctrmcompute(IX, IY, SPEC, CG, DT, U10ABS, U10DIR, S, D)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |          Lorenzo Mentaschi        |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         01-Oct-2018 |

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

    !/

    !/    Aug-2018 : Origination.                        ( version 6.07 )

    !/

    !  1. Purpose : estimates the UOST source term for a give spectrum

    !  2. Parameters :

    !

    !     Parameter list

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

    !     IGRID, integer:        id of the actual grid

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

    !

    !  3. Called by :

    !

    !      Name      Type  Module     Description

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

    !      W3SRCE    Subr. W3SRCEMD   Computation of the source terms

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

    !

    !  4. Source code :

    !

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

    IMPLICIT NONE


    INTEGER, INTENT(IN) :: ix, iy

    REAL, INTENT(IN) :: dt

    REAL, INTENT(IN) :: spec(srctrm%sgd%nspec), cg(srctrm%sgd%nk)

    REAL, INTENT(IN) :: u10abs, u10dir

    REAL, INTENT(OUT) :: s(srctrm%sgd%nspec), d(srctrm%sgd%nspec)

#ifdef W3_S

    INTEGER, SAVE           :: ient   = 0

#endif


#ifdef W3_S

    CALL strace (ient, 'UOST_SRCTRMCOMPUTE')

#endif


    CALL srctrm%COMPUTE(ix, iy, spec, cg, dt, u10abs, u10dir, s, d)

  END SUBROUTINE uost_srctrmcompute


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


  SUBROUTINE load_alphabeta(GRD, SGD, FILEUNIT)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |          Lorenzo Mentaschi        |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         01-Oct-2018 |

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

    !/

    !/    Aug-2018 : Origination.                        ( version 6.07 )

    !/

    !  1. Purpose : loads local and shadow alpha and beta from files

    !  2. Parameters :

    !

    !     Parameter list

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

    !     GRD, GRID type:        object representing the spatial grid to

    !                            be loaded

    !     SGD, SGRD type:        object representing the current spectral grid

    !     FILEUNIT, Integer:     unit id of the input files

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

    !

    !  3. Called by :

    !

    !      Name            Type  Module     Description

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

    !      UOST_INITGRID   Subr. W3UOSTMD   Initialization of the UOST grid

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

    !

    !  4. Source code :

    !

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

    IMPLICIT NONE


    TYPE(grid), INTENT(INOUT) :: grd

    TYPE(sgrd), INTENT(IN) :: sgd

    INTEGER, INTENT(IN) :: fileunit

    CHARACTER(256) :: filename

    LOGICAL :: fileexists

    INTEGER :: jg, j, l, i, ix, iy

#ifdef W3_S

    INTEGER, SAVE           :: ient   = 0

#endif


#ifdef W3_S

    CALL strace (ient, 'LOAD_ALPHABETA')

#endif


    ! LOADING LOCAL ALPHA/BETA

    filename = grd%UOSTFILELOCAL

    INQUIRE(file=filename, exist=fileexists)


    j = len_trim(filename)

    IF (.NOT. fileexists) THEN

      WRITE(ndse,*)'*** WAVEWATCH III ERROR IN W3UOST: '// &

           'FILE '//filename(:j)//' NOT FOUND. QUITTING'

      CALL extcde (9999)

    ENDIF

    WRITE(ndso,*)'FILE '//filename(:j)//' FOUND.'// &

         'LOADING UOST SETTINGS FOR GRID '//grd%GNAME


    CALL load_alphabeta_fromfile(fileunit, filename(:j), grd%NX, grd%NY, sgd%NK, sgd%NTH,&

         grd%UOSTABMULTFACTOR, grd%UOSTLOCALALPHA, grd%UOSTLOCALBETA,&

         grd%UOSTCELLSIZE, grd%UOST_LCL_OBSTRUCTED)


    ! LOADING SHADOW ALPHA/BETA

    filename = grd%UOSTFILESHADOW

    INQUIRE(file=filename, exist=fileexists)


    j = len_trim(filename)

    IF (.NOT. fileexists) THEN

      WRITE(ndse,*)'*** WAVEWATCH III ERROR IN W3UOST: '// &

           'FILE '//filename(:j)//' NOT FOUND. QUITTING'

      CALL extcde (9999)

    ENDIF

    WRITE(ndso,*)'FILE '//filename(:j)//' FOUND.'//&

         'LOADING UOST SETTINGS FOR GRID '//grd%GNAME


    CALL load_alphabeta_fromfile(fileunit, filename(:j), grd%NX, grd%NY, sgd%NK, sgd%NTH,&

         grd%UOSTABMULTFACTOR, grd%UOSTSHADOWALPHA, grd%UOSTSHADOWBETA,&

         grd%UOSTCELLSIZE, grd%UOST_SHD_OBSTRUCTED)


  END SUBROUTINE load_alphabeta


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


  SUBROUTINE load_alphabeta_fromfile(FILEUNIT, FILENAME, NX, NY, NK, NTH,&

       MULTFACTOR, ALPHAMTX, BETAMTX, CELLSIZE, ISOBSTRUCTED)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |          Lorenzo Mentaschi        |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         01-Oct-2018 |

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

    !/

    !/    Aug-2018 : Origination.                        ( version 6.07 )

    !/

    !  1. Purpose : loads alpha and beta from a single obstructions file

    !  2. Parameters :

    !

    !     Parameter list

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

    !     FILEUNIT, Integer:             unit of the file to be opened

    !     FILENAME, string:              name of the file

    !     NX, NY, NK, NTH, Integer:      size of the spatial/spectral grid

    !     MULTFACTOR, REAL:              multiplication factor for alpha and beta:

    !                                    alpha and beta should be real in [0,1]

    !                                    but to save memory the are stored in Integer*1

    !     ALPHAMTX, BETAMTX, Integer*1:  loaded alpha and beta spatial/spectral matrices

    !     CELLSIZE, REAL, REAL:          cell size for each spectral direction,

    !                                    also loaded from the file

    !     ISOBSTRUCTED, LOGICAL:         matrix of logicals, indicating for each cell

    !                                    if it is obstructed or not

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

    !

    !  3. Called by :

    !

    !      Name            Type  Module     Description

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

    !      LOAD_ALPHABETA  Subr. W3UOSTMD   Initialization of the UOST grid

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

    !

    !  4. Source code :

    !

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

    IMPLICIT NONE


    CHARACTER(*), INTENT(IN) :: filename

    REAL, INTENT(IN) :: multfactor

    INTEGER, INTENT(IN) :: fileunit, nx, ny, nk, nth

    INTEGER*1, INTENT(INOUT) :: alphamtx(:,:,:,:), betamtx(:,:,:,:)

    real*4, INTENT(INOUT) ::  cellsize(:,:,:)

    LOGICAL, INTENT(INOUT) ::  isobstructed(:,:)

    CHARACTER(LEN=600) :: line

    INTEGER :: fiostat

    LOGICAL :: header, filestart, readingcellsize, readingalpha

    INTEGER :: ix, iy, spgrds_size, ik

    REAL, ALLOCATABLE :: trans(:)

#ifdef W3_S

    INTEGER, SAVE           :: ient   = 0

#endif


#ifdef W3_S

    CALL strace (ient, 'LOAD_ALPHABETA_FROMFILE')

#endif


    !  INITIALIZING LOGICALS REPRESENTING THE DIFFERENT PHASES OF THE LOAD

    filestart = .true.

    header = .true.;

    readingcellsize = .false.

    readingalpha = .false.

    ik = 0


    ALLOCATE(trans(nth))


    OPEN(fileunit, file=filename, status='OLD', action='READ')

    read_loop: DO

      READ(fileunit, '(A)', iostat=fiostat) line


      IF (fiostat .NE. 0) EXIT read_loop


      IF (line(1:1) .EQ. '$') cycle


      IF (filestart) THEN

        !  reading the first line

        READ(line, '(I5)') spgrds_size

        filestart = .false.

      ELSEIF (header) THEN

        !  reading the position of an obstructed cell

        READ(line, *) ix, iy

        isobstructed(ix, iy) = .true.

        IF ((ix .GT. nx) .OR. (iy .GT. ny)) THEN

          WRITE(ndse,*) '*** WAVEWATCH III ERROR IN W3UOST: '// &

               'GRID INDICES OUT OF RANGE.'// &

               'CHECK FILE '//filename

          CALL extcde (9999)

        ENDIF

        !  marking the end of the reading of the header

        header = .false.

        ik = 1

        readingcellsize = .true.

      ELSEIF (readingcellsize) THEN

        !  reading the sizes of the cell

        READ(line, *) cellsize(ix, iy, :)

        readingcellsize = .false.

        readingalpha = .true.

      ELSE

        READ(line, *) trans

        IF (readingalpha) THEN

          !  reading alpha for frequency IK

          alphamtx(ix, iy, ik, :) = nint(trans*multfactor)

        ELSE

          !  reading beta for frequency IK

          betamtx(ix, iy, ik, :) = nint(trans*multfactor)

        ENDIF

        IF (ik .LT. nk) THEN

          ik = ik + 1

        ELSE IF (readingalpha) THEN

          !  preparing to read the next cell

          readingalpha = .false.

          ik = 1

        ELSE

          header = .true.

          ik = 1

        ENDIF

      ENDIF

    ENDDO read_loop

    CLOSE(fileunit)


    DEALLOCATE(trans)


  END SUBROUTINE load_alphabeta_fromfile


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

  SUBROUTINE uost_sourceterm_setgrid(THIS, GRD, SGD)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |          Lorenzo Mentaschi        |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         01-Oct-2018 |

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

    !/

    !/    Aug-2018 : Origination.                        ( version 6.07 )

    !/

    !  1. Purpose : method of the class UOST_SOURCETERM,

    !               to set the actual spatial and spectral grid

    !  2. Parameters :

    !

    !     Parameter list

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

    !     GRD, GRID type:        object representing the spatial grid to

    !                            be loaded

    !     SGD, SGRD type:        object representing the current spectral grid

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

    !

    !  3. Called by :

    !

    !      Name            Type  Module     Description

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

    !      UOST_SETGRID  Subr. W3UOSTMD   Setting the actual computation grid

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

    !

    !  4. Source code :

    !

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

    IMPLICIT NONE


    CLASS(uost_sourceterm), INTENT(INOUT) :: this

    TYPE(grid), TARGET, INTENT(IN) :: grd

    TYPE(sgrd), TARGET, INTENT(IN) :: sgd

    INTEGER :: ith, nth

#ifdef W3_S

    INTEGER, SAVE           :: ient   = 0

#endif


#ifdef W3_S

    CALL strace (ient, 'UOST_SOURCETERM_SETGRID')

#endif


    this%GRD => grd

    this%SGD => sgd


    IF (ALLOCATED(this%COSTH)) THEN

      DEALLOCATE(this%COSTH)

      DEALLOCATE(this%SINTH)

    ENDIF


    nth = this%SGD%NTH

    ALLOCATE(this%COSTH(nth))

    ALLOCATE(this%SINTH(nth))

    DO ith=1,nth

      this%COSTH(ith) = cos(sgd%TH(ith))

      this%SINTH(ith) = sin(sgd%TH(ith))

    ENDDO

  END SUBROUTINE uost_sourceterm_setgrid


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


  SUBROUTINE compute_reduction_psi(U10ABS, U10DIR, CGABS, CGDIR, DT, PSI)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |          Lorenzo Mentaschi        |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         01-Oct-2018 |

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

    !/

    !/    Aug-2018 : Origination.                        ( version 6.07 )

    !/

    !  1. Purpose : In conditions of wind sea, the effect

    !               of the unresolved obstacles is reduced.

    !               Here a reduction psi is computed, as a function of the

    !               wave age.

    !  2. Parameters :

    !

    !     Parameter list

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

    !     U10ABS, real:     absolute value of U10

    !     U10DIR, real:     direction of U10

    !     CGABS, real:      absolute value of the group velocity

    !     CGDIR, real:      direction of the group velocity

    !     DT, real:         time step

    !     PSI, real:        output psi factor

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

    !

    !  3. Called by :

    !

    !      Name                   Type  Module     Description

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

    !      UOST_SOURCETERM_COMPUTE_LD  Subr. W3UOSTMD   Computing the local dissipation

    !      UOST_SOURCETERM_COMPUTE_SE  Subr. W3UOSTMD   Computing the shadow effect

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

    !

    !  4. Source code :

    !

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

    USE constants, ONLY: pi


    IMPLICIT NONE


    REAL, PARAMETER :: tolerance = 0.000001

    REAL, PARAMETER :: whthr1 = .5, whthr2 = 1.5


    REAL, INTENT(IN) :: u10abs, u10dir, cgabs, cgdir, dt

    REAL, INTENT(OUT) :: psi

    REAL :: thdelta, cp, wa

#ifdef W3_S

    INTEGER, SAVE           :: ient   = 0

#endif


#ifdef W3_S

    CALL strace (ient, 'COMPUTE_REDUCTION_PSI')

#endif


    ! computing the wave age

    thdelta = abs(u10dir - cgdir)

    DO WHILE (thdelta .GT. pi)

      thdelta = thdelta - 2*pi

    ENDDO

    thdelta = abs(thdelta)

    IF (pi/2 - thdelta .GT. tolerance) THEN

      cp = cgabs*2 ! this is scrictly valid only in deep water

      wa = cp/u10abs/cos(thdelta)

    ELSE

      wa = 9999999 ! a very high number

    ENDIF


    IF (wa .LE. whthr1) THEN

      !  if the wave age is less that 0.5, psi = 0, i.e.

      !  no unresolved obstacle is considered

      psi = 0

    ELSEIF ((wa .GT. whthr1) .AND. (wa .LT. whthr2)) THEN

      !  if the wave age is between 0.5 and 1.5

      !  psi scales linearly with WA

      psi = (wa - whthr1)/(whthr2 - whthr1)

    ELSE

      !  if the wave age is greater than 1.5 psi = 1

      psi = 1

    ENDIF


  END SUBROUTINE compute_reduction_psi


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


  SUBROUTINE uost_sourceterm_compute_ld(THIS, IX, IY, SPEC, CG, DT, U10ABS, U10DIR, S, D)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |          Lorenzo Mentaschi        |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         01-Oct-2018 |

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

    !/

    !/    Aug-2018 : Origination.                        ( version 6.07 )

    !/

    !  1. Purpose : Method of the class UOST_SOURCETERM.

    !               Computation of the local dissipation of the spectrum

    !  2. Parameters :

    !

    !     Parameter list

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

    !     THIS: UOST_SOURCETERM      instance of UOST_SOURCETERM passed to the method

    !                           (compulsory in oo programming)

    !     IX, IY: Integer       coordinates of the actual cell

    !     SPEC: real            input spectrum

    !     CG: real              group velocity

    !     DT: real              time step

    !     U10ABS: real          absolute value of U10

    !     U10DIR: real          direction of U10

    !     S: real               source term

    !     D: real               differential of the source term over the spectrum

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

    !

    !  3. Called by :

    !

    !      Name                   Type  Module     Description

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

    !      UOST_SOURCETERM_COMPUTE     Subr. W3UOSTMD   Computing the source term

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

    !

    !  4. Source code :

    !

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

    IMPLICIT NONE


    CLASS(uost_sourceterm), INTENT(INOUT) :: THIS

    INTEGER, INTENT(IN) :: IX, IY

    REAL, INTENT(IN) :: SPEC(THIS%SGD%NSPEC), CG(THIS%SGD%NK)

    REAL, INTENT(OUT) :: S(THIS%SGD%NSPEC), D(THIS%SGD%NSPEC)

    REAL, INTENT(IN) :: U10ABS, U10DIR

    REAL, INTENT(IN) :: DT


    INTEGER :: IK, ITH, ISP, NK, NTH

    REAL :: ALPHA, BETA, CGI, CELLSIZE, SPECI, SFC

    LOGICAL :: CELLOBSTRUCTED

    REAL :: TH, PSI

#ifdef W3_S

    INTEGER, SAVE           :: IENT   = 0

#endif


#ifdef W3_S

    CALL strace (ient, 'UOST_SOURCETERM_COMPUTE_LD')

#endif


    s = 0

    d = 0


    cellobstructed =  this%GRD%UOST_LCL_OBSTRUCTED(ix, iy)

    IF (.NOT. cellobstructed) RETURN


    nk = this%SGD%NK

    nth = this%SGD%NTH


    DO ik = 1,nk

      cgi = cg(ik)

      DO ith = 1,nth


        !  Getting alpha and beta for local dissipation

        alpha = this%GRD%UOSTLOCALALPHA(ix, iy, ik, ith)/this%GRD%UOSTABMULTFACTOR

        alpha = max(min(alpha*this%GRD%UOSTLOCALFACTOR, 1.), 0.)

        beta = this%GRD%UOSTLOCALBETA(ix, iy, ik, ith)/this%GRD%UOSTABMULTFACTOR

        beta = max(min(beta*this%GRD%UOSTLOCALFACTOR, 1.), 0.)


        IF (alpha .EQ. 1) cycle


        !  Getting the size of the cell along direction ith

        cellsize = this%GRD%UOSTCELLSIZE(ix, iy, ith)*this%GRD%UOSTCELLSIZEFACTOR


        isp = ith + (ik-1)*nth

        speci = spec(isp)


        th = this%SGD%TH(ith)

        CALL compute_reduction_psi(u10abs, u10dir, cg(ik), th, dt, psi)


        IF (beta > 0.09) THEN

          !  Computing the local dissipation for a partially obstructed cell

          sfc = - cgi/cellsize * (1 - beta)/beta

        ELSE

          !  The cell is almost completely obstructed.

          !  Dissipating the energy almost completely.

          sfc = - cgi/cellsize * this%GAMMAUP

        ENDIF


        s(isp) = sfc * speci * psi

        d(isp) = sfc * psi

      ENDDO

    ENDDO


  END SUBROUTINE uost_sourceterm_compute_ld


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

  SUBROUTINE uost_sourceterm_compute_se(THIS, IX, IY, SPEC, CG, DT, U10ABS, U10DIR, S, D)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |          Lorenzo Mentaschi        |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         01-Oct-2018 |

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

    !/

    !/    Aug-2018 : Origination.                        ( version 6.07 )

    !/

    !  1. Purpose : Method of the class UOST_SOURCETERM.

    !               Computation of the shadow dissipation of the spectrum

    !  2. Parameters :

    !

    !     Parameter list

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

    !     THIS: UOST_SOURCETERM      instance of UOST_SOURCETERM passed to the method

    !                           (compulsory in oo programming)

    !     IX, IY: Integer       coordinates of the actual cell

    !     SPEC: real            input spectrum

    !     CG: real              group velocity

    !     DT: real              time step

    !     U10ABS: real          absolute value of U10

    !     U10DIR: real          direction of U10

    !     S: real               source term

    !     D: real               differential of the source term over the spectrum

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

    !

    !  3. Called by :

    !

    !      Name                   Type  Module     Description

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

    !      UOST_SOURCETERM_COMPUTE     Subr. W3UOSTMD   Computing the source term

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

    !

    !  4. Source code :

    !

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

    IMPLICIT NONE


    CLASS(uost_sourceterm), INTENT(INOUT), TARGET :: THIS

    INTEGER, INTENT(IN) :: IX, IY

    REAL, INTENT(IN) :: SPEC(THIS%SGD%NSPEC), CG(THIS%SGD%NK)

    REAL, INTENT(OUT) :: S(THIS%SGD%NSPEC), D(THIS%SGD%NSPEC)

    REAL, INTENT(IN) :: U10ABS, U10DIR

    REAL, INTENT(IN) :: DT


    INTEGER :: IK, ITH, IS

    REAL :: CGI, SPECI, SFC, CELLSIZE, &

         SFCLEFT, SFCRIGHT, SFCCENTER, THDIAG, CGDIAG, &

         ALPHASH, BETASH, GAMMMA, GG

    INTEGER :: N = 8, ithdiag, isp, nk, nth, nx, ny

    LOGICAL :: CELLOBSTRUCTED

    REAL :: TH, PSI

#ifdef W3_S

    INTEGER, SAVE           :: IENT   = 0

#endif


#ifdef W3_S

    CALL strace (ient, 'UOST_SOURCETERM_COMPUTE_SE')

#endif


    s = 0

    d = 0


    nk = this%SGD%NK

    nth = this%SGD%NTH

    nx = this%GRD%NX

    ny = this%GRD%NY


    IF ((ix .EQ. 1) .OR. (ix .EQ. nx) .OR. (iy .EQ. 1) .OR. (iy .EQ. ny)) RETURN


    cellobstructed = this%GRD%UOST_SHD_OBSTRUCTED(ix, iy)

    IF (.NOT. cellobstructed) RETURN


    DO ik=1,nk

      DO ith=1,nth


        !  Getting alpha and beta of the shadow

        alphash = this%GRD%UOSTSHADOWALPHA(ix, iy, ik, ith)/this%GRD%UOSTABMULTFACTOR

        alphash = max(min(alphash*this%GRD%UOSTSHADOWFACTOR, 1.), 0.)

        betash = this%GRD%UOSTSHADOWBETA(ix, iy, ik, ith)/this%GRD%UOSTABMULTFACTOR

        betash = max(min(betash*this%GRD%UOSTSHADOWFACTOR, 1.), 0.)


        IF (alphash .EQ. 1) cycle


        !  Getting the size of the cell along direction ith

        cellsize = this%GRD%UOSTCELLSIZE(ix, iy, ith)*this%GRD%UOSTCELLSIZEFACTOR


        cgi = cg(ik)


        gg = cgi/cellsize


        IF (alphash > 0.2) THEN

          !  Computing the shadow gamma coefficient for a partially obstructed cell

          gammma = (betash/alphash - 1)

        ELSE

          !  Alpha is small. The shadow dissipates the energy almost completely

          gammma = this%GAMMADOWN

        ENDIF


        th = this%SGD%TH(ith)

        !  Computing the reduction psi related with the wind component of the spectrum

        CALL compute_reduction_psi(u10abs, u10dir, cg(ik), th, dt, psi)


        sfc = - gg*gammma


        isp = ith + (ik-1)*nth

        speci = spec(isp)

        s(isp) = sfc * speci * psi

        d(isp) = sfc * psi

      ENDDO

    ENDDO

  END SUBROUTINE uost_sourceterm_compute_se


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

  SUBROUTINE uost_sourceterm_compute(THIS, IX, IY, SPEC, CG, DT, U10ABS, U10DIR, S, D)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |          Lorenzo Mentaschi        |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         01-Oct-2018 |

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

    !/

    !/    Aug-2018 : Origination.                        ( version 6.07 )

    !/

    !  1. Purpose : Method of the class UOST_SOURCETERM.

    !               Computation of the source term

    !  2. Parameters :

    !

    !     Parameter list

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

    !     THIS: UOST_SOURCETERM      instance of UOST_SOURCETERM passed to the method

    !                           (compulsory in oo programming)

    !     IX, IY: Integer       coordinates of the actual cell

    !     SPEC: real            input spectrum

    !     CG: real              group velocity

    !     DT: real              time step

    !     U10ABS: real          absolute value of U10

    !     U10DIR: real          direction of U10

    !     S: real               source term

    !     D: real               differential of the source term over the spectrum

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

    !

    !  3. Called by :

    !

    !      Name                   Type  Module     Description

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

    !      UOST_SRCTRMCOMPUTE     Subr. W3UOSTMD   Computing the source term

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

    !

    !  4. Source code :

    !

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

    IMPLICIT NONE


    CLASS(uost_sourceterm), INTENT(INOUT) :: THIS

    INTEGER, INTENT(IN) :: IX, IY

    REAL, INTENT(IN) :: SPEC(THIS%SGD%NSPEC), CG(THIS%SGD%NK)

    REAL, INTENT(IN) :: DT

    REAL, INTENT(IN) :: U10ABS, U10DIR

    REAL, INTENT(OUT) :: S(THIS%SGD%NSPEC), D(THIS%SGD%NSPEC)


    REAL :: S_LD(THIS%SGD%NSPEC), S_SE(THIS%SGD%NSPEC)

    REAL :: D_LD(THIS%SGD%NSPEC), D_SE(THIS%SGD%NSPEC)

#ifdef W3_S

    INTEGER, SAVE           :: IENT   = 0

#endif


#ifdef W3_S

    CALL strace (ient, 'UOST_SOURCETERM_COMPUTE')

#endif


    IF (.NOT. this%GRD%UOSTENABLED) THEN

      s = 0

      RETURN

    ENDIF


    !  Initializing the LD and SE components

    s_ld = 0

    s_se = 0

    !  Local dissipation

    CALL this%COMPUTE_LD(ix, iy, spec, cg, dt, u10abs, u10dir, s_ld, d_ld)

    !  Shadow effect

    CALL this%COMPUTE_SE(ix, iy, spec, cg, dt, u10abs, u10dir, s_se, d_se)

    s = s_ld + s_se

    d = d_ld + d_se

  END SUBROUTINE uost_sourceterm_compute


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


END MODULE w3uostmd


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