!> @file
!> @brief Contains the program W3UPRSTR.
!>
!> @author Stelios Flampouris @date 16-Feb-2017

#include "w3macros.h"
!/ ------------------------------------------------------------------- /
!>
!> @brief Update restart files based on Hs from DA.
!>
!> @details Update the WAVEWATCH III restart files based on the significant
!>  wave height analysis from any data assimilation system.
!>
!>  The W3UPRSTR is the intermediator between the background WW3
!>  and the analysis of the wave field, it modifies the original restart
!>  file according to the analysis.
!>  For the wave modeling and DA, the ww3_uprstr program applies the
!>  operator from the diagnostic to the prognostic variable.
!>
!> @author Stelios Flampouris @date 16-Feb-2017
!>
PROGRAM W3UPRSTR
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |         Stelios Flampouris        |
  !/                  |                        FORTRAN 90 |
  !/                  | First version:        16-Feb-2017 |
  !/                  +-----------------------------------+
  !/
  !/    16-Feb-2017 :  Original Code                       ( version 6.03 )
  !/    07-Jun-2017 :  Change of the Core
  !/    07-Jul-2017 :  Clean the code, add some flexibility, etc
  !/    04-Sep-2017 :  Simplified the code, take out a significant part of the
  !/                   flexibility (The code is still available at SVN/UpRest)
  !/    15-Sep-2017 :  Version 0.65                        ( version 6.03 )
  !/    01-Oct-2018 :  Fixes to preserve spectral energy correctly
  !/                   (Andy Saulter)                      ( version 6.06 )
  !/    17-Oct-2018 :  Version 0.95                        ( version 6.06 )
  !/                   Simplified the code, remove some user unfriendly
  !/                   options, add reg test ta1, add logical checks,
  !/                   unified the operator, add/update the documentation.
  !/    05-Oct-2019 :  Added UPD5 and UPD6 options, plus logic for running
  !/                   with SMC grids (Andy Saulter)       ( version 6.07 )
  !/    01-Nov-2019 :  UPD5 and UPD6 use wind data either from anl.XXX file
  !/                   or from restart under WRST switch (Andy Saulter)
  !/    06-Oct-2020 :  Added namelist input options        ( version 7.11 )
  !/    06-May-2021 :  Use SMCTYPE and FSWND for SMC grid. ( version 7.13 )
  !/
  !/    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 :
  !
  !     Update the WAVEWATCH III restart files based on the significant
  !     wave height analysis from any data assimilation system.
  !
  !  2. Method :
  !
  !  2.1. General:
  !     The W3UPRSTR is the intermediator between the background WW3
  !     and the analysis of the wave field, it modifies the original restart
  !     file according to the analysis.
  !     For the wave modeling and DA, the ww3_uprstr program applies the
  !     operator from the diagnostic to the prognostic variable.
  !
  !     See the chart below:
  !
  !                       +-------------------+
  !                       | WW3 Background Run|
  !                       +-------------------+
  !    +--------------+              |                   +-----+
  !    | Restart File | <------------|-----------------> |  Hs |
  !    +--------------+                                  +-----+
  !            |                                            |
  !            |                                            |
  !            |               +---------------+         +-----+
  !            |               |Hs Observations|-------> | D.A.|
  !            |               +---------------+         +-----+
  !            |                                            |
  !            |                  +----------+              |
  !            +----------------> | W3UPRSTR |<-/Analysis/--+
  !                               +----------+
  !                                     |
  !                          +----------------------+
  !                          | Updated Restart File |
  !                          +----------------------+
  !
  !     A. The WW3 Background Run has to provide two files:
  !         i. The field of Hs (for NCEP at grib2 format) and
  !        ii. the restart.ww3, at the WW3 format for restart files.
  !     Both of them, at the moment of the assimilation (Nevertheless, the WW3
  !     restart reader will fail when the timestamps are not identical).
  !
  !     B. The DA module produces the analysis and/or the difference (%) of
  !     the analysis from the first guess of Hs in the space of model and
  !     exports the results.
  !
  !     C. The algorithm
  !     The Hs correction is redistributed to each frequency and direction.
  !
  !     1. The W3UPRSTR imports: i. the restart.ww3,
  !                             ii. the analysis file and
  !                            iii. the input file: ww3_uprstr.inp, details at 2.2.2.i.
  !
  !     2. The W3UPRSTR updates the restart file according to the option at
  !     ww3_uprstr UPD[N]
  !     Note: With the version 6.06 some options have been removed, but the naming
  !     is consistent with the original version.
  !
  !     3. W3UPRSTR exports the updated spectrum in the same format as the
  !     restart.ww3. The name of the output file is: restart001.ww3 and it has to
  !     be renamed "restart.ww3" for the initialization of the next prediction
  !     cycle.
  !
  !     E. The user runs WW3 with the analysis restart file.
  !
  !     2.2. How to use ww3_uprstr
  !     The ww3_uperstr is one of the WW3 auxilary programs, therefore it works in
  !     a very similar way as the other auxilary programs.
  !
  !     A. To compile:
  !
  !     ww3_uprstr is included in the make_makefile.sh, to compile:
  !     $ ./w3_make ww3_uprstr
  !     or
  !     $ ./w3_make
  !
  !     And the executable "ww3_uprstr" will appear at [...]/model/exe/
  !
  !     B. To run:
  !     At the computational path:
  !        > ${EXE}/ww3_uprstr
  !     And it should run if the input files are at ./
  !
  !     C. Input Files:
  !
  !        i. ww3_uprstr.inp
  !     It includes some limited information for running the program:
  !
  ! -------------------------------------------------------------------- $
  ! WAVEWATCH III Update Restart input file                              $
  ! -------------------------------------------------------------------- $
  !
  ! Time of Assimilation ----------------------------------------------- $
  ! - Starting time in yyyymmdd hhmmss format.
  !
  ! This is the assimilation starting time and has to be the same with
  ! the time at the restart.ww3.
  !
  !   19680607 120000
  !
  ! Choose algorithm to update restart file
  !  UPDN for the Nth approach
  !  The UPDN*, with N<2 the same correction factor is applied at all the grid points
  !   UPD0C:: ELIMINATED
  !   UPDOF:: Option 0F  All the spectra are updated with a constant
  !           fac=HsAnl/HsBckg.
  !           Expected input: PRCNTG, as defined at fac
  !   UPD1 :: ELIMINATED
  !   UPDN, with N>1 each gridpoint has its own update factor.
  !   UPD2 :: Option 2   The fac(x,y,frq,theta), is calculated at each grid point
  !           according to the ratio of HsBckg and HsAnl (squared to preseve energy)
  !           Expected input: the Analysis field, grbtxt format
  !   UPD3 :: Option 3   The update factor is a surface with the shape of
  !           the background spectrum.
  !           Expected input: the Analysis field, grbtxt format
  !   UPD4 :: [NOT INCLUDED in this Version, Just keeping the spot]
  !           Option 4  The generalization of the UPD3. The update factor
  !           is the sum of surfaces which are applied on the background
  !           spectrum.
  !           The algorithm requires the mapping of each partition on the
  !           individual spectra; the map is used to determine the weighting
  !           surfaces.
  !           Expected input: the Analysis field, grbtxt format and the
  !           functions(frq,theta) of the update to be applied.
  !   UPD5 :: Option 5   Corrections are calculated as per UPD2 but are
  !           applied to wind-sea parts of the spectrum only when wind-sea
  !           is the dominant component, otherwise the whole spectrum is
  !           corrected
  !           Expected input: the Analysis Hs field plus background wind speed
  !                           and direction
  !   UPD6 :: Option 6   Corrections are calculated as per UPD5 but wind-sea
  !           components are also shifted in frequency space using Toba (1973)
  !           Expected input: the Analysis Hs field plus background wind speed
  !                           and direction
  !
  ! PRCNTG is input for option UPD0F and is the correction factor
  ! applied  to all the gridpoints (e.g. 1.)
  !
  !   0.475
  !
  ! PRCNTG_CAP is global input for option UPD2 and UPD3 and it is a cap on
  ! the maximum SWH correction factor applied to all the gridpoints, as
  ! both a multiple or divisor (e.g. cap at 5.0 means SWHANL/SWHBKG<=5.0
  ! and SWHANL/SWHBKG>=0.2). The value given should not be less than 1.0
  !
  !  5.0
  !
  ! Name of the file with the SWH analysis from the DA system            $
  ! suffix .grbtxt for text out of grib2 file.                           $
  !
  !   anl.grbtxt
  !
  ! -------------------------------------------------------------------- $
  ! WAVEWATCH III EoF ww3_uprstr.inp
  ! -------------------------------------------------------------------- $
  !
  !   ii. Data files anl.XXX
  !
  ! FOR UPD2,3 and UPD5,6 with WRST switch
  ! USE THE grbtxt FORMAT, See Format E.
  !
  ! Format E.
  !        Text file created by wgrib2. This format is tested more extensively
  !        and currently the only format supported for anl.grbtxt.
  !
  !        NX NY
  !        VAL0001
  !        VAL0002
  !        ...
  !        VALNX*NY
  !
  ! IMPORTANT : All the regtests are with the format E. strongly recommended.
  ! The order of the values in .grbtxt, is assumed the same by
  ! default as the order of spectral data in the restart file.
  !
  ! NOTE: It is recommended to use UPD5,6 with the WRST switch enabled and
  ! using SWH analysis data only as per Format E. However, the code includes
  ! an option to run using a text file in which case:
  ! USE THE grbtxtws format below
  !
  ! Text file with following lines:
  !        NX NY
  !        SWH0001 WSPD0001 WDIR0001
  !        SWH0002 WSPD0002 WDIR0002
  !        ...
  !        SWHNX*NY WSPDNX*NY WDIRNX*NY
  !
  ! The order of the values in .grbtxt, is assumed the same by
  ! default as the order of spectral data in the restart file.
  ! Wind speeds and directions in the anl.XXX file are assumed to be
  ! in CARTESIAN (GRID U,V) CONVENTION
  !
  ! NOTE About Format: if you prefer a different format; there are several
  ! I/O subroutines ready, not included in the current version of the code,
  ! contact the prgmr to get access to the source code.
  !
  !        iii. restart.ww3
  ! The restart file as came out of the background run, the name has to be
  ! restart.ww3, but the name of the output depends on the mod_def.ww3, the
  ! ww3_uprstr follows its content (be careful with ovewriting).
  !
  !  3. Example
  !     Use the regression tests ww3_ta1
  !
  !  4. Parameters :
  !
  !     Local parameters.
  !     ----------------------------------------------------------------
  !
  !     ----------------------------------------------------------------
  !
  !  5. Subroutines used :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !      W3NMOD    Subr. W3GDATMD Set number of model.
  !      W3SETG    Subr.   Id.    Point to selected model.
  !      W3NDAT    Subr. W3WDATMD Set number of model for wave data.
  !      W3SETW    Subr.   Id.    Point to selected model for wave data.
  !      W3NINP    Subr. W3IDATMD Set number of grids/models.
  !      W3SETI    Subr.   Id.    Point to data structure.
  !      W3DIMI    Subr.   Id.    Set array sizes in data structure.
  !      W2NAUX    Subr. W3ADATMD Set number of model for aux data.
  !      W3SETA    Subr.   Id.    Point to selected model for aux data.
  !      ITRACE    Subr. W3SERVMD Subroutine tracing initialization.
  !      NEXTLN    Subr.   Id.    Get next line from input file.
  !      EXTCDE    Subr.   Id.    Abort program as graceful as possible.
  !      W3IOGR    Subr. W3IOGRMD Reading/writing model definition file.
  !      WAVNU1    Subr. W3DISPMD

  !     ----------------------------------------------------------------
  !  Internal Subroutines:
  !
  !   READ_GRBTXT
  !   WORKER
  !   SWH_RSRT_1p
  !   WRITEMATRIX
  !
  !  6. Called by :
  !
  !     None, stand-alone program.
  !
  !  7. Error messages :
  !
  !  8. Remarks:
  !
  !     7.1 Use the grbtxt format for correction and analysis files.
  !
  !     7.2 There are some variables not used but declared, it's for future
  !      development.
  !
  !  9. Structure :
  !
  !     ----------------------------------------------------
  !     1.    Set up data structures.
  !     2.    Read model defintion file with base model ( W3IOGR )
  !     3.    Import restart file                       ( W3IORS )
  !     4.    Import correction percentage              (  )
  !        OR Import the analysis field                 (  )
  !     5.    Apply correction to the restart file      (  )
  !     6.    Export the updated restart file           ( W3IORS )
  !     ----------------------------------------------------
  !
  !  10. Switches :
  !
  !     !/SHRD  Switch for shared / distributed memory architecture.
  !     !/T
  !     !/S     Enable subroutine tracing.
  !
  ! 11. Known Bugs
  !
  !     1. Fix the format for the output (NSDO) of non strings, e.g. for
  !     TIME.
  !
  ! 12. Source code :
  !
  !/
  USE W3GDATMD, ONLY: W3NMOD, W3SETG
  USE W3WDATMD, ONLY: W3NDAT, W3SETW
  USE W3ADATMD, ONLY: W3NAUX, W3SETA
  USE W3ODATMD, ONLY: W3NOUT, W3SETO
  USE W3IORSMD, ONLY: W3IORS
  USE W3SERVMD, ONLY: ITRACE, NEXTLN, EXTCDE
  USE W3IOGRMD, ONLY: W3IOGR
  USE W3DISPMD, ONLY: WAVNU1
  !
  USE W3GDATMD, ONLY: GNAME, NX, NY, MAPSTA, SIG, NK, NTH, NSEA,  &
       NSEAL, MAPSF, DMIN, ZB, DSIP, DTH, RSTYPE,  &
       GTYPE, SMCTYPE
