w3wavemd.F90

Go to the documentation of this file.

 
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
MODULE w3wavemd
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         13-Sep-2022 |
  !/                  +-----------------------------------+
  !/
  !/    04-Feb-2000 : Origination.                        ( version 2.00 )
  !/                  For upgrades see subroutines.
  !/    14-Feb-2000 : Exact-NL added.                     ( version 2.01 )
  !/    05-Jan-2001 : Bug fix to allow model to run       ( version 2.05 )
  !/                  without output.
  !/    24-Jan-2001 : Flat grid version.                  ( version 2.06 )
  !/    09-Feb-2001 : Third propagation scheme added.     ( version 2.08 )
  !/    23-Feb-2001 : Check for barrier after source
  !/                  terms added ( W3NMIN ).     ( delayed version 2.07 )
  !/    16-Mar-2001 : Fourth propagation scheme added.    ( version 2.09 )
  !/    30-Mar-2001 : Sub-grid obstacles added.           ( version 2.10 )
  !/    23-May-2001 : Clean up and bug fixes.             ( version 2.11 )
  !/    10-Dec-2001 : Sub-grid obstacles for UQ schemes.  ( version 2.14 )
  !/    11-Jan-2002 : Sub-grid ice.                       ( version 2.15 )
  !/    24-Jan-2002 : Zero time step dor data ass.        ( version 2.17 )
  !/    18-Feb-2002 : Point output diagnostics added.     ( version 2.18 )
  !/    30-Apr-2002 : Add field output types 17-18.       ( version 2.20 )
  !/    09-May-2002 : Switch clean up.                    ( version 2.21 )
  !/    13-Nov-2002 : Add stress vector.                  ( version 3.00 )
  !/    26-Dec-2002 : Moving grid version.                ( version 3.02 )
  !/    01-Aug-2003 : Moving grid GSE correction.         ( version 3.03 )
  !/    20-Aug-2003 : Output server options added.        ( version 3.04 )
  !/    07-Oct-2003 : Output options for NN training.     ( version 3.05 )
  !/    29-Dec-2004 : Multiple grid version.              ( version 3.06 )
  !/                  W3INIT, W3MPII-O and WWVER moved to w3initmd.ftn
  !/    04-Feb-2005 : Add STAMP to par list of W3WAVE.    ( version 3.07 )
  !/    04-May-2005 : Change to MPI_COMM_WAVE.            ( version 3.07 )
  !/    28-Jun-2005 : Adding map recalc for W3ULEV call.  ( version 3.07 )
  !/    07-Sep-2005 : Updated boundary conditions.        ( version 3.08 )
  !/                  Fix NRQSG1/2 = 0 array bound issue.
  !/    13-Jun-2006 : Split STORE in G/SSTORE             ( version 3.09 )
  !/    26-Jun-2006 : Add output type 6.                  ( version 3.09 )
  !/    04-Jul-2006 : Consolidate stress arrays.          ( version 3.09 )
  !/    18-Oct-2006 : Partitioned spectral data output.   ( version 3.10 )
  !/    02-Feb-2007 : Add FLAGST test.                    ( version 3.10 )
  !/    02-Apr-2007 : Add partitioned field data.         ( version 3.11 )
  !/    07-May-2007 : Bug fix SKIP_O treatment.           ( version 3.11 )
  !/    17-May-2007 : Adding NTPROC/NAPROC separation.    ( version 3.11 )
  !/    08-Oct-2007 : Adding AS CX-Y to W3SRCE par. list. ( version 3.13 )
  !/    22-Feb-2008 : Initialize VGX-Y properly.          ( version 3.13 )
  !/    10-Apr-2008 : Bug fix writing log file (MPI).     ( version 3.13 )
  !/    29-May-2009 : Preparing distribution version.     ( version 3.14 )
  !/    30-Oct-2009 : Implement run-time grid selection.  ( version 3.14 )
  !/                  (W. E. Rogers & T. J. Campbell, NRL)
  !/    30-Oct-2009 : Implement curvilinear grid type.    ( version 3.14 )
  !/                  (W. E. Rogers & T. J. Campbell, NRL)
  !/    29-Mar-2010 : Adding coupling, ice in W3SRCE.     ( version 3.14_SHOM )
  !/    16-May-2010 : Adding transparencies in W3SCRE     ( version 3.14_SHOM )
  !/    23-Jun-2011 : Movable bed bottom friction BT4     ( version 4.04 )
  !/    03-Nov-2011 : Shoreline reflection on unst. grids ( version 4.04 )
  !/    02-Jul-2011 : Update for PALM coupling            ( version 4.07 )
  !/    06-Mar-2012 : Initializing ITEST as needed.       ( version 4.07 )
  !/    02-Jul-2012 : Update for PALM coupling            ( version 4.07 )
  !/    02-Sep-2012 : Clean up of open BC for UG grids    ( version 4.08 )
  !/    03-Sep-2012 : Fix format 902.                     ( version 4.10 )
  !/    07-Dec-2012 : Wrap W3SRCE with TMPn to limit WARN ( version 4.OF )
  !/    10-Dec-2012 : Modify field output MPI for new     ( version 4.OF )
  !/                  structure and smaller memory footprint.
  !/    12-Dec-2012 : Adding SMC grid.  JG_Li             ( version 4.08 )
  !/    26-Dec-2012 : Move FIELD init. to W3GATH.         ( version 4.OF )
  !/    16-Sep-2013 : Add Arctic part for SMC grid.       ( version 4.11 )
  !/    11-Nov-2013 : SMC and rotated grid incorporated in the main
  !/                  trunk                               ( version 4.13 )
  !/    14-Nov-2013 : Remove orphaned work arrays.        ( version 4.13 )
  !/    27-Nov-2013 : Fixes for OpenMP versions.          ( version 4.15 )
  !/    23-May-2014 : Adding ice fluxes to W3SRCE         ( version 5.01 )
  !/    27-May-2014 : Move to OMPG/X switch.              ( version 5.02 )
  !/    24-Apr-2015 : Adding OASIS coupling calls         ( version 5.07 )
  !/                  (M. Accensi & F. Ardhuin, IFREMER)
  !/    27-Aug-2015 : Update for ICEH, ICEF               ( version 5.08 )
  !/    14-Sep-2018 : Remove PALM implementation          ( version 6.06 )
  !/    15-Sep-2020 : Bugfix FIELD allocation. Remove     ( version 7.11 )
  !/                  defunct OMPX switches.
  !/    22-Mar-2021 : Update TAUA, RHOA                   ( version 7.13 )
  !/    06-May-2021 : Use ARCTC and SMCTYPE options. JGLi ( version 7.13 )
  !/    19-Jul-2021 : Momentum and air density support    ( version 7.14 )
  !/    11-Nov-2021 : Remove XYB since it is obsolete     ( version 7.xx )
  !/    13-Sep-2022 : Add OMP for W3NMIN loops. Hide
  !/                  W3NMIN in W3_DEBUGRUN for scaling.  ( version 7.xx )
  !/
  !/    Copyright 2009-2014 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 :
  !
  !  2. Variables and types :
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3WAVE    Subr. Public   Actual wave model.
  !      W3GATH    Subr. Public   Data transpose before propagation.
  !      W3SCAT    Subr. Public   Data transpose after propagation.
  !      W3NMIN    Subr. Public   Calculate minimum number of sea
  !                               points per processor.
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !      W3SETx    Subr. W3xDATMD Point to data structure.
  !
  !      W3UCUR    Subr. W3UPDTMD Interpolate current fields in time.
  !      W3UWND    Subr. W3UPDTMD Interpolate wind fields in time.
  !      W3UINI    Subr. W3UPDTMD Update initial conditions if init.
  !                               with initial wind conditions.
  !      W3UBPT    Subr. W3UPDTMD Update boundary points.
  !      W3UICE    Subr. W3UPDTMD Update ice coverage.
  !      W3ULEV    Subr. W3UPDTMD Transform the wavenumber grid.
  !      W3DDXY    Subr. W3UPDTMD Calculate dirivatives of the depth.
  !      W3DCXY    Subr. W3UPDTMD Calculate dirivatives of the current.
  !
  !      W3MAPn    Subr. W3PROnMD Preparation for  ropagation schemes.
  !      W3XYPn    Subr. W3PROnMD Longitude-latitude ("XY") propagation.
  !      W3KTPn    Subr. W3PROnMD Intra-spectral ("k-theta") propagation.
  !
  !      W3SRCE    Subr. W3SRCEMD Source term integration and calculation.
  !
  !      W3IOGR    Subr. W3IOGRMD Reading/writing model definition file.
  !      W3OUTG    Subr. W3IOGOMD Generate gridded output fields.
  !      W3IOGO    Subr. W3IOGOMD Read/write gridded output.
  !      W3IOPE    Subr. W3IOPOMD Extract point output.
  !      W3IOPO    Subr. W3IOPOMD Read/write point output.
  !      W3IOTR    Subr. W3IOTRMD Process spectral output along tracks.
  !      W3IORS    Subr. W3IORSMD Read/write restart files.
  !      W3IOBC    Subr. W3IOBCMD Read/write boundary conditions.
  !      W3CPRT    Subr. W3IOSFMD Partition spectra.
  !      W3IOSF    Subr.   Id.    Write partitioned spectral data.
  !
  !      STRACE    Subr. W3SERVMD Subroutine tracing.
  !      WWTIME    Subr.   Id.    System time in readable format.
  !      EXTCDE    Subr.   Id.    Program abort.
  !
  !      TICK21    Subr. W3TIMEMD Advance the clock.
  !      DSEC21    Func.   Id.    Difference between times.
  !      STME21    Subr.   Id.    Time in readable format.
  !
  !      MPI_BARRIER, MPI_STARTALL, MPI_WAITALL
  !                Subr.          Basic MPI routines.
  !     ----------------------------------------------------------------
  !
  !  5. Remarks : Call to W3NMIN hidden behind W3_DEBUGRUN. This call
  !               currently only serves to warn when one or more procs
  !               have no active seapoints. It has been hid as this
  !               dramatically increases runtime performance.
  !
  !  6. Switches :
  !
  !       !/SHRD  Switch for shared / distributed memory architecture.
  !       !/DIST  Id.
  !       !/MPI   Id.
  !       !/OMPG  Id.
  !
  !       !/PR1   First order propagation schemes.
  !       !/PR2   ULTIMATE QUICKEST scheme.
  !       !/PR3   Averaged ULTIMATE QUICKEST scheme.
  !       !/SMC   UNO2 scheme on SMC grid.
  !
  !       !/S     Enable subroutine tracing.
  !       !/T     Test output.
  !       !/MPIT  Test output for MPI specific code.
  !
  !  7. Source code :
  !
  !/ ------------------------------------------------------------------- /
  use w3servmd, only : print_memcheck
#ifdef W3_MPI
  USE w3adatmd, ONLY: mpibuf
#endif
  !module default
  implicit none
  !
  PUBLIC
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
 
  SUBROUTINE w3wave ( IMOD, ODAT, TEND, STAMP, NO_OUT &
#ifdef W3_OASIS
       ,ID_LCOMM, TIMEN                 &
#endif
       )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         22-Mar-2021 |
    !/                  +-----------------------------------+
    !/
    !/    17-Mar-1999 : Distributed FORTRAN 77 version.     ( version 1.18 )
    !/    04-Feb-2000 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/                  Major changes to logistics.
    !/    05-Jan-2001 : Bug fix to allow model to run       ( version 2.05 )
    !/                  without output.
    !/    24-Jan-2001 : Flat grid version.                  ( version 2.06 )
    !/    09-Feb-2001 : Third propagation scheme added.     ( version 2.08 )
    !/    23-Feb-2001 : Check for barrier after source
    !/                  terms added ( W3NMIN ).     ( delayed version 2.07 )
    !/    16-Mar-2001 : Fourth propagation scheme added.    ( version 2.09 )
    !/    30-Mar-2001 : Sub-grid obstacles added.           ( version 2.10 )
    !/    23-May-2001 : Barrier added for dry run, changed  ( version 2.10 )
    !/                  declaration of FLIWND.
    !/    10-Dec-2001 : Sub-grid obstacles for UQ schemes.  ( version 2.14 )
    !/    11-Jan-2002 : Sub-grid ice.                       ( version 2.15 )
    !/    24-Jan-2002 : Zero time step dor data ass.        ( version 2.17 )
    !/    09-May-2002 : Switch clean up.                    ( version 2.21 )
    !/    13-Nov-2002 : Add stress vector.                  ( version 3.00 )
    !/    26-Dec-2002 : Moving grid version.                ( version 3.02 )
    !/    01-Aug-2003 : Moving grid GSE correction.         ( version 3.03 )
    !/    07-Oct-2003 : Output options for NN training.     ( version 3.05 )
    !/    29-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    04-Feb-2005 : Add STAMP to par list.              ( version 3.07 )
    !/    04-May-2005 : Change to MPI_COMM_WAVE.            ( version 3.07 )
    !/    28-Jun-2005 : Adding map recalc for W3ULEV call.  ( version 3.07 )
    !/    07-Sep-2005 : Updated boundary conditions.        ( version 3.08 )
    !/    26-Jun-2006 : Add output type 6.                  ( version 3.09 )
    !/    04-Jul-2006 : Consolidate stress arrays.          ( version 3.09 )
    !/    18-Oct-2006 : Partitioned spectral data output.   ( version 3.10 )
    !/    02-Feb-2007 : Add FLAGST test.                    ( version 3.10 )
    !/    02-Apr-2007 : Add partitioned field data.         ( version 3.11 )
    !/                  Improve MPI_WAITALL call tests/allocations.
    !/    07-May-2007 : Bug fix SKIP_O treatment.           ( version 3.11 )
    !/    17-May-2007 : Adding NTPROC/NAPROC separation.    ( version 3.11 )
    !/    08-Oct-2007 : Adding AS CX-Y to W3SRCE par. list. ( version 3.13 )
    !/    22-Feb-2008 : Initialize VGX-Y properly.          ( version 3.13 )
    !/    10-Apr-2008 : Bug fix writing log file (MPI).     ( version 3.13 )
    !/    30-Oct-2009 : Implement run-time grid selection.  ( version 3.14 )
    !/                  (W. E. Rogers & T. J. Campbell, NRL)
    !/    30-Oct-2009 : Implement curvilinear grid type.    ( version 3.14 )
    !/                  (W. E. Rogers & T. J. Campbell, NRL)
    !/    31-Mar-2010 : Add reflections                     ( version 3.14.4 )
    !/    29-Oct-2010 : Implement unstructured grids        ( version 3.14.4 )
    !/                  (A. Roland and F. Ardhuin)
    !/    06-Mar-2011 : Output of max. CFL (F.Ardhuin)      ( version 3.14.4 )
    !/    05-Apr-2011 : Implement iteration for DTMAX <1s   ( version 3.14.4 )
    !/    02-Jul-2012 : Update for PALM coupling            ( version 4.07 )
    !/    02-Sep-2012 : Clean up of open BC for UG grids    ( version 4.08 )
    !/    03-Sep-2012 : Fix format 902.                     ( version 4.10 )
    !/    10-Dec-2012 : Modify field output MPI for new     ( version 4.OF )
    !/                  structure and smaller memory footprint.
    !/    16-Nov-2013 : Allows reflection on curvi. grids   ( version 4.13 )
    !/    27-Nov-2013 : Fixes for OpenMP versions.          ( version 4.15 )
    !/    23-May-2014 : Adding ice fluxes to W3SRCE         ( version 5.01 )
    !/    27-May-2014 : Move to OMPG/X switch.              ( version 5.02 )
    !/    24-Apr-2015 : Adding OASIS coupling calls         ( version 5.07 )
    !/                  (M. Accensi & F. Ardhuin, IFREMER)
    !/    27-Aug-2015 : Update for ICEH, ICEF               ( version 5.10 )
    !/    31-Mar-2016 : Current option for smc grid.        ( version 5.18 )
    !/    06-Jun-2018 : Add PDLIB/MEMCHECK/SETUP/NETCDF_QAD/TIMING
    !/                  OASIS/DEBUGINIT/DEBUGSRC/DEBUGRUN/DEBUGCOH
    !/                  DEBUGIOBP/DEBUGIOBC                 ( version 6.04 )
    !/    14-Sep-2018 : Remove PALM implementation          ( version 6.06 )
    !/    25-Sep-2020 : Oasis coupling at T+0               ( version 7.10 )
    !/    22-Mar-2021 : Update TAUA, RHOA                   ( version 7.13 )
    !/    06-May-2021 : Use ARCTC and SMCTYPE options. JGLi ( version 7.13 )
    !/
    !  1. Purpose :
    !
    !     Run WAVEWATCH III for a given time interval.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       IMOD    Int.   I   Model number.
    !       TEND    I.A.   I   Ending time of integration.
    !       STAMP   Log.   I   WRITE(*,*)time stamp (optional, defaults to T).
    !       NO_OUT  Log.   I   Skip output (optional, defaults to F).
    !                          Skip at ending time only!
    !     ----------------------------------------------------------------
    !
    !     Local parameters : Flags
    !     ----------------------------------------------------------------
    !       FLOUTG  Log.  Flag for running W3OUTG.
    !       FLPART  Log.  Flag for running W3CPRT.
    !       FLZERO  Log.  Flag for zero time interval.
    !       FLAG0   Log.  Flag for processors without tasks.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !     Any program shell or integrated model which uses WAVEWATCH III.
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !     - Currents are updated before winds as currents are used in wind
    !       and USTAR processing.
    !     - Ice and water levels can be updated only once per call.
    !     - If ice or water level time are undefined, the update
    !       takes place asap, otherwise around the "half-way point"
    !       betweem the old and new times.
    !     - To increase accuracy, the calculation of the intra-spectral
    !       propagation is performed in two parts around the spatial propagation.
    !
    !  8. Structure :
    !
    !     -----------------------------------------------------------
    !       0.  Initializations
    !         a Point to data structures
    !         b Subroutine tracing
    !         c Local parameter initialization
    !         d Test output
    !       1.  Check the consistency of the input.
    !         a Ending time versus initial time.
    !         b Water level time.
    !         c Current time interval.
    !         d Wind time interval.
    !         e Ice time.
    !       2.  Determine next time from ending and output
    !           time and get corresponding time step.
    !       3.  Loop over time steps (see below).
    !       4.  Perform output to file if requested.
    !         a Check if time is output time.
    !         b Processing and MPP preparations.  ( W3CPRT, W3OUTG )
    !         c Reset next output time.
    !        -------------- loop over output types ------------------
    !         d Perform output.                           ( W3IOxx )
    !         e Update next output time.
    !        -------------------- end loop --------------------------
    !       5.  Update log file.
    !       6.  If time is not ending time, branch back to 2.
    !     -----------------------------------------------------------
    !
    !      Section 3.
    !     ----------------------------------------------------------
    !       3.1  Interpolate winds and currents. ( W3UCUR, W3DCXY )
    !                                                    ( W3UWND )
    !                                                    ( W3UINI )
    !       3.2  Update boundary conditions.     ( W3IOBC, W3UBPT )
    !       3.3  Update ice coverage (if new ice map).   ( W3UICE )
    !       3.4  Transform grid (if new water level).    ( W3ULEV )
    !       3.5  Update maps and dirivatives.    ( W3MAPn, W3DDXY )
    !                                            ( W3NMIN, W3UTRN )
    !            Update grid advection vector.
    !       3.6  Perform propagation
    !          a Preparations.
    !          b Intra spectral part 1.                  ( W3KTPn )
    !          c Longitude-latitude       ( W3GATH, W3XYPn W3SCAT )
    !          b Intra spectral part 2.                  ( W3KTPn )
    !       3.7  Calculate and integrate source terms.   ( W3SRCE )
    !       3.8  Update global time step.
    !     ----------------------------------------------------------
    !
    !  9. Switches :
    !
    !     See module documentation.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    !/
    USE w3gdatmd
    USE w3wdatmd
    USE w3adatmd
    USE w3idatmd
    USE w3odatmd
    !/
    USE w3updtmd
    USE w3srcemd
