gx_outp.F90

Go to the documentation of this file.

 
 
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
PROGRAM gxoutp
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |            J.H. Alves             |
  !/                  |            F. Ardhuin             |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         27-Aug-2015 |
  !/                  +-----------------------------------+
  !/
  !/    30-Jun-1999 : Final FORTRAN 77                    ( version 1.18 )
  !/    24-Jan-2000 : Upgrade to FORTRAN 90               ( version 2.00 )
  !/    14-Feb-2000 : Exact nonlinear interactions        ( version 2.01 )
  !/    25-Jan-2001 : Cartesian grid version              ( version 2.06 )
  !/    02-Feb-2001 : Xnl version 3.0                     ( version 2.07 )
  !/    13-Nov-2002 : Add stress vector.                  ( version 3.00 )
  !/    27-Nov-2002 : First version of VDIA and MDIA.     ( version 3.01 )
  !/    01-Aug-2003 : Fix format for SH output points.    ( version 3.03 )
  !/    24-Dec-2004 : Multiple grid version.              ( version 3.06 )
  !/    23-Jun-2006 : Linear input added.                 ( version 3.09 )
  !/    29-Jun-2006 : Adding file name preamble.          ( version 3.09 )
  !/    03-Jul-2006 : Separate flux modules.              ( version 3.09 )
  !/    25-Jul-2006 : Grid ID for each point.             ( version 3.10 )
  !/    25-Apr-2007 : EMEAN in W3SPR2 par list.           ( version 3.11 )
  !/    09-Oct-2007 : WAM 4+ Sin and Sds added.           ( version 3.13 )
  !/                  (F. Ardhuin)
  !/    09-Oct-2007 : Experimental Sbs (BS1) added.       ( version 3.13 )
  !/                  (F. Ardhuin)
  !/    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-Aug-2010 : Adding ST4                          ( version 3.14 )
  !/    20-Apr-2010 : Fix initialization of USTAR.      ( version 3.14.1 )
  !/    23-Aug-2012 : Adding movable bed friction BT4     ( version 4.07 )
  !/    16-Jul-2012 : Move GMD (SNL3) and nonlinear filter (SNLS)
  !/                  from 3.15 (HLT).                    ( version 4.08 )
  !/    26-Dec-2012 : Modified obsolete declarations.     ( version 4.11 )
  !/    27-Aug-2015 : Sice add as additional output       ( version 5.10 )
  !/                  (in source terms)
  !/    19-Jul-2021 : Momentum and air density support    ( version 7.14 )
  !/
  !/    Copyright 2009-2012 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 :
  !
  !     Post-processing of point output for GrADS post-processing.
  !
  !  2. Method :
  !
  !     In order to be able to plot spectra and source terms as
  !     fields, spectral data is written as if it is fields data.
  !     The spectral direction becomes the longitude, 90.-FREQ
  !     become the latitude. This way, polar plots can be made
  !     using the GrADS 'NPS' map option. The level or z coordinate
  !     is used to store spectra and source terms for separate
  !     output points. The name of the output point is stored in
  !     the control file as the 'description' of the field.
  !     Also written is a separate file with mean input and wave
  !     parameters. This file contains per level and per time a
  !     single line containing :
  !
  !        Station ID, Longitude, Latitude, Depth, , Wind speed.
  !           U and V components, Air-Sea Temperature difference,
  !           Current velocity, U and V components, Significant
  !           wave height.
  !
  !     The files generated are :
  !
  !       ww3.spec.ctl     GrADS control file.
  !       ww3.spec.grads   GrADS data file.
  !       ww3.mean.grads   File with additional input and wave
  !                        parameters.
  !
  !     The first direction set to 90 degr. Grads NPS plot should
  !     therefore have 'set lon -180 180' for oceanographic directional
  !     convention.
  !
  !     Examples of using the three files can be found in spec.gs and
  !     source.gs.
  !
  !  3. Parameters :
  !
  !  4. Subroutines used :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !      W3NMOD    Subr. W3GDATMD Set number of model.
  !      W3SETG    Subr.   Id.    Point to selected model.
  !      W3NDAT    Subr. W3WDATMD Set number of model for wave data.
  !      W3SETW    Subr.   Id.    Point to selected model for wave data.
  !      W3NAUX    Subr. W3ADATMD Set number of model for aux data.
  !      W3SETA    Subr.   Id.    Point to selected model for aux data.
  !      W3NOUT    Subr. W3ODATMD Set number of model for output.
  !      W3SETO    Subr.   Id.    Point to selected model for output.
  !      ITRACE    Subr. W3SERVMD Subroutine tracing initialization.
  !      STRACE    Subr.   Id.    Subroutine tracing.
  !      NEXTLN    Subr.   Id.    Get next line from input filw
  !      EXTCDE    Subr.   Id.    Abort program as graceful as possible.
  !      STME21    Subr. W3TIMEMD Convert time to string.
  !      TICK21    Subr.   Id.    Advance time.
  !      DSEC21    Func.   Id.    Difference between times.
  !      W3IOGR    Subr. W3IOGRMD Reading/writing model definition file.
  !      W3IOPO    Subr. W3IOPOMD Reading/writing raw point output file.
  !      GXEXPO    Subr. Internal Execute point output.
  !     ----------------------------------------------------------------
  !
  !  5. Called by :
  !
  !     None, stand-alone program.
  !
  !  6. Error messages :
  !
  !     Checks on input, checks in W3IOxx.
  !     Check on grid type.
  !
  !  7. Remarks :
  !
  !     - Curvilinear grids currently not supported.
  !
  !  8. Structure :
  !
  !     See source code.
  !
  !  9. Switches :
  !
  !       !/S    Enable subroutine tracing.
  !
  ! 10. Source code :
  !
  !/ ------------------------------------------------------------------- /
  USE constants
  !/
  !     USE W3GDATMD, ONLY: W3NMOD, W3SETG
  USE w3wdatmd, ONLY: w3setw, w3ndat
#ifdef W3_NL1
  USE w3adatmd, ONLY: w3seta, w3naux
#endif
  USE w3odatmd, ONLY: w3seto, w3nout
  USE w3iogrmd, ONLY: w3iogr
#ifdef W3_BIN2NC
  USE w3iopomd, ONLY: w3iopon
#else 
  USE w3iopomd, ONLY: w3iopo
#endif
  USE w3servmd, ONLY : itrace, nextln, extcde
#ifdef W3_S
  USE w3servmd, ONLY : strace
