ww3_strt.F90

Go to the documentation of this file.

 
!
#include "w3macros.h"
 
!/ ------------------------------------------------------------------- /
PROGRAM w3strt
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         06-Jun-2018 |
  !/                  +-----------------------------------+
  !/
  !/    15-Jan-1999 : Final FORTRAN 77                    ( version 1.18 )
  !/    18-Jan-2000 : Upgrade to FORTRAN 90               ( version 2.00 )
  !/    11-Jan-2001 : Flat grid version                   ( version 2.06 )
  !/    11-Jun-2001 : Clean up.                           ( version 2.11 )
  !/    30-Apr-2002 : Updated W3IORS.                     ( version 2.20 )
  !/    13-Nov-2002 : Updated W3IORS.                     ( version 3.00 )
  !/    24-Dec-2004 : Multiple grid version.              ( version 3.06 )
  !/    28-Jun-2006 : Adding file name preamble.          ( version 3.09 )
  !/    08-May-2007 : Starting from calm as an option.    ( version 3.11 )
  !/    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)
  !/    31-Oct-2010 : Implement unstructured grid         ( version 3.14 )
  !/                  (A. Roland and F. Ardhuin)
  !/    05-Jul-2011 : Revert to X-Y gaussian shape        ( version 4.01 )
  !/    06-Mar-2012 : Hardening output.                   ( version 4.07 )
  !/    06-Jun-2018 : Add DEBUGINIT/EXPORTWWM             ( version 6.04 )
  !/
  !/
  !/    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 :
  !
  !     Generation of initial conditions for a "cold start" of
  !     WAVEWATCH III.
  !
  !  2. Method :
  !
  !     General model information is obtained from the model definition
  !     file using W3IOGR. The type of the initial field is read
  !     from the input file WW3_strt.inp (NDSI). Three types of initial
  !     conditions are available.
  !       1) Gaussian distribution in longitude, latitude and frequency,
  !          cos power in directions. Can default to single spectral
  !          bin.
  !       2) Predefined JONSWAP spectrum, Gaussian height distribution
  !          in space.
  !       3) Fetch-limited JONSWAP spectrum based on the actual wind
  !          speed. To avoid the need of reading a wind field, the
  !          restart file is a "dummy", and the actual initial field
  !          is constructed in the initialization routine W3INIT.
  !       4) User defined spectrum throughout the model.
  !       5) Starting from rest.
  !     The initial conditions are written to the restart.WW3 using the
  !     subroutine W3IORS. Note that the name of the restart file is set
  !     in W3IORS.
  !
  !  3. Parameters :
  !
  !     Local parameters.
  !     ----------------------------------------------------------------
  !       NDSI    Int.  Input unit number ("ww3_strt.inp").
  !       ITYPE   Int.  Type of field (see section 2).
  !       FP,SIP  Real  Peak frequency (Hz) and spread.  \
  !       XM,SIX  Real  Id. X (degr.).                   |
  !       YM,SIY  Real  Id. Y (degr.).                   |   ITYPE = 1
  !       HMAX    Real  Maximum wave height.             |
  !       NCOS    Real  Cosine power in dir. distr.      |
  !       THM     Real  Mean direction (cart. degr.)     / \
  !       ALFA    Real  Energy level of PM spectrum.       |
  !       FP      Real  Peak frequency (Hz).               | ITYPE = 2
  !       GAMMA   Real  Peak enhancement factor            |
  !       SIGA/B  Real  Spread with GAMA.                  /
  !     ----------------------------------------------------------------
  !
  !  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.
  !      W3DIMW    Subr.   Id.    Set array dims 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.
  !      EJ5P      Func.   Id.    Five parameter JONSWAP spectrum.
  !      PRT1DS    Subr. W3ARRYMD Print plot of 1-D spectrum.
  !      PRT2DS    Subr.   Id.    Print plot of 2-D spectrum.
  !      PRTBLK    Subr.   Id.    Print plot of array.
  !      WAVNU1    Subr. W3DISPMD Solve dispersion relation.
  !      W3IOGR    Subr. W3IOGRMD Reading/writing model definition file.
  !      W3IORS    Subr. W3IORSMD Reading/writing restart files.
  !      W3DIST    Subr. W3GSRUMD Compute distance between two points.
  !      MPI_xxx   Subr. mpif.h   Standard MPI routines.
  !     ----------------------------------------------------------------
  !
  !  5. Called by :
  !
  !     None, stand-alone program.
  !
  !  6. Error messages :
  !
  !  7. Remarks :
  !
  !     - While reading the restart file W3IORS will recognize the
  !       need for checking the time, as the restart file contains
  !       information on the origine of the file ("cold" or "hot").
  !     - User input for x-wise gaussian spread control, SIX, is
  !       now available again (option for SIX.NE.SIY available.)
  !       If user desires a distribution that is circular in real
  !       distances, user should input a negative number for SIX.
  !
  !  8. Structure :
  !
  !     ----------------------------------------------------
  !        1.a  Set up data structures.
  !                            ( W3NMOD , W3NDAT , W3NOUT
  !                              W3SETG , W3SETW , W3SETO )
  !          b  I-O setup.
  !          b  Print heading(s).
  !        2.a  Read model defintion file with base model
  !             data.                            ( W3IOGR )
  !          b  MPP initializations.
  !        3.   Get field type from the input file.
  !        4.   ITYPE = 1, Gaussian, cosine.
  !          a  Read parameters.
  !          b  Set-up 1-D spectrum.
  !          c  Set-up directional distribution.
  !          d  Normalize spectrum with Hmax.
  !          e  Distribute over grid.
  !        5.   ITYPE = 2, pre-defined JONSWAP.
  !          a  Read parameters.
  !          b  Set-up 1-D spectrum.
  !          c  2-D energy spectrum.
  !          d  Distribute over grid.
  !        6.   ITYPE = 3, fetch limited JONSWAP.
  !        7.   ITYPE = 4, user-defined spectrum.
  !          a  Read scale factor.
  !          b  Read and rescale spectrum.
  !          c  Distribute over grid.
  !        8.   ITYPE = 5, start from calm conditions.
  !        9.   Convert energy to action
  !       10.   Write restart file.              ( W3IORS )
  !     ----------------------------------------------------
  !
  !  9. Switches :
  !
  !     !/SHRD  Switch for shared / distributed memory architecture.
  !     !/DIST  Id.
  !
  !     !/SHRD  Switch for message passing method.
  !     !/MPI   Id.
  !
  !     !/S     Enable subroutine tracing.
  !
  !     !/O4    Output normalized 1-D energy spectrum.
  !     !/O5    Output normalized 2-D energy spectrum.
  !     !/O6    Output normalized wave heights (not MPP adapted).
  !
  ! 10. Source code :
  !
  !/ ------------------------------------------------------------------- /
  USE constants
  !/
  !     USE W3GDATMD, ONLY: W3NMOD, W3SETG
  !     USE W3WDATMD, ONLY: W3NDAT, W3SETW, W3DIMW