#ifdef W3_PR1
    USE w3pro1md
#endif
#ifdef W3_PR2
    USE w3pro2md
#endif
#ifdef W3_PR3
    USE w3pro3md
#endif
#ifdef W3_SMC
    USE w3psmcmd
#endif
    !
#ifdef W3_PR1
    USE w3profsmd
#endif
#ifdef W3_PR2
    USE w3profsmd
#endif
#ifdef W3_PR3
    USE w3profsmd
#endif
    !/
    USE w3triamd
    USE w3iogrmd
    USE w3iogomd
    USE w3iopomd
    USE w3iotrmd
    USE w3iorsmd
    USE w3iobcmd
    USE w3iosfmd
#ifdef W3_PDLIB
    USE pdlib_w3profsmd, only : apply_boundary_condition_va
    USE pdlib_w3profsmd, only : pdlib_w3xypug, pdlib_w3xypug_block_implicit, pdlib_w3xypug_block_explicit
    USE pdlib_w3profsmd, only : all_va_integral_print, all_vaold_integral_print, all_field_integral_print
    USE w3parall, only : pdlib_nseal, pdlib_nsealm
    USE yownodepool, only: npa, iplg, np
#endif
    !/
    USE w3servmd
    USE w3timemd
#ifdef W3_IC3
    USE w3sic3md
#endif
#ifdef W3_IS2
    USE w3sis2md
#endif
#ifdef W3_UOST
    USE w3uostmd, ONLY: uost_setgrid
#endif
    USE w3parall, ONLY : init_get_isea
#ifdef W3_SETUP
    USE w3wavset, only : wave_setup_computation
#endif
 
#ifdef W3_OASIS
    USE w3oacpmd, ONLY: id_oasis_time, cplt0
#endif
#ifdef W3_OASOCM
    USE w3ogcmmd, ONLY: snd_fields_to_ocean
#endif
#ifdef W3_OASACM
    USE w3agcmmd, ONLY: snd_fields_to_atmos
#endif
#ifdef W3_OASICM
    USE w3igcmmd, ONLY: snd_fields_to_ice
#endif
 
#ifdef W3_PDLIB
    USE pdlib_field_vec, only : do_output_exchanges
    USE pdlib_w3profsmd, ONLY: aspar_jac, aspar_diag_all, b_jac
    USE w3parall, only : lsloc
#endif
#ifdef W3_TIMINGS
    USE w3parall, only : print_my_time
#endif
    !
#ifdef W3_MPI
    include "mpif.h"
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)           :: IMOD, TEND(2),ODAT(35)
    LOGICAL, INTENT(IN), OPTIONAL :: STAMP, NO_OUT
#ifdef W3_OASIS
    INTEGER, INTENT(IN), OPTIONAL :: ID_LCOMM
    INTEGER, INTENT(IN), OPTIONAL :: TIMEN(2)
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters :
    !/
#ifdef W3_T
    INTEGER                 :: ILEN
#endif
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    INTEGER                 :: IP
    INTEGER                 :: TCALC(2), IT, IT0, NT, ITEST,        &
         ITLOC, ITLOCH, NTLOC, ISEA, JSEA,    &
         IX, IY, ISPEC, J, TOUT(2), TLST(2),  &
         REFLED(6), IK, ITH, IS, NKCFL
    INTEGER                 :: ISP, IP_glob
    INTEGER                 :: TTEST(2),DTTEST
    REAL                    :: ICEDAVE
    !
#ifdef W3_MPI
    LOGICAL                 :: SBSED
#endif
#ifdef W3_SEC1
    INTEGER                 :: ISEC1
#endif
#ifdef W3_SBS
    INTEGER                 :: JJ, NDSOFLG
#endif
#ifdef W3_MPI
    INTEGER                 :: IERR_MPI, NRQMAX
    INTEGER, ALLOCATABLE    :: STATCO(:,:), STATIO(:,:)
#endif
    INTEGER                 :: IXrel
    REAL                    :: DTTST, DTTST1, DTTST2, DTTST3,       &
         DTL0, DTI0, DTR0, DTI10, DTI50,      &
         DTGA, DTG, DTGpre, DTRES,            &
         FAC, VGX, VGY, FACK, FACTH,          &
         FACX, XXX, REFLEC(4),                &
         DELX, DELY, DELA, DEPTH, D50, PSIC
    REAL                     :: VSioDummy(NSPEC), VDioDummy(NSPEC), VAoldDummy(NSPEC)
    LOGICAL                  :: SHAVETOTioDummy
#ifdef W3_SEC1
    REAL                    :: DTGTEMP
#endif
    !
    REAL, ALLOCATABLE       :: FIELD(:)
    REAL                    :: TMP1(4), TMP2(3), TMP3(2), TMP4(2)
#ifdef W3_IC3
    REAL, ALLOCATABLE       :: WN_I(:)
#endif
#ifdef W3_REFRX
    REAL, ALLOCATABLE       :: CIK(:)
#endif
    !
    ! Orphaned arrays from old data structure
    !
    REAL, ALLOCATABLE       :: TAUWX(:), TAUWY(:)
    !
    LOGICAL                 :: FLACT, FLZERO, FLFRST, FLMAP, TSTAMP,&
         SKIP_O, FLAG_O, FLDDIR, READBC,      &
         FLAG0 = .false., floutg, flpfld,     &
         flpart, local, floutg2
    !
#ifdef W3_MPI
    LOGICAL                 :: FLGMPI(0:8)
#endif
#ifdef W3_IC3
    REAL                    :: FIXEDVISC,FIXEDDENS,FIXEDELAS
    REAL                    :: USE_CHENG, USE_CGICE, HICE
#endif
    LOGICAL                 :: UGDTUPDATE    ! true if time step should be updated for UG schemes
    CHARACTER(LEN=8)        :: STTIME
    CHARACTER(LEN=21)       :: IDACT
    CHARACTER(LEN=16)       :: OUTID
    CHARACTER(LEN=23)       :: IDTIME
    INTEGER eIOBP
    INTEGER ITH_F
#ifdef W3_PDLIB
    REAL ::             VS_SPEC(NSPEC)
    REAL ::             VD_SPEC(NSPEC)
#endif
    !
#ifdef W3_SBS
    CHARACTER(LEN=30)       :: FOUTNAME
#endif
    !
#ifdef W3_T
    REAL             :: INDSORT(NSEA), DTCFL1(NSEA)
#endif
    !/
#ifdef W3_SMC
    !Li   Temperature spectra for Arctic boundary update.
    REAL, ALLOCATABLE       :: BACSPEC(:)
    REAL                    :: BACANGL
#endif
    integer :: memunit
    !/ ------------------------------------------------------------------- /
    ! 0.  Initializations
    !
    xxx = undef
    memunit = 40000+iaproc
    ! 0.a Set pointers to data structure
    !
#ifdef W3_COU
    screen   =  333
#endif
    !
    IF ( ioutp  .NE. imod ) CALL w3seto ( imod, ndse, ndst )
    IF ( igrid  .NE. imod ) CALL w3setg ( imod, ndse, ndst )
    IF ( iwdata .NE. imod ) CALL w3setw ( imod, ndse, ndst )
    IF ( iadata .NE. imod ) CALL w3seta ( imod, ndse, ndst )
    IF ( iidata .NE. imod ) CALL w3seti ( imod, ndse, ndst )
#ifdef W3_UOST
    CALL uost_setgrid(imod)
#endif
 
#ifdef W3_DEBUGCOH
    CALL all_va_integral_print(imod, "W3WAVEMD, step 1", 1)
#endif
 
    !
    ALLOCATE(tauwx(nseal), tauwy(nseal))
#ifdef W3_REFRX
    ALLOCATE(cik(nseal))
#endif
    !
    IF ( PRESENT(stamp) ) THEN
      tstamp = stamp
    ELSE
      tstamp = .true.
    END IF
    !
    IF ( PRESENT(no_out) ) THEN
      skip_o = no_out
    ELSE
      skip_o = .false.
    END IF
#ifdef W3_DEBUGCOH
    CALL all_va_integral_print(imod, "W3WAVEMD, step 2", 1)
#endif
    !
    ! 0.b Subroutine tracing
    !
#ifdef W3_S
    CALL strace (ient, 'W3WAVE')
#endif
    !
    !
    ! 0.c Local parameter initialization
    !
    ipass  = ipass + 1
    idact  = '                 '
    outid  = '           '
    flact  = itime .EQ. 0
    flmap  = itime .EQ. 0
    flddir = itime .EQ. 0 .AND. ( flcth .OR. fsrefraction .OR. flck .OR. fsfreqshift )
    !
    flpfld = .false.
    DO j=1,noge(4)
      flpfld = flpfld .OR. flogrd(4,j) .OR. flogr2(4,j)
    END DO
    !
    IF ( iaproc .EQ. naplog ) backspace( ndso )
    !
    IF ( flcold ) THEN
      dtdyn = 0.
      fcut  = sig(nk) * tpiinv
    END IF
    !
    IF( gtype .EQ. smctype ) THEN
      j = 1
#ifdef W3_SMC
      !!Li   Use sea point only field for SMC grid.
      ALLOCATE ( field(ncel) )
#endif
    ELSE
      ALLOCATE ( field(1-ny:ny*(nx+2)) )
    ENDIF
    !
    local   = iaproc .LE. naproc
    ugdtupdate = .false.
    IF (flagll) THEN
      facx   =  1./(dera * radius)
    ELSE
      facx   =  1.
    END IF
    !
#ifdef W3_SBS
    ndsoflg = 99
#endif
#ifdef W3_MPI
    sbsed = .false.
#endif
#ifdef W3_SBS
    sbsed = .true.
#endif
    !
    tauwx  = 0.
    tauwy  = 0.
    !
    ! 0.d Test output
    !
#ifdef W3_T
    ilen   = len_trim(filext)
    WRITE (ndst,9000) imod, filext(:ilen), tend
#endif
    !
    ! 1.  Check the consistency of the input ----------------------------- /
    ! 1.a Ending time versus initial time
    !
    dttst  = dsec21( time , tend )
    flzero = dttst .EQ. 0.
#ifdef W3_T
    WRITE (ndst,9010) dttst, flzero
#endif
    IF ( dttst .LT. 0. ) THEN
      IF ( iaproc .EQ. naperr ) WRITE (ndse,1000)
      CALL extcde ( 1 )
    END IF
    !
    ! 1.b Water level time
    !
    IF ( fllev ) THEN
      IF ( tlev(1) .GE. 0. ) THEN
        dtl0   = dsec21( tlev , tln )
      ELSE
        dtl0   = 1.
      END IF
#ifdef W3_T
      WRITE (ndst,9011) dtl0
#endif
      IF ( dtl0 .LT. 0. ) THEN
        IF ( iaproc .EQ. naperr ) WRITE (ndse,1001)
        CALL extcde ( 2 )
      END IF
    ELSE
      dtl0   = 0.
    END IF
#ifdef W3_DEBUGCOH
    CALL all_va_integral_print(imod, "W3WAVEMD, step 4", 1)
#endif
    !
    ! 1.c Current interval
    !
    IF ( flcur ) THEN
      dttst1 = dsec21( tc0 , tcn )
      dttst2 = dsec21( tc0 , time )
      dttst3 = dsec21( tend , tcn )
#ifdef W3_T
      WRITE (ndst,9012) dttst1, dttst2, dttst3
#endif
      IF ( dttst1.LT.0. .OR. dttst2.LT.0. .OR. dttst3.LT.0. ) THEN
        IF ( iaproc .EQ. naperr ) WRITE (ndse,1002)
        CALL extcde ( 3 )
      END IF
      IF ( dttst2.EQ.0..AND. itime.EQ.0 ) THEN
        idact(7:7) = 'F'
        tofrst = time
      END IF
    END IF
    !
    ! 1.d Wind interval
    !
    IF ( flwind ) THEN
      dttst1 = dsec21( tw0 , twn )
      dttst2 = dsec21( tw0 , time )
      dttst3 = dsec21( tend , twn )
#ifdef W3_T
      WRITE (ndst,9013) dttst1, dttst2, dttst3
#endif
      IF ( dttst1.LT.0. .OR. dttst2.LT.0. .OR. dttst3.LT.0. ) THEN
        IF ( iaproc .EQ. naperr ) WRITE (ndse,1003)
        CALL extcde ( 4 )
      END IF
      IF ( dttst2.EQ.0..AND. itime.EQ.0 ) THEN
        idact(3:3) = 'F'
        tofrst = time
      END IF
    END IF
#ifdef W3_DEBUGCOH
    CALL all_va_integral_print(imod, "W3WAVEMD, step 5", 1)
#endif
    !
    ! 1.e Ice concentration interval
    !
    IF ( flice ) THEN
      IF ( tice(1) .GE. 0 ) THEN
        dti0   = dsec21( tice , tin )
      ELSE
        dti0   = 1.
      END IF