#endif
  USE w3timemd, ONLY: stme21, tick21, dsec21
  !/
  USE w3gdatmd
  USE w3wdatmd, ONLY: time
  USE w3odatmd, ONLY: ndse, ndst, ndso, nopts, ptloc, ptnme,      &
       dpo, wao, wdo, aso, cao, cdo, spco, fnmpre, &
       grdid, iceo, iceho, icefo
#ifdef W3_FLX5
  USE w3odatmd, ONLY: tauao, taudo, dairo
#endif
  !
  IMPLICIT NONE
  !/
  !/ ------------------------------------------------------------------- /
  !/ Local parameters
  !/
  INTEGER                 :: ndsi, ndsm, ndsop, ndsgrd, ndspnt,   &
       ndscgr, ndstrc, ntrace, ierr,        &
       iotest, i, tout(2), nout, tdum(2),   &
       nreq, ipoint, nlev, iout, time0(2),  &
       ih0, im0, id0, iid, ij0, iinc, ik,   &
       ireq, timen(2), j
#ifdef W3_S
  INTEGER, SAVE           :: ient   = 0
#endif
  REAL                    :: dtreq, dtest
  REAL                    :: undefp = -99.e20
  REAL                    :: fact
  LOGICAL                 :: flsrce(7)
  LOGICAL, ALLOCATABLE    :: flreq(:)
  CHARACTER               :: comstr*1, idtime*23, iddday*11,      &
       cinc*2
  CHARACTER(LEN=3)        :: mnth(12)
  CHARACTER(LEN=25)       :: idsrce(7)
  !/
  !/ ------------------------------------------------------------------- /
  !/
  DATA idsrce / 'Spectrum                 ' ,                     &
       'Wind-wave interactions   ' ,                     &
       'Nonlinear interactions   ' ,                     &
       'Dissipation              ' ,                     &
       'Wave-bottom interactions ' ,                     &
       'Wave-ice interactions ' ,                     &
       'Sum of selected sources  ' /
  DATA flsrce / .false. , .false. , .false. ,                     &
       .false. , .false. , .false., .false. /
  DATA time0  / -1, 0 /
  DATA mnth   / 'JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN',         &
       'JUL', 'AUG', 'SEP', 'OCT', 'NOV', 'DEC' /
  !
  !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  ! 1.  IO set-up.
  !
  CALL w3nmod ( 1, 6, 6 )
  CALL w3setg ( 1, 6, 6 )
  CALL w3ndat (    6, 6 )
  CALL w3setw ( 1, 6, 6 )
#ifdef W3_NL1
  CALL w3naux (    6, 6 )
  CALL w3seta ( 1, 6, 6 )
#endif
  CALL w3nout (    6, 6 )
  CALL w3seto ( 1, 6, 6 )
  !
  ndsi   = 10
  ndsm   = 20
  ndsop  = 20
  ndsgrd = 30
  ndspnt = 31
  ndscgr = 32
  !
  ndstrc =  6
  ntrace =  0
  !
  WRITE (ndso,900)
  !
  CALL itrace ( ndstrc, ntrace )
#ifdef W3_S
  CALL strace (ient, 'GXOUTP')