#ifdef W3_SMC
  USE W3GDATMD, ONLY: FSWND
#endif
  USE W3WDATMD, ONLY: VA, TIME
  USE W3ADATMD, ONLY: NSEALM
  USE W3ODATMD, ONLY: IAPROC, NAPERR, NAPLOG, NDS, NAPOUT
  USE W3ODATMD, ONLY: NDSE, NDSO, NDST, IDOUT, FNMPRE
#ifdef W3_WRST
  USE W3IDATMD
#endif
  !
  USE W3NMLUPRSTRMD
  !
  IMPLICIT NONE
  !/
  !/ ------------------------------------------------------------------- /
  !  Local variables
  !/
  INTEGER                 :: NDSI, NDSM, NDSTRC, NTRACE, IERR, I, J
  CHARACTER               :: COMSTR*1
  !
  TYPE(NML_RESTART_T)     :: NML_RESTART
  TYPE(NML_UPDATE_T)      :: NML_UPDATE
  !
  !      REAL, ALLOCATABLE       :: BETAW(:)
  !      LOGICAL, ALLOCATABLE    :: MASK(:)
  LOGICAL                 :: anl_exists, CORWSEA, FLGNML
  INTEGER                 :: IMOD,  NDSEN, IX, IY, IK, ITH, &
       IXW, IYW
  REAL, ALLOCATABLE       :: UPDPRCNT(:,:),VATMP(:), HSIG(:,:),     &
       A(:), HS_ANAL(:,:), gues(:,:),         &
       HS_DIF(:,:),SWHANL(:,:), SWHBCKG(:,:), &
       SWHUPRSTR(:,:),VATMP_NORM(:),          &
       WSBCKG(:,:),WDRBCKG(:,:)
  INTEGER, ALLOCATABLE    :: VAMAPWS(:)
  REAL                    :: PRCNTG, PRCNTG_CAP, THRWSEA
  INTEGER                 :: ROWS, COLS, ISEA
  CHARACTER(128)          :: FLNMCOR, FLNMANL
  CHARACTER(16)           :: UPDPROC
  !     for howv
  REAL                    :: SWHTMP,SWHBCKG_1, SWHANL_1,            &
       DEPTH, WN, CG, ETOT, E1I,              &
       SWHTMP1,SUMVATMP, SWHBCKG_W, SWHBCKG_S
  REAL                    :: K
  CHARACTER(8), PARAMETER :: MYNAME='W3UPRSTR'
  LOGICAL                 :: SMCGRD = .FALSE.
  LOGICAL                 :: SMCWND = .FALSE.
  LOGICAL                 :: WRSTON = .FALSE.
  !/
  !/ ------------------------------------------------------------------- /
  !/
  !  1.  IO set-up.
  CALL W3NMOD ( 1, 6, 6 )
  CALL W3SETG ( 1, 6, 6 )
  CALL W3NDAT (    6, 6 )
  CALL W3SETW ( 1, 6, 6 )
  CALL W3NAUX (    6, 6 )
  CALL W3SETA ( 1, 6, 6 )
  CALL W3NOUT (    6, 6 )
  CALL W3SETO ( 1, 6, 6 )
#ifdef W3_WRST
  CALL W3NINP (    6, 6 )
  CALL W3SETI ( 1, 6, 6 )
#endif
  !
  NDSE   = 6
  NDSI   = 10
  NDSM   = 20
  !
  IAPROC = 1
  NAPOUT = 1
  NAPERR = 1
  IMOD   = 1
  NAPLOG = 1
  !
  NDSTRC =  6
  NTRACE = 10
  CALL ITRACE ( NDSTRC, NTRACE )
  !
  IF ( IAPROC .EQ. NAPERR ) THEN
    NDSEN  = NDSE
  ELSE
    NDSEN  = -1
  END IF
  !
  WRITE (NDSO,900)
  !
#ifdef W3_WRST
  !Compiling with WRST will allow access to options UPD5/6
  WRSTON = .TRUE.
  WRITE (NDSO,*) '*** UPRSTR will read wind from restart files'