#ifdef W3_T
      WRITE (ndst,9014) dti0
#endif
      IF ( dti0 .LT. 0. ) THEN
        IF ( iaproc .EQ. naperr ) WRITE (ndse,1004)
        CALL extcde ( 5 )
      END IF
    ELSE
      dti0   = 0.
    END IF
#ifdef W3_DEBUGCOH
    CALL all_va_integral_print(imod, "W3WAVEMD, step 6", 1)
#endif
    !
    ! 1.f Momentum interval
    !
    IF ( fltaua ) THEN
      dttst1 = dsec21( tu0 , tun )
      dttst2 = dsec21( tu0 , time )
      dttst3 = dsec21( tend , tun )
#ifdef W3_T
      WRITE (ndst,9017) dttst1, dttst2, dttst3
#endif
      IF ( dttst1.LT.0. .OR. dttst2.LT.0. .OR. dttst3.LT.0. ) THEN
        IF ( iaproc .EQ. naperr ) WRITE (ndse,1007)
        CALL extcde ( 3 )
      END IF
      IF ( dttst2.EQ.0..AND. itime.EQ.0 ) THEN
        idact(9:9) = 'F'
        tofrst = time
      END IF
    END IF
    !
    ! 1.g Air density time
    !
    IF ( flrhoa ) THEN
      dttst1 = dsec21( tr0 , trn )
      dttst2 = dsec21( tr0 , time )
      dttst3 = dsec21( tend , trn )
#ifdef W3_T
      WRITE (ndst,9018) dttst1, dttst2, dttst3
#endif
      IF ( dttst1.LT.0. .OR. dttst2.LT.0. .OR. dttst3.LT.0. ) THEN
        IF ( iaproc .EQ. naperr ) WRITE (ndse,1008)
        CALL extcde ( 2 )
      END IF
      IF ( dttst2.EQ.0..AND. itime.EQ.0 ) THEN
        idact(11:11) = 'F'
        tofrst = time
      END IF
    END IF
    !
    ! 1.e Ice thickness interval
    !
    IF ( flic1 ) THEN
      IF ( tic1(1) .GE. 0 ) THEN
        dti10   = dsec21( tic1 , ti1 )
      ELSE
        dti10   = 1.
      END IF
#ifdef W3_T
      WRITE (ndst,9015) dti10
#endif
      IF ( dti10 .LT. 0. ) THEN
        IF ( iaproc .EQ. naperr ) WRITE (ndse,1005)
        CALL extcde ( 5 )
      END IF
    ELSE
      dti10   = 0.
    END IF
    !
    ! 1.e Ice floe interval
    !
#ifdef W3_IS2
    IF ( flic5 ) THEN
      IF ( tic5(1) .GE. 0 ) THEN
        dti50   = dsec21( tic5 , ti5 )
      ELSE
        dti50   = 1.
      END IF
#ifdef W3_T
      WRITE (ndst,9016) dti50
#endif
      IF ( dti50 .LT. 0. ) THEN
        IF ( iaproc .EQ. naperr ) WRITE (ndse,1006)
        CALL extcde ( 5 )
      END IF
    ELSE
      dti50   = 0.
    END IF
#endif
    !
    ! 2.  Determine next time from ending and output --------------------- /
    !     time and get corresponding time step.
    !
    flfrst = .true.
    DO
#ifdef W3_TIMINGS
      CALL print_my_time("First entry in the TIME LOOP")
#endif
 
#ifdef W3_DEBUGCOH
      CALL all_va_integral_print(imod, "W3WAVEMD, step 6.1", 1)
#endif
      !
      !
      ! 2.a Pre-calculate table for IC3 ------------------------------------ /
#ifdef W3_IC3
      use_cheng=ic3pars(9)
      IF( use_cheng==1.0 )THEN
        fixedvisc=ic3pars(14)
        fixeddens=ic3pars(15)
        fixedelas=ic3pars(16)
        IF ( (fixedvisc.LT.0.0).OR.(fixeddens.LT.0.0) .OR. (fixedelas.LT.0.0) ) THEN
          IF ( iaproc .EQ. naperr )                          &
               WRITE(ndse,*)'Cheng method requires stationary',   &
               ' and uniform rheology from namelist.'
          CALL extcde(2)
        END IF
        IF (calledic3table==0) THEN
          CALL ic3table_cheng(fixedvisc,fixeddens,fixedelas)
          calledic3table = 1
        ENDIF
      ENDIF
#endif
 
      ! 2.b Update group velocity and wavenumber from ice parameters ------- /
      !     from W3SIC3MD module. ------------------------------------------ /
      !     Note: "IF FLFRST" can be added for efficiency, but testing req'd
 
      jsea=1 ! no switch (intentional)
 
#ifdef W3_IC3
      use_cgice=ic3pars(12)
      IF ( use_cgice==1.0 ) THEN
        IF ( iaproc .EQ. naperr ) WRITE(screen,920)
#endif
 
#ifdef W3_IC3
        DO jsea=1,nseal
#endif
#ifdef W3_DIST
          isea   = iaproc + (jsea-1)*naproc
#endif
#ifdef W3_SHRD
          isea   = jsea
#endif
#ifdef W3_IC3
          ALLOCATE(wn_i(SIZE(wn(:,isea))))
          wn_i(:) = 0.
          depth  = max( dmin , dw(isea) )
          ix     = mapsf(isea,1)
          iy     = mapsf(isea,2)
#endif
 
          ! 2.b.1 Using Cheng method: requires stationary/uniform rheology.
          !       However, ice thickness may be input by either method
 
#ifdef W3_IC3
          IF ( use_cheng==1.0 ) THEN
            IF (flic1) THEN
              hice=icep1(ix,iy)
            ELSEIF (ic3pars(13).GE.0.0)THEN
              hice=ic3pars(13)
            ELSE
              IF ( iaproc .EQ. naperr ) then
                WRITE(ndse,*)'ICE THICKNESS NOT AVAILABLE FOR CG CALC'
              end if
              CALL extcde(2)
            ENDIF
            IF (hice > 0.0) THEN ! non-zero ice
              CALL w3ic3wncg_cheng(wn(:,isea),wn_i(:), cg(:,isea),hice,fixedvisc, &
                   fixeddens, fixedelas, depth)
            END IF ! non-zero ice
#endif
 
#ifdef W3_IC3
          ELSE ! not using Cheng method
#endif
            ! 2.b.2 If not using Cheng method: require FLIC1 to FLIC4 (not strictly
            !       necesssary, but makes code simpler)
 
#ifdef W3_IC3
            IF (flic1.AND.flic2.AND.flic3.AND.flic4) THEN
              IF (icep1(ix,iy)>0.0) THEN ! non-zero ice
                CALL w3ic3wncg_v1(wn(:,isea),wn_i(:), cg(:,isea),icep1(ix,iy), &
                     icep2(ix,iy), icep3(ix,iy),icep4(ix,iy),depth)
              END IF ! non-zero ice
            ELSE
              IF ( iaproc .EQ. naperr ) then
                WRITE(ndse,*)'ICE PARAMETERS NOT AVAILABLE FOR CG CALC'
              end if
              CALL extcde(2)
            END IF
          ENDIF ! IF USE_CHENG...
#endif
 
#ifdef W3_IC3
          DEALLOCATE(wn_i)
        END DO ! DO JSEA=1,NSEAL
      END IF !  IF USE_CGICE ...
#endif
      !
      IF ( tofrst(1) .GT. 0 ) THEN
        dttst  = dsec21( tend , tofrst )
      ELSE
        dttst  = 0.
      ENDIF
      !
      IF ( dttst.GE.0. ) THEN
        tcalc = tend
      ELSE
        tcalc = tofrst
      END IF
      !
      dttst  = dsec21( time , tcalc )
      nt     = 1 + int( dttst / dtmax - 0.001 )
      dtga   = dttst / real(nt)
      IF ( dttst .EQ. 0. ) THEN
        it0    = 0
        IF ( .NOT.flzero ) itime  = itime - 1
        nt     = 0
      ELSE
        it0    = 1
      END IF
      call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE')
      !
#ifdef W3_T
      WRITE (ndst,9020) it0, nt, dtga
#endif
      !
      ! ==================================================================== /
      !
      ! 3.  Loop over time steps
      !
      dtres  = 0.
 
      !
      DO it = it0, nt
#ifdef W3_TIMINGS
        CALL print_my_time("Begin of IT loop")
#endif
#ifdef W3_SETUP
        CALL wave_setup_computation
#endif
        ! copy old values
#ifdef W3_PDLIB
        DO ip=1,nseal
          DO ispec=1,nspec
            vaold(ispec,ip)=va(ispec,ip)
          END DO
        END DO
#endif
        !
#ifdef W3_DEBUGCOH
        CALL all_va_integral_print(imod, "Beginning time loop", 1)
#endif
#ifdef W3_TIMINGS
        CALL print_my_time("After assigning VAOLD")
#endif
        !
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 0')
        !
        itime  = itime + 1
        !
        dtg    = real(nint(dtga+dtres+0.0001))
        dtres  = dtres + dtga - dtg
        IF ( abs(dtres) .LT. 0.001 ) dtres  = 0.
        CALL tick21 ( time , dtg )
        !
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 1')
 
        IF ( tstamp .AND. screen.NE.ndso .AND. iaproc.EQ.napout ) THEN
          CALL wwtime ( sttime )
          CALL stme21 ( time , idtime )
          WRITE (screen,950) idtime, sttime
        END IF
        !
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 2')
 
        vgx = 0.
        vgy = 0.
        IF(inflags1(10)) THEN
          dttst1 = dsec21( time, tgn )
          dttst2 = dsec21( tg0, tgn )
          fac    = dttst1 / max( 1. , dttst2 )
          vgx    = (fac*ga0+(1.-fac)*gan) * cos(fac*gd0+(1.-fac)*gdn)
          vgy    = (fac*ga0+(1.-fac)*gan) * sin(fac*gd0+(1.-fac)*gdn)
        END IF
#ifdef W3_TIMINGS
        CALL print_my_time("After VGX/VGY assignation")
#endif
        !
#ifdef W3_T
        WRITE (ndst,9021) itime, it, time, flmap, flddir, vgx, vgy, dtg, dtres
#endif
        !
        ! 3.1 Interpolate winds, currents, and momentum.
        !     (Initialize wave fields with winds)
        !
#ifdef W3_DEBUGDCXDX
        WRITE(740+iaproc,*) 'Debug DCXDX FLCUR=', flcur
#endif
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 3a')
 
        IF ( flcur  ) THEN
#ifdef W3_DEBUGCOH
          CALL all_va_integral_print(imod, "Before UCUR", 1)
#endif
#ifdef W3_TIMINGS
          CALL print_my_time("W3WAVE, step 6.4.1")
#endif
          CALL w3ucur ( flfrst )
 
          call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 3b')
 
          IF (gtype .EQ. smctype) THEN
            ix = 1
#ifdef W3_SMC
            !!Li  Use new sub for DCXDX/Y and DCYDX/Y assignment.
            CALL smcdcxy
#endif
          ELSE IF (gtype .EQ. ungtype) THEN
#ifdef W3_DEBUGDCXDX
            WRITE(740+iaproc,*) 'Before call to UG_GRADIENT for assigning DCXDX/DCXDY array'
#endif
            CALL ug_gradients(cx, dcxdx, dcxdy)
            CALL ug_gradients(cy, dcydx, dcydy)
            ugdtupdate=.true.
            cflxymax = 0.
          ELSE
            CALL w3dzxy(cx(1:ubound(cx,1)),'m/s',dcxdx, dcxdy) !CX GRADIENT
            CALL w3dzxy(cy(1:ubound(cy,1)),'m/s',dcydx, dcydy) !CY GRADIENT
          ENDIF  !! End GTYPE
          !
          call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 4')
          !
        ELSE IF ( flfrst ) THEN
          ugdtupdate=.true.
          cflxymax = 0.
          cx = 0.
          cy = 0.
        END IF ! FLCUR
#ifdef W3_TIMINGS
        CALL print_my_time("After CX/CY assignation")
#endif
        !
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 5')
 
        IF ( flwind ) THEN
          IF ( flfrst ) asf = 1.
          CALL w3uwnd ( flfrst, vgx, vgy )
        ELSE IF ( flfrst ) THEN
          u10    = 0.01
          u10d   = 0.
          ust    = 0.05
          ustdir = 0.05
        END IF
 
#ifdef W3_DEBUGRUN
        DO jsea = 1, nseal
          DO is = 1, nspec
            IF (va(is, jsea) .LT. 0.) THEN
              WRITE(740+iaproc,*) 'TEST W3WAVE 5', va(is,jsea)
              CALL flush(740+iaproc)
            ENDIF
          ENDDO
        ENDDO
        IF (sum(va) .NE. sum(va)) THEN
          WRITE(740+iaproc,*) 'NAN in ACTION 5', ix, iy, sum(va)
          CALL flush(740+iaproc)
          stop
        ENDIF
#endif
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 6')
 
#ifdef W3_TIMINGS
        CALL print_my_time("After U10, etc. assignation")
#endif
        !
#ifdef W3_DEBUGCOH
        CALL all_va_integral_print(imod, "Before call to W3UINI", 1)
#endif
#ifdef W3_TIMINGS
        CALL print_my_time("Before call W3UINI")
#endif
        IF ( fliwnd .AND. local ) CALL w3uini ( va )
        !
        IF ( fltaua ) THEN
          CALL w3utau ( flfrst )
        ELSE IF ( flfrst ) THEN
          taua    = 0.01
          tauadir = 0.
        END IF
        !
        IF ( flrhoa ) THEN
          CALL w3urho ( flfrst )
        ELSE IF ( flfrst ) THEN
          rhoair = dair
        END IF
        !
        ! 3.2 Update boundary conditions if boundary flag is true (FLBPI)
        !
#ifdef W3_DEBUGCOH
        CALL all_va_integral_print(imod, "Before boundary update", 1)
#endif
#ifdef W3_TIMINGS
        CALL print_my_time("Before boundary update")
#endif
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 7')
 
        IF ( flbpi .AND. local ) THEN
          !
          DO
            IF ( tbpin(1) .EQ. -1 ) THEN
              readbc = .true.
              idact(1:1) = 'F'
            ELSE
              readbc = dsec21(time,tbpin).LT.0.
              IF (readbc.AND.idact(1:1).EQ.' ') idact(1:1) = 'X'
            END IF
            flact  = readbc .OR. flact
 
            IF ( readbc ) THEN
              CALL w3iobc ( 'READ', nds(9), tbpi0, tbpin, itest, imod )
              IF ( itest .NE. 1 ) CALL w3ubpt
            ELSE
              itest  = 0
            END IF
            IF ( itest .LT. 0 ) idact(1:1) = 'L'
            IF ( itest .GT. 0 ) idact(1:1) = ' '
            IF ( .NOT. (readbc.AND.flbpi) ) EXIT
          END DO
 
        END IF
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 7')
 
#ifdef W3_PDLIB
        CALL apply_boundary_condition_va
#ifdef W3_DEBUGCOH
        CALL all_va_integral_print(imod, "After FLBPI and LOCAL", 1)
#endif
#endif
#ifdef W3_TIMINGS
        CALL print_my_time("After FLBPI and LOCAL")
#endif
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 8')
        !
        ! 3.3.1 Update ice coverage (if new ice map).
        !     Need to be run on output nodes too, to update MAPSTx
        !
        IF ( flice .AND. dti0.NE.0. ) THEN
          !
          IF ( tice(1).GE.0 ) THEN
            IF ( dti0 .LT. 0. ) THEN
              idact(13:13) = 'B'
            ELSE
              dttst  = dsec21( time, tin )
              IF ( dttst .LE. 0.5*dti0 ) idact(13:13) = 'U'
            END IF
          ELSE
            idact(13:13) = 'I'
          END IF
          !
          IF ( idact(13:13).NE.' ' ) THEN
            CALL w3uice ( va )
            dti0   = 0.
            flact  = .true.
            flmap  = .true.
          END IF
        END IF
#ifdef W3_DEBUGCOH
        CALL all_va_integral_print(imod, "After FLICE and DTI0", 1)
#endif
#ifdef W3_TIMINGS
        CALL print_my_time("After FLICE and DTI0")