#ifdef W3_NL1
  USE w3adatmd, ONLY: w3naux, w3seta
#endif
  USE w3odatmd, ONLY: w3nout, w3seto, flogrr
  USE w3servmd, ONLY: itrace, nextln, ej5p, extcde
#ifdef W3_S
  USE w3servmd, ONLY : strace
#endif
#ifdef W3_O4
  USE w3arrymd, ONLY : prt1ds
#endif
#ifdef W3_O5
  USE w3arrymd, ONLY : prt2ds
#endif
#ifdef W3_O6
  USE w3arrymd, ONLY : prtblk
#endif
  USE w3dispmd, ONLY : wavnu1
  USE w3iogrmd, ONLY: w3iogr
  USE w3iorsmd, ONLY: w3iors
  USE w3gsrumd, ONLY: w3dist
  !/
  USE w3gdatmd
  USE w3wdatmd
  USE w3odatmd, ONLY: ndse, ndst, ndso, naproc, iaproc,           &
       napout, naperr, fnmpre
#ifdef W3_WRST
  USE w3idatmd, ONLY: w3ninp
#endif
  !/
  IMPLICIT NONE
  !
#ifdef W3_MPI
  include "mpif.h"
#endif
  !/
  !/ ------------------------------------------------------------------- /
  !/ Local parameters
  !/
  INTEGER                 :: ndsi, ndsm, ndsr, ndstrc, ntrace,    &
       ndsen, ierr, itype, ncos, ikm, ik,   &
       ithm, ith, jsea, isea, ix, iy, j
#ifdef W3_MPI
  INTEGER                 :: ierr_mpi
#endif
#ifdef W3_S
  INTEGER, SAVE           :: ient = 0
#endif
#ifdef W3_O6
  INTEGER                 :: nsx, nsy
  INTEGER, ALLOCATABLE    :: mapo(:,:)
