ww3_uprstr.F90

Go to the documentation of this file.

 
 
#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)
  !/    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
  !/
  !/ ------------------------------------------------------------------- /
  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 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
  !/
  !/ ------------------------------------------------------------------- /
  !/
  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
  !/
  !/ ------------------------------------------------------------------- /
  !/
  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
  !/
  !/ ------------------------------------------------------------------- /
  !/
  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
  !/
  !/ ------------------------------------------------------------------- /
  !/
  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
  !/
  !/ ------------------------------------------------------------------- /
  !/
  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
  !/
  !/ ------------------------------------------------------------------- /
  !/
  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
  !/
  !/ ------------------------------------------------------------------- /
  !/
  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
  !/
  !/ ------------------------------------------------------------------- /
  !/
  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