#endif
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 9')
        !
        ! 3.3.2 Update ice thickness
        !
        IF ( flic1 .AND. dti10.NE.0. ) THEN
          !
          IF ( tic1(1).GE.0 ) THEN
            IF ( dti10 .LT. 0. ) THEN
              idact(15:15) = 'B'
            ELSE
              dttst  = dsec21( time, ti1 )
              IF ( dttst .LE. 0.5*dti10 ) idact(15:15) = 'U'
            END IF
          ELSE
            idact(15:15) = 'I'
          END IF
 
          !
          IF ( idact(15:15).NE.' ' ) THEN
            CALL w3uic1 ( flfrst )
            dti10   = 0.
            flact  = .true.
            flmap  = .true.
          END IF
          !
        END IF
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 10')
        !
        ! 3.3.3 Update ice floe diameter
        !
#ifdef W3_IS2
        IF ( flic5 .AND. dti50.NE.0. ) THEN
          !
          IF ( tic5(1).GE.0 ) THEN
            IF ( dti50 .LT. 0. ) THEN
              idact(18:18) = 'B'
            ELSE
              dttst  = dsec21( time, ti5 )
              IF ( dttst .LE. 0.5*dti50 ) idact(18:18) = 'U'
            END IF
          ELSE
            idact(18:18) = 'I'
          END IF
          !
          IF ( idact(18:18).NE.' ' ) THEN
            CALL w3uic5( flfrst )
            dti50   = 0.
            flact  = .true.
            flmap  = .true.
          END IF
          !
        END IF
#endif
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 11a')
        !
        ! 3.4 Transform grid (if new water level).
        !
        IF ( fllev .AND. dtl0 .NE.0. ) THEN
          !
          IF ( tlev(1) .GE. 0 ) THEN
            IF ( dtl0 .LT. 0. ) THEN
              idact(5:5) = 'B'
            ELSE
              dttst  = dsec21( time, tln )
              IF ( dttst .LE. 0.5*dtl0 ) idact(5:5) = 'U'
            END IF
          ELSE
            idact(5:5) = 'I'
          END IF
          !
          IF ( idact(5:5).NE.' ' ) THEN
 
            CALL w3ulev ( va, va )
 
            ugdtupdate=.true.
            cflxymax = 0.
            dtl0   = 0.
            flact  = .true.
            flmap  = .true.
            flddir = flddir .OR. flcth .OR. fsrefraction .OR. flck .OR. fsfreqshift
          END IF
        END IF
#ifdef W3_DEBUGCOH
        CALL all_va_integral_print(imod, "After FFLEV and DTL0", 1)
#endif
#ifdef W3_TIMINGS
        CALL print_my_time("After FFLEV and DTL0")
#endif
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 11b')
        !
        ! 3.5 Update maps and derivatives.
        !
        IF ( flmap ) THEN
          IF ( gtype .NE. smctype ) THEN
#ifdef W3_PR1
            CALL w3map1 ( mapsta )
#endif
#ifdef W3_PR2
            CALL w3map2
#endif
#ifdef W3_PR3
            CALL w3map3
#endif
            CALL w3utrn ( trnx, trny )
#ifdef W3_PR3
            CALL w3mapt
#endif
          END IF  !! GTYPE
 
          !! Hides call to W3NMIN, which currently only serves to warn when
          !! one or more procs have zero active seapoints.
#ifdef W3_DEBUGRUN
          CALL w3nmin ( mapsta, flag0 )
          IF ( flag0 .AND. iaproc.EQ.naperr ) WRITE (ndse,1030) imod
#endif
          flmap  = .false.
        END IF
        !
        !
        IF ( flddir ) THEN
          IF (gtype .EQ. smctype) THEN
            ix = 1
#ifdef W3_SMC
            !!Li  Use new sub for DDDX and DDDY assignment.
            CALL smcdhxy
#endif
          ELSE IF (gtype .EQ. ungtype) THEN
            CALL ug_gradients(dw, dddx, dddy)
          ELSE
            CALL w3dzxy(dw(1:ubound(dw,1)),'m',dddx,dddy)
          END IF
          flddir = .false.
        END IF
 
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 12')
        !
        ! Calculate PHASE SPEED GRADIENT.
        dcdx = 0.
        dcdy = 0.
#ifdef W3_REFRX
        cik  = 0.
        !
        IF (gtype .NE. ungtype) THEN
          DO ik=0,nk+1
            cik = sig(ik) / wn(ik,1:nsea)
            CALL w3dzxy(cik,'m/s',dcdx(ik,:,:),dcdy(ik,:,:))
          END DO
        ELSE
          WRITE (ndse,1040)
          CALL extcde(2)
          ! CALL UG_GRADIENTS(CMN, DCDX, DCDY) !/ Stefan, to be confirmed!
        END IF
#endif
        !
        !
        fliwnd = .false.
        flfrst = .false.
        !
#ifdef W3_PDLIB
#ifdef W3_DEBUGSRC
        WRITE(740+iaproc,*) 'ITIME=', itime, ' IT=', it
        CALL all_va_integral_print(imod, "VA before W3SRCE_IMP_PRE", 1)
        CALL all_field_integral_print(vstot, "VSTOT before W3SRCE_IMP_PRE")
        CALL all_field_integral_print(vdtot, "VDTOT before W3SRCE_IMP_PRE")
        IF (debug_node .le. nseal) THEN
          WRITE(740+iaproc,*) '     Values for DEBUG_NODE=', debug_node
          WRITE(740+iaproc,*) '     sum(VA)=', sum(va(:,debug_node))
          WRITE(740+iaproc,*) '     sum(VSTOT)=', sum(vstot(:,debug_node))
          WRITE(740+iaproc,*) '     sum(VDTOT)=', sum(vdtot(:,debug_node))
        END IF
#endif
        IF (it .eq. 0) THEN
          dtgpre = 1.
        ELSE
          dtgpre = dtg
        END IF
#endif
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 13')
        !
#ifdef W3_PDLIB
        IF (lpdlib .and. flsou .and. fssource) THEN
#endif
 
#ifdef W3_OMP0
          !$OMP PARALLEL DO PRIVATE (JSEA,ISEA,IX,IY) SCHEDULE (DYNAMIC,1)
#endif
 
#ifdef W3_PDLIB
          d50=0.0002
          reflec(:)=0.
          refled(:)=0
          psic=0.
#endif
 
#ifdef W3_PDLIB
          IF (lsloc) THEN
            b_jac     = 0.
            aspar_jac = 0.
          ELSE
            vstot = 0.
            vdtot = 0.
          ENDIF
#endif
 
 
#ifdef W3_PDLIB
 
          DO jsea = 1, np
 
            CALL init_get_isea(isea, jsea)
 
            ix     = mapsf(isea,1)
            iy     = mapsf(isea,2)
            dela=1.
            delx=1.
            dely=1.
 
#ifdef W3_REF1
            IF (gtype.EQ.rlgtype) THEN
              delx=sx*clats(isea)/facx
              dely=sy/facx
              dela=delx*dely
            END IF
            IF (gtype.EQ.clgtype) THEN
              ! Maybe what follows works also for RLGTYPE ... to be verified
              delx=hpfac(iy,ix)/ facx
              dely=hqfac(iy,ix)/ facx
              dela=delx*dely
            END IF
            reflec=reflc(:,isea)
            reflec(4)=berg(isea)*reflec(4)
            refled=refld(:,isea)
#endif
 
#ifdef W3_BT4
            d50=sed_d50(isea)
            psic=sed_psic(isea)
#endif
            !
#ifdef W3_DEBUGSRC
            IF (ix .eq. debug_node) THEN
              WRITE(740+iaproc,*) 'NODE_SRCE_IMP_PRE : IX=', ix, ' JSEA=', jsea
            END IF
            WRITE(740+iaproc,*) 'IT/IX/IY/IMOD=', it, ix, iy, imod
            WRITE(740+iaproc,*) 'ISEA/JSEA=', isea, jsea
            WRITE(740+iaproc,*) 'Before sum(VA)=', sum(va(:,jsea))
            FLUSH(740+iaproc)
#endif
            CALL w3srce(srce_imp_pre, it, isea, jsea, ix, iy, imod, &
                 vaold(:,jsea), va(:,jsea),                         &
                 vsiodummy, vdiodummy, shavetot(jsea),              &
                 alpha(1:nk,jsea), wn(1:nk,isea),                   &
                 cg(1:nk,isea), clats(isea), dw(isea), u10(isea),   &
                 u10d(isea),                                        &
#ifdef W3_FLX5
                 taua(isea), tauadir(isea),                         &
#endif
                 as(isea), ust(isea),                               &
                 ustdir(isea), cx(isea), cy(isea),                  &
                 ice(isea), iceh(isea), icef(isea),                 &
                 icedmax(isea),                                     &
                 reflec, refled, delx, dely, dela,                  &
                 trnx(iy,ix), trny(iy,ix), berg(isea),              &
                 fpis(isea), dtdyn(jsea),                           &
                 fcut(jsea), dtgpre, tauwx(jsea), tauwy(jsea),      &
                 tauox(jsea), tauoy(jsea), tauwix(jsea),            &
                 tauwiy(jsea), tauwnx(jsea),                        &
                 tauwny(jsea),  phiaw(jsea), charn(jsea),           &
                 tws(jsea), phioc(jsea), tmp1, d50, psic, tmp2,     &
                 phibbl(jsea), tmp3, tmp4, phice(jsea),             &
                 tauocx(jsea), tauocy(jsea), wnmean(jsea),          &
                 rhoair(isea), asf(isea))
            IF (.not. lsloc) THEN
              vstot(:,jsea) = vsiodummy
              vdtot(:,jsea) = vdiodummy
            ENDIF
#ifdef W3_DEBUGSRC
            WRITE(740+iaproc,*) 'After sum(VA)=', sum(va(:,jsea))
            WRITE(740+iaproc,*) '   sum(VSTOT)=', sum(vstot(:,jsea))
            WRITE(740+iaproc,*) '   sum(VDTOT)=', sum(vdtot(:,jsea))
            WRITE(740+iaproc,*) '     SHAVETOT=', shavetot(jsea)
            FLUSH(740+iaproc)
#endif
          END DO ! JSEA
        END IF ! PDLIB
#endif
 
 
#ifdef W3_PDLIB
#ifdef W3_DEBUGSRC
        WRITE(740+iaproc,*) 'ITIME=', itime, ' IT=', it
        CALL all_va_integral_print(imod, "VA after W3SRCE_IMP_PRE", 1)
        CALL all_field_integral_print(vstot, "VSTOT after W3SRCE_IMP_PRE")
        CALL all_field_integral_print(vdtot, "VDTOT after W3SRCE_IMP_PRE")
        IF (debug_node .le. nseal) THEN
          WRITE(740+iaproc,*) '     Values for DEBUG_NODE=', debug_node
          WRITE(740+iaproc,*) '     sum(VA)=', sum(va(:,debug_node))
          WRITE(740+iaproc,*) '     sum(VSTOT)=', sum(vstot(:,debug_node))
          WRITE(740+iaproc,*) '     sum(VDTOT)=', sum(vdtot(:,debug_node))
        END IF
#endif
#endif
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 14')
 
        IF ( flzero ) THEN
#ifdef W3_T
          WRITE (ndst,9022)
#endif
          GOTO 400
        END IF
        IF ( it.EQ.0 ) THEN
          dtg = 1.
          !            DTG = 60.
          GOTO 370
        END IF
        IF ( fldry .OR. iaproc.GT.naproc ) THEN
#ifdef W3_T
          WRITE (ndst,9023)
#endif
          GOTO 380
        END IF
        !
        ! Estimation of the local maximum CFL for XY propagation
        !
#ifdef W3_T
        WRITE(ndse,*) 'Computing CFLs .... ',nseal
#endif
        IF ( flogrd(9,3).AND. ugdtupdate ) THEN
          IF (fstotalimp .eqv. .false.) THEN
            nkcfl=nk
#ifdef W3_T
            nkcfl=1
#endif
            !
#ifdef W3_OMPG
            !$OMP PARALLEL DO PRIVATE (JSEA,ISEA) SCHEDULE (DYNAMIC,1)
#endif
            !
            DO jsea=1, nseal
              CALL init_get_isea(isea, jsea)
#ifdef W3_PR3
              IF (gtype .EQ. ungtype) THEN
                IF ( flogrd(9,3) ) THEN
#endif
#ifdef W3_T
                  IF (mod(isea,100).EQ.0) WRITE(ndse,*) 'COMPUTING CFL FOR NODE:',isea
#endif
#ifdef W3_PDLIB
                  IF (.NOT. lpdlib) THEN
#endif
#ifdef W3_PR3
                    CALL w3cflug ( isea, nkcfl, facx, facx, dtg, mapfs, cflxymax(jsea), &
                         vgx, vgy )
#endif
#ifdef W3_PDLIB
                  ENDIF
#endif
#ifdef W3_PR3
                END IF
              ELSE
                CALL w3cflxy ( isea, dtg, mapsta, mapfs, cflxymax(jsea), vgx, vgy )
              END IF
#endif
            END DO
            !
#ifdef W3_OMPG
            !$OMP END PARALLEL DO
#endif
            !
          END IF
        END IF
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 15')
        !
 
        !
#ifdef W3_T
        IF (gtype .EQ. ungtype) THEN
          IF ( flogrd(9,3) ) THEN
            dtcfl1(:)=1.
            DO jsea=1,nseal
              indsort(jsea)=float(jsea)
              dtcfl1(jsea)=dtg/cflxymax(jsea)
            END DO
            CALL ssort1 (dtcfl1, indsort, nseal, 2)
            IF ( iaproc .EQ. naperr ) WRITE(ndse,*) 'Nodes requesting smallest timesteps:'
            IF ( iaproc .EQ. naperr ) WRITE(ndse,'(A,10I10)')   'Nodes      ',nint(indsort(1:10))
            IF ( iaproc .EQ. naperr ) WRITE(ndse,'(A,10F10.2)') 'time steps ',dtcfl1(1:10)
            DO jsea = 1, min(nseal,200)
              isea   = nint(indsort(jsea))            ! will not work with MPI
              ix     = mapsf(isea,1)
              IF (jsea.EQ.1) then
                WRITE(995,*) '       IP  dtmax_exp(ip)        x-coord        y-coord        z-coord'
              end if
              WRITE(995,'(I10,F10.2,3F10.4)') ix,  dtcfl1(jsea), xgrd(1,ix), ygrd(2,ix), zb(ix)
            END DO ! JSEA
            CLOSE(995)
          END IF
        END IF
#endif
 
        !
        ! 3.6 Perform Propagation = = = = = = = = = = = = = = = = = = = = = = =
        ! 3.6.1 Preparations
        !
#ifdef W3_SEC1
        dtgtemp=dtg
        dtg=dtg/nitersec1
        DO isec1=1,nitersec1
#endif
          ntloc  = 1 + int( dtg/dtcfli - 0.001 )
#ifdef W3_SEC1
          IF ( iaproc .EQ. napout ) then
            WRITE(ndse,'(A,I4,A,I4)') '   SUBSECOND STEP:',isec1,' out of ',nitersec1
          end if
#endif
          !
          facth  = dtg / (dth*real(ntloc))
          fack   = dtg / real(ntloc)
 
          ttest(1) = time(1)
          ttest(2) = 0
          dttest = dsec21(ttest,time)
          itloch = ( ntloc + 1 - mod(nint(dttest/dtg),2) ) / 2
          !
          ! 3.6.2 Intra-spectral part 1
          !
#ifdef W3_DEBUGCOH
          CALL all_va_integral_print(imod, "Before intraspectral part 1", 1)
#endif
#ifdef W3_TIMINGS
          CALL print_my_time("Before intraspectral")
#endif
          IF ( flcth .OR. flck ) THEN
            DO itloc=1, itloch
              !
#ifdef W3_OMPG
              !$OMP PARALLEL PRIVATE (JSEA,ISEA,IX,IY,DEPTH,IXrel)
              !$OMP DO SCHEDULE (DYNAMIC,1)
#endif
              !
              DO jsea=1, nseal
                CALL init_get_isea(isea, jsea)
                ix     = mapsf(isea,1)
                iy     = mapsf(isea,2)
 
 
                IF ( gtype .EQ. ungtype ) THEN
                  IF (lpdlib) THEN
#ifdef W3_PDLIB
                    IF (iobp_loc(jsea) .NE. 1) cycle