#endif
  REAL                    :: fp, sip, thm, xm, six, ym, siy, hmax,&
       chsip, frrel, etot, e1i, factor, x,  &
       y, rdsqr, alfa, gamma, siga, sigb,   &
       yln, fr, beta, frr, s, sumd, ang,    &
       arg, facs, depth, wn, cg, hpqmax
  REAL, ALLOCATABLE       :: e1(:), dd(:), e2(:,:), e21(:), finp(:,:)
#ifdef W3_O5
  REAL, ALLOCATABLE       :: e2out(:,:)
#endif
#ifdef W3_O6
  REAL, ALLOCATABLE       :: hsig(:,:)
#endif
  CHARACTER               :: comstr*1, inxout*4
#ifdef W3_EXPORTWWM
  INTEGER                 :: ispec
#endif
  LOGICAL                 :: flone,nosix
  !/
  !/ ------------------------------------------------------------------- /
  !
  ! 1.a Initialize data structure
  !
  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 )
#ifdef W3_WRST
  CALL w3ninp(    6, 6 )
#endif
  !
  ! 1.b IO set-up.
  !
  ndsi   = 10
  ndsm   = 20
  ndsr   = 20
  !
  flogrr(:,:) = .false.
  !
  ndstrc =  6
  ntrace = 10
  CALL itrace ( ndstrc, ntrace )
  !
#ifdef W3_S
  CALL strace (ient, 'W3STRT')
#endif
  !
  ! 1.c MPP initializations
  !
#ifdef W3_SHRD
  naproc = 1
  iaproc = 1
#endif
  !
#ifdef W3_MPI
  CALL mpi_init      ( ierr_mpi )
  CALL mpi_comm_size ( mpi_comm_world, naproc, ierr_mpi )
  CALL mpi_comm_rank ( mpi_comm_world, iaproc, ierr_mpi )
  iaproc = iaproc + 1
#endif
  !
  IF ( iaproc .EQ. naperr ) THEN
    ndsen  = ndse
  ELSE
    ndsen  = -1
  END IF
  !
  IF ( iaproc .EQ. napout ) WRITE (ndso,900)
  !
  j      = len_trim(fnmpre)
  OPEN (ndsi,file=fnmpre(:j)//'ww3_strt.inp',status='OLD',        &
       err=800,iostat=ierr)
  rewind(ndsi)
  READ (ndsi,'(A)',END=801,ERR=802) comstr
  IF (comstr.EQ.' ') comstr = '$'
  IF ( iaproc .EQ. napout ) WRITE (ndso,901) comstr
  !
  !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  ! 2.  Read model definition file and mpp initializations.
  ! 2.a Reading file
  !
  CALL w3iogr ( 'READ', ndsm )
  !
  IF ( iaproc .EQ. napout ) WRITE (ndso,902) gname
  !
  ! 2.b MPP initializations
  !
#ifdef W3_SHRD
  nseal  = nsea
#endif
  !
#ifdef W3_DIST
  nseal  = 1 + (nsea-iaproc)/naproc
  IF ( nsea .LT. naproc ) GOTO 803
#endif
  !
  CALL w3dimw ( 1, ndse, ndst )
  ALLOCATE ( e1(nk), dd(nth), e2(nth,nk), e21(nspec),           &
       finp(nk,nth) )
  !
  !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  ! 3.  Read type from input file.
  !
  CALL nextln ( comstr , ndsi , ndsen )
  READ (ndsi,*,END=801,ERR=802) itype
  IF ( itype.LT.1 .OR. itype.GT.5 ) THEN
    IF ( iaproc .EQ. naperr ) WRITE (ndse,1010) itype
    CALL extcde ( 1 )
  END IF
  IF ( iaproc .EQ. napout ) WRITE (ndso,930) itype
  !
  !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  ! 4.  ITYPE = 1, Gaussian, cosine.
  !
  IF ( itype .EQ. 1 ) THEN
    inxout = 'COLD'
    !
    ! 4.a Read parameters.
    !
    CALL nextln ( comstr , ndsi , ndsen )
    READ (ndsi,*,END=801,ERR=802)                            &
         fp, sip, thm, ncos, xm, six, ym, siy, hmax
    fp     = max( 0.5 * tpiinv * sig(1) , fp )
    sip    = max( 0. , sip )
    DO
      IF ( thm .LT. 0. ) THEN
        thm    = thm + 360.
      ELSE
        EXIT
      END IF
    END DO
    thm    = mod( thm , 360. )
    ncos   = max( 0 , 2*(ncos/2) )
 
    nosix=.false.
    IF(six.LT.0.0)THEN
      IF ( iaproc .EQ. napout ) WRITE (ndso,903)
      nosix=.true.
    END IF
 
    hpqmax=-999.0
    DO jsea=1, nseal
