#include "w3macros.h"
!/ ------------------------------------------------------------------- /
   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)
!/
!/    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
      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
!
      IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!  Local variables
!/
      INTEGER                 :: NDSI, NDSM, NDSTRC, NTRACE, IERR, I, J
      CHARACTER               :: COMSTR*1
!
!      REAL, ALLOCATABLE       :: BETAW(:)
!      LOGICAL, ALLOCATABLE    :: MASK(:)
      LOGICAL                 :: anl_exists, CORWSEA
      INTEGER                 :: IMOD, INTYPE, 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 )
!
      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)
!
!/
!/ ------------------------------------------------------------------- /
!  2. Read the ww3_uprstr.inp
!/
      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 )
      READ (NDSI,*,END=2001,ERR=2002) PRCNTG
      CALL NEXTLN ( COMSTR , NDSI , NDSEN )
      READ (NDSI,*,END=2001,ERR=2002) PRCNTG_CAP
      CALL NEXTLN ( COMSTR , NDSI , NDSEN )
      READ (NDSI,*,END=2001,ERR=2002) FLNMANL
!/
!/ ------------------------------------------------------------------- /
!  3.  Read model definition file.
!/
      CALL W3IOGR ( 'READ', NDSM )
      NSEAL  = NSEA
      WRITE (NDSO,920) GNAME
!/
!/ ------------------------------------------------------------------- /
!  4. Read restart file
!/
      CALL W3IORS ( 'READ', NDS(6), SIG(NK), INTYPE, IMOD )!
      IF ( IAPROC .EQ. NAPLOG ) THEN
         IF (INTYPE.EQ.0.OR.INTYPE.EQ.1.OR.INTYPE.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'
         END IF
      END IF
!/
!/ ------------------------------------------------------------------- /
! 5. Update restart spectra array according to the selected option
!/
      SELECT CASE (UPDPROC)
!/
!/ ------------------------------------------------------------------- /
! UPD0F
!/
      CASE ('UPD0F')
         WRITE (NDSO,902) 'UPD0F'
         DO ISEA=1, NSEA, 1
            CALL UPDATE_VA(PRCNTG, VA(:,ISEA))
         END DO
!/
!/ ------------------------------------------------------------------- /
! 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'
!
!        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)) )
         ENDIF
!
!        Read additional Input: Analysis Field
         INQUIRE(FILE=FLNMANL, EXIST=anl_exists)
         IF (anl_exists) THEN
            CALL READ_GRBTXT(SWHANL, FLNMANL, SMCGRD)
         ELSE
            WRITE (NDSO,*) trim(FLNMANL), ' does not exist, stopping...'
            DEALLOCATE( SWHANL,VATMP,SWHBCKG )
            STOP
         END IF
!
!        Calculation
         DO ISEA=1, NSEA, 1
            IF (.NOT. SMCGRD) THEN
               IX = MAPSF(ISEA,1)
               IY = MAPSF(ISEA,2)
            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)
               CALL CHECK_PRCNTG (PRCNTG,PRCNTG_CAP)
               CALL UPDATE_VA(PRCNTG, VA(:,ISEA))
            END IF
         END DO
!
         DEALLOCATE( SWHANL,VATMP,SWHBCKG )
!/
!/ ------------------------------------------------------------------- /
! UPD3
! As per UPD2, but the update factor is a surface with the shape of the
! background spectrum
!/
      CASE ('UPD3')
         WRITE (NDSO,902) 'UPD3'
!
!        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)) )
         ENDIF
!
!        Read additional Input: Analysis Field
         INQUIRE(FILE=FLNMANL, EXIST=anl_exists)
         IF (anl_exists) THEN
            CALL READ_GRBTXT(SWHANL, FLNMANL, SMCGRD)
         ELSE
            WRITE (NDSO,*) trim(FLNMANL), ' does not exist, stopping...'
            DEALLOCATE( SWHANL,VATMP,SWHBCKG,VATMP_NORM,A )
            STOP
         END IF
!
!        Calculation
         DO ISEA=1, NSEA, 1
            IF (.NOT. SMCGRD) THEN
               IX = MAPSF(ISEA,1)
               IY = MAPSF(ISEA,2)
            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)
               CALL CHECK_PRCNTG(PRCNTG,PRCNTG_CAP)
               VATMP_NORM=VATMP/SUM(VATMP)
               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)
               CALL CHECK_PRCNTG (PRCNTG,PRCNTG_CAP)
               CALL UPDATE_VA(PRCNTG, VATMP)
               VA(:,ISEA)=VATMP
            END IF
         END DO
!
         DEALLOCATE( SWHANL,VATMP,SWHBCKG,VATMP_NORM,A )