#endif
                  ELSE
                    IF (iobp(isea) .NE. 1) cycle
                  ENDIF
                ENDIF
 
                IF ( mapsta(iy,ix) .EQ. 1 ) THEN
                  depth  = max( dmin , dw(isea) )
                  IF (lpdlib) THEN
                    ixrel = jsea
                  ELSE
                    ixrel = ix
                  END IF
                  !
                  IF( gtype .EQ. smctype ) THEN
                    j = 1
#ifdef W3_SMC
                    !!Li    Refraction and GCT in theta direction is done by rotation.
                    CALL w3krtn ( isea, facth, fack, cthg0s(isea),       &
                         cg(:,isea), wn(:,isea), depth,                  &
                         dhdx(isea), dhdy(isea), dhlmt(:,isea),          &
                         cx(isea), cy(isea), dcxdx(iy,ix),               &
                         dcxdy(iy,ix), dcydx(iy,ix), dcydy(iy,ix),       &
                         dcdx(:,iy,ix), dcdy(:,iy,ix), va(:,jsea) )
#endif
                    !
                  ELSE
                    j = 1
                    !
#ifdef W3_PR1
                    CALL w3ktp1 ( isea, facth, fack, cthg0s(isea),       &
                         cg(:,isea), wn(:,isea), depth,                  &
                         dddx(iy,ixrel), dddy(iy,ixrel), cx(isea),       &
                         cy(isea), dcxdx(iy,ixrel), dcxdy(iy,ixrel),     &
                         dcydx(iy,ixrel), dcydy(iy,ixrel),               &
                         dcdx(:,iy,ixrel), dcdy(:,iy,ixrel), va(:,jsea))
#endif
#ifdef W3_PR2
                    CALL w3ktp2 ( isea, facth, fack, cthg0s(isea),       &
                         cg(:,isea), wn(:,isea), depth,                  &
                         dddx(iy,ixrel), dddy(iy,ixrel), cx(isea),       &
                         cy(isea), dcxdx(iy,ixrel), dcxdy(iy,ixrel),     &
                         dcydx(iy,ixrel), dcydy(iy,ixrel),               &
                         dcdx(:,iy,ixrel), dcdy(:,iy,ixrel), va(:,jsea))
#endif
#ifdef W3_PR3
                    CALL w3ktp3 ( isea, facth, fack, cthg0s(isea),       &
                         cg(:,isea), wn(:,isea), depth,                  &
                         dddx(iy,ixrel), dddy(iy,ixrel), cx(isea),       &
                         cy(isea), dcxdx(iy,ixrel), dcxdy(iy,ixrel),     &
                         dcydx(iy,ixrel), dcydy(iy,ixrel),               &
                         dcdx(:,iy,ixrel), dcdy(:,iy,ixrel), va(:,jsea), &
                         cflthmax(jsea), cflkmax(jsea) )
#endif
                    !
                  END IF  !!  GTYPE
                  !
                END IF
              END DO
              !
#ifdef W3_OMPG
              !$OMP END DO
              !$OMP END PARALLEL
#endif
              !
            END DO
          END IF
 
          call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 16')
 
#ifdef W3_DEBUGCOH
          CALL all_va_integral_print(imod, "Before spatial advection", 1)
#endif
#ifdef W3_TIMINGS
          CALL print_my_time("Before spatial advection")
#endif
          !
          ! 3.6.3 Longitude-latitude
          !       (time step correction in routine)
          !
          IF (gtype .EQ. ungtype) THEN
            IF (flagll) THEN
              facx   =  1./(dera * radius)
            ELSE
              facx   =  1.
            END IF
          END IF
 
          IF (lpdlib) THEN
            !
#ifdef W3_PDLIB
            IF (flcx .or. flcy) THEN
              IF (.NOT. fstotalimp .AND. .NOT. fstotalexp) THEN
                DO ispec=1,nspec
                  CALL pdlib_w3xypug ( ispec, facx, facx, dtg, vgx, vgy, ugdtupdate )
                END DO
              END IF
            END IF
#endif
            !
#ifdef W3_PDLIB
            IF (fstotalimp .and. (it .ne. 0)) THEN
#endif
#ifdef W3_DEBUGCOH
              CALL all_va_integral_print(imod, "Before Block implicit", 1)
#endif
#ifdef W3_PDLIB
              CALL pdlib_w3xypug_block_implicit(imod, facx, facx, dtg, vgx, vgy, ugdtupdate )
#endif
#ifdef W3_PDLIB
            ELSE IF(fstotalexp .and. (it .ne. 0)) THEN
#endif
#ifdef W3_PDLIB
              CALL pdlib_w3xypug_block_explicit(imod, facx, facx, dtg, vgx, vgy, ugdtupdate )
#endif
#ifdef W3_PDLIB
            ENDIF
#endif
          ELSE
            IF (flcx .or. flcy) THEN
              !
#ifdef W3_MPI
              IF ( nrqsg1 .GT. 0 ) THEN
                CALL mpi_startall (nrqsg1, irqsg1(1,1), ierr_mpi)
                CALL mpi_startall (nrqsg1, irqsg1(1,2), ierr_mpi)
              END IF
#endif
              !
              !
              ! Initialize FIELD variable
              field = 0.
              !
              DO ispec=1, nspec
                IF ( iappro(ispec) .EQ. iaproc ) THEN
                  !
                  IF( gtype .EQ. smctype ) THEN
                    ix = 1
#ifdef W3_SMC
                    !!Li   Use SMC sub to gether field
                    CALL w3gathsmc ( ispec, field )
#endif
                  ELSE IF (.NOT.lpdlib ) THEN
                    CALL w3gath ( ispec, field )
                  END IF   !! GTYPE
                  !
                  IF (gtype .EQ. smctype) THEN
                    ix = 1
#ifdef W3_SMC
                    !!Li   Propagation on SMC grid uses UNO2 scheme.
                    CALL w3psmc ( ispec, dtg, field )
#endif
                    !
                  ELSE IF (gtype .EQ. ungtype) THEN
                    ix = 1
#ifdef W3_MPI
                    IF (.NOT. lpdlib) THEN
#endif
#ifdef W3_PR1
                      CALL w3xypug ( ispec, facx, facx, dtg, field, vgx, vgy, ugdtupdate )
#endif
#ifdef W3_PR2
                      CALL w3xypug ( ispec, facx, facx, dtg, field, vgx, vgy, ugdtupdate )
#endif
#ifdef W3_PR3
                      CALL w3xypug ( ispec, facx, facx, dtg, field, vgx, vgy, ugdtupdate )
#endif
#ifdef W3_MPI
                    END IF
#endif
                    !
                  ELSE
                    ix = 1
#ifdef W3_PR1
                    CALL w3xyp1 ( ispec, dtg, mapsta, field, vgx, vgy )
#endif
#ifdef W3_PR2
                    CALL w3xyp2 ( ispec, dtg, mapsta, mapfs, field, vgx, vgy )
#endif
#ifdef W3_PR3
                    CALL w3xyp3 ( ispec, dtg, mapsta, mapfs, field, vgx, vgy )
#endif
                    !
                  END IF   !! GTYPE
                  !
                  IF( gtype .EQ. smctype ) THEN
                    ix = 1
#ifdef W3_SMC
                    !!Li   Use SMC sub to scatter field
                    CALL w3scatsmc ( ispec, mapsta, field )
#endif
                  ELSE IF (.NOT.lpdlib ) THEN
                    CALL w3scat ( ispec, mapsta, field )
                  END IF   !! GTYPE
 
                END IF
              END DO
              !
#ifdef W3_MPI
              IF ( nrqsg1 .GT. 0 ) THEN
                ALLOCATE ( statco(mpi_status_size,nrqsg1) )
                CALL mpi_waitall (nrqsg1, irqsg1(1,1), statco, ierr_mpi)
                CALL mpi_waitall (nrqsg1, irqsg1(1,2), statco, ierr_mpi)
                DEALLOCATE ( statco )
              END IF
#endif
              call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 17')
              !
              !Li   Initialise IK IX IY in case ARC option is not used to avoid warnings.
              ik=1
              ix=1
              iy=1
#ifdef W3_SMC
              !Li    Find source boundary spectra and assign to SPCBAC
              IF( arctc ) THEN
 
                DO ik = 1, nbac
                  IF( ik .LE. (nbac-nbgl) ) THEN
                    iy = iclbac(ik)
                  ELSE
                    iy = nglo + ik
                  ENDIF
 
                  !Li    Work out root PE (ISPEC) and JSEA numbers for IY
#ifdef W3_DIST
                  ispec = mod( iy-1, naproc )
                  jsea = 1 + (iy - ispec - 1)/naproc
#endif
#ifdef W3_SHRD
                  ispec = 0
                  jsea = iy
#endif
#endif
                  ! W3_SMC ...
                  !
#ifdef W3_SMC
                  !!Li   Assign boundary cell spectra.
                  IF( iaproc .EQ. ispec+1 ) THEN
                    spcbac(:,ik)=va(:,jsea)
                  ENDIF
#endif
                  !
#ifdef W3_SMC
                  !!Li   Broadcast local SPCBAC(:,IK) to all other PEs.
#ifdef W3_MPI
                  CALL mpi_bcast(spcbac(1,ik),nspec,mpi_real,ispec,mpi_comm_wave,ierr_mpi)
                  CALL mpi_barrier (mpi_comm_wave,ierr_mpi)
#endif
#endif
                  !
#ifdef W3_SMC
                END DO   !! Loop IK ends.
#endif
                !
#ifdef W3_SMC
                !!Li    Update Arctic boundary cell spectra if within local range
                ALLOCATE ( bacspec(nspec) )
                DO ik = 1, nbac
                  IF( ik .LE. (nbac-nbgl) ) THEN
                    ix = nglo + ik
                    bacangl = angarc(ik)
                  ELSE
                    ix = iclbac(ik)
                    bacangl = - angarc(ik)
                  ENDIF
 
                  !!Li    Work out boundary PE (ISPEC) and JSEA numbers for IX
#ifdef W3_DIST
                  ispec = mod( ix-1, naproc )
                  jsea = 1 + (ix - ispec - 1)/naproc
#endif
#ifdef W3_SHRD
                  ispec = 0
                  jsea = ix
#endif
#endif
                  !
#ifdef W3_SMC
                  IF( iaproc .EQ. ispec+1 ) THEN
                    bacspec = spcbac(:,ik)
 
                    CALL w3acturn( nth, nk, bacangl, bacspec )
 
                    va(:,jsea) = bacspec
                    !!Li              WRITE(NDSE,*) "IAPROC, IX, JSEAx, IK=", IAPROC, IX, JSEA, IK
                  ENDIF
 
                END DO  !! Loop IK ends.
                DEALLOCATE ( bacspec )
 
              ENDIF  !! ARCTC
#endif
              !
              ! End of test FLCX.OR.FLCY
            END IF
            !
          END IF
 
#ifdef W3_DEBUGCOH
          CALL all_va_integral_print(imod, "After spatial advection", 1)
#endif
#ifdef W3_TIMINGS
          CALL print_my_time("After spatial advection")
#endif
          !
          ! 3.6.4 Intra-spectral part 2
          !
          IF ( flcth .OR. flck ) THEN
            DO itloc=itloch+1, ntloc
              !
#ifdef W3_OMPG
              !$OMP PARALLEL PRIVATE (JSEA,ISEA,IX,IY,DEPTH,IXrel)
              !$OMP DO SCHEDULE (DYNAMIC,1)
#endif
              !
              DO jsea = 1, nseal
 
                CALL init_get_isea(isea, jsea)
                ix     = mapsf(isea,1)
                iy     = mapsf(isea,2)
                depth  = max( dmin , dw(isea) )
 
                IF ( gtype .EQ. ungtype ) THEN
                  IF (lpdlib) THEN
#ifdef W3_PDLIB
                    IF (iobp_loc(jsea) .NE. 1) cycle
#endif
                  ELSE
                    IF (iobp(isea) .NE. 1) cycle
                  ENDIF
                ENDIF
 
                IF ( mapsta(iy,ix) .EQ. 1 ) THEN
                  IF (lpdlib) THEN
                    ixrel = jsea
                  ELSE
                    ixrel = ix
                  END IF
                  !
                  IF( gtype .EQ. smctype ) THEN
                    j = 1
#ifdef W3_SMC
                    !!Li    Refraction and GCT in theta direction is done by rotation.
                    CALL w3krtn ( isea, facth, fack, cthg0s(isea),       &
                         cg(:,isea), wn(:,isea), depth,                  &
                         dhdx(isea), dhdy(isea), dhlmt(:,isea),          &
                         cx(isea), cy(isea), dcxdx(iy,ix),               &
                         dcxdy(iy,ix), dcydx(iy,ix), dcydy(iy,ix),       &
                         dcdx(:,iy,ix), dcdy(:,iy,ix), va(:,jsea) )
#endif
                    !
                  ELSE
                    j = 1
#ifdef W3_PR1
                    CALL w3ktp1 ( isea, facth, fack, cthg0s(isea),       &
                         cg(:,isea), wn(:,isea), depth,                  &
                         dddx(iy,ixrel), dddy(iy,ixrel), cx(isea),       &
                         cy(isea), dcxdx(iy,ixrel), dcxdy(iy,ixrel),     &
                         dcydx(iy,ixrel), dcydy(iy,ixrel),               &
                         dcdx(:,iy,ixrel), dcdy(:,iy,ixrel), va(:,jsea))
#endif
#ifdef W3_PR2
                    CALL w3ktp2 ( isea, facth, fack, cthg0s(isea),       &
                         cg(:,isea), wn(:,isea), depth,                  &
                         dddx(iy,ixrel), dddy(iy,ixrel), cx(isea),       &
                         cy(isea), dcxdx(iy,ixrel), dcxdy(iy,ixrel),     &
                         dcydx(iy,ixrel), dcydy(iy,ixrel),               &
                         dcdx(:,iy,ixrel), dcdy(:,iy,ixrel), va(:,jsea))
#endif
#ifdef W3_PR3
                    CALL w3ktp3 ( isea, facth, fack, cthg0s(isea),       &
                         cg(:,isea), wn(:,isea), depth,                  &
                         dddx(iy,ixrel), dddy(iy,ixrel), cx(isea),       &
                         cy(isea), dcxdx(iy,ixrel), dcxdy(iy,ixrel),     &
                         dcydx(iy,ixrel), dcydy(iy,ixrel),               &
                         dcdx(:,iy,ixrel), dcdy(:,iy,ixrel), va(:,jsea), &
                         cflthmax(jsea), cflkmax(jsea) )
#endif
                    !
                  END IF  !! GTYPE
                  !
                END IF
              END DO
              !
#ifdef W3_OMPG
              !$OMP END DO
              !$OMP END PARALLEL
#endif
              !
            END DO
          END IF
#ifdef W3_DEBUGCOH
          CALL all_va_integral_print(imod, "After intraspectral adv.", 1)
#endif
#ifdef W3_TIMINGS
          CALL print_my_time("fter intraspectral adv.")
#endif
 
          !
          ugdtupdate = .false.
          !
          ! 3.6 End propapgation  = = = = = = = = = = = = = = = = = = = = = = = =
 
          ! 3.7 Calculate and integrate source terms.
          !
370       CONTINUE
          IF ( flsou ) THEN
            !
            d50=0.0002
            reflec(:)=0.
            refled(:)=0
            psic=0.
#ifdef W3_PDLIB
#ifdef W3_DEBUGSRC
            WRITE(740+iaproc,*) 'ITIME=', itime, ' IT=', it
            CALL all_vaold_integral_print("VAOLD before W3SRCE_IMP_POST", 1)
            CALL all_va_integral_print(imod, "VA before W3SRCE_IMP_POST", 1)
            IF (debug_node .le. nseal) THEN
              WRITE(740+iaproc,*) '     Values for DEBUG_NODE=', debug_node
              WRITE(740+iaproc,*) '     sum(VA)=', sum(va(:,debug_node))
              WRITE(740+iaproc,*) '     sum(VAOLD)=', sum(vaold(:,debug_node))
              WRITE(740+iaproc,*) '     sum(VSTOT)=', sum(vstot(:,debug_node))
              WRITE(740+iaproc,*) '     sum(VDTOT)=', sum(vdtot(:,debug_node))
            END IF
#endif
#endif
            !
#ifdef W3_OMPG
            !$OMP PARALLEL PRIVATE (JSEA,ISEA,IX,IY,DELA,DELX,DELY,        &
            !$OMP&                  REFLEC,REFLED,D50,PSIC,TMP1,TMP2,TMP3,TMP4)
            !$OMP DO SCHEDULE (DYNAMIC,1)