#ifdef W3_DIST
      isea   = iaproc + (jsea-1)*naproc
#endif
#ifdef W3_SHRD
      isea   = jsea
#endif
      ix     = mapsf(isea,1)
      iy     = mapsf(isea,2)
      IF(hpfac(iy,ix).GT.hpqmax)THEN
        hpqmax=hpfac(iy,ix)
      ENDIF
    END DO
    six = max(0.01*hpqmax,six)
 
    hpqmax=-999.0
    DO jsea=1, nseal
#ifdef W3_DIST
      isea   = iaproc + (jsea-1)*naproc
#endif
#ifdef W3_SHRD
      isea   = jsea
#endif
      ix     = mapsf(isea,1)
      iy     = mapsf(isea,2)
      IF(hqfac(iy,ix).GT.hpqmax)THEN
        hpqmax=hqfac(iy,ix)
      ENDIF
    END DO
    siy = max(0.01*hpqmax,siy)
 
    hmax   = max( 0. , hmax )
    !
    IF ( iaproc .EQ. napout ) THEN
      IF ( flagll ) THEN
        factor = 1.
        WRITE (ndso,940) fp, sip, thm, ncos, &
             factor*xm, min(9999.99,factor*six), factor*ym,  &
             min(9999.99,factor*siy), hmax
      ELSE
        factor = 1.e-3
        WRITE (ndso,941) fp, sip, thm, ncos, &
             factor*xm, min(9999.99,factor*six), factor*ym,  &
             min(9999.99,factor*siy), hmax
      END IF
    END IF
    !
    fp     = fp  * tpi
    sip    = sip * tpi
    thm    = mod( 630. - thm , 360. ) * dera
    !
    ! 4.b Make 1-D spectrum.
    !
    chsip  = 0.1 * dsip(1)
    flone  = sip .LT. chsip
    ikm    = nint( 1. + (log(fp)-log(fr1*tpi))/log(xfr) )
    ikm    = max( 1 , min( nk , ikm ) )
    !
    DO ik=1, nk
      IF ( flone ) THEN
        IF (ik.EQ.ikm) THEN
          e1(ik) = 1.
        ELSE
          e1(ik) = 0.
        END IF
      ELSE
        frrel  = (sig(ik)-fp)/sip
        IF (abs(frrel).LT.10) THEN
          e1(ik) = exp( -0.125 * frrel**2 )
        ELSE
          e1(ik) = 0.
        END IF
      END IF
    END DO
    !
#ifdef W3_O4
    IF ( iaproc .EQ. napout ) CALL prt1ds                      &
         (ndso, nk, e1, sig(1:), '  ', 10, 0.,                &
         'Unscaled 1-D', ' ', 'TEST E(f)')
#endif
    !
    ! 4.c Make directional distribution.
    !
    flone  = ncos .GT. 20
    ithm   = 1 + nint( thm / dth )
    DO ith=1, nth
      IF (flone) THEN
        IF ( ith .EQ. ithm ) THEN
          dd(ith) = 1.
        ELSE
          dd(ith) = 0.
        END IF
      ELSE
        dd(ith) = max( cos(th(ith)-thm) , 0. )**ncos
      END IF
    END DO
    !
    ! 4.d 2-D energy spectrum.
    !
    etot   = 0.
    DO ik=1, nk
      e1i    = 0.
      DO ith=1, nth
        e2(ith,ik) = e1(ik) * dd(ith)
        e1i        = e1i + e2(ith,ik)
      END DO
      etot   = etot + e1i * dsip(ik)
    END DO
    etot   = etot * dth
    factor = hmax**2 / ( 16. * etot )
    !
    e2     = factor * e2
    !
#ifdef W3_O5
    ALLOCATE ( e2out(nk,nth) )
    DO ith=1, nth
      DO ik=1, nk
        e2out(ik,ith) = tpi * e2(ith,ik)
      END DO
    END DO
#endif
    !