!/
!/ ------------------------------------------------------------------- /
! 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'
         ! 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)) )
         ENDIF
!
!        Read additional Input: Analysis Field
         INQUIRE(FILE=FLNMANL, EXIST=anl_exists)
         IF (anl_exists) THEN
            CALL READ_GRBTXTWS(SWHANL,WSBCKG,WDRBCKG,FLNMANL,SMCGRD)
         ELSE
            WRITE (NDSO,*) trim(FLNMANL), ' does not exist, stopping...'
            DEALLOCATE( SWHANL,VATMP,SWHBCKG,VAMAPWS,WSBCKG,WDRBCKG )
            STOP
         END IF
!
!        Calculation
         DO ISEA=1, NSEA, 1
            IF (.NOT. SMCGRD) THEN
               IX = MAPSF(ISEA,1)
               IY = MAPSF(ISEA,2)
               IXW = IX
               IYW = IY
            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
               VA(:,ISEA)=VATMP
            END IF
         END DO
!
         DEALLOCATE( SWHANL,VATMP,SWHBCKG,VAMAPWS,WSBCKG,WDRBCKG )
!/
!/ ------------------------------------------------------------------- /
! 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'
         ! 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)) )
         ENDIF
!
!        Read additional Input: Analysis Field
         INQUIRE(FILE=FLNMANL, EXIST=anl_exists)
         IF (anl_exists) THEN
            CALL READ_GRBTXTWS(SWHANL,WSBCKG,WDRBCKG,FLNMANL,SMCGRD)
         ELSE
            WRITE (NDSO,*) trim(FLNMANL), ' does not exist, stopping...'
            DEALLOCATE( SWHANL,VATMP,SWHBCKG,VAMAPWS,WSBCKG,WDRBCKG )
            STOP
         END IF
!
!        Calculation
         DO ISEA=1, NSEA, 1
            IF (.NOT. SMCGRD) THEN
               IX = MAPSF(ISEA,1)
               IY = MAPSF(ISEA,2)
               IXW = IX
               IYW = IY
            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
               VA(:,ISEA)=VATMP
            END IF
         END DO
!
         DEALLOCATE( SWHANL,VATMP,SWHBCKG,VAMAPWS,WSBCKG,WDRBCKG )
!/
!/ ------------------------------------------------------------------- /
!     End of update options
!/
      END SELECT
!/
!/ ------------------------------------------------------------------- /
! 6. Write updated restart file
!/
      WRITE (NDSO,903)
      INTYPE = 3
      CALL W3IORS ( 'HOT', NDS(6), SIG(NK), INTYPE, 1 )
!
!/
!/ ------------------------------------------------------------------- /
! 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_uprsrt         '/)
!
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, I5, I5/)
!
!/
   CONTAINS
!/
!/ ------------------------------------------------------------------- /
!/
   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
!/
!/ ---------------------------------------------------------------------
!/
   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 x.xx )
!/
!/    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'
      IF ( PRCNTG > 1. ) THEN
         PRCNTG = MIN(PRCNTG, 1. * PRCNTG_CAP)
      ELSE IF ( PRCNTG < 1. ) THEN
         PRCNTG = MAX(PRCNTG, 1. / PRCNTG_CAP)
      END IF
!
   END SUBROUTINE CHECK_PRCNTG
!/
!/ ------------------------------------------------------------------- /
!/
      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'
!
         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
         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
         END IF
!
         CLOSE(IP_FID)
!
      END SUBROUTINE READ_GRBTXT
!/
!/ ------------------------------------------------------------------- /
!/
      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'
!
         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
         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
         ENDIF
!
         CLOSE(IP_FID)
!
      END SUBROUTINE READ_GRBTXTWS
!/
!/ ------------------------------------------------------------------- /
!/
      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'
!
         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 )
!
      END SUBROUTINE SWH_RSRT_1p
!/
!/ ------------------------------------------------------------------- /
!/
      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'
!
         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 )
!
      END SUBROUTINE SWH_RSRT_1pw
!/
!/ ------------------------------------------------------------------- /
!/
      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
!
         DO ISEA=1, NSEA, 1
            IF (.NOT. SMCGRD) THEN
               IX = MAPSF(ISEA,1)
               IY = MAPSF(ISEA,2)
            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
!
      END SUBROUTINE UVTOCART
!/
!/ ------------------------------------------------------------------- /
!/
      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'
!
         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
!
      END SUBROUTINE UPDTWSPEC
!/
!/ ------------------------------------------------------------------- /
!/
      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(:)
!
         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)
!
      END SUBROUTINE UPDTWSPECF
!/
!/ ------------------------------------------------------------------- /
!/
      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