#endif
 
            !
            DO jsea=1, nseal
              CALL init_get_isea(isea, jsea)
              ix     = mapsf(isea,1)
              iy     = mapsf(isea,2)
              dela=1.
              delx=1.
              dely=1.
#ifdef W3_REF1
              IF (gtype.EQ.rlgtype) THEN
                delx=sx*clats(isea)/facx
                dely=sy/facx
                dela=delx*dely
              END IF
              IF (gtype.EQ.clgtype) THEN
                ! Maybe what follows works also for RLGTYPE ... to be verified
                delx=hpfac(iy,ix)/ facx
                dely=hqfac(iy,ix)/ facx
                dela=delx*dely
              END IF
#endif
              !
#ifdef W3_REF1
              reflec=reflc(:,isea)
              reflec(4)=berg(isea)*reflec(4)
              refled=refld(:,isea)
#endif
#ifdef W3_BT4
              d50=sed_d50(isea)
              psic=sed_psic(isea)
#endif
 
 
              IF ( mapsta(iy,ix) .EQ. 1 .AND. flagst(isea)) THEN
                tmp1   = whitecap(jsea,1:4)
                tmp2   = bedforms(jsea,1:3)
                tmp3   = taubbl(jsea,1:2)
                tmp4   = tauice(jsea,1:2)
#ifdef W3_PDLIB
                IF (fssource) THEN
                  CALL w3srce(srce_imp_post,it,isea,jsea,ix,iy,imod,     &
                       vaold(:,jsea), va(:,jsea),                        &
                       vsiodummy,vdiodummy,shavetot(jsea),               &
                       alpha(1:nk,jsea), wn(1:nk,isea),                  &
                       cg(1:nk,isea), clats(isea), dw(isea), u10(isea),  &
                       u10d(isea),                                       &
#ifdef W3_FLX5
                       taua(isea), tauadir(isea),                        &
#endif
                       as(isea), ust(isea),                              &
                       ustdir(isea), cx(isea), cy(isea),                 &
                       ice(isea), iceh(isea), icef(isea),                &
                       icedmax(isea),                                    &
                       reflec, refled, delx, dely, dela,                 &
                       trnx(iy,ix), trny(iy,ix), berg(isea),             &
                       fpis(isea), dtdyn(jsea),                          &
                       fcut(jsea), dtg, tauwx(jsea), tauwy(jsea),        &
                       tauox(jsea), tauoy(jsea), tauwix(jsea),           &
                       tauwiy(jsea), tauwnx(jsea),                       &
                       tauwny(jsea),  phiaw(jsea), charn(jsea),          &
                       tws(jsea),phioc(jsea), tmp1, d50, psic, tmp2,     &
                       phibbl(jsea), tmp3, tmp4, phice(jsea),            &
                       tauocx(jsea), tauocy(jsea), wnmean(jsea),         &
                       rhoair(isea), asf(isea))
                ELSE
#endif
                  CALL w3srce(srce_direct, it, isea, jsea, ix, iy, imod, &
                       vaolddummy, va(:,jsea),                           &
                       vsiodummy, vdiodummy, shavetotiodummy,            &
                       alpha(1:nk,jsea), wn(1:nk,isea),                  &
                       cg(1:nk,isea), clats(isea), dw(isea), u10(isea),  &
                       u10d(isea),                                       &
#ifdef W3_FLX5
                       taua(isea), tauadir(isea),                        &
#endif
                       as(isea), ust(isea),                              &
                       ustdir(isea), cx(isea), cy(isea),                 &
                       ice(isea), iceh(isea), icef(isea),                &
                       icedmax(isea),                                    &
                       reflec, refled, delx, dely, dela,                 &
                       trnx(iy,ix), trny(iy,ix), berg(isea),             &
                       fpis(isea), dtdyn(jsea),                          &
                       fcut(jsea), dtg, tauwx(jsea), tauwy(jsea),        &
                       tauox(jsea), tauoy(jsea), tauwix(jsea),           &
                       tauwiy(jsea), tauwnx(jsea),                       &
                       tauwny(jsea),  phiaw(jsea), charn(jsea),          &
                       tws(jsea), phioc(jsea), tmp1, d50, psic,tmp2,     &
                       phibbl(jsea), tmp3, tmp4 , phice(jsea),           &
                       tauocx(jsea), tauocy(jsea), wnmean(jsea),         &
                       rhoair(isea), asf(isea))
#ifdef W3_PDLIB
                END IF
#endif
                whitecap(jsea,1:4) = tmp1
                bedforms(jsea,1:3) = tmp2
                taubbl(jsea,1:2) = tmp3
                tauice(jsea,1:2) = tmp4
              ELSE
                ust(isea) = undef
                ustdir(isea) = undef
                dtdyn(jsea) = undef
                fcut(jsea) = undef
                !                    VA(:,JSEA)  = 0.
              END IF
            END DO
 
            !
#ifdef W3_OMPG
            !$OMP END DO
            !$OMP END PARALLEL
#endif
            !
#ifdef W3_PDLIB
#ifdef W3_DEBUGSRC
            WRITE(740+iaproc,*) 'ITIME=', itime, ' IT=', it
            CALL all_vaold_integral_print("VAOLD after W3SRCE_IMP_PRE_POST", 1)
            CALL all_va_integral_print(imod, "VA after W3SRCE_IMP_PRE_POST", 1)
            IF (debug_node .le. nseal) THEN
              WRITE(740+iaproc,*) '     Values for DEBUG_NODE=', debug_node
              WRITE(740+iaproc,*) '     sum(VA)=', sum(va(:,debug_node))
              WRITE(740+iaproc,*) '     min/max(VA)=', minval(va(:,debug_node)), maxval(va(:,debug_node))
            END IF
#endif
#endif
          END IF
#ifdef W3_DEBUGCOH
          CALL all_va_integral_print(imod, "After source terms", 1)
#endif
#ifdef W3_TIMINGS
          CALL print_my_time("After source terms")
#endif
          !
          ! End of interations for DTMAX < 1s
          !
#ifdef W3_SEC1
          IF (it.EQ.0) EXIT
        END DO
        IF (it.GT.0) dtg=dtgtemp
#endif
 
 
 
 
        !
        !
        ! 3.8 Update global time step.
        !     (Branch point FLDRY, IT=0)
        !
380     CONTINUE
        !
        IF (it.NE.nt) THEN
          dttst  = dsec21( time , tcalc )
          dtg    = dttst / real(nt-it)
        END IF
        !
        IF ( flact .AND. it.NE.nt .AND. iaproc.EQ.naplog ) THEN
          CALL stme21 ( time , idtime )
          IF ( idlast .NE. time(1) ) THEN
            WRITE (ndso,900) itime, ipass, idtime(01:19), idact, outid
            idlast = time(1)
          ELSE
            WRITE (ndso,901) itime, ipass, idtime(12:19), idact, outid
          END IF
          flact  = .false.
          idact  = '         '
        END IF
        !
#ifdef W3_DEBUGCOH
        CALL all_va_integral_print(imod, "end of time loop", 1)
#endif
#ifdef W3_TIMINGS
        CALL print_my_time("end of time loop")
#endif
        !
        !
      END DO
 
#ifdef W3_TIMINGS
      CALL print_my_time("W3WAVE, step 6.21.1")
#endif
      !
#ifdef W3_T
      WRITE (ndst,9030)
#endif
      call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE END TIME LOOP')
      !
      !     End of loop over time steps
      ! ==================================================================== /
      !
400   CONTINUE
      !
      ! 4.  Perform output to file if requested ---------------------------- /
      ! 4.a Check if time is output time
      !     Delay if data assimilation time.
      !
      !
      IF ( tofrst(1)  .EQ. -1 ) THEN
        dttst  = 1.
      ELSE
        dttst   = dsec21( time, tofrst )
      END IF
      !
      IF ( tdn(1)  .EQ. -1 ) THEN
        dttst1 = 1.
      ELSE
        dttst1  = dsec21( time, tdn )
      END IF
      !
      dttst2 = dsec21( time, tend )
      flag_o = .NOT.skip_o .OR. ( skip_o .AND. dttst2.NE.0. )
      !
#ifdef W3_T
      WRITE (ndst,9040) tofrst, tdn, dttst, dttst1, flag_o
#endif
      !
      IF ( dttst.LE.0. .AND. dttst1.NE.0. .AND. flag_o ) THEN
        !
#ifdef W3_T
        WRITE (ndst,9041)
#endif
        !
        ! 4.b Processing and MPP preparations
        !
        IF ( flout(1) ) THEN
          floutg = dsec21(time,tonext(:,1)).EQ.0.
        ELSE
          floutg = .false.
        END IF
        !
        IF ( flout(7) ) THEN
          floutg2 = dsec21(time,tonext(:,7)).EQ.0.
        ELSE
          floutg2 = .false.
        END IF
        !
        flpart = .false.
        IF ( flout(1) .AND. flpfld ) flpart = flpart .OR. dsec21(time,tonext(:,1)).EQ.0.
        IF ( flout(6) ) flpart = flpart .OR. dsec21(time,tonext(:,6)).EQ.0.
        !
#ifdef W3_T
        WRITE (ndst,9042) local, flpart, floutg
#endif
        !
        IF ( local .AND. flpart ) CALL w3cprt ( imod )
        IF ( local .AND. (floutg .OR. floutg2) ) then
          CALL w3outg ( va, flpfld, floutg, floutg2 )
        end if
        !
#ifdef W3_MPI
        flgmpi = .false.
        nrqmax = 0
#endif
        !
#ifdef W3_MPI
        IF ( (floutg) .OR. (floutg2 .AND. sbsed) ) THEN
          IF (.NOT. lpdlib) THEN
            IF (nrqgo.NE.0 ) THEN
#endif
#ifdef W3_MPI
              CALL mpi_startall ( nrqgo, irqgo , ierr_mpi )
#endif
 
#ifdef W3_MPI
              flgmpi(0) = .true.
              nrqmax    = max( nrqmax , nrqgo )
#endif
#ifdef W3_MPIT
              WRITE (ndst,9043) '1a', nrqgo, nrqmax, napfld
#endif
#ifdef W3_MPI
            END IF
#endif
            !
#ifdef W3_MPI
            IF (nrqgo2.NE.0 ) THEN
#endif
#ifdef W3_MPI
              CALL mpi_startall ( nrqgo2, irqgo2, ierr_mpi )
#endif
#ifdef W3_MPI
              flgmpi(1) = .true.
              nrqmax    = max( nrqmax , nrqgo2 )
#endif
#ifdef W3_MPIT
              WRITE (ndst,9043) '1b', nrqgo2, nrqmax, napfld
#endif
#ifdef W3_MPI
            END IF
          ELSE
#endif
#ifdef W3_PDLIB
            CALL do_output_exchanges(imod)
#endif
#ifdef W3_MPI
          END IF ! IF (.NOT. LPDLIB) THEN
        END IF
#endif
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE AFTER TIME LOOP 1')
        !
#ifdef W3_MPI
        IF ( flout(2) .AND. nrqpo.NE.0 ) THEN
          IF ( dsec21(time,tonext(:,2)).EQ.0. ) THEN
            CALL mpi_startall ( nrqpo, irqpo1, ierr_mpi )
            flgmpi(2) = .true.
            nrqmax    = max( nrqmax , nrqpo )
#endif
#ifdef W3_MPIT
            WRITE (ndst,9043) '2 ', nrqpo, nrqmax, nappnt
#endif
#ifdef W3_MPI
          END IF
        END IF
#endif
        !
#ifdef W3_MPI
        IF ( flout(4) .AND. nrqrs.NE.0 ) THEN
          IF ( dsec21(time,tonext(:,4)).EQ.0. ) THEN
            CALL mpi_startall ( nrqrs, irqrs , ierr_mpi )
            flgmpi(4) = .true.
            nrqmax    = max( nrqmax , nrqrs )
#endif
#ifdef W3_MPIT
            WRITE (ndst,9043) '4 ', nrqrs, nrqmax, naprst
#endif
#ifdef W3_MPI
          END IF
        END IF
#endif
        !
#ifdef W3_MPI
        IF ( flout(8) .AND. nrqrs.NE.0 ) THEN
          IF ( dsec21(time,tonext(:,8)).EQ.0. ) THEN
            CALL mpi_startall ( nrqrs, irqrs , ierr_mpi )
            flgmpi(8) = .true.
            nrqmax    = max( nrqmax , nrqrs )
#endif
#ifdef W3_MPIT
            WRITE (ndst,9043) '8 ', nrqrs, nrqmax, naprst
#endif
#ifdef W3_MPI
          END IF
        END IF
#endif
        !
#ifdef W3_MPI
        IF ( flout(5) .AND. nrqbp.NE.0 ) THEN
          IF ( dsec21(time,tonext(:,5)).EQ.0. ) THEN
            CALL mpi_startall ( nrqbp , irqbp1, ierr_mpi )
            flgmpi(5) = .true.
            nrqmax    = max( nrqmax , nrqbp )
#endif
#ifdef W3_MPIT
            WRITE (ndst,9043) '5a', nrqbp, nrqmax, napbpt
#endif
#ifdef W3_MPI
          END IF
        END IF
#endif
        !
#ifdef W3_MPI
        IF ( flout(5) .AND. nrqbp2.NE.0 .AND. iaproc.EQ.napbpt) THEN
          IF ( dsec21(time,tonext(:,5)).EQ.0. ) THEN
            CALL mpi_startall (nrqbp2,irqbp2,ierr_mpi)
            nrqmax    = max( nrqmax , nrqbp2 )
#endif
#ifdef W3_MPIT
            WRITE (ndst,9043) '5b', nrqbp2, nrqmax, napbpt
#endif
#ifdef W3_MPI
          END IF
        END IF
#endif
        !
#ifdef W3_MPI
        IF ( nrqmax .NE. 0 ) ALLOCATE ( statio(mpi_status_size,nrqmax) )
#endif
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE AFTER TIME LOOP 2')
        !
        ! 4.c Reset next output time
 
        !
        tofrst(1) = -1
        tofrst(2) =  0
        !
        DO j=1, notype
 
          IF ( flout(j) ) THEN
            !
            ! 4.d Perform output
            !
#ifdef W3_NL5
            IF (j .EQ. 2) tosnl5  = tonext(:, 2)
#endif
            tout(:) = tonext(:,j)
            dttst   = dsec21( time, tout )
            !
            IF ( dttst .EQ. 0. ) THEN
              IF ( ( j .EQ. 1 )              &
#ifdef W3_SBS
                   .OR. ( j .EQ. 7 )         &
#endif
                   ) THEN
                IF ( iaproc .EQ. napfld ) THEN
#ifdef W3_MPI
                  IF ( flgmpi(1) ) CALL mpi_waitall ( nrqgo2, irqgo2, statio, ierr_mpi )
                  flgmpi(1) = .false.
#endif
                  !
#ifdef W3_SBS
                  IF ( j .EQ. 1 ) THEN
#endif
                    CALL w3iogo( 'WRITE', nds(7), itest, imod &
#ifdef W3_ASCII
                            ,nds(14)                          &
#endif
                            )
#ifdef W3_SBS
                  ENDIF
#endif
                  !
#ifdef W3_SBS
                  !
                  !     Generate output flag file for fields and SBS coupling.
                  !
                  jj = len_trim( filext )
                  CALL stme21 ( time, idtime )
                  foutname = 'Field_done.' // idtime(1:4) &
                       // idtime(6:7) // idtime(9:10) &
                       // idtime(12:13) // '.' // filext(1:jj)
#endif
                  !
#ifdef W3_SBS
                  OPEN( unit=ndsoflg, file=foutname)
                  CLOSE( ndsoflg )
#endif
                END IF
                !
              ELSE IF ( j .EQ. 2 ) THEN
                !
                !   Point output
                !
                IF ( iaproc .EQ. nappnt ) THEN
                  !
                  !   Gets the necessary spectral data
                  !
                  CALL w3iope ( va )
#ifdef W3_BIN2NC
                  CALL w3iopon ( 'WRITE', nds(8), itest, imod )
#else
                  CALL w3iopo ( 'WRITE', nds(8), itest, imod &                          
#ifdef W3_ASCII
                          ,nds(15)                           &
#endif
                          )