#ifdef W3_O5
    IF ( iaproc .EQ. napout ) CALL prt2ds                       &
         ( ndso, nk, nk, nth, e2out, sig(1:), ' ', dera*tpi, &
         0., 0.0001, 'Energy', 'm2s', 'TEST 2-D')
    DEALLOCATE ( e2out )
#endif
    !
    ! 4.e Distribute over grid.
    !
 
    DO ik=1, nk
      e21(1+(ik-1)*nth:ik*nth) = e2(:,ik)
    END DO
    !
    DO jsea=1, nseal
      !
#ifdef W3_DIST
      isea   = iaproc + (jsea-1)*naproc
#endif
#ifdef W3_SHRD
      isea   = jsea
#endif
      IF (gtype .EQ. ungtype) THEN
        ix     = mapsf(isea,1)
        x      = xgrd(1,ix)
        y      = ygrd(1,ix)
      ELSE
        ix     = mapsf(isea,1)
        iy     = mapsf(isea,2)
        x      = xgrd(iy,ix)
        y      = ygrd(iy,ix)
      ENDIF
      IF(nosix)THEN
        rdsqr  =(w3dist(flagll,x,y,xm,ym)/siy)**2
      ELSE
        rdsqr  =((x-xm)/six)**2 + ((y-ym)/siy)**2
      ENDIF
      IF ( rdsqr .GT. 40. ) THEN
        factor = 0.
      ELSE
        factor = exp( -0.5 * rdsqr )
      END IF
      !
#ifdef W3_EXPORTWWM
      factor = 1.
#endif
      va(:,jsea) = factor * e21
      !
 
      !
    END DO
    !
    !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ! 5.  ITYPE = 2, pre-defined JONSWAP.
    !
  ELSE IF ( itype .EQ. 2 ) THEN
    inxout = 'COLD'
    !
    ! 5.a Read parameters.
    !
    CALL nextln ( comstr , ndsi , ndsen )
    READ (ndsi,*,END=801,ERR=802)                               &
         alfa, fp, thm, gamma, siga, sigb, xm, six, ym, siy
    !
    IF (alfa.LE.0.) alfa = 0.0081
    IF (fp  .LE.0.) fp   = 0.10
    IF (siga.LE.0.) siga = 0.07
    IF (sigb.LE.0.) sigb = 0.09
    fp     = max( 0.5 * tpiinv * sig(1) , fp )
    fp     = min( tpiinv * sig(nk) , fp )
 
    nosix=.false.
    IF(six.LT.0.0)THEN
      IF ( iaproc .EQ. napout ) WRITE (ndso,903)
      nosix=.true.
    END IF
 
    hpqmax=-999.0
    DO jsea=1, nseal
#ifdef W3_DIST
      isea   = iaproc + (jsea-1)*naproc
#endif
#ifdef W3_SHRD
      isea   = jsea
#endif
      ix     = mapsf(isea,1)
      iy     = mapsf(isea,2)
      IF(hpfac(iy,ix).GT.hpqmax)THEN
        hpqmax=hpfac(iy,ix)
      ENDIF
    END DO
    six = max(0.01*hpqmax,six)
 
    hpqmax=-999.0
    DO jsea=1, nseal
#ifdef W3_DIST
      isea   = iaproc + (jsea-1)*naproc
#endif
#ifdef W3_SHRD
      isea   = jsea
#endif
      ix     = mapsf(isea,1)
      iy     = mapsf(isea,2)
      IF(hqfac(iy,ix).GT.hpqmax)THEN
        hpqmax=hqfac(iy,ix)
      ENDIF
    END DO
    siy = max(0.01*hpqmax,siy)
 
    DO
      IF ( thm .LT. 0. ) THEN
        thm    = thm + 360.
      ELSE
        EXIT
      END IF
    END DO
    thm    = mod( thm , 360. )
    gamma  = max(gamma,1.)
    yln    = log(gamma)
    !
    IF ( iaproc .EQ. napout ) THEN
      IF ( flagll ) THEN
        factor = 1.
        WRITE (ndso,950) alfa, fp, thm, gamma, siga, sigb,  &
             factor*xm, factor*six, factor*ym, factor*siy
      ELSE
        factor = 1.e-3
        WRITE (ndso,951) alfa, fp, thm, gamma, siga, sigb,  &
             factor*xm, factor*six, factor*ym, factor*siy
      END IF
    END IF
    thm    = mod( 630. - thm , 360. ) * dera
    !
    ! 5.b Make 1-D spectrum.
    !
    DO ik=1, nk
      fr     = sig(ik) * tpiinv
      e1(ik) = ej5p(fr, alfa, fp, yln, siga, sigb )
    END DO
    !