#endif
  !
  j      = len_trim(fnmpre)
  OPEN (ndsi,file=fnmpre(:j)//'gx_outp.inp',status='OLD',         &
       err=800,iostat=ierr)
  READ (ndsi,'(A)',END=801,ERR=802) comstr
  IF (comstr.EQ.' ') comstr = '$'
  WRITE (ndso,901) comstr
  !
  OPEN (ndsgrd,file=fnmpre(:j)//'ww3.spec.grads',                 &
       form='UNFORMATTED', convert=file_endian)
  OPEN (ndspnt,file=fnmpre(:j)//'ww3.mean.grads',form='FORMATTED')
  OPEN (ndscgr,file=fnmpre(:j)//'ww3.spec.ctl',form='FORMATTED')
  !
  !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  ! 2.  Read model definition file.
  !
  CALL w3iogr ( 'READ', ndsm )
  WRITE (ndso,920) gname
  IF ( flagll ) THEN
    fact   = 1.
  ELSE
    fact   = 1.e-3
  END IF
  !
  !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  ! 3.  Read general data and first fields from file
  !
#ifdef W3_BIN2NC
  CALL w3iopon ( 'READ', ndsop, iotest )
#else
  CALL w3iopo ( 'READ', ndsop, iotest )
#endif
  ALLOCATE ( flreq(nopts) )
  !
  WRITE (ndso,930)
  DO i=1, nopts
    IF ( flagll ) THEN
      WRITE (ndso,931) ptnme(i), fact*ptloc(1,i), fact*ptloc(2,i)
    ELSE
      WRITE (ndso,932) ptnme(i), fact*ptloc(1,i), fact*ptloc(2,i)
    END IF
  END DO
  !
  !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  ! 4.  Read requests from input file.
  !     Output times
  !
  CALL nextln ( comstr , ndsi , ndse )
  READ (ndsi,*,END=801,ERR=802) TOUT, DTREQ, nout
  dtreq  = max( 0. , dtreq )
  IF ( dtreq.EQ.0 ) nout = 1
  nout   = max( 1 , nout )
  !
  CALL stme21 ( tout , idtime )
  WRITE (ndso,940) idtime
  !
  tdum = 0
  CALL tick21 ( tdum , dtreq )
  CALL stme21 ( tdum , idtime )
  IF ( dtreq .GE. 86400. ) THEN
    WRITE (iddday,'(I10,1X)') int(dtreq/86400.)
  ELSE
    iddday = '           '
  END IF
  idtime(1:11) = iddday
  idtime(21:23) = '   '
  WRITE (ndso,941) idtime, nout
  !
  ! ... Output points
  !
  flreq = .false.
  nreq   = 0
  !
  DO
    CALL nextln ( comstr , ndsi , ndse )
    READ (ndsi,*,END=801,ERR=802) ipoint
    IF ( ipoint .GT. 0 ) THEN
      IF ( ipoint .LE. nopts ) THEN
        IF ( .NOT. flreq(ipoint) ) nreq   = nreq + 1
        flreq(ipoint) = .true.
      END IF
    ELSE
      EXIT
    END IF
  END DO
  !
  ! ... Output of output points
  !
  WRITE (ndso,950) nreq
  DO i=1, nopts
    IF (flreq(i)) THEN
      IF ( flagll ) THEN
        WRITE (ndso,951) ptnme(i), fact*ptloc(1,i), &
             fact*ptloc(2,i)
      ELSE
        WRITE (ndso,956) ptnme(i), fact*ptloc(1,i), &
             fact*ptloc(2,i)
      END IF
    END IF
  END DO
  !
  ! ... Output of output points
  !
  CALL nextln ( comstr , ndsi , ndse )
  READ (ndsi,*,END=801,ERR=802) flsrce
  WRITE (ndso,952)
  nlev   = 0
  DO i=1, 7
    IF ( flsrce(i) ) THEN
      WRITE (ndst,953) idsrce(i)
      nlev   = nlev + 1
    END IF
  END DO
  !
  WRITE (ndso,955)
  !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  ! 5.  Time management.
  !
  iout   = 0
  !
  DO
    dtest  = dsec21( time , tout )
    IF ( dtest .GT. 0. ) THEN
#ifdef W3_BIN2NC
      CALL w3iopon ( 'READ', ndsop, iotest )
#else 
      CALL w3iopo ( 'READ', ndsop, iotest )
#endif
      IF ( iotest .EQ. -1 ) THEN
        WRITE (ndso,998)
        EXIT
      END IF
      cycle
    END IF
    IF ( dtest .LT. 0. ) THEN
      CALL tick21 ( tout , dtreq )
      cycle
    END IF
    !
    iout   = iout + 1
    CALL stme21 ( tout , idtime )
    !
    CALL gxexpo
    timen  = tout
    !
    IF ( time0(1) .EQ. -1 ) time0 = time
    !
    CALL tick21 ( tout , dtreq )
    IF ( iout .GE. nout ) EXIT
  END DO
  !
  !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  ! 6.  Close data file and write control file
  ! 6.a Close data sets
  !
  WRITE (ndso,960)
  !
  WRITE (ndso,961)
  CLOSE (ndsgrd)
  CLOSE (ndspnt)
  !
  WRITE (ndso,962)
  !
  ! 6.b Set up timing info
  !
  ih0    = time0(2)/10000
  im0    = mod(time0(2)/100,100)
  id0    = mod(time0(1),100)
  iid    = mod(time0(1)/100,100)
  ij0    = time0(1)/10000
  !
  IF ( iout .GT. 1 ) dtreq  = dsec21( time0, timen ) / real(iout-1)
  IF ( iout .EQ. 1 ) dtreq  = 3600.
  IF ( dtreq .GT. 3599. ) THEN
    cinc   = 'HR'
    iinc   = nint(dtreq/3600.)
    IF ( mod(nint(dtreq),3600) .NE. 0 ) GOTO 820
  ELSE
    cinc   = 'MN'
    iinc   = nint(dtreq/60.)
  END IF
  !
  WRITE (ndso,963) iout, ih0, im0, id0, mnth(iid), ij0, iinc, cinc
  !
  ! 6.c Write control file for spectral data
  !
  WRITE (ndso,964)
  !
  WRITE (ndscgr,970) undefp, nth, 90.+th(1)*rade, dth*rade,       &
       nk, (90.-tpiinv*sig(ik),ik=nk,max(1,nk-4),-1)
  WRITE (ndscgr,971) (90.-tpiinv*sig(ik),ik=nk-5,1,-1)
  WRITE (ndscgr,972) nlev, 1., 1.,                                &
       iout, ih0, im0, id0, mnth(iid), ij0,         &
       iinc, cinc, nreq
  !
  ireq    = 0
  DO i=1, nopts
    IF ( flreq(i) ) THEN
      ireq   = ireq + 1
      WRITE (ndscgr,973) ireq, nlev, 99, ptnme(i)
    END IF
  END DO
  !
  WRITE (ndscgr,974)
  !
  GOTO 888
  !
  ! Escape locations read errors :
  !
800 CONTINUE
  WRITE (ndse,1000) ierr
  CALL extcde ( 10 )
  !
801 CONTINUE
  WRITE (ndse,1001)
  CALL extcde ( 11 )
  !
802 CONTINUE
  WRITE (ndse,1002) ierr
  CALL extcde ( 12 )
  !
820 CONTINUE
  WRITE (ndse,1020) dtreq
  CALL extcde ( 20 )
  !
821 CONTINUE
  WRITE (ndse,1021)
  CALL extcde ( 21 )
  !
888 CONTINUE
  !
  WRITE (ndso,999)
  !
  ! Formats
  !
900 FORMAT (/12x,'    *** WAVEWATCH III GrADS point output post.***    '/ &
       12x,'====================================================='/)
901 FORMAT ( '  Comment character is ''',a,''''/)
  !
920 FORMAT ( '  Grid name : ',a/)
  !
930 FORMAT ( '  Points in file : '/                                 &
       ' ------------------------------------')
 
931 FORMAT ( '      ',a,2f10.2)
 
932 FORMAT ( '      ',a,2(f8.1,'E3'))
  !
940 FORMAT (/'  Output time data : '/                               &
       ' --------------------------------------------------'/ &
       '      First time         : ',a)
941 FORMAT ( '      Interval           : ',a/                       &
       '      Number of requests : ',i4)
  !
950 FORMAT (/'  Requested output for',i3,' points : '/              &
       ' --------------------------------------------------')
 
951 FORMAT ( '      ',a,2f10.2)
 
956 FORMAT ( '      ',a,2(f8.1,'E3'))
 
952 FORMAT (/'  Requested output fields :'/                         &
       ' --------------------------------------------------')
953 FORMAT ( '      ',a)
955 FORMAT (/'  Output times :'/                                    &
       ' --------------------------------------------------')
  !
960 FORMAT (//'  Final file management '/                           &
       ' -----------------------------------------------------')
961 FORMAT ( '      Closing file ww3.spec.grads'/                   &
       '      Closing file ww3.mean.grads')
962 FORMAT ( '      Preparing control files :')
963 FORMAT ( '         Number of times : ',i6/                      &
       '         Initial time ID : ',i2.2,':',i2.2,'Z',i2.2,a3,i4/ &
       '         Time step ID    : ',i2,a2)
964 FORMAT ( '      Writing ww3.spec.ctl'/)
  !
970 FORMAT ('DSET      ww3.spec.grads'/                             &
       'TITLE     WAVEWATCH III spectra and source terms'/     &
       'OPTIONS   sequential'/                                 &
       'OPTIONS   big_endian'/                                 &
       'UNDEF    ',e10.2/                                      &
       'XDEF     ',i4,'  LINEAR ',2f8.2/                       &
       'YDEF     ',i4,'  LEVELS ',5f8.4)
971 FORMAT (22x,5f8.4)
972 FORMAT ('ZDEF     ',i4,'  LINEAR ',2f8.2/                       &
       'TDEF     ',i4,'  LINEAR ',i6.2,':',i2.2,'Z',i2.2,a3,i4, &
       2x,i2,a2/                                              &
       'VARS     ',i4)
973 FORMAT ('LOC',i3.3,2i4,2x,a)
974 FORMAT ('ENDVARS')
  !
998 FORMAT (/'      End of file reached '/)
  !
999 FORMAT (/'  End of program '/                                   &
       ' ========================================='/          &
       '         WAVEWATCH III GrADS point output '/)
  !
1000 FORMAT (/' *** WAVEWATCH III ERROR IN GXOUTP : '/               &
       '     ERROR IN OPENING INPUT FILE'/                    &
       '     IOSTAT =',i5/)
  !
1001 FORMAT (/' *** WAVEWATCH III ERROR IN GXOUTP : '/               &
       '     PREMATURE END OF INPUT FILE'/)
  !
1002 FORMAT (/' *** WAVEWATCH III ERROR IN GXOUTP : '/               &
       '     ERROR IN READING FROM INPUT FILE'/               &
       '     IOSTAT =',i5/)
  !
1020 FORMAT (/' *** WAVEWATCH III ERROR IN GXOUTF : '/               &
       '     FIELD INCREMENT > 1HR BUT NOT MULTIPLE',f10.0/)
  !
1021 FORMAT (/' *** WAVEWATCH III ERROR IN GXOUTF : '/               &
       '     UPDATE PARS IN LOOP 610 !!!'/)
  !/
  !/ Internal subroutine GXEXPO ---------------------------------------- /
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
  SUBROUTINE gxexpo
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         16-Jul-2012 |
    !/                  +-----------------------------------+
    !/
    !/    30-Jun-1999 : Final FORTRAN 77                    ( version 1.18 )
    !/    24-Jan-2000 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/                  Massive changes to logistics
    !/    25-Jan-2001 : Cartesian grid version              ( version 2.06 )
    !/    02-Feb-2001 : Xnl version 5                       ( version 2.07 )
    !/    01-Aug-2003 : Fix format for SH output points.    ( version 3.03 )
    !/    24-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    23-Jun-2006 : Linear input added.                 ( version 3.09 )
    !/    03-Jul-2006 : Separate flux modules.              ( version 3.09 )
    !/    25-Jul-2006 : Grid ID for each point.             ( version 3.10 )
    !/    25-Apr-2007 : EMEAN in W3SPR2 par list.           ( version 3.11 )
    !/    09-Oct-2007 : WAM 4+ Sin and Sds added.           ( version 3.13 )
    !/                  (F. Ardhuin)
    !/    09-Oct-2007 : Experimental Sbs (BS1) added.       ( version 3.13 )
    !/                  (F. Ardhuin)
    !/    16-Jul-2012 : Move GMD (SNL3) and nonlinear filter (SNLS)
    !/                  from 3.15 (HLT).                    ( version 4.08 )
    !/    18-Aug-2018 : S_{ice} IC5 (Q. Liu)                ( version 6.06 )
    !/
    !  1. Purpose :
    !
    !     Perform actual point output.
    !
    !  3. Parameters :
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3SPRn    Subr. W3SRCnMD Mean wave parameters for use in
    !                               source terms.
    !      W3FLXn    Subr. W3FLXnMD Flux/stress computation.
    !      W3SLNn    Subr. W3SLNnMD Linear input.
    !      W3SINn    Subr. W3SRCnMD Input source term.
    !      W3SDSn    Subr. W3SRCnMD Whitecapping source term
    !      W3SNLn    Subr. W3SNLnMD Nonlinear interactions.
    !      W3SBTn    Subr. W3SBTnMD Bottom friction source term.
    !      W3SDBn    Subr. W3SBTnMD Depth induced breaking source term.
    !      W3STRn    Subr. W3STRnMD Triad interaction source term.
    !      W3SBSn    Subr. W3SBSnMD Bottom scattering source term.
    !      W3SXXn    Subr. W3SXXnMD Unclassified source term.
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !      STME21    Subr. W3TIMEMD Convert time to string.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !     Program in which it is contained.
    !
    !  6. Error messages :
    !
    !     None.
    !
    !  7. Remarks :
    !
    !     - Spectra are relative frequency energy spectra.
    !     - Note that arrays CX and CY of the main program now contain
    !       the absolute current speed and direction respectively.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !       !/S      Enable subroutine tracing.
    !       !/T      Enable test output.
    !
    !       !/FLXx   Flux/stress computation.
    !       !/LNx    Linear input package
    !       !/STx    Source term package
    !       !/NLx    Nonlinear interaction package
    !       !/BTx    Bottom friction package
    !       !/ICx    Ice source term package
    !       !/DBx    Depth-induced breaking package
    !       !/TRx    Triad interaction package
    !       !/BSx    Bottom scattering package
    !
    !       !/STAB2  Stability correction for !/ST2
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
#ifdef W3_FLX1
    USE w3flx1md
#endif
#ifdef W3_FLX2
    USE w3flx2md
#endif
#ifdef W3_FLX3
    USE w3flx3md
#endif
#ifdef W3_FLX4
    USE w3flx4md
#endif
#ifdef W3_FLX5
    USE w3flx5md
#endif
#ifdef W3_LN1
    USE w3sln1md
#endif
#ifdef W3_ST1
    USE w3src1md
#endif
#ifdef W3_ST2
    USE w3src2md
#endif
#ifdef W3_ST3
    USE w3src3md
#endif
#ifdef W3_ST4
    USE w3src4md, ONLY : w3spr4, w3sin4, w3sds4
#endif
#ifdef W3_ST6
    USE w3src6md
    USE w3swldmd, ONLY : w3swl6
    USE w3gdatmd, ONLY : swl6s6
#endif
#ifdef W3_NL1
    USE w3snl1md
#endif
#ifdef W3_NL2
    USE w3snl2md
#endif
#ifdef W3_NL3
    USE w3snl3md
#endif
#ifdef W3_NL4
    USE w3snl4md
#endif
#ifdef W3_NLS
    USE w3snlsmd
#endif
#ifdef W3_BT1
    USE w3sbt1md
#endif
#ifdef W3_BT4
    USE w3sbt4md
#endif
#ifdef W3_BT8
    USE w3sbt8md
#endif
#ifdef W3_IC1
    USE w3sic1md
#endif
#ifdef W3_IC2
    USE w3sic2md
#endif
#ifdef W3_IC3
    USE w3sic3md
#endif
#ifdef W3_IC4
    USE w3sic4md
#endif
#ifdef W3_IC5
    USE w3sic5md
#endif
#ifdef W3_DB1
    USE w3sdb1md
#endif
#ifdef W3_BS1
    USE w3sbs1md
#endif
#ifdef W3_IS2
    USE w3sis2md
#endif
    !/
    USE w3dispmd, ONLY: liu_forward_dispersion, nar1d, dfac, n1max, &
         ecg1, ewn1, dsie
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: J, I1, I2, IK, ITH, ISPEC, IKM, IKL, &
         IKH, ITT, IX, IY, ISEA
#ifdef W3_S
    INTEGER, SAVE           :: IENT   = 0
#endif
    REAL                    :: XL, XH, XL2, XH2, DEPTH, SQRTH, UDIR,&
         UDIRR, UABS, CDIR, SIX, R1, R2, ET,  &
         EWN, ETR, ETX, ETY, EBND, EBX, EBY,  &
         HSIG, WLEN, TMEAN, THMEAN, THSPRD,   &
         EMAX, EL, EH, DENOM, FP, THP, SPP,   &
         FACTOR, CD, USTAR, FHIGH, ZWND, ICE, &
         USTD, Z0, CHARN, EMEAN, FMEAN, WNMEAN,&
         ICETHICK, ICECON, ICEF
#ifdef W3_FLX5
    REAL                     ::TAUA, TAUADIR, RHOAIR
#endif
#ifdef W3_IS2
    REAL                   :: ICEDMAX
#endif
#ifdef W3_ST1
    REAL                    :: AMAX, FH1, FH2
#endif
#ifdef W3_ST2
    REAL                    :: AMAX, ALPHA(NK), FPI
#endif
#ifdef W3_ST3
    REAL                    :: FMEANS, FMEANWS, TAUWX, TAUWY, AMAX, &
         TAUWNX, TAUWNY
#endif
#ifdef W3_ST4
    REAL                    :: FMEANWS, TAUWX, TAUWY, AMAX, &
         TAUWNX, TAUWNY,  FMEAN1, WHITECAP(1:4), DLWMEAN
#endif
#ifdef W3_ST6
    REAL                    :: AMAX, TAUWX, TAUWY, TAUWNX, TAUWNY
#endif
#ifdef W3_BS1
    REAL                    :: TAUSCX, TAUSCY
#endif
#ifdef W3_BT3
    REAL                    :: D50
#endif
#ifdef W3_BT4
    REAL                    :: D50, PSIC, BEDFORM(3), TAUBBL(2)
#endif
#ifdef W3_STAB2
    REAL                    :: STAB0, STAB, THARG1, THARG2, COR1,   &
         COR2, ASFAC
#endif
    REAL                    :: HSMIN = 0.05
    REAL                    :: WN(NK), CG(NK), E(NK,NTH), E1(NK),   &
         APM(NK), THBND(NK), SPBND(NK),       &
         A(NTH,NK), WN2(NTH,NK),WN_R(NK),     &
         ALPHA_LIU(NK), CG_ICE(NK), R(NK)
    REAL                    :: DIA(NTH,NK), SWI(NK,NTH), SNL(NK,NTH),&
         SDS(NK,NTH), SBT(NK,NTH), SIS(NK,NTH),&
         STT(NK,NTH), DIA2(NK,NTH)
    REAL                    :: XLN(NTH,NK), XWI(NTH,NK), XNL(NTH,NK),&
         XTR(NTH,NK), XDS(NTH,NK), XDB(NTH,NK),&
         XBT(NTH,NK), XBS(NTH,NK), XXX(NTH,NK),&
         XWL(NTH,NK), XIS(NTH,NK)
    LOGICAL                 :: LBREAK
#ifdef W3_ST3
    LOGICAL                :: LLWS(NTH,NK)
#endif
#ifdef W3_ST4
    LOGICAL                :: LLWS(NTH,NK)
    REAL                   :: LAMBDA(NSPEC)
#endif
    CHARACTER               :: DTME21*23
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'GXEXPO')
#endif
    !
    xl     = 1./xfr - 1.
    xh     =  xfr - 1.
    xl2    = xl**2
    xh2    = xh**2
    ice = 0.
    !
    xln = 0.
    xwi = 0.
    xnl = 0.
    xtr = 0.
    xds = 0.
    xdb = 0.
    xbt = 0.
    xbs = 0.
    xwl = 0.
    xis = 0.
    xxx = 0.
    !
#ifdef W3_T
    WRITE (ndst,9000) (flreq(j),j=1,nopts)
    WRITE (ndst,9001) flsrce
#endif
    !
    !     Output of time
    !
    CALL stme21 ( time , dtme21 )
    WRITE (ndso,905) dtme21
    !
    !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !     Loop over output points.
    !
    DO j=1, nopts
      IF ( flreq(j) ) THEN
        !
#ifdef W3_T
        WRITE (ndst,9002) ptnme(j)
#endif
        !
        ! 2. Calculate grid parameters using and inlined version of WAVNU1.
        !
        depth  = max( dmin, dpo(j) )
        sqrth  = sqrt( depth )
        udir   = mod( 270. - wdo(j)*rade , 360. )
        udirr  = wdo(j)
        uabs   = max( 0.001 , wao(j) )
#ifdef W3_FLX5
        taua     = max( 0.001 , tauao(j))
        tauadir  = mod( 270. - taudo(j)*rade , 360. )
        rhoair   = max( 0. , dairo(j))
#endif
        cdir   = mod( 270. - cdo(j)*rade , 360. )
#ifdef W3_IS2
        icedmax = max( 0., icefo(j))
#endif
#ifdef W3_IC2
        icef     = 0.
#endif
#ifdef W3_IS2
        icef     = icedmax
#endif
        icethick = max(0., iceho(j))
        icecon = max(0., iceo(j))
        !
#ifdef W3_STAB2
        stab0  = zwind * grav / 273.
        stab   = stab0 * aso(j) / max(5.,wao(j))**2
        stab   = max( -1. , min( 1. , stab ) )
        tharg1 = max( 0. , ffng*(stab-ofstab))
        tharg2 = max( 0. , ffps*(stab-ofstab))
        cor1   = ccng * tanh(tharg1)
        cor2   = ccps * tanh(tharg2)
        asfac  = sqrt( (1.+cor1+cor2)/shstab )
#endif
        !
#ifdef W3_T
        WRITE (ndst,9010) depth
#endif
        DO ik=1, nk
          six    = sig(ik) * sqrth
          i1     = int(six/dsie)
          IF (i1.LE.n1max) THEN
            i2 = i1 + 1
            r1 = six/dsie - real(i1)
            r2 = 1. - r1
            wn(ik) = ( r2*ewn1(i1) + r1*ewn1(i2) ) / depth
            cg(ik) = ( r2*ecg1(i1) + r1*ecg1(i2) ) * sqrth
          ELSE
            wn(ik) = sig(ik)*sig(ik)/grav
            cg(ik) = 0.5 * grav / sig(ik)
          END IF
#ifdef W3_T
          WRITE (ndst,9011) ik, tpi/sig(ik), wn(ik), cg(ik)
#endif
          !
        END DO
 
        IF (iicedisp) THEN
          CALL liu_forward_dispersion  (icethick,0.,depth, &
               sig,wn_r,cg_ice,alpha_liu)
        ELSE
          wn_r=wn
          cg_ice=cg
        END IF
        r(:)=1 ! In case IC2 is defined but not IS2
 
        !
        ! 3.  Prepare spectra etc.
        ! 3.a Mean wave parameters.
        !
        et     = 0.
        ewn    = 0.
        etr    = 0.
        etx    = 0.
        ety    = 0.
        DO ik=1, nk
          ebnd   = 0.
          ebx    = 0.
          eby    = 0.
          DO ith=1, nth
            ispec  = ith + (ik-1)*nth
            e(ik,ith) = spco(ispec,j)
            ebnd   = ebnd + spco(ispec,j)
            ebx    = ebx  + spco(ispec,j)*ecos(ith)
            eby    = eby  + spco(ispec,j)*esin(ith)
          END DO
          e1(ik) = ebnd * dth
          apm(ik)= e1(ik) / ( tpi * grav**2 / sig(ik)**5  )
          IF ( e1(ik) .GT. 1.e-5) THEN
            thbnd(ik) = mod(630.- rade*atan2(eby,ebx),360.)
            spbnd(ik) = rade * sqrt( max( 0. , 2.*( 1. -      &
                 sqrt( max(0.,(ebx**2+eby**2)/ebnd**2) ) ) ) )
          ELSE
            thbnd(ik) = -999.9
            spbnd(ik) = -999.9
          END IF
          ebnd   = e1(ik) * dsii(ik) * tpiinv
          et     = et  + ebnd
          ewn    = ewn + ebnd / wn(ik)
          etr    = etr + ebnd / sig(ik)
          etx    = etx + ebx * dsii(ik)
          ety    = ety + eby * dsii(ik)
        END DO
        !
        ! tail factors for radian action etc ...!
        !
        ebnd   = e1(nk) * tpiinv / ( sig(nk) * dth )
        et     = et  + fte *ebnd
        ewn    = ewn + ftwl*ebnd
        etr    = etr + fttr*ebnd
        etx    = dth*etx*tpiinv + fte*ebx*tpiinv/sig(nk)
        ety    = dth*ety*tpiinv + fte*eby*tpiinv/sig(nk)
        !
        hsig   = 4. * sqrt( et )
        IF ( hsig .GT. hsmin ) THEN
          wlen   = ewn / et * tpi
          tmean  = etr / et * tpi
          thmean = mod( 630. - rade*atan2(ety,etx) , 360. )
          thsprd = rade * sqrt( max( 0. , 2.*( 1. - sqrt(     &
               max(0.,(etx**2+ety**2)/et**2) ) ) ) )
        ELSE
          wlen   = 0.
          tmean  = 0.
          thmean = 0.
          thsprd = 0.
          DO ik=1, nk
            e1(ik) = 0.
            DO ith=1, nth
              e(ik,ith) = 0.
            END DO
          END DO
        END IF
        !
        ! peak frequency
        !
        emax   = e1(nk)
        ikm    = nk
        !
        DO ik=nk-1, 1, -1
          IF ( e1(ik) .GT. emax ) THEN
            emax   = e1(ik)
            ikm    = ik
          END IF
        END DO
        !
        ikl    = max(  1 , ikm-1 )
        ikh    = min( nk , ikm+1 )
        el     = e1(ikl) - e1(ikm)
        eh     = e1(ikh) - e1(ikm)
        denom  = xl*eh - xh*el
        !
        IF ( hsig .GE. hsmin ) THEN
          fp     = sig(ikm) * ( 1. + 0.5 * ( xl2*eh - xh2*el )  &
               / sign( max(abs(denom),1.e-15) , denom ) )
          thp    = thbnd(ikm)
          spp    = spbnd(ikm)
        ELSE
          fp     = 0.
          thp    = 0.
          spp    = 0.
        END IF
        !
        ! 3.4 source terms
        !
        DO ik=1, nk
          factor = tpiinv * cg(ik) / sig(ik)
          DO ith=1, nth
            ispec  = ith + (ik-1)*nth
            a(ith,ik)   = factor * spco(ispec,j)
            wn2(ith,ik) = wn(ik)
          END DO
        END DO
        !
#ifdef W3_STAB2
        uabs   = uabs / asfac
#endif
        !
#ifdef W3_ST0
        zwnd   = 10.
#endif
#ifdef W3_ST1
        zwnd   = 10.
#endif
#ifdef W3_ST2
        zwnd   = zwind
#endif
#ifdef W3_ST3
        zwnd   = zzwnd
        tauwx  = 0.
        tauwy  = 0.
        llws(:,:)  = .true.
#endif
        ustar  = 1.
#ifdef W3_ST4
        zwnd   = zzwnd
        tauwx  = 0.
        tauwy  = 0.
#endif
#ifdef W3_ST6
        zwnd   = 10.
#endif
        !
#ifdef W3_ST0
        fhigh  = sig(nk)
#endif
#ifdef W3_ST1
        CALL w3spr1 (a, cg, wn, emean, fmean, wnmean, amax)
        fp     = 0.85 * fmean
        fh1    = fxfm * fmean
        fh2    = fxpm / ustar
        fhigh  = max( fh1 , fh2 )
#endif
#ifdef W3_ST2
        CALL w3spr2 (a, cg, wn, depth, fp , uabs, ustar,     &
             emean, fmean, wnmean, amax, alpha, fp )
#endif
#ifdef W3_ST3
        CALL w3spr3 (a, cg, wn, emean, fmean, fmeans,        &
             wnmean, amax, uabs, udirr, ustar, ustd, &
             tauwx, tauwy, cd, z0, charn, llws, fmeanws)
#endif
#ifdef W3_ST4
        CALL w3spr4 (a, cg, wn, emean, fmean,  fmean1,      &
             wnmean, amax, uabs, udirr,             &
#ifdef W3_FLX5
             taua, tauadir, rhoair,           &
#endif
             ustar, ustd, tauwx, tauwy, cd, z0,     &
             charn, llws, fmeanws, dlwmean )
#endif
#ifdef W3_ST6
        CALL w3spr6 (a, cg, wn, emean, fmean, wnmean, amax, fp)
#endif
        !
#ifdef W3_FLX1
        CALL w3flx1 ( zwnd, uabs, udirr,                     &
             ustar, ustd, z0, cd )
#endif
#ifdef W3_FLX2
        CALL w3flx2 ( zwnd, depth, fp, uabs, udirr,          &
             ustar, ustd, z0, cd )
#endif
#ifdef W3_FLX3
        CALL w3flx3 ( zwnd, depth, fp, uabs, udirr,          &
             ustar, ustd, z0, cd )
#endif
#ifdef W3_FLX4
        CALL w3flx4 ( zwnd, uabs, udirr, ustar, ustd, z0, cd )
#endif
#ifdef W3_FLX5
        CALL w3flx5 ( zwnd, uabs, udirr, taua, tauadir,  &
             rhoair, ustar, ustd, z0, cd, charn )
#endif
        !
        DO itt=1, 3
#ifdef W3_ST2
          CALL w3sin2 (a, cg, wn2, uabs, udirr, cd, z0,      &
               fpi, xwi, dia )
          CALL w3spr2 (a, cg, wn, depth, fpi, uabs, ustar,   &
               emean, fmean, wnmean, amax, alpha, fp )
#endif
#ifdef W3_ST3
          CALL w3sin3 (a, cg, wn2, uabs, ustar, dair/dwat,   &
               aso(j), udirr, z0, cd, tauwx, tauwy,  &
               tauwnx, tauwny,                       &
               ice, xwi, dia, llws, ix, iy )
          CALL w3spr3 (a, cg, wn, emean, fmean, fmeans,        &
               wnmean, amax, uabs, udirr, ustar, ustd, &
               tauwx, tauwy, cd, z0, charn, llws, fmeanws)
#endif
#ifdef W3_ST4
          CALL w3sin4 (a, cg, wn2, uabs, ustar, dair/dwat,   &
               aso(j), udirr, z0, cd, tauwx, tauwy,  &
               tauwnx, tauwny, xwi, dia, llws, ix, iy, lambda )
          CALL w3spr4 (a, cg, wn, emean, fmean,  fmean1,      &
               wnmean, amax, uabs, udirr,               &
#ifdef W3_FLX5
               taua, tauadir, rhoair,             &
#endif
               ustar, ustd, tauwx, tauwy, cd, z0,       &
               charn, llws, fmeanws, dlwmean )
#endif
#ifdef W3_FLX2
          CALL w3flx2 ( zwnd, depth, fp, uabs, udirr,        &
               ustar, ustd, z0, cd )
#endif
#ifdef W3_FLX3
          CALL w3flx3 ( zwnd, depth, fp, uabs, udirr,        &
               ustar, ustd, z0, cd )
#endif
        END DO
        !
#ifdef W3_ST2
        fhigh  = xfc * fpi
#endif
        !
        IF ( flsrce(2) ) THEN
#ifdef W3_LN1
          CALL w3sln1 ( wn, fhigh, ustar, udirr, xln )
#endif
          !
#ifdef W3_ST1
          CALL w3sin1 (a, wn2, ustar, udirr,      xwi, dia )
#endif
#ifdef W3_ST2
          CALL w3sin2 (a, cg, wn2, uabs, udirr, cd, z0,    &
               fpi, xwi, dia )
#endif
#ifdef W3_ST3
          CALL w3sin3 (a, cg, wn2, uabs, ustar, dair/dwat, &
               aso(j), udirr, z0, cd,              &
               tauwx, tauwy, tauwnx, tauwny,       &
               ice, xwi, dia, llws, ix, iy )
#endif
#ifdef W3_ST4
          CALL w3sin4 (a, cg, wn2, uabs, ustar, dair/dwat, &
               aso(j), udirr, z0, cd,              &
               tauwx, tauwy, tauwnx, tauwny,       &
               xwi, dia,  llws, ix, iy, lambda )
#endif
#ifdef W3_ST6
          CALL w3sin6 (a, cg, wn2, uabs, ustar, udirr, cd, &
               dair, tauwx, tauwy, tauwnx, tauwny, xwi, dia )
#endif
        END IF
        IF ( flsrce(3) ) THEN
#ifdef W3_NL1
          CALL w3snl1 ( a, cg, wnmean*depth,      xnl, dia )
#endif
#ifdef W3_NL2
          CALL w3snl2 ( a, cg, depth,             xnl, dia )
#endif
#ifdef W3_NL3
          CALL w3snl3 ( a, cg, wn, depth,         xnl, dia )
#endif
#ifdef W3_NL4
          CALL w3snl4 ( a, cg, wn, depth,         xnl, dia )
#endif
        END IF
        IF ( flsrce(4) ) THEN
#ifdef W3_ST1
          CALL w3sds1 ( a, wn2, emean, fmean, wnmean, xds, dia )
#endif
#ifdef W3_ST2
          CALL w3sds2 ( a, cg, wn, fpi, ustar, alpha, xds, dia )
#endif
#ifdef W3_ST3
          CALL w3sds3 ( a, wn, cg, emean, fmeans, wnmean,         &
               ustar, ustd, depth, xds, dia, ix, iy )
#endif
#ifdef W3_ST4
          CALL w3sds4 ( a, wn, cg,          &
               ustar, ustd, depth, dair, xds, dia, ix, iy, lambda, whitecap , dlwmean)
#endif
#ifdef W3_ST6
          CALL w3sds6 ( a, cg, wn,            xds, dia )
          IF (swl6s6) CALL w3swl6 ( a, cg, wn,  xwl, dia )
#endif
          !
#ifdef W3_DB1
          CALL w3sdb1 ( j, a, depth, emean, fmean, wnmean, cg, &
               lbreak, xdb, dia )
#endif
          !
        END IF
        IF ( flsrce(5) ) THEN
 
#ifdef W3_BT1
          CALL w3sbt1 ( a, cg, wn, depth,   xbt, dia )
#endif
 
#ifdef W3_IC1
          CALL w3sic1 ( a, depth, cg,      ix, iy, xbt, dia )
#endif
#ifdef W3_IC2
          CALL w3sic2 ( a, depth, icethick, icef ,cg,  wn, ix, iy, xbt, dia, wn_r, &
               cg_ice, alpha_liu, r )
#endif
#ifdef W3_IC3
          CALL w3sic3 ( a, depth, cg,  wn, ix, iy, xbt, dia )
#endif
#ifdef W3_IC4
          CALL w3sic4 ( a, depth, cg,      ix, iy, xbt, dia )
#endif
#ifdef W3_IC5
          CALL w3sic5 ( a, depth, cg,  wn, ix, iy, xbt, dia )
#endif
 
#ifdef W3_BT4
          ix=1    ! to be fixed later
          iy=1    ! to be fixed later
          isea=1  ! to be fixed later
          d50 = sed_d50(isea)
          psic= sed_psic(isea)
#endif
 
#ifdef W3_BT4
          CALL w3sbt4 ( a, cg, wn, depth, d50, psic, taubbl,   &
               bedform, xbt, dia, ix, iy )
#endif
          !
 
#ifdef W3_BT8
          CALL w3sbt8 ( a, depth, xbt, dia, ix, iy )
#endif
 
#ifdef W3_BS1
          CALL w3sbs1 ( a, cg, wn, depth, cao(j)*cos(cdo(j)), &
               cao(j)*sin(cdo(j)),     &
               tauscx, tauscy, xbs, dia )
#endif
        END IF
 
        IF ( flsrce(6) ) THEN
 
#ifdef W3_IS2
          CALL w3sis2(a, depth, icecon, icethick, icef, icedmax, ix, iy, &
               xis, dia, dia2, wn, cg, wn_r, cg_ice, r)
#endif
        END IF
        !
#ifdef W3_STAB2
        uabs   = uabs * asfac
#endif
        !
        DO ik=1, nk
          factor = tpi / cg(ik) * sig(ik)
          DO ith=1, nth
            ispec       = ith + (ik-1)*nth
            e(ik,ith) = spco(ispec,j)
            swi(ik,ith) = ( xwi(ith,ik) + xln(ith,ik) ) * factor
            snl(ik,ith) = ( xnl(ith,ik) + xtr(ith,ik) ) * factor
            sds(ik,ith) = ( xds(ith,ik) + xdb(ith,ik) ) * factor
#ifdef W3_ST6
            sds(ik,ith) =   sds(ik,ith) +(xwl(ith,ik)   * factor)
#endif
            sbt(ik,ith) = ( xbt(ith,ik) + xbs(ith,ik) ) * factor
            sis(ik,ith) = xis(ith,ik) * factor
            stt(ik,ith) = xxx(ith,ik) * factor
          END DO
        END DO
        stt    = stt + swi + snl + sds + sbt + sis
 
        !
        ! 4.a Perform output
        !
        IF ( flsrce(1) ) WRITE (ndsgrd)                           &
             ((e(ik,ith),ith=1,nth),ik=nk,1,-1)
        IF ( flsrce(2) ) WRITE (ndsgrd)                           &
             ((swi(ik,ith),ith=1,nth),ik=nk,1,-1)
        IF ( flsrce(3) ) WRITE (ndsgrd)                           &
             ((snl(ik,ith),ith=1,nth),ik=nk,1,-1)
        IF ( flsrce(4) ) WRITE (ndsgrd)                           &
             ((sds(ik,ith),ith=1,nth),ik=nk,1,-1)
        IF ( flsrce(5) ) WRITE (ndsgrd)                           &
             ((sbt(ik,ith),ith=1,nth),ik=nk,1,-1)
        IF ( flsrce(6) ) WRITE (ndsgrd)                           &
             ((sis(ik,ith),ith=1,nth),ik=nk,1,-1)
        IF ( flsrce(7) ) WRITE (ndsgrd)                           &
             ((stt(ik,ith),ith=1,nth),ik=nk,1,-1)
        !
        IF ( flagll ) THEN
          WRITE (ndspnt,940) ptnme(j),                          &
               fact*ptloc(1,j), fact*ptloc(2,j), dpo(j), wao(j),  &
               wao(j)*cos(wdo(j)), wao(j)*sin(wdo(j)), aso(j),    &
               cao(j), cao(j)*cos(cdo(j)), cao(j)*sin(cdo(j)),    &
               hsig, grdid(j)
        ELSE
          WRITE (ndspnt,941) ptnme(j),                          &
               fact*ptloc(1,j), fact*ptloc(2,j), dpo(j), wao(j),  &
               wao(j)*cos(wdo(j)), wao(j)*sin(wdo(j)), aso(j),    &
               cao(j), cao(j)*cos(cdo(j)), cao(j)*sin(cdo(j)),    &
               hsig, grdid(j)
        END IF
        !
        ! ... End of points loop
        !
      END IF
    END DO
    !
    RETURN
    !
    ! Formats
    !
905 FORMAT (9x,a)
 
940 FORMAT (a10,1x,2f6.1,f7.1,3f7.1,f8.2,3f7.2,f6.2,2x,a)
 
941 FORMAT (a10,1x,2f8.1,f7.1,3f7.1,f8.2,3f7.2,f6.2,2x,a)
 
    !
#ifdef W3_T
9000 FORMAT (' TEST GXEXPO : FLAGS  :',40l2)
9001 FORMAT (' TEST GXEXPO : FLSRCE :',6l2)
9002 FORMAT (' TEST GXEXPO : OUTPUT POINT : ',a)
9010 FORMAT (' TEST GXEXPO : DEPTH =',f7.1,'  IK, T, K, CG :')
9011 FORMAT ('               ',i3,f8.2,f8.4,f8.2)
#endif
    !/
    !/ End of GXEXPO ----------------------------------------------------- /
    !/
  END SUBROUTINE gxexpo
  !/
  !/ End of GXOUTP ----------------------------------------------------- /
  !/
END PROGRAM gxoutp