!> @file
!> @brief Parmeterization of the unresolved obstacles.
!>
!> @author Lorenzo Mentaschi
!> @date 08-Oct-2018
!>
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
!>
!> @brief Parmeterization of the unresolved obstacles.
!>
!> @author Lorenzo Mentaschi
!> @date 08-Oct-2018
!>
!> @copyright Copyright 2009-2022 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.
!>
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
!/ ------------------------------------------------------------------- /
!>
!> @brief Allocate the UOST variables for a given grid, and load
!> them from file.
!>
!> @param[in] IGRID Id of the grid being initialized.
!> @param[in] FILELOCAL File from where the alpha/beta coefficient
!> for the local dissipation are loaded.
!> @param[in] FILESHADOW File from where the alpha/beta coefficient
!> for the shadow effect are loaded.
!> @param[in] LOCALFACTOR Adjustment parameter for the local
!> dissipation alpha and beta.
!> @param[in] SHADOWFACTOR Adjustment parameter for the shadow
!> dissipation alpha and beta.
!> @author Lorenzo Mentaschi
!> @date 01-Oct-2018
!>
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
!/ ------------------------------------------------------------------- /
!>
!> @brief Sets the current grid in the sourceterm object.
!>
!> @param[in] IGRID Id of the actual grid.
!>
!> @author Lorenzo Mentaschi
!> @date 01-Oct-2018
!>
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
!/ ------------------------------------------------------------------- /
!/ ------------------------------------------------------------------- /
!>
!> @brief Estimates the UOST source term for a give spectrum.
!>
!> @param[in] IX
!> @param[in] IY
!> @param[in] SPEC
!> @param[in] CG
!> @param[in] DT
!> @param[in] U10ABS
!> @param[in] U10DIR
!> @param[out] S
!> @param[out] D
!>
!> @author Lorenzo Mentaschi
!> @date 01-Oct-2018
!>
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
!/ ------------------------------------------------------------------- /
!>
!> @brief Loads local and shadow alpha and beta from files.
!>
!> @param[inout] GRD Object representing the spatial grid to be loaded.
!> @param[inout] SRD Object representing the current spectral grid.
!> @param[inout] FILEUNIT Unit id of the input files.
!>
!> @author Lorenzo Mentaschi
!> @date 01-Oct-2018
!>
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
!/ ------------------------------------------------------------------- /
!>
!> @brief Loads alpha and beta from a single obstructions file.
!>
!> @param[in] FILEUNIT Unit of the file to be opened.
!> @param[in] FILENAME Name of the file.
!> @param[in] NX Size of spatial grid X dim.
!> @param[in] NY Size of spatial grid Y dim.
!> @param[in] NK Size of spectral grid K dim.
!> @param[in] NTH Size of spectral grid TH dim.
!> @param[in] MULTFACTOR Multiplication factor for alpha and beta.
!> @param[inout] ALPHAMTX Loaded alpha spatial/spectral matrices.
!> @param[inout] BETAMTX Loaded beta spatial/spectral matrices.
!> @param[inout] CELLSIZE Cell size for each spectral direction,
!> also loaded from the file.
!> @param[inout] ISOBSTRUCTED Matrix of logicals, indicating for each cell
!> if it is obstructed or not.
!>
!> @author Lorenzo Mentaschi
!> @date 01-Oct-2018
!>
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
!/ ------------------------------------------------------------------- /
!>
!> @brief Method of the class UOST_SOURCETERM, to set the
!> actual spatial and spectral grid.
!>
!> @param[inout] THIS
!> @param[in] GRD Object representing the spatial grid to be loaded.
!> @param[in] SGD Object representing the current spectral grid.
!>
!> @author Lorenzo Mentaschi
!> @date 01-Oct-2018
!>
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
!/ ------------------------------------------------------------------- /
!>
!> @brief In conditions of wind sea, the effect of the unresolved
!> obstacles is reduced.
!>
!> @details Here a reduction psi is computed, as a function of the
!> wave age.
!>
!> @param[in] U10ABS Absolute value of U10.
!> @param[in] U10DIR Direction of U10.
!> @param[in] CGABS Absolute value of the group velocity.
!> @param[in] CGDIR Direction of the group velocity.
!> @param[in] DT Time step.
!> @param[out] PSI Output psi factor.
!>
!> @author Lorenzo Mentaschi
!> @date 01-Oct-2018
!>
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
!/ ------------------------------------------------------------------- /
!>
!> @brief Method of the class UOST_SOURCETERM.
!>
!> @details Computation of the local dissipation of the spectrum.
!>
!> @param[inout] THIS Instance of UOST_SOURCETERM passed to the method.
!> @param[in] IX X coordinates of actual call.
!> @param[in] IY Y coordinates of actual call.
!> @param[in] SPEC Input spectrum.
!> @param[in] CG Group velocity.
!> @param[in] DT Time step.
!> @param[in] U10ABS Absolute value of U10.
!> @param[in] U10DIR Direction of U10.
!> @param[out] S Source term.
!> @param[out] D Differential of the source term over the spectrum.
!>
!> @author Lorenzo Mentaschi
!> @date 01-Oct-2018
!>
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
!/ ------------------------------------------------------------------- /
!>
!> @brief Method of the class UOST_SOURCETERM.
!>
!> @details Computation of the shadow dissipation of the spectrum.
!>
!> @param[inout] THIS Instance of UOST_SOURCETERM passed to the method.
!> @param[in] IX X coordinates of actual call.
!> @param[in] IY Y coordinates of actual call.
!> @param[in] SPEC Input spectrum.
!> @param[in] CG Group velocity.
!> @param[in] DT Time step.
!> @param[in] U10ABS Absolute value of U10.
!> @param[in] U10DIR Direction of U10.
!> @param[out] S Source term.
!> @param[out] D Differential of the source term over the spectrum.
!>
!> @author Lorenzo Mentaschi
!> @date 01-Oct-2018
!>
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
!/ ------------------------------------------------------------------- /
!>
!> @brief Method of the class UOST_SOURCETERM.
!>
!> @details Computation of the source term.
!>
!> @param[inout] THIS Instance of UOST_SOURCETERM passed to the method.
!> @param[in] IX X coordinates of actual call.
!> @param[in] IY Y coordinates of actual call.
!> @param[in] SPEC Input spectrum.
!> @param[in] CG Group velocity.
!> @param[in] DT Time step.
!> @param[in] U10ABS Absolute value of U10.
!> @param[in] U10DIR Direction of U10.
!> @param[out] S Source term.
!> @param[out] D Differential of the source term over the spectrum.
!>
!> @author Lorenzo Mentaschi
!> @date 01-Oct-2018
!>
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
!/ ------------------------------------------------------------------- /