#ifdef W3_O4
    IF ( iaproc .EQ. napout ) CALL prt1ds                   &
         (ndso, nk, e1, sig(1:), '  ', 18, 0.,      &
         'E(f)', ' ', 'TEST 1-D')
#endif
    !
    ! 5.c 2-D energy spectrum.
    !     Factor 2pi to go to E(sigma,theta)
    !
    DO ik = 1,nk
      fr     = sig(ik) * tpiinv
      IF (fr.LT.fp) THEN
        beta =  4.06
      ELSE
        beta = -2.34
      END IF
      frr    = min( 2.5 , fr/fp )
      s      = 9.77 * frr**beta
      sumd   = 0.
      DO ith = 1,nth
        ang    = cos( 0.5 * ( thm - th(ith) ) )**2
        dd(ith) = 0.
        IF(ang.GT.1.e-20) THEN
          arg    = s * log(ang)
          IF(arg.GT.-170) dd(ith) = exp(arg)
        END IF
        sumd    = sumd + dd(ith)
      END DO
      factor = 1. / (tpi*sumd*dth)
      DO ith = 1,nth
        e2(ith,ik) = factor * e1(ik) * dd(ith)
      END DO
    END DO
    !
#ifdef W3_O5
    ALLOCATE ( e2out(nk,nth) )
    DO ith=1, nth
      DO ik=1, nk
        e2out(ik,ith) = tpi * e2(ith,ik)
      END DO
    END DO
#endif
    !
#ifdef W3_O5
    IF ( iaproc .EQ. napout ) CALL prt2ds                   &
         (ndso, nk, nk, nth, e2out, sig(1:), ' ', 1.,   &
         0., 0.0001, 'E(f,theta)', 'm2s', 'TEST 2-D')
    DEALLOCATE ( e2out )
#endif
    !
    ! 5.d Distribute over grid.
    !
 
    DO ik=1, nk
      e21(1+(ik-1)*nth:ik*nth) = e2(:,ik)
    END DO
    !
    !
    DO jsea=1, nseal
      !
#ifdef W3_DIST
      isea   = iaproc + (jsea-1)*naproc
#endif
#ifdef W3_SHRD
      isea   = jsea
#endif
      IF (gtype .EQ. ungtype) THEN
        ix     = mapsf(isea,1)
        x      = xgrd(1,ix)
        y      = ygrd(1,ix)
      ELSE
        ix     = mapsf(isea,1)
        iy     = mapsf(isea,2)
        x      = xgrd(iy,ix)
        y      = ygrd(iy,ix)
      ENDIF
      IF(nosix)THEN
        rdsqr  =(w3dist(flagll,x,y,xm,ym)/siy)**2
      ELSE
        rdsqr  =((x-xm)/six)**2 + ((y-ym)/siy)**2
      ENDIF
      IF ( rdsqr .GT. 40. ) THEN
        factor = 0.
      ELSE
        factor = exp( -0.5 * rdsqr )
      END IF
      !
      va(:,jsea) = factor * e21
      !
    END DO
    !
    !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ! 6.  ITYPE = 3, fetch limited JONSWAP.
    !
  ELSE IF ( itype .EQ. 3 ) THEN
    inxout = 'WIND'
    IF ( iaproc .EQ. napout ) WRITE (ndso,960)
    !
    !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ! 7.  ITYPE = 4, User defined.
    !
  ELSE IF ( itype .EQ. 4 ) THEN
    inxout = 'COLD'
    !
    ! 7.a Read parameters.
    !
    CALL nextln ( comstr , ndsi , ndsen )
    READ (ndsi,*,END=801,ERR=802) facs
    IF ( facs .LE. 0. ) facs = 1.
    IF ( iaproc .EQ. napout ) WRITE (ndso,970) facs
    !
    ! 7.b Read and rescale spectrum.
    !
    CALL nextln ( comstr , ndsi , ndsen )
    READ (ndsi,*,END=801,ERR=802)                               &
         ((finp(ik,ith),ik=1,nk),ith=1,nth)
    !
    finp = finp * facs / tpi
    !
#ifdef W3_O5
    IF ( iaproc .EQ. napout ) CALL prt2ds                   &
         (ndso, nk, nk, nth, finp, sig(1:), ' ', tpi,    &
         0., 0.0001, 'Energy', 'm2s', 'TEST 2-D')
