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