#endif 
                  END IF
                !
              ELSE IF ( j .EQ. 3 ) THEN
                !
                ! Track output
                !
                CALL w3iotr ( nds(11), nds(12), va, imod )
              ELSE IF ( j .EQ. 4 ) THEN
                CALL w3iors ('HOT', nds(6), xxx, imod, flout(8) )
                itest = rstype
              ELSE IF ( j .EQ. 5 ) THEN
                IF ( iaproc .EQ. napbpt ) THEN
#ifdef W3_MPI
                  IF (nrqbp2.NE.0) CALL mpi_waitall ( nrqbp2, irqbp2,statio, ierr_mpi )
#endif
                  CALL w3iobc ( 'WRITE', nds(10),         &
                       time, time, itest, imod )
                END IF
              ELSE IF ( j .EQ. 6 ) THEN
                CALL w3iosf ( nds(13), imod )
#ifdef W3_OASIS
              ELSE IF ( j .EQ. 7 ) THEN
                !
                ! Send variables to atmospheric or ocean circulation or ice model
                !
                IF (dtout(7).NE.0) THEN
                  IF ( (mod(id_oasis_time,nint(dtout(7))) .EQ. 0 ) .AND. &
                       (dsec21(time00, time) .GT. 0.0) ) THEN
                    IF ( (cplt0 .AND. (dsec21(time, timen) .GT. 0.0)) .OR. &
                         .NOT. cplt0 ) THEN
                      IF (cplt0) id_oasis_time = nint(dsec21( time00 , time ))
 
#endif
#ifdef W3_OASACM
                      CALL snd_fields_to_atmos()
#endif
#ifdef W3_OASOCM
                      CALL snd_fields_to_ocean()
#endif
#ifdef W3_OASICM
                      CALL snd_fields_to_ice()
#endif
#ifdef W3_OASIS
                      IF (.NOT. cplt0) id_oasis_time = nint(dsec21( time00 , time ))
                    ENDIF
                  ENDIF
                ENDIF
#endif
              END IF
              !
              CALL tick21 ( tout, dtout(j) )
              tonext(:,j) = tout
              tlst        = tolast(:,j)
              dttst       = dsec21( tout , tlst )
              flout(j)    = dttst.GE.0.
              IF ( flout(j) ) THEN
                outid(2*j-1:2*j-1) = 'X'
#ifdef W3_OASIS
                IF ( (dtout(7).NE.0) .AND.           &
                     (dsec21(time,time00).EQ.0 .OR.  &
                     dsec21(time,timeend).EQ.0) ) outid(13:13) = ' '
#endif
              ELSE
                outid(2*j-1:2*j-1) = 'L'
              END IF
            END IF
            !
            ! 4.e Update next output time
            !
            IF ( flout(j) ) THEN
              IF ( tofrst(1).EQ.-1 ) THEN
                tofrst = tout
              ELSE
                dttst  = dsec21( tout , tofrst )
                IF ( dttst.GT.0.) THEN
                  tofrst = tout
                END IF
              END IF
            END IF
            !
          END IF
          !
        END DO
        !
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE AFTER TIME LOOP 3')
 
        ! If there is a second stream of restart files then J=8 and FLOUT(8)=.TRUE.
        j=8
        IF ( flout(j) ) THEN
          !
          ! 4.d Perform output
          !
          tout(:) = tonext(:,j)
          dttst   = dsec21( time, tout )
          IF ( dttst .EQ. 0. ) THEN
            CALL w3iors ('HOT', nds(6), xxx, imod, flout(8) )
            itest = rstype
            CALL tick21 ( tout, dtout(j) )
            tonext(:,j) = tout
            tlst        = tolast(:,j)
            dttst       = dsec21( tout , tlst )
            flout(j)    = dttst.GE.0.
            IF ( flout(j) ) THEN
              outid(2*j-1:2*j-1) = 'X'
#ifdef W3_OASIS
              IF ( (dtout(7).NE.0) .AND.           &
                   (dsec21(time,time00).EQ.0 .OR.  &
                   dsec21(time,timeend).EQ.0) ) outid(13:13) = ' '
#endif
            ELSE
              outid(2*j-1:2*j-1) = 'L'
            END IF
          END IF
          !
          ! 4.e Update next output time
          !
          IF ( flout(j) ) THEN
            IF ( tofrst(1).EQ.-1 ) THEN
              tofrst = tout
            ELSE
              dttst  = dsec21( tout , tofrst )
              IF ( dttst.GT.0.) THEN
                tofrst = tout
              END IF
            END IF
          END IF
        END IF
        !        END OF CHECKPOINT
        !
        call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE AFTER TIME LOOP 3')
        !
#ifdef W3_MPI
        IF ( flgmpi(0) ) CALL mpi_waitall ( nrqgo, irqgo , statio, ierr_mpi )
        IF ( flgmpi(2) ) CALL mpi_waitall ( nrqpo, irqpo1, statio, ierr_mpi )
        IF ( flgmpi(4) ) CALL mpi_waitall ( nrqrs, irqrs , statio, ierr_mpi )
        IF ( flgmpi(8) ) CALL mpi_waitall ( nrqrs, irqrs , statio, ierr_mpi )
        IF ( flgmpi(5) ) CALL mpi_waitall ( nrqbp, irqbp1, statio, ierr_mpi )
        IF ( nrqmax .NE. 0 ) DEALLOCATE ( statio )
#endif
        !
#ifdef W3_T
        WRITE (ndst,9044)
#endif
      END IF
#ifdef W3_TIMINGS
      CALL print_my_time("Before update log file")
#endif
      call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE AFTER TIME LOOP 4')
      !
      ! 5.  Update log file ------------------------------------------------ /
      !
      IF ( iaproc.EQ.naplog ) THEN
        !
        CALL stme21 ( time , idtime )
        IF ( flcur ) THEN
          dttst  = dsec21( time , tcn )
          IF ( dttst .EQ. 0. ) idact(7:7) = 'X'
        END IF
        IF ( flwind ) THEN
          dttst  = dsec21( time , twn )
          IF ( dttst .EQ. 0. ) idact(3:3) = 'X'
        END IF
        IF ( fltaua ) THEN
          dttst  = dsec21( time , tun )
          IF ( dttst .EQ. 0. ) idact(9:9) = 'X'
        END IF
        IF ( flrhoa ) THEN
          dttst  = dsec21( time , trn )
          IF ( dttst .EQ. 0. ) idact(11:11) = 'X'
        END IF
        IF ( tdn(1) .GT. 0  ) THEN
          dttst  = dsec21( time , tdn )
          IF ( dttst .EQ. 0. ) idact(21:21) = 'X'
        END IF
        !
        IF ( idlast.NE.time(1) ) THEN
          WRITE (ndso,900) itime, ipass, idtime(1:19), idact, outid
          idlast = time(1)
        ELSE
          WRITE (ndso,901) itime, ipass, idtime(12:19), idact, outid
        END IF
        !
      END IF
      !
      idact  = '         '
      outid  = '           '
      flact  = .false.
      !
      ! 6.  If time is not ending time, branch back to 2 ------------------- /
      !
      dttst  = dsec21( time, tend )
      IF ( dttst .EQ. 0. ) EXIT
#ifdef W3_TIMINGS
      CALL print_my_time("Continuing the loop")
#endif
    END DO
    call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE AFTER TIME LOOP 5')
    !
 
    IF ( tstamp .AND. screen.NE.ndso .AND. iaproc.EQ.napout ) THEN
      CALL wwtime ( sttime )
      WRITE (screen,951) sttime
    END IF
 
    IF ( iaproc .EQ. naplog ) WRITE (ndso,902)
    !
    DEALLOCATE(field)
    DEALLOCATE(tauwx, tauwy)
    !
    call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE END W3WAVE')
    !
    RETURN
    !
    ! Formats
    !
900 FORMAT (4x,i6,'|',i6,'| ', a19  ,' | ',a,' | ',a,' |')
901 FORMAT (4x,i6,'|',i6,'| ',11x,a8,' | ',a,' | ',a,' |')
902 FORMAT (2x,'--------+------+---------------------+'                  &
         ,'-----------------------+------------------+')
    !
#ifdef W3_IC3
920 FORMAT ('     Updating k and Cg from ice param. 1,2,3,4.'/)
#endif
950 FORMAT ('  WAVEWATCH III calculating for ',a,' at ',a)
951 FORMAT ('  WAVEWATCH III reached the end of a computation',          &
         ' loop at ',a)
1000 FORMAT (/' *** WAVEWATCH III ERROR IN W3WAVE :'/                    &
         '     ENDING TIME BEFORE STARTING TIME '/)
1001 FORMAT (/' *** WAVEWATCH III ERROR IN W3WAVE :'/                    &
         '     NEW WATER LEVEL BEFORE OLD WATER LEVEL '/)
1002 FORMAT (/' *** WAVEWATCH III ERROR IN W3WAVE :'/                    &
         '     ILLEGAL CURRENT INTERVAL '/)
1003 FORMAT (/' *** WAVEWATCH III ERROR IN W3WAVE :'/                    &
         '     ILLEGAL WIND INTERVAL '/)
1004 FORMAT (/' *** WAVEWATCH III ERROR IN W3WAVE :'/                    &
         '     NEW ICE FIELD BEFORE OLD ICE FIELD '/)
1005 FORMAT (/' *** WAVEWATCH III ERROR IN W3WAVE :'/                    &
         '     NEW IC1 FIELD BEFORE OLD IC1 FIELD '/)
1007 FORMAT (/' *** WAVEWATCH III ERROR IN W3WAVE :'/                    &
         '     NEW ATM MOMENTUM BEFORE OLD ATM MOMENTUM '/)
1008 FORMAT (/' *** WAVEWATCH III ERROR IN W3WAVE :'/                    &
         '     NEW AIR DENSITY BEFORE OLD AIR DENSITY '/)
#ifdef W3_IS2
1006 FORMAT (/' *** WAVEWATCH III ERROR IN W3WAVE :'/                    &
         '     NEW IC5 FIELD BEFORE OLD IC5 FIELD '/)
#endif
1030 FORMAT (/' *** WAVEWATCH III WARING IN W3WAVE :'/                   &
         '     AT LEAST ONE PROCESSOR HAS 0 ACTIVE POINTS',              &
         ' IN GRID',i3)
#ifdef W3_REFRX
1040 FORMAT (/' *** WAVEWATCH III ERROR IN W3WAVE :'/                    &
         '     EXPERIMENTAL FEATURE !/REFRX NOT FULLY IMPLEMENTED.'/)
#endif
    !
#ifdef W3_T
9000 FORMAT (                                                            &
         '============================================================', &
         '===================='/                                         &
         ' TEST W3WAVE : RUN MODEL',i3,' FILEXT [',a,                    &
         '] UP TO ',i8.8,i7.6 /                                          &
         '====================',                                         &
         '============================================================')
9010 FORMAT (' TEST W3WAVE : DT INT. =',f12.1,'   FLZERO = ',l1)
9011 FORMAT (' TEST W3WAVE : DT LEV. =',f12.1)
9012 FORMAT (' TEST W3WAVE : DT CUR. =',f12.1/                           &
         '                        ',f12.1/                               &
         '                        ',f12.1)
9013 FORMAT (' TEST W3WAVE : DT WIND =',f12.1/                           &
         '                        ',f12.1/                               &
         '                        ',f12.1)
9014 FORMAT (' TEST W3WAVE : DT ICE  =',f12.1)
9015 FORMAT (' TEST W3WAVE : DT IC1  =',f12.1)
9016 FORMAT (' TEST W3WAVE : DT IC5  =',f12.1)
9017 FORMAT (' TEST W3WAVE : DT TAU  =',f12.1)
9018 FORMAT (' TEST W3WAVE : DT RHO  =',f12.1)
9020 FORMAT (' TEST W3WAVE : IT0, NT, DTG :',2i4,f8.1)
9021 FORMAT (' TEST W3WAVE : ITIME etc',i6,i4,i10.8,i7.6,1x,2l1,         &
         2f6.2,f7.1,f6.2)
9022 FORMAT (' TEST W3WAVE : SKIP TO 400 IN 3.5')
9023 FORMAT (' TEST W3WAVE : SKIP TO 380 IN 3.5')
9030 FORMAT (' TEST W3WAVE : END OF COMPUTATION LOOP')
9040 FORMAT (' TEST W3WAVE : CHECKING FOR OUTPUT'/                       &
         '               TOFRST           :',i9.8,i7.6/                  &
         '               TND              :',i9.8,i7.6/                  &
         '               DTTST[1], FLAG_O :',2f8.1,l4)
9041 FORMAT (' TEST W3WAVE : PERFORMING OUTPUT')
9042 FORMAT (' TEST W3WAVE : OUTPUT COMPUTATION FLAGS: ',3l2)
#endif
#ifdef W3_MPIT
9043 FORMAT (' TEST W3WAVE : TYPE, NRQ, NRQMAX, NA : ',a2,3i6)
#endif
#ifdef W3_T
9044 FORMAT (' TEST W3WAVE : END OF OUTPUT')
#endif
    !/
    !/ End of W3WAVE ----------------------------------------------------- /
    !/
  END SUBROUTINE w3wave
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3gath ( ISPEC, FIELD )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         26-Dec-2012 |
    !/                  +-----------------------------------+
    !/
    !/    04-Jan-1999 : Distributed FORTRAN 77 version.     ( version 1.18 )
    !/    13-Jan-2000 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/                  Major changes to logistics.
    !/    29-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    13-Jun-2006 : Split STORE in G/SSTORE             ( version 3.09 )
    !/    26-Dec-2012 : Move FIELD init. to W3GATH.         ( version 4.OF )
    !/
    !  1. Purpose :
    !
    !     Gather spectral bin information into a propagation field array.
    !
    !  2. Method :
    !
    !     Direct copy or communication calls (MPP version).
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       ISPEC   Int.   I   Spectral bin considered.
    !       FIELD   R.A.   O   Full field to be propagated.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !
    !      MPI_STARTALL, MPI_WAITALL
    !                Subr. mpif.h   MPI persistent comm. routines (!/MPI).
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !     - The field is extracted but not converted.
    !     - MPI version requires posing of send and receive calls in
    !       W3WAVE to match local calls.
    !     - MPI version does not require an MPI_TESTALL call for the
    !       posted gather operation as MPI_WAITALL is mandatory to
    !       reset persistent communication for next time step.
    !     - MPI version allows only two new pre-fetch postings per
    !       call to minimize chances to be slowed down by gathers that
    !       are not yet needed, while maximizing the pre-loading
    !       during the early (low-frequency) calls to the routine
    !       where the amount of calculation needed for proagation is
    !       the largest.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/SHRD  Switch for message passing method.
    !     !/MPI   Id.
    !
    !     !/S     Enable subroutine tracing.
    !     !/MPIT  MPI test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    USE w3gdatmd, ONLY: nspec, nx, ny, nsea, nseal, mapsf, dmin
    USE w3parall, ONLY: init_get_isea
    USE w3wdatmd, ONLY: a => va
#ifdef W3_MPI
    USE w3adatmd, ONLY: mpibuf, bstat, ibfloc, isploc, bispl, &
         nsploc, nrqsg2, irqsg2, gstore
    USE w3odatmd, ONLY: ndst, iaproc, naproc, notype
#endif
    !/
    !
#ifdef W3_MPI
    include "mpif.h"
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: ISPEC
    REAL, INTENT(OUT)       :: FIELD(1-NY:NY*(NX+2))
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
#ifdef W3_SHRD
    INTEGER                 :: ISEA, IXY
#endif
#ifdef W3_MPI
    INTEGER                 :: STATUS(MPI_STATUS_SIZE,NSPEC),  &
         IOFF, IERR_MPI, JSEA, ISEA,     &
         IXY, IS0, IB0, NPST, J
#endif
#ifdef W3_S
    INTEGER, SAVE           :: IENT
#endif
#ifdef W3_MPIT
    CHARACTER(LEN=15)       :: STR(MPIBUF), STRT
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3GATH')
#endif
    !
    field  = 0.
    !
    ! 1.  Shared memory version ------------------------------------------ /
    !
#ifdef W3_SHRD
    DO isea=1, nsea
      ixy        = mapsf(isea,3)
      field(ixy) = a(ispec,isea)
    END DO
#endif
    !
#ifdef W3_SHRD
    RETURN
#endif
    !
    ! 2.  Distributed memory version ( MPI ) ----------------------------- /
    ! 2.a Update counters
    !