#endif
    !
    ! 7.c Distribute over grid.
    !
    DO jsea=1, nseal
      !
#ifdef W3_DIST
      isea   = iaproc + (jsea-1)*naproc
#endif
#ifdef W3_SHRD
      isea   = jsea
#endif
      DO ik=1, nk
        DO ith=1, nth
          va(ith+(ik-1)*nth,jsea) = finp(ik,ith)
        END DO
      END DO
    END DO
    !
    !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ! 8.  ITYPE = 5, fetch limited JONSWAP.
    !
  ELSE
    inxout = 'CALM'
    IF ( iaproc .EQ. napout ) WRITE (ndso,980)
    !
  END IF
  !
  !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  ! 9.  Convert E(sigma) to N(k)
  !
  IF ( itype.NE.3 .AND. itype.NE.5 ) THEN
    IF ( iaproc .EQ. napout ) WRITE (ndso,990)
    !
#ifdef W3_O6
    ALLOCATE ( hsig(nx,ny) )
    hsig   = 0.
#endif
    !
    DO jsea=1, nseal
#ifdef W3_DIST
      isea   = iaproc + (jsea-1)*naproc
#endif
#ifdef W3_SHRD
      isea   = jsea
#endif
      depth  = max( dmin , -zb(isea) )
#ifdef W3_O6
      etot   = 0.
#endif
      DO ik=1, nk
        CALL wavnu1 ( sig(ik), depth, wn, cg )
#ifdef W3_O6
        e1i    = 0.
#endif
        DO ith=1, nth
#ifdef W3_O6
          e1i    = e1i + va(ith+(ik-1)*nth,jsea)
#endif
          va(ith+(ik-1)*nth,jsea) = va(ith+(ik-1)*nth,jsea) *   &
               cg / sig(ik)
        END DO
#ifdef W3_O6
        etot   = etot + e1i*dsip(ik)
#endif
      END DO
#ifdef W3_O6
      ix     = mapsf(isea,1)
      iy     = mapsf(isea,2)
      hsig(ix,iy) = 4. * sqrt( etot * dth )
#endif
#ifdef W3_EXPORTWWM
      IF (jsea .eq. 1) THEN
        DO ith=1,nth
          DO ik=1,nk
            ispec = ith + nth * (ik-1)
            WRITE(10003) ith, ik, va(ispec,jsea)
          END DO
        END DO
        WRITE(740+iaproc,*) 'FINAL : sum(VA)=', sum(va(:,jsea))
      END IF
#endif
    END DO
    !
#ifdef W3_O6
    ALLOCATE ( mapo(nx,ny) )
    DO ix=1, nx
      DO iy=1, ny
        mapo(ix,iy) = mapsta(iy,ix)
      END DO
    END DO
#endif
    !
#ifdef W3_MPI
    IF ( naproc .EQ. 1 ) THEN
#endif
#ifdef W3_O6
      nsx    = 1 + nx/35
      nsy    = 1 + ny/35
      IF ( iaproc .EQ. napout ) CALL prtblk                   &
           (ndso, nx, ny, nx, hsig, mapo, 0, 0.,       &
           1, nx, nsx, 1, ny, nsy, 'Hs', 'm')
#endif
#ifdef W3_MPI
    END IF
#endif
    !
  END IF
  !
  !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !10.  Write restart file.
  !
  IF ( iaproc .EQ. napout ) WRITE (ndso,995)
  CALL w3iors ( inxout, ndsr, sig(nk) )
  !
  GOTO 888
  !
  ! Escape locations read errors :
  !
800 CONTINUE
  IF ( iaproc .EQ. naperr ) WRITE (ndse,1000) ierr
  CALL extcde ( 10 )
  !
801 CONTINUE
  IF ( iaproc .EQ. naperr ) WRITE (ndse,1001)
  CALL extcde ( 11 )
  !
802 CONTINUE
  IF ( iaproc .EQ. naperr ) WRITE (ndse,1002) ierr
  CALL extcde ( 12 )
  !
#ifdef W3_DIST
803 CONTINUE
  IF ( iaproc .EQ. naperr ) WRITE (ndse,1003) nsea, naproc
  CALL extcde ( 13 )
#endif
  !
888 CONTINUE
  IF ( iaproc .EQ. napout ) WRITE (ndso,999)
#ifdef W3_MPI
  CALL mpi_finalize  ( ierr_mpi )
#endif
  !
  ! Formats
  !
