#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, IAPROC, NAPOUT
      USE W3SERVMD, ONLY: EXTCDE
!/S      USE W3SERVMD, ONLY: STRACE

      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
!/S      INTEGER, SAVE           :: IENT   = 0

!/S      CALL STRACE (IENT, 'UOST_INITGRID') 

      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
        IF (IAPROC .EQ. NAPOUT) THEN
          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
      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
!/S      INTEGER, SAVE           :: IENT   = 0

!/S      CALL STRACE (IENT, 'UOST_SETGRID') 

      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)
!/S      INTEGER, SAVE           :: IENT   = 0

!/S      CALL STRACE (IENT, 'UOST_SRCTRMCOMPUTE') 

      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
!/S      INTEGER, SAVE           :: IENT   = 0

!/S      CALL STRACE (IENT, 'LOAD_ALPHABETA') 
    

      ! 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  
      IF (IAPROC .EQ. NAPOUT) THEN
        WRITE(NDSO,*)'FILE '//FILENAME(:J)//' FOUND.'// &  
               'LOADING UOST SETTINGS FOR GRID '//GRD%GNAME
      ENDIF
  
      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  
      IF (IAPROC .EQ. NAPOUT) THEN
        WRITE(NDSO,*)'FILE '//FILENAME(:J)//' FOUND.'//& 
                'LOADING UOST SETTINGS FOR GRID '//GRD%GNAME
      ENDIF
  
      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(:)
!/S      INTEGER, SAVE           :: IENT   = 0

!/S      CALL STRACE (IENT, 'LOAD_ALPHABETA_FROMFILE') 
    
      !  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
!/S      INTEGER, SAVE           :: IENT   = 0

!/S      CALL STRACE (IENT, 'UOST_SOURCETERM_SETGRID') 
    
      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
!/S      INTEGER, SAVE           :: IENT   = 0

!/S      CALL STRACE (IENT, 'COMPUTE_REDUCTION_PSI') 

      ! 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
!/S      INTEGER, SAVE           :: IENT   = 0

!/S      CALL STRACE (IENT, 'UOST_SOURCETERM_COMPUTE_LD') 

      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
!/S      INTEGER, SAVE           :: IENT   = 0

!/S      CALL STRACE (IENT, 'UOST_SOURCETERM_COMPUTE_SE') 
 
      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)
!/S      INTEGER, SAVE           :: IENT   = 0

!/S      CALL STRACE (IENT, 'UOST_SOURCETERM_COMPUTE') 
      
      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  

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