#ifdef W3_MPI
    isploc = isploc + 1
    ibfloc = ibfloc + 1
    IF ( ibfloc .GT. mpibuf ) ibfloc = 1
#endif
    !
#ifdef W3_MPIT
    IF ( isploc .EQ. 1 ) THEN
      str = '--------------+'
      WRITE (ndst,9000) str
    END IF
    str    = '              |'
    strt   = str(ibfloc)
    strt(9:9) = 'A'
#endif
    !
    ! 2.b Check status of present buffer
    ! 2.b.1 Scatter (send) still in progress, wait to end
    !
#ifdef W3_MPI
    IF ( bstat(ibfloc) .EQ. 2 ) THEN
      ioff =  1 + (bispl(ibfloc)-1)*nrqsg2
      IF ( nrqsg2 .GT. 0 ) CALL mpi_waitall ( nrqsg2, irqsg2(ioff,2), status, ierr_mpi )
      bstat(ibfloc) = 0
#endif
#ifdef W3_MPIT
      strt(13:13) = 'S'
#endif
#ifdef W3_MPI
    END IF
#endif
    !
    ! 2.b.2 Gather (recv) not yet posted, post now
    !
#ifdef W3_MPI
    IF ( bstat(ibfloc) .EQ. 0 ) THEN
      bstat(ibfloc) = 1
      bispl(ibfloc) = isploc
      ioff =  1 + (isploc-1)*nrqsg2
      IF ( nrqsg2 .GT. 0 ) CALL mpi_startall ( nrqsg2, irqsg2(ioff,1), ierr_mpi )
#endif
#ifdef W3_MPIT
      strt(10:10) = 'g'
#endif
#ifdef W3_MPI
    END IF
#endif
    !
    ! 2.c Put local spectral densities in store
    !
#ifdef W3_MPI
    DO jsea=1, nseal
      CALL init_get_isea(isea, jsea)
      gstore(isea,ibfloc) = a(ispec,jsea)
    END DO
#endif
    !
    ! 2.d Wait for remote spectral densities
    !
#ifdef W3_MPI
    ioff =  1 + (bispl(ibfloc)-1)*nrqsg2
    IF ( nrqsg2 .GT. 0 ) CALL mpi_waitall ( nrqsg2, irqsg2(ioff,1), status, ierr_mpi )
#endif
    !
#ifdef W3_MPIT
    strt(11:11) = 'G'
    WRITE (strt(1:7),'(I2,I5)') bstat(ibfloc), isploc
    str(ibfloc) = strt
#endif
    !
    ! 2.e Convert storage array to field.
    !
#ifdef W3_MPI
    DO isea=1, nsea
      ixy        = mapsf(isea,3)
      field(ixy) = gstore(isea,ibfloc)
    END DO
#endif
    !
    ! 2.f Pre-fetch data in available buffers
    !
#ifdef W3_MPI
    is0    = isploc
    ib0    = ibfloc
    npst   = 0
#endif
    !
#ifdef W3_MPI
    DO j=1, mpibuf-1
      is0    = is0 + 1
      IF ( is0 .GT. nsploc ) EXIT
      ib0    = 1 + mod(ib0,mpibuf)
      IF ( bstat(ib0) .EQ. 0 ) THEN
        bstat(ib0) = 1
        bispl(ib0) = is0
        ioff       = 1 + (is0-1)*nrqsg2
        IF ( nrqsg2 .GT. 0 ) CALL mpi_startall ( nrqsg2, irqsg2(ioff,1), ierr_mpi )
        npst       = npst + 1
#endif
#ifdef W3_MPIT
        strt        = str(ib0)
        strt(10:10) = 'g'
        WRITE (strt(1:7),'(I2,I5)') bstat(ib0), bispl(ib0)
        str(ib0)    = strt
#endif
#ifdef W3_MPI
      END IF
      IF ( npst .GE. 2 ) EXIT
    END DO
#endif
    !
    ! 2.g Test output
    !
#ifdef W3_MPIT
    DO ib0=1, mpibuf
      strt   = str(ib0)
      IF ( strt(2:2) .EQ. ' ' ) THEN
        IF ( bstat(ib0) .EQ. 0 ) THEN
          WRITE (strt(1:2),'(I2)') bstat(ib0)
        ELSE
          WRITE (strt(1:7),'(I2,I5)') bstat(ib0), bispl(ib0)
        END IF
        str(ib0) = strt
      END IF
    END DO
    WRITE (ndst,9010) isploc, str
#endif
    !
#ifdef W3_MPI
    RETURN
#endif
    !
    ! Formats
    !
#ifdef W3_MPIT
9000 FORMAT ( ' TEST OF BUFFER MANAGEMENT MPI :'/              &
         ' -------------------------------'/              &
         '      RECORDS ALTERNATELY WRITTEN BY W3GATH AND W3SCAT'/ &
         '      FRIST COLLUMN  : LOCAL ISPEC'/                     &
         '      OTHER COLLUMNS : BUFFER STATUS INDICATOR '/        &
         '                        0 : INACTIVE'/                   &
         '                        1 : RECEIVING'/                  &
         '                        2 : SENDING'/                    &
         '                       LOCAL ISPEC FOR BUFFER'/          &
         '                       A  : ACTIVE BUFFER'/              &
         '                       g/G: START/FINISH RECIEVE'/       &
         '                       s/S: START/FINISH SEND'/          &
         ' +-----+',8a15)
9010 FORMAT ( ' |',i4,' |',8a15)
#endif
    !/
    !/ End of W3GATH ----------------------------------------------------- /
    !/
  END SUBROUTINE w3gath
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3scat ( ISPEC, MAPSTA, FIELD )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         13-Jun-2006 |
    !/                  +-----------------------------------+
    !/
    !/    04-Jan-1999 : Distributed FORTRAN 77 version.     ( version 1.18 )
    !/    13-Jan-2000 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/                  Major changes to logistics.
    !/    28-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    07-Sep-2005 : Updated boundary conditions.        ( version 3.08 )
    !/    13-Jun-2006 : Split STORE in G/SSTORE             ( version 3.09 )
    !/
    !  1. Purpose :
    !
    !     'Scatter' data back to spectral storage after propagation.
    !
    !  2. Method :
    !
    !     Direct copy or communication calls (MPP version).
    !     See also W3GATH.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       ISPEC   Int.   I   Spectral bin considered.
    !       MAPSTA  I.A.   I   Status map for spatial grid.
    !       FIELD   R.A.   I   Full field to be propagated.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !
    !      MPI_STARTALL, MPI_WAITALL, MPI_TESTALL
    !                Subr. mpif.h   MPI persistent comm. routines (!/MPI).
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !     None.
    !
    !  7. Remarks :
    !
    !     - The field is put back but not converted !
    !     - MPI persistent communication calls initialize in W3MPII.
    !     - See W3GATH and W3MPII for additional comments on data
    !       buffering.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/SHRD  Switch for message passing method.
    !     !/MPI   Id.
    !
    !     !/S     Enable subroutine tracing.
    !     !/MPIT  MPI test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3gdatmd, ONLY: nsea, nseal, mapsf, nspec, nx, ny
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    USE w3wdatmd, ONLY: a => va
#ifdef W3_MPI
    USE w3adatmd, ONLY: mpibuf, bstat, ibfloc, isploc, bispl, &
         nsploc, nrqsg2, irqsg2, sstore
#endif
    USE w3odatmd, ONLY: ndst
#ifdef W3_MPI
    USE w3odatmd, ONLY: iaproc, naproc
#endif
    USE constants, ONLY : lpdlib
    USE w3parall, only: init_get_isea
    !/
    !
#ifdef W3_MPI
    include "mpif.h"
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: ISPEC, MAPSTA(NY*NX)
    REAL, INTENT(IN)        :: FIELD(1-NY:NY*(NX+2))
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
#ifdef W3_SHRD
    INTEGER                 :: ISEA, IXY
#endif
#ifdef W3_MPI
    INTEGER                 :: ISEA, IXY, IOFF, IERR_MPI, J,   &
         STATUS(MPI_STATUS_SIZE,NSPEC),  &
         JSEA, IB0
#endif
#ifdef W3_S
    INTEGER, SAVE           :: IENT
#endif
#ifdef W3_MPIT
    CHARACTER(LEN=15)       :: STR(MPIBUF), STRT
#endif
#ifdef W3_MPI
    LOGICAL                 :: DONE
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3SCAT')
#endif
    !
    ! 1.  Shared memory version ------------------------------------------ *
    !
#ifdef W3_SHRD
    DO isea=1, nsea
      ixy           = mapsf(isea,3)
      IF ( mapsta(ixy) .NE. 0 ) a(ispec,isea) = field(ixy)
    END DO
#endif
    !
#ifdef W3_SHRD
    RETURN
#endif
    !
    ! 2.  Distributed memory version ( MPI ) ----------------------------- *
    ! 2.a Initializations
    !
#ifdef W3_MPIT
    DO ib0=1, mpibuf
      str(ib0) = '              |'
    END DO
#endif
    !
#ifdef W3_MPIT
    strt   = str(ibfloc)
    strt(9:9) = 'A'
#endif
    !
    ! 2.b Convert full grid to sea grid, active points only
    !
#ifdef W3_MPI
    DO isea=1, nsea
      ixy    = mapsf(isea,3)
      IF ( mapsta(ixy) .NE. 0 ) sstore(isea,ibfloc) = field(ixy)
    END DO
#endif
    !
    ! 2.c Send spectral densities to appropriate remote
    !
#ifdef W3_MPI
    ioff   = 1 + (isploc-1)*nrqsg2
    IF ( nrqsg2 .GT. 0 ) CALL mpi_startall ( nrqsg2, irqsg2(ioff,2), ierr_mpi )
    bstat(ibfloc) = 2
#endif
#ifdef W3_MPIT
    strt(12:12) = 's'
    WRITE (strt(1:7),'(I2,I5)') bstat(ibfloc), isploc
    str(ibfloc) = strt
#endif
    !
    ! 2.d Save locally stored results
    !
#ifdef W3_MPI
    DO jsea=1, nseal
      CALL init_get_isea(isea, jsea)
      ixy    = mapsf(isea,3)
      IF (mapsta(ixy) .NE. 0) a(ispec,jsea) = sstore(isea,ibfloc)
    END DO
#endif
    !
    ! 2.e Check if any sends have finished
    !
#ifdef W3_MPI
    ib0    = ibfloc
#endif
    !
#ifdef W3_MPI
    DO j=1, mpibuf
      ib0    = 1 + mod(ib0,mpibuf)
      IF ( bstat(ib0) .EQ. 2 ) THEN
        ioff   = 1 + (bispl(ib0)-1)*nrqsg2
        IF ( nrqsg2 .GT. 0 ) THEN
          CALL mpi_testall ( nrqsg2, irqsg2(ioff,2), done, status, ierr_mpi )
        ELSE
          done   = .true.
        END IF
        IF ( done .AND. nrqsg2.GT.0 ) THEN
          CALL mpi_waitall ( nrqsg2, irqsg2(ioff,2), status, ierr_mpi )
        END IF
        IF ( done ) THEN
          bstat(ib0) = 0
#endif
#ifdef W3_MPIT
          strt        = str(ib0)
          WRITE (strt(1:7),'(I2,I5)') bstat(ib0), bispl(ib0)
          strt(13:13) = 'S'
          str(ib0)    = strt
#endif
#ifdef W3_MPI
        END IF
      END IF
    END DO
#endif
    !
    ! 2.f Last component, finish message passing, reset buffer control
    !
#ifdef W3_MPI
    IF ( isploc .EQ. nsploc ) THEN
#endif
      !
#ifdef W3_MPI
      DO ib0=1, mpibuf
        IF ( bstat(ib0) .EQ. 2 ) THEN
          ioff   = 1 + (bispl(ib0)-1)*nrqsg2
          IF ( nrqsg2 .GT. 0 ) CALL mpi_waitall ( nrqsg2, irqsg2(ioff,2), status, ierr_mpi )
          bstat(ib0) = 0
#endif
#ifdef W3_MPIT
          strt        = str(ib0)
          WRITE (strt(1:7),'(I2,I5)') bstat(ib0), bispl(ib0)
          strt(13:13) = 'S'
          str(ib0)    = strt
#endif
#ifdef W3_MPI
        END IF
      END DO
#endif
      !
#ifdef W3_MPI
      isploc = 0
      ibfloc = 0
#endif
      !
#ifdef W3_MPI
    END IF
#endif
    !
    ! 2.g Test output
    !
#ifdef W3_MPIT
    DO ib0=1, mpibuf
      strt   = str(ib0)
      IF ( strt(2:2) .EQ. ' ' ) THEN
        IF ( bstat(ib0) .EQ. 0 ) THEN
          WRITE (strt(1:2),'(I2)') bstat(ib0)
        ELSE
          WRITE (strt(1:7),'(I2,I5)') bstat(ib0), bispl(ib0)
        END IF
        str(ib0) = strt
      END IF
    END DO
#endif
    !
#ifdef W3_MPIT
    WRITE (ndst,9000) str
#endif
    !
#ifdef W3_MPIT
    IF ( isploc .EQ. 0 ) THEN
      DO ib0=1, mpibuf
        str(ib0) = '--------------+'
      END DO
      WRITE (ndst,9010) str
      WRITE (ndst,*)
    END IF
#endif
    !
#ifdef W3_MPI
    RETURN
#endif
    !
    ! Formats
    !
#ifdef W3_MPIT
9000 FORMAT ( ' |     |',8a15)
9010 FORMAT ( ' +-----+',8a15)
#endif
    !/
    !/ End of W3SCAT ----------------------------------------------------- /
    !/
  END SUBROUTINE w3scat
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3nmin ( MAPSTA, FLAG0 )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         28-Dec-2004 |
    !/                  +-----------------------------------+
    !/
    !/    23-Feb-2001 : Origination.                        ( version 2.07 )
    !/    28-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/
    !  1. Purpose :
    !
    !     Check minimum number of active sea points at given processor to
    !     evaluate the need for a MPI_BARRIER call.
    !
    !  2. Method :
    !
    !     Evaluate mapsta.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       MAPSTA  I.A.   I   Status map for spatial grid.
    !       FLAG0   log.   O   Flag to identify 0 as minimum.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !     None.
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S     Enable subroutine tracing.
    !     !/T     Test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    USE w3gdatmd, ONLY: nsea, mapsf, nx, ny
    USE w3odatmd, ONLY: ndst, naproc
    USE w3parall, ONLY: init_get_jsea_isproc
    !/
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: MAPSTA(NY*NX)
    LOGICAL, INTENT(OUT)    :: FLAG0
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: NMIN, IPROC, NLOC, ISEA, IXY
    INTEGER                 :: JSEA, ISPROC
#ifdef W3_S
    INTEGER, SAVE           :: IENT
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3NMIN')
#endif
    !
    nmin   = nsea
    !
#ifdef W3_OMPG
    !$OMP PARALLEL PRIVATE (IPROC,NLOC,ISEA,JSEA,ISPROC,IXY,NMIN)
    !$OMP DO SCHEDULE (DYNAMIC,1)
#endif
    DO iproc=1, naproc
      nloc   = 0
      DO isea=1, nsea
        CALL init_get_jsea_isproc(isea, jsea, isproc)
        IF (isproc .eq. iproc) THEN
          ixy    = mapsf(isea,3)
          IF ( mapsta(ixy) .EQ. 1 ) nloc = nloc + 1
        END IF
      END DO
#ifdef W3_SMC
      !!Li   For SMC grid, local sea points are equally NSEA/NAPROC
      !!Li   so the NLOC is overwirte by a constant.
      nloc = nsea/naproc
#endif
      !
#ifdef W3_T
      WRITE (ndst,9000) iproc, nloc
#endif
      nmin   = min( nmin , nloc )
    END DO
#ifdef W3_OMPG
    !$OMP END DO
    !$OMP END PARALLEL
#endif
    !
    flag0  = nmin .EQ. 0
#ifdef W3_T
    WRITE (ndst,9001) nmin, flag0
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT ( ' TEST W3NMIN : IPROC =',i3,'  NLOC =',i5)
9001 FORMAT ( ' TEST W3NMIN : NMIN =',i5,'  FLAG0 =',l2)
#endif
    !/
    !/ End of W3NMIN ----------------------------------------------------- /
    !/
  END SUBROUTINE w3nmin
  !/
  !/ End of module W3WAVEMD -------------------------------------------- /
  !/
END MODULE w3wavemd