900 FORMAT (/15x,'   *** WAVEWATCH III  Initial conditions ***   '/ &
       15x,'==============================================='/)
901 FORMAT ( '  Comment character is ''',a,''''/)
902 FORMAT ( '  Grid name : ',a/)
903 FORMAT ( '  Negative SIX was provided by user.         '/       &
       '  WW3 will create a gaussian distribution    '/       &
       '  that is circular in real space. ')
  !
930 FORMAT (/'  Initial field ITYPE =',i2/                          &
       ' --------------------------------------------------')
  !
940 FORMAT ( '       Gaussian / cosine power spectrum '//           &
       '       Peak frequency and spread (Hz)    :',2x,2f8.4/ &
       '       Mean direction (Naut., degr.)     :',f7.1/     &
       '       Cosine power of dir. distribution :',i5/       &
       '       Mean longitude and spread (degr.) :',2f8.2/    &
       '       Mean latitude and spread (degr.)  :',2f8.2/    &
       '       Maximum wave height               :',f8.2/)
  !
950 FORMAT ( '       JONSWAP spectrum'//                            &
       '       alfa                          (-) : ',f12.5/   &
       '       Peak frequecy                (Hz) : ',f11.4/   &
       '       Mean direction       (Naut.,deg.) : ',f 8.1/   &
       '       gamma                         (-) : ',f 9.2/   &
       '       sigma-A                       (-) : ',f11.4/   &
       '       sigma-B                       (-) : ',f11.4/   &
       '       Mean longitude and spread (degr.) : ',2f9.2/   &
       '       Mean latitude and spread  (degr.) : ',2f9.2)
941 FORMAT ( '       Gaussian / cosine power spectrum '//           &
       '       Peak frequency and spread (Hz)    :',2x,2f8.4/ &
       '       Mean direction (Naut., degr.)     :',f7.1/     &
       '       Cosine power of dir. distribution :',i5/       &
       '       Mean X and spread (km)            :',2f8.2/    &
       '       Mean Y and spread (km)            :',2f8.2/    &
       '       Maximum wave height               :',f8.2/)
  !
951 FORMAT ( '       JONSWAP spectrum'//                            &
       '       alfa                      (-) : ',f12.5/   &
       '       Peak frequecy            (Hz) : ',f11.4/   &
       '       Mean direction   (Naut.,deg.) : ',f 8.1/   &
       '       gamma                     (-) : ',f 9.2/   &
       '       sigma-A                   (-) : ',f11.4/   &
       '       sigma-B                   (-) : ',f11.4/   &
       '       Mean X and spread        (km) : ',2f9.2/   &
       '       Mean Y and spread        (km) : ',2f9.2)
  !
960 FORMAT ( '       Fetch-limited JONSWAP spectra based on local '/ &
       '       wind speed (fetch related to grid increment).')
  !
970 FORMAT ( '       User-defined energy spectrum F(f,theta).'//    &
       '       Scale factor             (-) : ',e12.4/)
  !
980 FORMAT ( '       Starting from calm conditions (Hs = 0)')
  !
990 FORMAT (/'       Converting energy to action ... ')
995 FORMAT (/'       Writing restart file  ... '/)
  !
999 FORMAT (/'  End of program '/                                   &
       ' ========================================='/          &
       '         WAVEWATCH III Initial conditions '/)
  !
1000 FORMAT (/' *** WAVEWATCH III ERROR IN W3STRT : '/               &
       '     ERROR IN OPENING INPUT FILE'/                    &
       '     IOSTAT =',i5/)
  !
1001 FORMAT (/' *** WAVEWATCH III ERROR IN W3STRT : '/               &
       '     PREMATURE END OF INPUT FILE'/)
  !
1002 FORMAT (/' *** WAVEWATCH III ERROR IN W3STRT : '/               &
       '     ERROR IN READING FROM INPUT FILE'/               &
       '     IOSTAT =',i5/)
  !
1010 FORMAT (/' *** WAVEWATCH III ERROR IN W3STRT : '/               &
       '     ILLEGAL TYPE, ITYPE =',i4/)
  !
#ifdef W3_DIST
1003 FORMAT (/' *** WAVEWATCH III ERROR IN W3STRT : '/               &
       '     NUMBER OF SEA POINTS LESS THAN NUMBER OF PROC.'/ &
       '     NSEA, NAPROC =',2i8/)
#endif
  !/
  !/ End of W3STRT ----------------------------------------------------- /
  !/
END PROGRAM w3strt