#endif
  !/
  !/ ------------------------------------------------------------------- /
  !  2. Read the ww3_uprstr input data
  !
  ! process ww3_uprstr namelist
  !
  INQUIRE(FILE=TRIM(FNMPRE)//"ww3_uprstr.nml", EXIST=FLGNML)
  IF (FLGNML) THEN
    ! Read namelist
    CALL W3NMLUPRSTR (NDSI, TRIM(FNMPRE)//'ww3_uprstr.nml', NML_RESTART, &
         NML_UPDATE, IERR)
    READ(NML_RESTART%RESTARTTIME, *) TIME
    UPDPROC = NML_UPDATE%UPDPROC
    PRCNTG = NML_UPDATE%PRCNTG
    PRCNTG_CAP = NML_UPDATE%PRCNTGCAP
    THRWSEA = NML_UPDATE%THRWSEA
    FLNMANL = NML_UPDATE%FILE
  END IF
  !/
  ! otherwise read from the .inp file
  IF (.NOT. FLGNML) THEN
    J      = LEN_TRIM(FNMPRE)
    OPEN (NDSI,FILE=FNMPRE(:J)//'ww3_uprstr.inp',STATUS='OLD',       &
         ERR=800,IOSTAT=IERR)
    READ (NDSI,'(A)',END=801,ERR=802) COMSTR
    IF (COMSTR.EQ.' ') COMSTR = '$'
    WRITE (NDSO,901) COMSTR
    !
    CALL NEXTLN ( COMSTR , NDSI , NDSEN )
    READ (NDSI,*,END=2001,ERR=2002) TIME
    CALL NEXTLN ( COMSTR , NDSI , NDSEN )
    READ (NDSI,*,END=2001,ERR=2002) UPDPROC
    CALL NEXTLN ( COMSTR , NDSI , NDSEN )
    IF (UPDPROC .EQ. 'UPD0F') THEN
      READ (NDSI,*,END=2001,ERR=2002) PRCNTG
    ELSE
      IF ((UPDPROC .EQ. 'UPD2') .OR. (UPDPROC .EQ. 'UPD3')) THEN
        !         CALL NEXTLN ( COMSTR , NDSI , NDSEN )
        READ (NDSI,*,END=2001,ERR=2002) PRCNTG_CAP
#ifdef W3_F
        CALL NEXTLN ( COMSTR , NDSI , NDSEN )
        READ (NDSI,*,END=2001,ERR=2002) FLNMCOR
#endif
      ELSE
        READ (NDSI,*,END=2001,ERR=2002) PRCNTG_CAP, THRWSEA
      END IF
      CALL NEXTLN ( COMSTR , NDSI , NDSEN )
      READ (NDSI,*,END=2001,ERR=2002) FLNMANL
    END IF
  ENDIF
#ifdef W3_T
  WRITE (NDSO,*)' TIME: ',TIME
#endif
  !/
  !/ ------------------------------------------------------------------- /
  !  3.  Read model definition file.
  !/
  CALL W3IOGR ( 'READ', NDSM )
  NSEAL  = NSEA
  WRITE (NDSO,920) GNAME
  !/
#ifdef W3_SMC
  !! SMC grid option is activated if GTYPE .EQ. SMCTYPE.  JGLi06May2021
  IF( GTYPE .EQ. SMCTYPE ) SMCGRD = .TRUE.
  !! SMC sea-point wind option is activated if FSWND=.TRUE.  JGLi06May2021
  IF( FSWND )              SMCWND = .TRUE.
#endif
#ifdef W3_WRST
  ! Override SMCWND - at present restarts only store wind on
  ! a regular grid
  SMCWND = .FALSE.
#endif
#ifdef W3_SMC
  WRITE (NDSO,*) '*** UPRSTR set to work with SMC grid model'
#endif
  !/
  !/ ------------------------------------------------------------------- /
  !  4. Read restart file
  !/
#ifdef W3_WRST
  ! Set the wind flag to true when reading restart wind
  INFLAGS1(3) = .TRUE.
  CALL W3DIMI ( 1, 6, 6 ) !Needs to be called after w3iogr to have correct dimensions?
#endif
  CALL W3IORS ( 'READ', NDS(6), SIG(NK), IMOD )!
  IF ( IAPROC .EQ. NAPLOG ) THEN
    IF (RSTYPE.EQ.0.OR.RSTYPE.EQ.1.OR.RSTYPE.EQ.4) THEN
      WRITE (NDSO,1004) 'Terminating ww3_uprstr: The restart ' // &
           'file is not read'
      CALL EXTCDE ( 1 )
    ELSE
      WRITE (NDSO,1004) '  Updating Restart File'
      WRITE (NDSO,*) ' TIME: ',TIME
    END IF
  END IF
#ifdef W3_T
  WRITE (NDST,*), MYNAME,' : Exporting VA as imported to VA01.txt'
  CALL writeMatrix('VA01.txt', REAL(VA))
#endif
  !/
  !/ ------------------------------------------------------------------- /
  ! 5. Update restart spectra array according to the selected option
  !/
  SELECT CASE (UPDPROC)
    !/
    !/ ------------------------------------------------------------------- /
    ! UPD0F
    !/
  CASE ('UPD0F')
    WRITE (NDSO,902) 'UPD0F'
    WRITE (NDSO,1005) ' PRCNTG = ',PRCNTG
#ifdef W3_T
    ALLOCATE( VATMP  (SIZE(VA    ,1)                ))
    ALLOCATE( SWHANL (SIZE(MAPSTA,1), SIZE(MAPSTA,2)))
    ALLOCATE( SWHBCKG(SIZE(MAPSTA,1), SIZE(MAPSTA,2)))
#endif
    DO ISEA=1, NSEA, 1
#ifdef W3_T
      IX = MAPSF(ISEA,1)
      IY = MAPSF(ISEA,2)
      VATMP = VA(:,ISEA)
      CALL SWH_RSRT_1p (VATMP, ISEA, SWHBCKG_1)
      SWHBCKG(IY,IX)=SWHBCKG_1
#endif
      CALL UPDATE_VA(PRCNTG, VA(:,ISEA))
#ifdef W3_T
      VATMP = VA(:,ISEA)
      CALL SWH_RSRT_1p (VATMP, ISEA, SWHANL_1)
      SWHANL(IY,IX)=SWHANL_1
      WRITE (NDSO,*) ' =========== UPD0F Output ==========='
      WRITE (NDSO,*)'ISEA = ', ISEA,' PRCNTG = ',PRCNTG,    &
           ' SWHBCKG = ',SWHBCKG(IY,IX),           &
           ' SWHANL= ', SWHANL(IY,IX)
#endif
    END DO
#ifdef W3_T
    CALL writeMatrix('SWHBCKG_UPD0F.txt', REAL(SWHBCKG))
    CALL writeMatrix('SWHANL_UPD0F.txt' , REAL(SWHANL ))
    CALL writeMatrix('SWHRSTR_UPD0F.txt', REAL(SWHANL ))

    DEALLOCATE ( VATMP, SWHBCKG, SWHANL )
#endif
    !/
    !/ ------------------------------------------------------------------- /
    ! UPD2
    ! Apply a bulk correction to the wave spectrum at each grid cell based
    ! on the ratio of HsBckg and HsAnl
    !/
  CASE ('UPD2')
    WRITE (NDSO,902) 'UPD2'
    WRITE (NDSO,1005) ' PRCNTG_CAP = ',PRCNTG_CAP
    WRITE (NDSO,1006) ' Reading updated SWH from: ',trim(FLNMANL)
    !
    !        Array allocation
    ALLOCATE ( VATMP(SIZE(VA,1)))
    IF (.NOT. SMCGRD) THEN
      ALLOCATE( SWHBCKG(SIZE(MAPSTA,1), SIZE(MAPSTA,2)) )
      ALLOCATE( SWHANL(SIZE(MAPSTA,1), SIZE(MAPSTA,2)) )
#ifdef W3_SMC
    ELSE
      ALLOCATE( SWHBCKG(NSEA,1) )
      ALLOCATE( SWHANL(NSEA,1) )
#endif
    ENDIF
#ifdef W3_T
    IF (.NOT. SMCGRD) THEN
      ALLOCATE( SWHUPRSTR(SIZE(MAPSTA,1), SIZE(MAPSTA,2)) )
    ELSE
      ALLOCATE( SWHUPRSTR(NSEA,1) )
    ENDIF
#endif
    !
    !        Read additional Input: Analysis Field
    INQUIRE(FILE=FLNMANL, EXIST=anl_exists)
    IF (anl_exists) THEN
#ifdef W3_T
      WRITE (NDSO,*) 'shape(SWHANL)', shape(SWHANL)
#endif
      CALL READ_GRBTXT(SWHANL, FLNMANL, SMCGRD)
#ifdef W3_T
      CALL writeMatrix('SWHANL_IN.txt',SWHANL)
#endif
    ELSE
      WRITE (NDSO,*) trim(FLNMANL), ' does not exist, stopping...'
      DEALLOCATE( SWHANL,VATMP,SWHBCKG )
#ifdef W3_T
      DEALLOCATE( SWHUPRSTR )
#endif
      STOP
    END IF
    !
    !        Calculation
    DO ISEA=1, NSEA, 1
      IF (.NOT. SMCGRD) THEN
        IX = MAPSF(ISEA,1)
        IY = MAPSF(ISEA,2)
#ifdef W3_SMC
      ELSE
        IX = 1
        IY = ISEA
#endif
      ENDIF
      VATMP = VA(:,ISEA)
      CALL SWH_RSRT_1p (VATMP, ISEA, SWHBCKG_1)
      SWHBCKG(IY,IX)=SWHBCKG_1
      !
      IF ( SWHBCKG(IY,IX) > 0.01 .AND. SWHANL(IY,IX) > 0.01 ) THEN
        PRCNTG=(SWHANL(IY,IX)/SWHBCKG_1)
#ifdef W3_T
        WRITE (NDSO,*) 'ISEA = ', ISEA,' IX = ',IX,' IY = ', IY,         &
             ' PRCNTG = ',PRCNTG,' SWHBCKG = ',SWHBCKG(IY,IX), &
             ' SWHANL = ', SWHANL(IY,IX)
#endif
        CALL CHECK_PRCNTG (PRCNTG,PRCNTG_CAP)
        CALL UPDATE_VA(PRCNTG, VA(:,ISEA))
#ifdef W3_T
        CALL SWH_RSRT_1p (VA(:,ISEA), ISEA, SWHUPRSTR(IY,IX))
        WRITE (NDSO,*) ' =========== UPD2 Output ==========='
        WRITE (NDSO,*)'ISEA = ',ISEA,                            &
             'SWH_BCKG = ', SWHBCKG(IY,IX),             &
             'SWH_ANL = ', SWHANL(IY,IX),               &
             'PRCNTG = ', PRCNTG,                       &
             'SWH_RSTR = ',SWHUPRSTR(IY,IX)
#endif
      END IF
    END DO
#ifdef W3_T
    CALL writeMatrix('SWHBCKG_UPD2.txt', REAL(SWHBCKG  ))
    CALL writeMatrix('SWHANL_UPD2.txt' , REAL(SWHANL   ))
    CALL writeMatrix('SWHRSTR_UPD2.txt', REAL(SWHUPRSTR))
#endif
    !
    DEALLOCATE( SWHANL,VATMP,SWHBCKG )
#ifdef W3_T
    DEALLOCATE( SWHUPRSTR )
#endif
    !/
    !/ ------------------------------------------------------------------- /
    ! UPD3
    ! As per UPD2, but the update factor is a surface with the shape of the
    ! background spectrum
    !/
  CASE ('UPD3')
    WRITE (NDSO,902) 'UPD3'
    WRITE (NDSO,1005) ' PRCNTG_CAP = ',PRCNTG_CAP
    WRITE (NDSO,1006) ' Reading updated SWH from: ',trim(FLNMANL)
    !
    !        Array allocation
    ALLOCATE ( VATMP(SIZE(VA,1)))
    ALLOCATE ( VATMP_NORM(SIZE(VA,1)))
    ALLOCATE ( A(SIZE(VA,1)))
    IF (.NOT. SMCGRD) THEN
      ALLOCATE( SWHBCKG(SIZE(MAPSTA,1), SIZE(MAPSTA,2)) )
      ALLOCATE( SWHANL(SIZE(MAPSTA,1), SIZE(MAPSTA,2)) )
#ifdef W3_SMC
    ELSE
      ALLOCATE( SWHBCKG(NSEA,1) )
      ALLOCATE( SWHANL(NSEA,1) )
#endif
    ENDIF
#ifdef W3_T
    IF (.NOT. SMCGRD) THEN
      ALLOCATE( SWHUPRSTR(SIZE(MAPSTA,1), SIZE(MAPSTA,2)) )
    ELSE
      ALLOCATE( SWHUPRSTR(NSEA,1) )
    ENDIF
#endif
    !
    !        Read additional Input: Analysis Field
    INQUIRE(FILE=FLNMANL, EXIST=anl_exists)
    IF (anl_exists) THEN
#ifdef W3_T
      WRITE (NDSO,*) 'shape(SWHANL)', shape(SWHANL)
#endif
      CALL READ_GRBTXT(SWHANL, FLNMANL, SMCGRD)
#ifdef W3_T
      CALL writeMatrix('SWHANL_IN.txt',SWHANL)
#endif
    ELSE
      WRITE (NDSO,*) trim(FLNMANL), ' does not exist, stopping...'
      DEALLOCATE( SWHANL,VATMP,SWHBCKG,VATMP_NORM,A )
#ifdef W3_T
      DEALLOCATE( SWHUPRSTR )
#endif
      STOP
    END IF
    !
    !        Calculation
    DO ISEA=1, NSEA, 1
      IF (.NOT. SMCGRD) THEN
        IX = MAPSF(ISEA,1)
        IY = MAPSF(ISEA,2)
#ifdef W3_SMC
      ELSE
        IX = 1
        IY = ISEA
#endif
      ENDIF
      VATMP = VA(:,ISEA)
      CALL SWH_RSRT_1p (VATMP, ISEA, SWHBCKG_1)
      SWHBCKG(IY,IX)=SWHBCKG_1
      !
      IF ( SWHBCKG(IY,IX) > 0.01 .AND. SWHANL(IY,IX) > 0.01 ) THEN
        !Step 1.
        PRCNTG=(SWHANL(IY,IX)/SWHBCKG_1)
#ifdef W3_T
        WRITE (NDSO,*) ' =========== Step 1. ==========='
        WRITE (NDSO,*) ' ISEA = ', ISEA,' IX = ',IX,' IY = ', IY,         &
             ' PRCNTG = ',PRCNTG,' SWHBCKG = ',SWHBCKG(IY,IX), &
             ' SWHANL = ', SWHANL(IY,IX)
#endif
        CALL CHECK_PRCNTG(PRCNTG,PRCNTG_CAP)
        VATMP_NORM=VATMP/SUM(VATMP)
#ifdef W3_T
        WRITE (NDSO,*)' ISEA =', ISEA,' IX = ',IX,' IY = ', IY, &
             ' PRCNTG = ',PRCNTG,   &
             ' SWHBCKG = ',SWHBCKG(IY,IX), ' SWHANL = ', SWHANL(IY,IX)
#endif
        IF (PRCNTG > 1.) THEN
          A=PRCNTG**2*(1 + VATMP_NORM)
        ELSE
          A=PRCNTG**2*(1 - VATMP_NORM)
        END IF
        VATMP=A*VATMP
        CALL SWH_RSRT_1p (VATMP, ISEA, SWHTMP)
        PRCNTG=(SWHANL(IY,IX)/SWHTMP)
#ifdef W3_T
        SWHUPRSTR(IY,IX)=SWHTMP
        WRITE (NDSO,*) ' =========== Step 2. ==========='
        WRITE (NDSO,*)'ISEA = ', ISEA, ' PRCNTG = ',PRCNTG,        &
             ' SWHANL= ', SWHANL(IY,IX),                  &
             ' SWHUPRSTR(IY,IX) = ', SWHUPRSTR(IY,IX)
#endif
        CALL CHECK_PRCNTG (PRCNTG,PRCNTG_CAP)
        CALL UPDATE_VA(PRCNTG, VATMP)
        VA(:,ISEA)=VATMP
#ifdef W3_T
        CALL SWH_RSRT_1p (VATMP, ISEA, SWHTMP)
        SWHUPRSTR(IY,IX)=SWHTMP
        WRITE (NDSO,*) ' =========== UPD3 Output ==========='
        WRITE (NDSO,*)'ISEA = ',ISEA,'SWH_BCKG = ', SWHBCKG(IY,IX), &
             'SWH_ANL = ', SWHANL(IY,IX),                 &
             'SWH_RSTR = ',SWHUPRSTR(IY,IX)
#endif
      END IF
    END DO
#ifdef W3_T
    CALL writeMatrix('SWHBCKG_UPD3.txt', REAL(SWHBCKG))
    CALL writeMatrix('SWHANL_UPD3.txt' , REAL(SWHANL ))
    CALL writeMatrix('SWHRSTR_UPD3.txt', REAL(SWHUPRSTR))
#endif
    !
    DEALLOCATE( SWHANL,VATMP,SWHBCKG,VATMP_NORM,A )
#ifdef W3_T
    DEALLOCATE( SWHUPRSTR )
#endif
    !/
    !/ ------------------------------------------------------------------- /
    ! UPD5
    ! Corrects wind-sea only in wind dominated conditions - bulk correction
    ! The fac(x,y,frq,theta), is calculated at each grid point according to
    ! HsBckg and HsAnl
    !/
  CASE ('UPD5')
    WRITE (NDSO,902) 'UPD5'
    WRITE (NDSO,1005) ' PRCNTG_CAP = ',PRCNTG_CAP
    WRITE (NDSO,1005) ' THRWSEA = ',THRWSEA
    WRITE (NDSO,1006) ' Reading updated SWH from: ',trim(FLNMANL)
    ! Presently set hardwired THRWSEA energy threshold here
    ! not user defined in input file
    ! THRWSEA = 0.7
    !
    !        Array allocation
    ALLOCATE ( VATMP(SIZE(VA,1)))
    ALLOCATE ( VAMAPWS(SIZE(VA,1)))
    IF (.NOT. SMCGRD) THEN
      ! SWH arrays allocated using Y,X convention as per wgrib write
      ALLOCATE( SWHBCKG(SIZE(MAPSTA,1), SIZE(MAPSTA,2)) )
      ALLOCATE( SWHANL(SIZE(MAPSTA,1), SIZE(MAPSTA,2)) )
      ! Wind arrays allocated using X,Y convention as in w3idatmd
      ALLOCATE( WSBCKG(SIZE(MAPSTA,2), SIZE(MAPSTA,1)) )
      ALLOCATE( WDRBCKG(SIZE(MAPSTA,2), SIZE(MAPSTA,1)) )
#ifdef W3_SMC
    ELSE
      ALLOCATE( SWHBCKG(NSEA,1) )
      ALLOCATE( SWHANL(NSEA,1) )
      ! Use SMCWND to determine if reading a seapoint aray for wind
      IF( SMCWND ) THEN
        ALLOCATE( WSBCKG(NSEA,1) )
        ALLOCATE( WDRBCKG(NSEA,1) )
      ELSE
        ALLOCATE(WSBCKG(SIZE(MAPSTA,2), SIZE(MAPSTA,1)))
        ALLOCATE(WDRBCKG(SIZE(MAPSTA,2), SIZE(MAPSTA,1)))
      ENDIF
#endif
    ENDIF
#ifdef W3_T
    IF (.NOT. SMCGRD) THEN
      ALLOCATE( SWHUPRSTR(SIZE(MAPSTA,1), SIZE(MAPSTA,2)) )
    ELSE
      ALLOCATE( SWHUPRSTR(NSEA,1) )
    ENDIF
#endif
    !
    !        Read additional Input: Analysis Field
    INQUIRE(FILE=FLNMANL, EXIST=anl_exists)
    IF (anl_exists) THEN
#ifdef W3_T
      WRITE (NDSO,*) 'shape(SWHANL)', shape(SWHANL)
#endif
#ifdef W3_WRST
      ! For WRST switch read only corrected SWH
      ! Wind will have been read from the restart
      IF (WRSTON) THEN
        CALL READ_GRBTXT(SWHANL, FLNMANL, SMCGRD)
      ELSE
#endif
        CALL READ_GRBTXTWS(SWHANL,WSBCKG,WDRBCKG,FLNMANL,SMCGRD)
#ifdef W3_WRST
      ENDIF
#endif
#ifdef W3_T
      CALL writeMatrix('SWHANL_IN.txt',SWHANL)
#endif
    ELSE
      WRITE (NDSO,*) trim(FLNMANL), ' does not exist, stopping...'
      DEALLOCATE( SWHANL,VATMP,SWHBCKG,VAMAPWS,WSBCKG,WDRBCKG )
#ifdef W3_T
      DEALLOCATE( SWHUPRSTR )
#endif
      STOP
    END IF
    !
#ifdef W3_WRST
    !Calculate wind speed and direction values from u,v..
    !..using cartesian direction convention
    !At present assume only needed for data read from restart
    CALL UVTOCART(WXNwrst,WYNwrst,WSBCKG,WDRBCKG,SMCWND)
#endif
    !
    !        Calculation
    DO ISEA=1, NSEA, 1
      IF (.NOT. SMCGRD) THEN
        IX = MAPSF(ISEA,1)
        IY = MAPSF(ISEA,2)
        IXW = IX
        IYW = IY
#ifdef W3_SMC
      ELSE
        IX = 1
        IY = ISEA
        IF( SMCWND ) THEN
          ! Wind arrays allocated using (X,Y) convention for regular grids
          ! but overriding here for the SMC grid which are always defined
          ! as (NSEA,1) by switching the IY and IX dimension values around
          IXW = IY
          IYW = IX
        ELSE
          IXW = MAPSF(ISEA,1)
          IYW= MAPSF(ISEA,2)
        ENDIF
#endif
      ENDIF
      VATMP = VA(:,ISEA)
      CALL SWH_RSRT_1pw (VATMP, WSBCKG(IXW,IYW), WDRBCKG(IXW,IYW), ISEA, &
           SWHBCKG_1, SWHBCKG_W, SWHBCKG_S, VAMAPWS)
      SWHBCKG(IY,IX)=SWHBCKG_1
      !
      IF ( SWHBCKG(IY,IX) > 0.01 .AND. SWHANL(IY,IX) > 0.01 ) THEN
        ! If wind-sea is dominant energy component apply correction to
        ! wind-sea part only
        IF ( (SWHBCKG_W / SWHBCKG_1)**2.0 > THRWSEA ) THEN
          ! Apply spectrum updates to wind-sea bins only
          PRCNTG=SQRT((SWHANL(IY,IX)**2.0-SWHBCKG_S**2.0)/SWHBCKG_W**2.0)
          CALL CHECK_PRCNTG(PRCNTG,PRCNTG_CAP)
          CALL UPDTWSPEC(VATMP, PRCNTG, VAMAPWS)
          ! else correct the whole spectrum as for UPD2
        ELSE
          PRCNTG=(SWHANL(IY,IX)/SWHBCKG_1)
          CALL CHECK_PRCNTG(PRCNTG,PRCNTG_CAP)
          CALL UPDATE_VA(PRCNTG,VATMP)
        END IF
#ifdef W3_T
        WRITE (NDSO,*) 'ISEA = ', ISEA,' IX = ',IX,' IY = ', IY,         &
             ' PRCNTG = ',PRCNTG,' SWHBCKG = ',SWHBCKG(IY,IX), &
             ' SWHANL = ', SWHANL(IY,IX)
#endif
        VA(:,ISEA)=VATMP
#ifdef W3_T
        CALL SWH_RSRT_1p (VATMP, ISEA, SWHTMP)
        SWHUPRSTR(IY,IX)=SWHTMP
        WRITE (NDSO,*) ' =========== UPD5 Output ==========='
        WRITE (NDSO,*)'ISEA = ',ISEA,'SWH_BCKG = ', SWHBCKG(IY,IX), &
             'SWH_ANL = ', SWHANL(IY,IX),                 &
             'SWH_RSTR = ',SWHUPRSTR(IY,IX)
#endif
      END IF
    END DO
#ifdef W3_T
    CALL writeMatrix('SWHBCKG_UPD5.txt', REAL(SWHBCKG  ))
    CALL writeMatrix('SWHANL_UPD5.txt' , REAL(SWHANL   ))
    CALL writeMatrix('SWHRSTR_UPD5.txt', REAL(SWHUPRSTR))
#endif
    !
    DEALLOCATE( SWHANL,VATMP,SWHBCKG,VAMAPWS,WSBCKG,WDRBCKG )
#ifdef W3_T
    DEALLOCATE( SWHUPRSTR )
#endif
    !/
    !/ ------------------------------------------------------------------- /
    ! UPD6
    ! Hybrid of Lionello et al. and Kohno methods
    ! Corrects wind-sea only in wind dominated conditions - including fp shift
    ! The fac(x,y,frq,theta), is calculated at each grid point according to
    ! HsBckg and HsAnl
    !/
  CASE ('UPD6')
    WRITE (NDSO,902) 'UPD6'
    WRITE (NDSO,1005) ' PRCNTG_CAP = ',PRCNTG_CAP
    WRITE (NDSO,1005) ' THRWSEA = ',THRWSEA
    WRITE (NDSO,1006) ' Reading updated SWH from: ',trim(FLNMANL)
    ! Presently set hardwired CORWSEA logical and THRWSEA energy
    ! thresholds here, not user defined in input file
    CORWSEA = .FALSE.
    !THRWSEA = 0.7
    !
    !        Array allocation
    ALLOCATE ( VATMP(SIZE(VA,1)))
    ALLOCATE ( VAMAPWS(SIZE(VA,1)))
    IF (.NOT. SMCGRD) THEN
      ! SWH arrays allocated using Y,X convention as per wgrib write
      ALLOCATE( SWHBCKG(SIZE(MAPSTA,1), SIZE(MAPSTA,2)) )
      ALLOCATE( SWHANL(SIZE(MAPSTA,1), SIZE(MAPSTA,2)) )
      ! Wind arrays allocated using X,Y convention as in w3idatmd
      ALLOCATE( WSBCKG(SIZE(MAPSTA,2), SIZE(MAPSTA,1)) )
      ALLOCATE( WDRBCKG(SIZE(MAPSTA,2), SIZE(MAPSTA,1)) )
#ifdef W3_SMC
    ELSE
      ALLOCATE( SWHBCKG(NSEA,1) )
      ALLOCATE( SWHANL(NSEA,1) )
      ! Use SMCWND to determine if reading a seapoint aray for wind
      IF( SMCWND ) THEN
        ALLOCATE( WSBCKG(NSEA,1) )
        ALLOCATE( WDRBCKG(NSEA,1) )
      ELSE
        ALLOCATE(WSBCKG(SIZE(MAPSTA,2), SIZE(MAPSTA,1)))
        ALLOCATE(WDRBCKG(SIZE(MAPSTA,2), SIZE(MAPSTA,1)))
      ENDIF
#endif
    ENDIF
#ifdef W3_T
    IF (.NOT. SMCGRD) THEN
      ALLOCATE( SWHUPRSTR(SIZE(MAPSTA,1), SIZE(MAPSTA,2)) )
    ELSE
      ALLOCATE( SWHUPRSTR(NSEA,1) )
    ENDIF
#endif
    !
    !        Read additional Input: Analysis Field
    INQUIRE(FILE=FLNMANL, EXIST=anl_exists)
    IF (anl_exists) THEN
#ifdef W3_T
      WRITE (NDSO,*) 'shape(SWHANL)', shape(SWHANL)
#endif
#ifdef W3_WRST
      ! For WRST switch read only corrected SWH
      ! Wind will have been read from the restart
      IF (WRSTON) THEN
        CALL READ_GRBTXT(SWHANL, FLNMANL, SMCGRD)
      ELSE
#endif
        CALL READ_GRBTXTWS(SWHANL,WSBCKG,WDRBCKG,FLNMANL,SMCGRD)
#ifdef W3_WRST
      ENDIF
#endif
#ifdef W3_T
      CALL writeMatrix('SWHANL_IN.txt',SWHANL)
#endif
    ELSE
      WRITE (NDSO,*) trim(FLNMANL), ' does not exist, stopping...'
      DEALLOCATE( SWHANL,VATMP,SWHBCKG,VAMAPWS,WSBCKG,WDRBCKG )
#ifdef W3_T
      DEALLOCATE( SWHUPRSTR )
#endif
      STOP
    END IF
    !
#ifdef W3_WRST
    !Calculate wind speed and direction values from u,v..
    !..using cartesian direction convention
    !At present assume only needed for data read from restart
    CALL UVTOCART(WXNwrst,WYNwrst,WSBCKG,WDRBCKG,SMCWND)
#endif
    !
    !        Calculation
    DO ISEA=1, NSEA, 1
      IF (.NOT. SMCGRD) THEN
        IX = MAPSF(ISEA,1)
        IY = MAPSF(ISEA,2)
        IXW = IX
        IYW = IY
#ifdef W3_SMC
      ELSE
        IX = 1
        IY = ISEA
        IF( SMCWND ) THEN
          ! Wind arrays allocated using (X,Y) convention for regular grids
          ! but overriding here for the SMC grid which are always defined
          ! as (NSEA,1) by switching the IY and IX dimension values around
          IXW = IY
          IYW = IX
        ELSE
          IXW = MAPSF(ISEA,1)
          IYW = MAPSF(ISEA,2)
        ENDIF
#endif
      ENDIF
      VATMP = VA(:,ISEA)
      CALL SWH_RSRT_1pw (VATMP, WSBCKG(IXW,IYW), WDRBCKG(IXW,IYW), ISEA, &
           SWHBCKG_1, SWHBCKG_W, SWHBCKG_S, VAMAPWS)
      SWHBCKG(IY,IX)=SWHBCKG_1
      !/
      IF ( SWHBCKG(IY,IX) > 0.01 .AND. SWHANL(IY,IX) > 0.01 ) THEN
        ! If wind-sea is dominant energy component apply correction to
        ! wind-sea part only
        IF ( (SWHBCKG_W / SWHBCKG_1)**2.0 > THRWSEA ) THEN
          ! Apply spectrum updates to wind-sea bins only
          PRCNTG=SQRT((SWHANL(IY,IX)**2.0-SWHBCKG_S**2.0)/SWHBCKG_W**2.0)
          CALL CHECK_PRCNTG(PRCNTG,PRCNTG_CAP)
          CALL UPDTWSPECF(VATMP, PRCNTG, VAMAPWS, ISEA, .FALSE.)
          ! else correct the whole spectrum
        ELSE
          PRCNTG=(SWHANL(IY,IX)/SWHBCKG_1)
          CALL CHECK_PRCNTG(PRCNTG,PRCNTG_CAP)
          IF (CORWSEA) THEN
            ! Include frequency shifts in wind-sea update
            CALL UPDTWSPECF(VATMP, PRCNTG, VAMAPWS, ISEA, .TRUE.)
          ELSE
            ! bulk correction only, as per UPD2
            CALL UPDATE_VA(PRCNTG,VATMP)
          END IF
        END IF
#ifdef W3_T
        WRITE (NDSO,*) 'ISEA = ', ISEA,' IX = ',IX,' IY = ', IY,         &
             ' PRCNTG = ',PRCNTG,' SWHBCKG = ',SWHBCKG(IY,IX), &
             ' SWHANL = ', SWHANL(IY,IX)
#endif
        VA(:,ISEA)=VATMP
#ifdef W3_T
        CALL SWH_RSRT_1p (VATMP, ISEA, SWHTMP)
        SWHUPRSTR(IY,IX)=SWHTMP
        WRITE (NDSO,*) ' =========== UPD6 Output ==========='
        WRITE (NDSO,*)'ISEA = ',ISEA,'SWH_BCKG = ', SWHBCKG(IY,IX), &
             'SWH_ANL = ', SWHANL(IY,IX),                 &
             'SWH_RSTR = ',SWHUPRSTR(IY,IX)
#endif
      END IF
    END DO
#ifdef W3_T
    CALL writeMatrix('SWHBCKG_UPD6.txt', REAL(SWHBCKG  ))
    CALL writeMatrix('SWHANL_UPD6.txt' , REAL(SWHANL   ))
    CALL writeMatrix('SWHRSTR_UPD6.txt', REAL(SWHUPRSTR))
#endif
    !
    DEALLOCATE( SWHANL,VATMP,SWHBCKG,VAMAPWS,WSBCKG,WDRBCKG )
#ifdef W3_T
    DEALLOCATE( SWHUPRSTR )
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !     End of update options
    !/
  END SELECT
  !/
  !/ ------------------------------------------------------------------- /
  ! 6. Write updated restart file
  !/
#ifdef W3_WRST
  ! Copy read wind values from restart for write out
  WXN = WXNwrst
  WYN = WYNwrst
#endif
  WRITE (NDSO,903)
  RSTYPE = 3
  CALL W3IORS ( 'HOT', NDS(6), SIG(NK), 1 )
#ifdef W3_T
  WRITE (NDST,*), MYNAME,' : Exporting VA at the end of the re-analysis'
  CALL writeMatrix('VA02.txt', REAL(VA))
#endif
  !
  !/
  !/ ------------------------------------------------------------------- /
  ! Escape locations read errors 08k:
  !/
  GOTO 888
  !
800 CONTINUE
  WRITE (NDSE,1000) IERR
  CALL EXTCDE ( 10 )
  !
801 CONTINUE
  WRITE (NDSE,1001)
  CALL EXTCDE ( 11 )
  !
802 CONTINUE
  WRITE (NDSE,1002) IERR
  CALL EXTCDE ( 12 )
  !
888 CONTINUE
  WRITE (NDSO,999)
  !/
  !/ ------------------------------------------------------------------- /
  ! Escape locations read errors 2k:
  !/
  GOTO 2222
  !
2001 CONTINUE
  IF ( IAPROC .EQ. NAPERR ) WRITE (NDSE,1001)
  GOTO 2222
  !
2002 CONTINUE
  IF ( IAPROC .EQ. NAPERR ) WRITE (NDSE,1002) IERR
  GOTO 2222
  !
2222 CONTINUE
  !/
  !/ ------------------------------------------------------------------- /
  !  Formats
  !/
900 FORMAT (/15X,'   *** WAVEWATCH III ww3_uprstr Initializing ***   '/ &
       15X,'  ==============================================='/)
901 FORMAT ( '  Comment character is ''',A,''''/)
  !
902 FORMAT ( '  The Option ''',A,''' is used.'/)
  !
903 FORMAT ( '  Exporting the Updated Restart file to "restart001.ww3"'/)
  !
920 FORMAT ( '  Grid name : ',A/)
  !
930 FORMAT (/'  Time interval : '/                                  &
       ' --------------------------------------------------')
  !
931 FORMAT ( '       Starting time : ',A)
  !
932 FORMAT ( '       Ending time   : ',A/)
  !
999 FORMAT (/'  End of program '/                                   &
       ' ========================================='/          &
       '         WAVEWATCH III ww3_uprstr         '/)
  !
1000 FORMAT (/' *** WAVEWATCH III ERROR IN W3UPRSTR : '/              &
       '     ERROR IN OPENING INPUT FILE'/                    &
       '     IOSTAT =',I5/)
  !
1001 FORMAT (/' *** WAVEWATCH III ERROR IN W3UPRSTR : '/               &
       '     PREMATURE END OF INPUT FILE'/)

1002 FORMAT (/' *** WAVEWATCH III ERROR IN W3UPRSTR : '/              &
       '     ERROR IN READING FROM INPUT FILE'/               &
       '     IOSTAT =',I5/)
1004 FORMAT (/' '/,A/)
1005 FORMAT (' ',A, F6.3/)
1006 FORMAT (' ',A, A/)
  !
  !/
CONTAINS
  !/
  !/ ------------------------------------------------------------------- /
  !>
  !> @brief Apply correction to the spectrum.
  !>
  !> @details The factor is (swh_anal/swh_bkg)**2 as applying to wave energy.
  !>
  !> @param[in] PRCNTG
  !> @param[inout] VATMP
  !> @author Stelios Flampouris  @date 16-Oct-2018
  !>
  SUBROUTINE UPDATE_VA (PRCNTG, VATMP)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |         Stelios Flampouris        |
    !/                  |                        FORTRAN 90 |
    !/                  | Created       :       16-Oct-2018 |
    !/                  +-----------------------------------+
    !/
    !/    16-Oct-2018 : Original Code                       ( version 6.06 )
    !/
    !/    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 :
    !     Apply correction to the spectrum
    !
    !  2. Method :
    !     The factor is (swh_anal/swh_bkg)**2 as applying to wave energy
    !  3. Parameters :
    !
    !     Local parameters.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     ----------------------------------------------------------------
    !  Internal Subroutines:
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    !     !/T
    !
    ! 10. Source code :
    !
    !/
    REAL, INTENT(IN) :: PRCNTG
    REAL, DIMENSION(:), INTENT(INOUT)  :: VATMP
    !
    VATMP = (PRCNTG**2)*VATMP
    !
  END SUBROUTINE UPDATE_VA
  !/
  !/ ---------------------------------------------------------------------
  !/
  !>
  !> @brief Last sanity check before the update of the spectrum.
  !>
  !> @details The percentage of change is compared against a user defined cap.
  !>
  !> @param[inout] PRCNTG
  !> @param[inout] PRCNTG_CAP
  !> @author Stelios Flampouris  @date 16-Oct-2018
  !>
  SUBROUTINE CHECK_PRCNTG (PRCNTG,PRCNTG_CAP)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |         Stelios Flampouris        |
    !/                  |                        FORTRAN 90 |
    !/                  | Created       :       16-Oct-2018 |
    !/                  +-----------------------------------+
    !/
    !/    16-Oct-2018 : Original Code                       ( version 6.06 )
    !/    24-Oct-2018 : Update by Andy Saulter              ( version 7.14 )
    !/
    !/    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 :
    !     Last sanity check before the update of the spectrum
    !  2. Method :
    !     The percentage of change is compared against a user defined cap.
    !  3. Parameters :
    !
    !     Local parameters.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     ----------------------------------------------------------------
    !  Internal Subroutines:
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    !     !/T
    !
    ! 10. Source code :
    !
    !/
    REAL, INTENT(INOUT)  ::    PRCNTG
    REAL, INTENT(IN   )  ::    PRCNTG_CAP
    ! local
    CHARACTER(12), PARAMETER :: MYNAME='CHECK_PRCNTG'
#ifdef W3_T

    WRITE (NDSO,*) trim(MYNAME)," The original correction is ",PRCNTG
    WRITE (NDSO,*) trim(MYNAME)," The cap is ",PRCNTG_CAP
#endif
    IF ( PRCNTG_CAP < 1. ) THEN
      WRITE (NDSO,*) trim(MYNAME)," WARNING: PRCNTG_CAP set < 1."
      WRITE (NDSO,*) trim(MYNAME),"          This may introduce spurious corrections"
    END IF
#ifdef W3_T
    WRITE (NDSO,*) trim(MYNAME)," The cap is ",PRCNTG_CAP
#endif
    IF ( PRCNTG > 1. ) THEN
#ifdef W3_T
      WRITE (NDSO,*) trim(MYNAME)," PRCNTG > 1."
#endif
      PRCNTG = MIN(PRCNTG, 1. * PRCNTG_CAP)
    ELSE IF ( PRCNTG < 1. ) THEN
#ifdef W3_T
      WRITE (NDSO,*) trim(MYNAME)," PRCNTG < 1."
#endif
      PRCNTG = MAX(PRCNTG, 1. / PRCNTG_CAP)
#ifdef W3_T

#endif
    END IF
#ifdef W3_T
    WRITE (NDSO,*) trim(MYNAME)," The updated correction is ",PRCNTG
#endif
    !
  END SUBROUTINE CHECK_PRCNTG
  !/
  !/ ------------------------------------------------------------------- /
  !/
  !> @brief Read gribtxt files.
  !>
  !> @param[inout] UPDPRCNT
  !> @param[inout] FLNMCOR
  !> @param[inout] SMCGRD
  !> @author Stelios Flampouris @date 16-Oct-2018
  !>
  SUBROUTINE READ_GRBTXT(UPDPRCNT,FLNMCOR,SMCGRD)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |         Stelios Flampouris        |
    !/                  |                        FORTRAN 90 |
    !/                  | Created       :       15-Mar-2017 |
    !/                  | Last Update   :       16-Oct-2018 |
    !/                  +-----------------------------------+
    !/
    !/    15-Mar-2017 : Original Code                       ( version 6.04 )
    !/    16-Oct-2018 : Generalization of the reader        ( version 6.06 )
    !/
    !/    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 :
    !     Read gribtxt files
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Local parameters.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     ----------------------------------------------------------------
    !  Internal Subroutines:
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    !     !/T
    !
    ! 10. Source code :
    !
    !/
    REAL, DIMENSION(:,:), INTENT(OUT) :: UPDPRCNT
    CHARACTER(*), INTENT(IN) :: FLNMCOR
    LOGICAL, INTENT(IN) :: SMCGRD
    ! Local Variables
    INTEGER            :: I, J, IERR
    INTEGER            :: K, L, M, N
    REAL :: A
    INTEGER, PARAMETER :: IP_FID = 123
    CHARACTER(25), PARAMETER::myname='read_grbtxt'
    !
#ifdef W3_T
    WRITE (NDSO,*) trim(myname), ' starts'
#endif
    J     = LEN_TRIM(FNMPRE)
    OPEN (IP_FID,FILE=FNMPRE(:J)//TRIM(FLNMCOR),STATUS='OLD' &
         ,ACTION='read',IOSTAT=IERR)
    !
    ! Read text header and check dimensions match expected values
    IF (.NOT. SMCGRD) THEN
      READ( IP_FID, *) M,N
      IF (( SIZE(UPDPRCNT,1) /= N) .OR. ( SIZE(UPDPRCNT,2) /= M )) THEN
        WRITE (NDSO,*) trim(myname),': These are not the grid ' // &
             'dimensions: M=',M,' N=',N
        STOP
      END IF
#ifdef W3_SMC
    ELSE
      READ( IP_FID, *) N
      IF ( SIZE(UPDPRCNT,1) /= N ) THEN
        WRITE (NDSO,*) trim(myname),': These are not the grid ' // &
             'dimensions: N=',N
        STOP
      END IF
#endif
    END IF
    UPDPRCNT=0
    !
    ! Read the data into its allocated array
    IF (.NOT. SMCGRD) THEN
      DO L=1,N
        DO K=1,M
          A=0.
          READ(IP_FID,*)A
          UPDPRCNT(N+1-L,K)=A
        END DO
      END DO
#ifdef W3_SMC
    ELSE
      DO L=1,N
        A=0.
        READ(IP_FID,*)A
        UPDPRCNT(L,1)=A
      END DO
#endif
    END IF
    !
    CLOSE(IP_FID)
    !
#ifdef W3_T
    WRITE (NDSO,*) trim(myname), ' ends'
#endif
  END SUBROUTINE READ_GRBTXT
  !/
  !/ ------------------------------------------------------------------- /
  !/
  !>
  !> @brief Read txt files that include wind data.
  !>
  !> @param[out] UPDPRCNT
  !> @param[out] WSPD
  !> @param[out] WDIR
  !> @param[in] FLNMCOR
  !> @param[in] SMCGRD
  !>
  !> @author Andy Saulter @date 24-Oct-2018
  !>
  SUBROUTINE READ_GRBTXTWS(UPDPRCNT,WSPD,WDIR,FLNMCOR,SMCGRD)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |         Andy Saulter              |
    !/                  |                        FORTRAN 90 |
    !/                  | Original code :       24-Oct-2018 |
    !/                  | Last update :         05-Oct-2019 |
    !/                  +-----------------------------------+
    !/
    !/    24-Oct-2018 : Original Code                       ( version 6.07 )
    !/
    !/    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 :
    !     Read txt files that include wind data
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Local parameters.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     ----------------------------------------------------------------
    !  Internal Subroutines:
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    !     !/T
    !
    ! 10. Source code :
    !
    !/
    REAL, DIMENSION(:,:), INTENT(OUT) :: UPDPRCNT, WSPD, WDIR
    CHARACTER(*), INTENT(IN) :: FLNMCOR
    LOGICAL, INTENT(IN) :: SMCGRD
    ! Local Variables
    INTEGER            :: I, J, IERR
    INTEGER            :: K, L, M, N
    REAL :: A, WS, WD
    INTEGER, PARAMETER :: IP_FID = 123
    CHARACTER(25), PARAMETER::myname='read_grbtxt'
    !
#ifdef W3_T
    WRITE (NDSO,*) trim(myname), ' starts'
#endif
    J     = LEN_TRIM(FNMPRE)
    OPEN (IP_FID,FILE=FNMPRE(:J)//TRIM(FLNMCOR),STATUS='OLD' &
         ,ACTION='read',IOSTAT=IERR)
    !
    ! Read text header and check dimensions match expected values
    IF (.NOT. SMCGRD) THEN
      READ( IP_FID, *) M,N
      IF (( SIZE(UPDPRCNT,1) /= N) .OR. ( SIZE(UPDPRCNT,2) /= M )) THEN
        WRITE (NDSO,*) trim(myname),': These are not the grid ' // &
             'dimensions: M=',M,' N=',N
        STOP
      END IF
#ifdef W3_SMC
    ELSE
      READ( IP_FID, *) N
      IF ( SIZE(UPDPRCNT,1) /= N ) THEN
        WRITE (NDSO,*) trim(myname),': These are not the grid ' // &
             'dimensions: N=',N
        STOP
      END IF
#endif
    END IF
    UPDPRCNT=0
    WSPD=0.
    WDIR=0.
    !
    ! Read the data into allocated arrays
    IF (.NOT. SMCGRD) THEN
      DO L=1,N
        DO K=1,M
          A=0.
          WS=0.
          WD=0.
          READ(IP_FID,*)A, WS, WD
          !SWH data read onto Y,X grid
          UPDPRCNT(N+1-L,K)=A
          !Wind data read onto X,Y grid
          WSPD(K,N+1-L)=WS
          WDIR(K,N+1-L)=WD
        END DO
      END DO
#ifdef W3_SMC
    ELSE
      DO L=1,N
        A=0.
        READ(IP_FID,*)A, WS, WD
        UPDPRCNT(L,1)=A
        WSPD(L,1)=WS
        WDIR(L,1)=WD
      END DO
#endif
    ENDIF
    !
    CLOSE(IP_FID)
    !
    !
#ifdef W3_T
    WRITE (NDSO,*) trim(myname), ' ends'
#endif
  END SUBROUTINE READ_GRBTXTWS
  !/
  !/ ------------------------------------------------------------------- /
  !/
  !> @brief Calculate the significant wave height from the restart file for 1 point.
  !>
  !> @param[in] VA1p
  !> @param[in] ISEA1p
  !> @param[out] HSIG1p
  !> @author Stelios Flampouris  @date 15-May-2017
  !>
  SUBROUTINE SWH_RSRT_1p (VA1p, ISEA1p, HSIG1p )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |         Stelios Flampouris        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         15-Mar-2017 |
    !/                  +-----------------------------------+
    !/
    !/    15-Mar-2017 : Original Code                       ( version 6.04 )
    !/
    !/    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 :
    !     Calculate the significant wave height from the restart file for 1 point
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Local parameters.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     ----------------------------------------------------------------
    !  Internal Subroutines:
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    !     !/T
    !
    ! 10. Source code :
    !
    !/
    REAL, INTENT(OUT) :: HSIG1p
    INTEGER, INTENT(IN) :: ISEA1p
    REAL, DIMENSION(:), INTENT(IN)  :: VA1p
    CHARACTER(25),PARAMETER :: myname='SWH_RSRT_1p'
    !
#ifdef W3_FT
    WRITE (NDSO,*)' '
    WRITE (NDSO,*) trim(myname), ' starts'
#endif
    HSIG1p = 0.
    DEPTH  = MAX ( DMIN , -ZB(ISEA1p) )
    ETOT   = 0.
    !
    DO IK=1, NK
      CALL WAVNU1 ( SIG(IK), DEPTH, WN, CG )
      E1I    = 0.
      DO ITH=1, NTH
        E1I    = E1I + VA1p(ITH+(IK-1)*NTH)  *  SIG(IK) / CG
      END DO
      ETOT   = ETOT + E1I*DSIP(IK)
    END DO
    !
    HSIG1p = 4. * SQRT ( ETOT * DTH )
    !
#ifdef W3_FT
    WRITE (NDSO,*) ' ', trim(myname), ' ends'
    WRITE (NDSO,*)' '
#endif
  END SUBROUTINE SWH_RSRT_1p
  !/
  !/ ------------------------------------------------------------------- /
  !/
  !> @brief Calculate Hs from restart for 1 point.
  !>
  !> @details Calculate the significant wave height for total, wind sea and
  !>  swell components from the restart file for 1 point
  !>
  !> @param[in] VA1p
  !> @param[in] WS
  !> @param[in] WD
  !> @param[in] ISEA1p
  !> @param[out] HSIG1p
  !> @param[out] HSIGwp
  !> @param[out] HSIGsp
  !> @param[out] VAMAPWS
  !>
  !> @author Andy Saulter  @date 05-0ct-2019
  !>
  SUBROUTINE SWH_RSRT_1pw (VA1p, WS, WD, ISEA1p, HSIG1p, HSIGwp, HSIGsp, VAMAPWS )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |         Andy Saulter              |
    !/                  |                        FORTRAN 90 |
    !/                  | Original code :       24-Oct-2018 |
    !/                  | Last update :         05-Oct-2019 |
    !/                  +-----------------------------------+
    !/
    !/    24-Oct-2018 : Original Code                       ( version 6.07 )
    !/
    !/    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 :
    !     Calculate the significant wave height for total, wind sea and
    !     swell components from the restart file for 1 point
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Local parameters.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     ----------------------------------------------------------------
    !  Internal Subroutines:
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    !     !/T
    !
    ! 10. Source code :
    !
    !/
    USE W3GDATMD, ONLY: TH
    USE W3ODATMD, ONLY: WSMULT !same wind sea cut-off factor for sea/swell outputs
    !
    REAL, INTENT(OUT) :: HSIG1p, HSIGwp, HSIGsp
    INTEGER, INTENT(IN) :: ISEA1p
    REAL, INTENT(IN) :: WS, WD
    REAL, DIMENSION(:), INTENT(IN)  :: VA1p
    INTEGER, DIMENSION(:), INTENT(OUT) :: VAMAPWS ! Wind-sea id for spectral bins
    REAL :: RELWS, ETOTw, ETOTs, EwI, EsI
    CHARACTER(25),PARAMETER :: myname='SWH_RSRT_1pw'
    !
#ifdef W3_T
    WRITE (NDSO,*) trim(myname), ' starts'
#endif
    HSIG1p = 0.
    HSIGwp = 0.
    HSIGsp = 0.
    DEPTH  = MAX ( DMIN , -ZB(ISEA1p) )
    ETOT   = 0.
    ETOTw  = 0.
    ETOTs  = 0.
    !
    DO IK=1, NK
      CALL WAVNU1 ( SIG(IK), DEPTH, WN, CG )
      E1I    = 0.
      EwI    = 0.
      EsI    = 0.
      DO ITH=1, NTH
        ! Relative wind-sea calc assumes input with in direction toward
        ! i.e. same as for the wave spectrum
        RELWS = WSMULT * WS * MAX(0.0, COS(WD - TH(ITH)))
        E1I    = E1I + VA1p(ITH+(IK-1)*NTH) * SIG(IK) / CG
        IF ( RELWS > (SIG(IK)/WN) ) THEN
          EwI = EwI + VA1p(ITH+(IK-1)*NTH) * SIG(IK) / CG
          VAMAPWS(ITH+(IK-1)*NTH) = 1
        ELSE
          EsI = EsI + VA1p(ITH+(IK-1)*NTH) * SIG(IK) / CG
          VAMAPWS(ITH+(IK-1)*NTH) = 0
        END IF
      END DO
      ETOT   = ETOT + E1I*DSIP(IK)
      ETOTw  = ETOTw + EwI*DSIP(IK)
      ETOTs  = ETOTs + EsI*DSIP(IK)
    END DO
    !
    HSIG1p = 4. * SQRT ( ETOT * DTH )
    HSIGwp = 4. * SQRT ( ETOTw * DTH )
    HSIGsp = 4. * SQRT ( ETOTs * DTH )
    !
#ifdef W3_T
    WRITE (NDSO,*) trim(myname), ' ends'
#endif
  END SUBROUTINE SWH_RSRT_1pw
  !/
  !/ ------------------------------------------------------------------- /
  !/
  !> @brief Calculate speed and cartesian convention directions from u,v
  !>  input vectors.
  !>
  !> @param[in] UVEC
  !> @param[in] VVEC
  !> @param[out] SPD
  !> @param[out] DCART
  !> @param[in] SMCGRD
  !>
  !> @author Andy Saulter  @date 05-Oct-2019
  !>
  SUBROUTINE UVTOCART (UVEC, VVEC, SPD, DCART, SMCGRD)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |         Andy Saulter              |
    !/                  |                        FORTRAN 90 |
    !/                  | Original code :       05-Oct-2019 |
    !/                  +-----------------------------------+
    !/
    !/    05-Oct-2019 : Original Code                       ( version 6.07 )
    !/
    !/    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 :
    !     Calculate speed and cartesian convention directions from u,v
    !     input vectors
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Local parameters.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     ----------------------------------------------------------------
    !  Internal Subroutines:
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    !     !/T
    !
    ! 10. Source code :
    !
    !/
    USE CONSTANTS, ONLY: TPI
    !
    REAL, DIMENSION(:,:), INTENT(OUT) :: SPD, DCART
    REAL, DIMENSION(:,:), INTENT(IN) :: UVEC, VVEC
    LOGICAL, INTENT(IN) :: SMCGRD
    !
#ifdef W3_T
    WRITE (NDSO,*) trim(myname), ' starts'
#endif
    !
    DO ISEA=1, NSEA, 1
      IF (.NOT. SMCGRD) THEN
        IX = MAPSF(ISEA,1)
        IY = MAPSF(ISEA,2)
#ifdef W3_SMC
      ELSE
        IX = 1
        IY = ISEA
#endif
      ENDIF
      !
      SPD(IY,IX) = SQRT( UVEC(IY,IX)**2 + VVEC(IY,IX)**2 )
      IF( SPD(IY,IX) .GT. 1.E-7) THEN
        DCART = MOD( TPI+ATAN2(UVEC(IY,IX),VVEC(IY,IX)) , TPI )
      ELSE
        DCART = 0
      END IF
      SPD(IY,IX) = MAX( SPD(IY,IX) , 0.001 )
    END DO
    !
#ifdef W3_T
    WRITE (NDSO,*) trim(myname), ' ends'
#endif
  END SUBROUTINE UVTOCART
  !/
  !/ ------------------------------------------------------------------- /
  !/
  !> @brief Updates the wind-sea part of the wave spectrum only.
  !>
  !> @param[inout] VATMP
  !> @param[in] PRCNTG
  !> @param[in] VAMAPWS
  !>
  !> @author Andy Saulter  @date 05-Oct-2019
  !>
  !>
  SUBROUTINE UPDTWSPEC(VATMP, PRCNTG, VAMAPWS)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |         Andy Saulter              |
    !/                  |                        FORTRAN 90 |
    !/                  | Original code :       24-Oct-2018 |
    !/                  | Last update :         05-Oct-2019 |
    !/                  +-----------------------------------+
    !/
    !/    24-Oct-2018 : Original Code                       ( version 6.07 )
    !/
    !/    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 :
    !     Updates the wind-sea part of the wave spectrum only
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Local parameters.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     ----------------------------------------------------------------
    !  Internal Subroutines:
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    !     !/T
    !
    ! 10. Source code :
    !
    !/
    REAL, DIMENSION(:), INTENT(INOUT)  :: VATMP
    INTEGER, DIMENSION(:), INTENT(IN) :: VAMAPWS
    REAL, INTENT(IN) :: PRCNTG
    CHARACTER(25),PARAMETER :: myname='UPDTWSPEC'
    !
#ifdef W3_T
    WRITE (NDSO,*) trim(myname), ' starts'
#endif
    DO IK=1, NK
      DO ITH=1, NTH
        IF ( VAMAPWS(ITH+(IK-1)*NTH) .EQ. 1 ) THEN
          VATMP(ITH+(IK-1)*NTH) = VATMP(ITH+(IK-1)*NTH) * PRCNTG**2
        END IF
      END DO
    END DO
    !
#ifdef W3_T
    WRITE (NDSO,*) trim(myname), ' ends'
#endif
  END SUBROUTINE UPDTWSPEC
  !/
  !/ ------------------------------------------------------------------- /
  !/
  !> @brief Updates the wind-sea part of the wave spectrum and shifts
  !>  in frequency space.
  !>
  !> @param[inout] VATMP
  !> @param[in] PRCNTG
  !> @param[in] VAMAPWS
  !> @param[in] ISEA1p
  !> @param[in] ADJALL
  !>
  !> @author Andy Saulter @date 05-Oct-2019
  !>
  SUBROUTINE UPDTWSPECF(VATMP, PRCNTG, VAMAPWS, ISEA1p, ADJALL)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |         Andy Saulter              |
    !/                  |                        FORTRAN 90 |
    !/                  | Original code :       24-Oct-2018 |
    !/                  | Last update :         05-Oct-2019 |
    !/                  +-----------------------------------+
    !/
    !/    24-Oct-2018 : Original Code                       ( version 6.07 )
    !/
    !/    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 :
    !     Updates the wind-sea part of the wave spectrum and shifts in frequency
    !     space
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Local parameters.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     ----------------------------------------------------------------
    !  Internal Subroutines:
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    !     !/T
    !
    ! 10. Source code :
    !
    !/
    REAL, DIMENSION(:), INTENT(INOUT)  :: VATMP
    INTEGER, DIMENSION(:), INTENT(IN) :: VAMAPWS
    REAL, INTENT(IN) :: PRCNTG
    INTEGER, INTENT(IN) :: ISEA1p
    LOGICAL, INTENT(IN) :: ADJALL
    CHARACTER(25),PARAMETER :: myname='UPDTWSPECF'
    REAL :: FFAC, SIGSHFT, FDM1, FDM2, WN1, CG1, WN2, CG2
    INTEGER :: LPF, M1, M2
    REAL, ALLOCATABLE :: VASHFT(:)
    !
#ifdef W3_T
    WRITE (NDSO,*) trim(myname), ' starts'
#endif
    DEPTH  = MAX( DMIN , -ZB(ISEA1p))
    ALLOCATE(VASHFT(SIZE(VATMP)))
    VASHFT(:) = 0.0
    !
    ! 1st iteration shifts wind-sea energy in freq space
    FFAC = (1. / PRCNTG**2)**(1.0/3.0) ! uses Toba's relationship
    DO IK=1, NK
      CALL WAVNU1(SIG(IK), DEPTH, WN, CG)
      SIGSHFT = FFAC * SIG(IK)
      DO ITH=1, NTH
        IF ( VAMAPWS(ITH+(IK-1)*NTH) .EQ. 1 ) THEN
          ! Interpolate frequency bin according to f-shift
          LPF = 1
          DO WHILE (LPF < NK)
            IF (SIG(LPF) >= SIGSHFT) THEN
              IF (LPF .EQ. 1) THEN
                CALL WAVNU1(SIG(LPF), DEPTH, WN1, CG1)
                VASHFT(ITH+(LPF-1)*NTH) = VASHFT(ITH+(LPF-1)*NTH) + &
                     VATMP(ITH+(IK-1)*NTH) *  &
                     (DSIP(IK)*SIG(IK)/CG) /  &
                     (DSIP(LPF)*SIG(LPF)/CG1)
              ELSE
                M2 = LPF
                M1 = LPF - 1
                FDM1 = SIGSHFT - SIG(M1)
                FDM2 = SIG(M2) - SIG(M1)
                CALL WAVNU1(SIG(M1), DEPTH, WN1, CG1)
                CALL WAVNU1(SIG(M2), DEPTH, WN2, CG2)
                VASHFT(ITH+(M1-1)*NTH) = VASHFT(ITH+(M1-1)*NTH) + &
                     (FDM1 / FDM2) *         &
                     VATMP(ITH+(IK-1)*NTH) * &
                     (DSIP(IK)*SIG(IK)/CG) / &
                     (DSIP(M1)*SIG(M1)/CG1)
                VASHFT(ITH+(M2-1)*NTH) = VASHFT(ITH+(M2-1)*NTH) + &
                     (1.0 - FDM1 / FDM2) *   &
                     VATMP(ITH+(IK-1)*NTH) * &
                     (DSIP(IK)*SIG(IK)/CG) / &
                     (DSIP(M2)*SIG(M2)/CG2)
              END IF
              LPF = NK + 1
            ENDIF
            LPF = LPF + 1
          END DO
          IF (LPF .EQ. NK) THEN
            CALL WAVNU1(SIG(LPF), DEPTH, WN1, CG1)
            VASHFT(ITH+(LPF-1)*NTH) = VASHFT(ITH+(LPF-1)*NTH) + &
                 VATMP(ITH+(IK-1)*NTH) *  &
                 (DSIP(IK)*SIG(IK)/CG) /  &
                 (DSIP(LPF)*SIG(LPF)/CG1)
          END IF
        END IF
      END DO
    END DO
    ! 2nd iteration scales wind-sea energy
    DO IK=1, NK
      DO ITH=1, NTH
        IF ( VAMAPWS(ITH+(IK-1)*NTH) .EQ. 1 ) THEN
          VASHFT(ITH+(IK-1)*NTH) = VASHFT(ITH+(IK-1)*NTH) * PRCNTG**2
        END IF
      END DO
    END DO
    ! 3rd iteration combines wind-sea and swell energy
    DO IK=1, NK
      DO ITH=1, NTH
        IF ( VAMAPWS(ITH+(IK-1)*NTH) .EQ. 1 ) THEN
          VATMP(ITH+(IK-1)*NTH) = VASHFT(ITH+(IK-1)*NTH)
        ELSE
          IF ( ADJALL ) THEN
            ! Swell components are also re-scaled
            VATMP(ITH+(IK-1)*NTH) = VATMP(ITH+(IK-1)*NTH) * &
                 PRCNTG**2            + &
                 VASHFT(ITH+(IK-1)*NTH)
          ELSE
            ! Re-scaling wind-sea only
            VATMP(ITH+(IK-1)*NTH) = VATMP(ITH+(IK-1)*NTH) + &
                 VASHFT(ITH+(IK-1)*NTH)
          END IF
        END IF
      END DO
    END DO
    !
    DEALLOCATE(VASHFT)
    !
#ifdef W3_T
    WRITE (NDSO,*) trim(myname), ' ends'
#endif
  END SUBROUTINE UPDTWSPECF
  !/
  !/ ------------------------------------------------------------------- /
  !/
  !> @brief Writes a 2D array to text file, column by column.
  !>
  !> @param[inout] FILENAME Path to the output file.
  !> @param[inout] RDA_A 2D array to write.
  !>
  !> @author Stelios Flampouris  @date 15-Mar-2017
  !>
  SUBROUTINE WRITEMATRIX(FILENAME, RDA_A)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |         Stelios Flampouris        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         15-Mar-2017 |
    !/                  +-----------------------------------+
    !/
    !/    15-Mar-2017 : Original Code                       ( version 6.04 )
    !/
    !/    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 :
    !     Writes a 2D array to text file, column by column
    !  2. Method :
    !
    !  3. Parameters :
    !     fileName path to the output file
    !     rda_A 2D array to write
    !
    !     Local parameters.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     ----------------------------------------------------------------
    !  Internal Subroutines:
    !
    !  5. Called by :
    !  Any routine that has to write 2d arrays !?!
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    !     !/T
    !
    ! 10. Source code :
    !
    !/
    REAL, DIMENSION(:, :), INTENT(IN) :: RDA_A
    CHARACTER(*)         , INTENT(IN) :: FILENAME
    INTEGER IB_I, IB_J, IL_IOS
    INTEGER, PARAMETER :: IP_FID = 123
    !
    OPEN( UNIT = IP_FID, FILE = FILENAME, STATUS = 'REPLACE', &
         FORM = 'FORMATTED', IOSTAT = IL_IOS)
    IF (IL_IOS /= 0) PRINT*,'In writeMatrix : Error creating file'//FILENAME
    DO IB_J = 1, SIZE(RDA_A,2)
      DO IB_I = 1, SIZE(RDA_A,1)
        !               write(unit=ip_fid, fmt='(I, $)') rda_A(ib_i,ib_j)
        WRITE(UNIT=IP_FID, FMT='(E18.8, $)') RDA_A(IB_I,IB_J)
      END DO
      WRITE(UNIT=IP_FID, FMT=*)''
    END DO
    CLOSE(IP_FID)
    !
  END SUBROUTINE WRITEMATRIX
  !/
  !/ ------------------------------------------------------------------- /
  !/
END PROGRAM W3UPRSTR