w3servmd.F90

Go to the documentation of this file.

#include "w3macros.h"
!/ ------------------------------------------------------------------- /
MODULE w3servmd
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         15-Jan-2021 |
  !/                  +-----------------------------------+
  !/
  !/    For update log see individual subroutines.
  !/    12-Jun-2012 : Add /RTD option or rotated grid option.
  !/                  (Jian-Guo Li)                       ( version 4.06 )
  !/    11-Nov-2013 : SMC and rotated grid incorporated in the main
  !/                  trunk                               ( version 4.13 )
  !/    18-Aug-2016 : Add dist_sphere: angular distance   ( version 5.11 )
  !/    01-Mar-2016 : Added W3THRTN and W3XYRTN for post  ( version 6.02 )
  !/                  processing rotated grid data
  !/    15-Jan-2021 : Added UV_TO_MAG_DIR routine         ( version 7.12 )
  !/
  !/    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 :
  !
  !     In this module all WAVEWATCH specific service routines have
  !     been gathered.
  !
  !  2. Variables and types :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      NDSTRC    Int.  Private  Data set number for output of STRACE
  !                               (set in ITRACE).
  !      NTRACE    Int.  Private  Maximum number of trace prints in
  !                               strace (set in ITRACE).
  !     ----------------------------------------------------------------
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      ITRACE    Subr. Public   (Re-) Initialization for STRACE.
  !      STRACE    Subr. Public   Enable subroutine tracing, usually
  !                               activated with the !/S switch.
  !      NEXTLN    Subr. Public   Get to next line in input command file.
  !      W3S2XY    Subr. Public   Grid conversion routine.
  !      EJ5P      R.F.  Public   Five parameter JONSWAP spectrum.
  !      WWDATE    Subr. Public   Get system date.
  !      WWTIME    Subr. Public   Get system time.
  !      EXTCDE    Subr. Public   Abort program with exit code.
  !     Four subs for rotated grid are appended to this module.  As they
  !     are shared with SMC grid, they are not quoted by option /RTD but
  !     are available for general use.     JGLi12Jun2012
  !     W3SPECTN       turns wave spectrum anti-clockwise by AnglD
  !     W3ACTURN       turns wave action(k,nth) anti-clockwise by AnglD.
  !     W3LLTOEQ       convert standard into rotated lat/lon, plus AnglD
  !     W3EQTOLL       revers of the LLTOEQ, but AnglD unchanged.
  !     W3THTRN        turns direction value anti-clockwise by AnglD
  !     W3XYTRN        turns 2D vectors anti-clockwise by AnglD
  !
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !     None.
  !
  !  5. Remarks :
  !
  !  6. Switches
  !
  !       !/S    Enable subroutine tracing using STRACE in this module.
  !
  !  7. Source code :
  !
  !/ ------------------------------------------------------------------- /
 
  ! module default
  implicit none
 
  PUBLIC
  !
  INTEGER, PRIVATE        :: NDSTRC = 6, ntrace = 0
  !
CONTAINS
  !/ ------------------------------------------------------------------- /
  SUBROUTINE itrace (NDS, NMAX)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         23-Nov-1999 |
    !/                  +-----------------------------------+
    !/
    !/    23-Nov-1999 : First version of routine.           ( version 2.00 )
    !/
    !  1. Purpose :
    !
    !     (Re-) initialization for module version of STRACE.
    !
    !  3. Parameter list
    !     ----------------------------------------------------------------
    !       NDS     Int.   I   Data set number ofr trace file.
    !       NMAX    Int.   I   Maximum number of traces per routine.
    !     ----------------------------------------------------------------
    !
    !     Private to module :
    !     ----------------------------------------------------------------
    !       NDSTRC  Int.  Output unit number for trace.     ( from NDS  )
    !       NTRACE  Int.  Maximum number of trace prints.   ( from NMAX )
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     None.
    !
    !  5. Called by :
    !
    !     Any program, multiple calls allowed.
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: NDS, NMAX
    !/
    !/ ------------------------------------------------------------------- /
    !/
    ntrace = max( 0 , nmax )
    ndstrc = nds
    !
    RETURN
    !/
    !/ End of ITRACE ----------------------------------------------------- /
    !/
  END SUBROUTINE itrace
  !/ ------------------------------------------------------------------- /
  SUBROUTINE strace (IENT, SNAME)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         25-Jan-2000 |
    !/                  +-----------------------------------+
    !/                                   Original version by N. Booij, DUT
    !/
    !/    30-Mar-1993 : Final FORTRAN 77                    ( version 1.18 )
    !/    23-Nov-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    25-Jan-2000 : Force flushing of uniit.            ( version 2.00 )
    !/                  This was taken out around version 3.01.
    !/
    !  1. Purpose :
    !
    !     Keep track of entered subroutines.
    !
    !  3. Parameter list
    !     ----------------------------------------------------------------
    !       IENT    Int.  I/O  Number of times that STRACE has been
    !                          called by the routine.
    !       SNAME   Char.  I   Name of the subroutine (max. 6 characters)
    !     ----------------------------------------------------------------
    !
    !     Private to module :
    !     ----------------------------------------------------------------
    !       NDSTRC  Int.  Output unit number for trace.
    !       NTRACE  Int.  Maximum number of trace prints.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     None.
    !
    !  5. Called by :
    !
    !     Any program, after private variables have been set by NTRACE.
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(INOUT)  :: IENT
    CHARACTER, INTENT(IN)   :: SNAME*(*)
    !/
    !/ ------------------------------------------------------------------- /
    !/
    IF (ntrace.EQ.0 .OR. ient.GE.ntrace) RETURN
    !
    ient = ient + 1
    IF (ient.EQ.1) THEN
      WRITE (ndstrc,10) sname
    ELSE
      WRITE (ndstrc,11) sname, ient
    END IF
    !
    RETURN
    !
    ! Formats
    !
10  FORMAT (' ---> TRACE SUBR : ',a6)
11  FORMAT (' ---> TRACE SUBR : ',a6,'  ENTRY: ',i6)
    !/
    !/ End of STRACE ----------------------------------------------------- /
    !/
  END SUBROUTINE strace
  !/ ------------------------------------------------------------------- /
  SUBROUTINE nextln ( CHCKC , NDSI , NDSE )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         10-Dec-2014 |
    !/                  +-----------------------------------+
    !/
    !/    15-Jan-1999 : Final FORTRAN 77                    ( version 1.18 )
    !/    18-Nov-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    10-Dec-2014 : Skip blank lines and leading blanks ( version 5.04 )
    !/
    !  1. Purpose :
    !
    !     Sets file pointer to next active line of input file, by skipping
    !     blank lines and lines starting with the character CHCKC. Leading
    !     white space is allowed before the character CHCKC.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       CHCKC   C*1   I  Check character for defining comment line.
    !       NDSI    Int.  I  Input dataset number.
    !       NDSE    Int.  I  Error output dataset number.
    !                        (No output if NDSE < 0).
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !       STRACE ( !/S switch )
    !
    !  5. Called by :
    !
    !       Any routine.
    !
    !  6. Error messages :
    !
    !     - On EOF or error in input file.
    !
    !  9. Switches :
    !
    !     !/S  Enable subroutine tracing.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: NDSI, NDSE
    CHARACTER, INTENT(IN)   :: CHCKC*1
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    INTEGER                 :: IERR
    CHARACTER(128)          :: MSG
    CHARACTER(256)          :: LINE, TEST
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'NEXTLN')
#endif
    !
100 CONTINUE
    ! read line
    READ ( ndsi, 900, END=800, ERR=801, IOSTAT=IERR, IOMSG=MSG ) line
    ! leading blanks removed and placed on the right
    test = adjustl( line )
    IF ( test(1:1).EQ.chckc .OR. len_trim(test).EQ.0 ) THEN
      ! if comment or blank line, then skip
      GOTO 100
    ELSE
      ! otherwise, backup to beginning of line
      backspace( ndsi, err=802, iostat=ierr, iomsg=msg )
    ENDIF
    RETURN
    !
800 CONTINUE
    IF ( ndse .GE. 0 ) WRITE (ndse,910)
    CALL extcde ( 1 )
    !
801 CONTINUE
    IF ( ndse .GE. 0 ) WRITE (ndse,911) ierr, trim(msg)
    CALL extcde ( 2 )
    !
802 CONTINUE
    IF ( ndse .GE. 0 ) WRITE (ndse,912) ierr, trim(msg)
    CALL extcde ( 3 )
    !
    ! Formats
    !
900 FORMAT (a)
910 FORMAT (/' *** WAVEWATCH III ERROR IN NEXTLN : '/ &
         '     PREMATURE END OF INPUT FILE'/)
911 FORMAT (/' *** WAVEWATCH III ERROR IN NEXTLN : '/ &
         '     ERROR IN READING FROM FILE'/           &
         '     IOSTAT =',i5,/                         &
         '     IOMSG = ',a/)
912 FORMAT (/' *** WAVEWATCH III ERROR IN NEXTLN : '/ &
         '     ERROR ON BACKSPACE'/                   &
         '     IOSTAT =',i5,/                         &
         '     IOMSG = ',a/)
    !/
    !/ End of NEXTLN ----------------------------------------------------- /
    !/
  END SUBROUTINE nextln
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3s2xy ( NSEA, MSEA, MX, MY, S, MAPSF, XY )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III            NOAA/NMC |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         23-Nov-1999 |
    !/                  +-----------------------------------+
    !/
    !/    11-Dec-1996 : Final FORTRAN 77                    ( version 1.18 )
    !/    23-Nov-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/
    !  1. Purpose :
    !
    !     Convert a data array on the storage grid to a data array on the
    !     full spatial grid. Land and ice points in the full grid are
    !     not touched. Output array of conventional type XY(IX,IY).
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       NSEA    Int.   I    Number of sea points.
    !       MSEA, MX, MY
    !               Int.   I    Array dimensions.
    !       S       R.A.   I    Data on storage grid.
    !       MAPSF   I.A.   I    Storage map for IX and IY, resp.
    !       XY      R.A.   O    Data on XY grid.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     None.
    !
    !  5. Called by :
    !
    !     Any WAVEWATCH III routine.
    !
    !  9. Switches :
    !
    !     None.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: MSEA, NSEA, MX, MY, MAPSF(MSEA,2)
    REAL, INTENT(IN)        :: S(MSEA)
    REAL, INTENT(OUT)       :: XY(MX,MY)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: ISEA, IX, IY
    !/
    !/ ------------------------------------------------------------------- /
    !/
    DO isea=1, nsea
      ix     = mapsf(isea,1)
      iy     = mapsf(isea,2)
      xy(ix,iy) = s(isea)
    end do
    !/
    !/ End of W3S2XY ----------------------------------------------------- /
    !/
  END SUBROUTINE w3s2xy
  !/ ------------------------------------------------------------------- /
  REAL FUNCTION EJ5P ( F, ALFA, FP, YLN, SIGA, SIGB )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         23-Nov-1999 |
    !/                  +-----------------------------------+
    !/
    !/    23-AMy-1985 : Original by G. Ph. van Vledder.
    !/    23-Nov-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/
    !  1. Purpose :
    !
    !     Computation of spectral density using a 5-parameter
    !     JONSWAP-spectrum
    !
    !  2. Method
    !
    !     EJ5P(F) = A.EXP(B + LN(Y).EXP(C))
    !
    !     where: A = ALFA * 0.06175 * F**(-5)
    !            B = -1.25*(FP/F)**4
    !            C = -0.5 * ((F - FP)/(SIG * FP))**2
    !     and
    !            GRAV**2/(2.PI)**4 = 0.06175
    !
    !  3. Parameters :
    !
    !     Parameter list
    !
    !     ----------------------------------------------------------------
    !       F       Real   I    Frequency in Hz
    !       ALFA    Real   I    Energy scaling factor
    !       FP      Real   I    Peak frequency in Hz
    !       YLN     Real   I    Peak overshoot factor, given by LN-value
    !       SIGA    Real   I    Spectral width, for F < FP
    !       SIGB    Real   I    Spectral width, FOR F > FP
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     None.
    !
    !  5. Called by :
    !
    !     Any.
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !     EXPMIN is a machine dependant constant such that
    !     EXP(EXPMIN) can be successfully evaluated without
    !     underflow by the compiler supllied EXP routine.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     None.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(IN)        :: f, alfa, fp, yln, siga, sigb
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    REAL                    :: sig, a, b, c
    REAL, SAVE              :: eps=1.e-4, expmin=-180.
    !/
    !/ ------------------------------------------------------------------- /
    !/
    IF(f.LT.eps) THEN
      ej5p = 0.0
      RETURN
    END IF
    !
    a = alfa * 0.06175 / f**5
    b = -1.25 * (fp/f)**4
    b = max(b,expmin)
    !
    IF (yln.LT.eps) THEN
      ej5p = a * exp(b)
    ELSE
      IF( f.LE.fp) THEN
        sig = siga
      ELSE
        sig = sigb
      END IF
      c = -0.5 * ((f - fp)/(sig * fp))**2
      c = max(c,expmin)
      ej5p = a * exp(b + exp(c) * yln)
    END IF
    !
    RETURN
    !/
    !/ End of NEXTLN ----------------------------------------------------- /
    !/
  END FUNCTION ej5p
  !/ ------------------------------------------------------------------- /
  REAL function dist_sphere ( lo1,la1,lo2,la2 )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           F. Ardhuin              |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         18-Aug-2016 |
    !/                  +-----------------------------------+
    !/
    !/    18-Aug-2016 :  Creation                           ( version 5.11 )
    !/
    !  1. Purpose :
    !
    !     Computes distance between two points on a sphere
    !
    !  2. Method
    !
    !
    !  3. Parameters :
    !
    !     Parameter list
    !
    !     ----------------------------------------------------------------
    !       LO1     Real   I    Longitude of 1st point
    !       LA1     Real   I    Latitude of 1st point
    !       LO2     Real   I    Longitude of 2nd point
    !       LA2     Real   I    Latitude of 2nd point
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     None.
    !
    !  5. Called by :
    !
    !     WW3_BOUNC
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !     None.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     None.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(IN)        :: lo1, la1, lo2, la2
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    !  None
    !/
    !/ ------------------------------------------------------------------- /
    !/
    dist_sphere=acos(sin(la2*dera)*sin(la1*dera)+ &
         cos(la2*dera)*cos(la1*dera)*cos((lo2-lo1)*dera))*rade
    !
    RETURN
    !/
    !/ End of NEXTLN ----------------------------------------------------- /
    !/
  END FUNCTION dist_sphere
  !/ ------------------------------------------------------------------- /
 
  !/ ------------------------------------------------------------------- /
  SUBROUTINE wwdate (STRNG)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         26-Dec-2012 |
    !/                  +-----------------------------------+
    !/
    !/    23-Dec-1998 : Final FORTRAN 77                    ( version 1.18 )
    !/    23-Nov-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    18-Sep-2000 : PGI switch added                    ( version 2.04 )
    !/    13-Mar-2001 : LF95 switch added                   ( version 2.09 )
    !/    08-May-2002 : Replace obsolete switches with F90  ( version 2.21 )
    !/    26-Dec-2012 : Modified obsolete declarations.     ( version 4.11 )
    !/
    !  1. Purpose :
    !
    !     Get date from machine dependent routine.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       STRNG   C*10   O   String with date in format YYYY/MM/DD
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     Machine dependent.
    !
    !  5. Called by :
    !
    !     Any routine.
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    CHARACTER, INTENT(OUT)  :: STRNG*10
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    CHARACTER(LEN=8)        :: DATE
    CHARACTER(LEN=10)       :: TIME
    CHARACTER(LEN=5)        :: ZONE
    INTEGER                 :: VALUES(8)
    !/
    !/ ------------------------------------------------------------------- /
    !/
    strng = '----/--/--'
    CALL date_and_time ( date, time, zone, values )
    strng(1:4) = date(1:4)
    strng(6:7) = date(5:6)
    strng(9:10) = date(7:8)
    !
    !
    RETURN
    !/
    !/ End of WWDATE ----------------------------------------------------- /
    !/
  END SUBROUTINE wwdate
  !/ ------------------------------------------------------------------- /
  SUBROUTINE wwtime (STRNG)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         26-Dec-2012 |
    !/                  +-----------------------------------+
    !/
    !/    23-Dec-1998 : Final FORTRAN 77                    ( version 1.18 )
    !/    23-Nov-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    18-Sep-2000 : PGI switch added                    ( version 2.04 )
    !/    13-Mar-2001 : LF95 switch added                   ( version 2.09 )
    !/    08-May-2002 : Replace obsolete switches with F90  ( version 2.21 )
    !/    26-Dec-2012 : Modified obsolete declarations.     ( version 4.11 )
    !/
    !  1. Purpose :
    !
    !     Get time from machine dependent routine.
    !
    !  2. Method :
    !
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       STRNG   C*8    O   String with time in format hh:mm:ss
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     Machine dependent.
    !
    !  5. Called by :
    !
    !     Any routine.
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    CHARACTER, INTENT(OUT)  :: STRNG*8
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    CHARACTER(LEN=8)        :: DATE
    CHARACTER(LEN=10)       :: TIME
    CHARACTER(LEN=5)        :: ZONE
    INTEGER                 :: VALUES(8)
    !/
    !/ ------------------------------------------------------------------- /
    !/
    !
    strng = '--:--:--'
    CALL date_and_time ( date, time, zone, values )
    strng(1:2) = time(1:2)
    strng(4:5) = time(3:4)
    strng(7:8) = time(5:6)
    !
    RETURN
    !/
    !/ End of WWTIME ----------------------------------------------------- /
    !/
  END SUBROUTINE wwtime
  !/ ------------------------------------------------------------------- /
  SUBROUTINE extcde ( IEXIT, UNIT, MSG, FILE, LINE, COMM )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         06-Jun-2018 |
    !/                  +-----------------------------------+
    !/
    !/    06-Jan-1998 : Final FORTRAN 77                    ( version 1.18 )
    !/    23-Nov-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    10-Dec-2014 : Add checks for allocate status      ( version 5.04 )
    !/    11-Mar-2015 : Allow non-error exit (iexit=0)      ( version 5.04 )
    !/    20-Jan-2017 : Add optional MPI communicator arg   ( version 6.02 )
    !/    06-Jun-2018 : Add optional MPI                    ( version 6.04 )
    !/
    !  1. Purpose :
    !
    !     Perform a program stop with an exit code.
    !
    !     If exit code IEXIT=0, then it is not an error, but
    !     a stop has been requested by the calling routine:
    !     wait for other processes in communicator to catch up.
    !
    !     If exit code IEXIT.ne.0, then abort program w/out
    !     waiting for other processes to catch up (important for example
    !     when not all processes are used by WW3).
    !
    !  2. Method :
    !
    !     Machine dependent.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       IEXIT   Int.   I   Exit code to be used.
    !       UNIT    Int.   I   (optional) file unit to write error message
    !       MSG     Str.   I   (optional) error message
    !       FILE    Str.   I   (optional) name of source code file
    !       LINE    Int.   I   (optional) line number in source code file
    !       COMM    Int.   I   (optional) MPI communicator
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !  5. Called by :
    !
    !     Any.
    !
    !  9. Switches :
    !
    !     !/MPI  MPI finalize interface if active
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !
#ifdef W3_MPI
    include "mpif.h"
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN) :: IEXIT
    INTEGER,      INTENT(IN), OPTIONAL :: UNIT
    CHARACTER(*), INTENT(IN), OPTIONAL :: MSG
    CHARACTER(*), INTENT(IN), OPTIONAL :: FILE
    INTEGER,      INTENT(IN), OPTIONAL :: LINE
    INTEGER,      INTENT(IN), OPTIONAL :: COMM
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_MPI
    INTEGER                 :: IERR_MPI
    LOGICAL                 :: RUN
#endif
    INTEGER                 :: IUN
    CHARACTER(256)          :: LMSG = ""
    CHARACTER(6)            :: LSTR
    CHARACTER(10)           :: PREFIX = "WW3 ERROR:"
    !/
    !/ Set file unit for error output
    !/
    iun = 0
    IF (PRESENT(unit)) iun = unit
    !/
    !/ Report error message
    !/
    IF (PRESENT(msg)) THEN
      WRITE (iun,"(A)") prefix//" "//trim(msg)
    END IF
    !/
    !/ Report context
    !/
    IF ( PRESENT(file) ) THEN
      lmsg = trim(lmsg)//" FILE="//trim(file)
    END IF
    IF ( PRESENT(line) ) THEN
      WRITE (lstr,'(I0)') line
      lmsg = trim(lmsg)//" LINE="//trim(lstr)
    END IF
    IF ( len_trim(lmsg).GT.0 ) THEN
      WRITE (iun,"(A)") prefix//trim(lmsg)
    END IF
    !/
    !/ Handle MPI exit
    !/
#ifdef W3_MPI
    CALL mpi_initialized ( run, ierr_mpi )
    IF ( run ) THEN
      IF ( iexit.EQ.0 ) THEN ! non-error state
        IF ( PRESENT(comm) ) CALL mpi_barrier ( comm, ierr_mpi )
        CALL mpi_finalize (ierr_mpi )
      ELSE ! error state
        WRITE(*,'(/A,I6/)') 'EXTCDE MPI_ABORT, IEXIT=', iexit
        IF (PRESENT(unit)) THEN
          WRITE(*,'(/A,I6/)') 'EXTCDE UNIT=', unit
        END IF
        IF (PRESENT(msg)) THEN
          WRITE(*,'(/2A/)') 'EXTCDE MSG=', msg
        END IF
        IF (PRESENT(file)) THEN
          WRITE(*,'(/2A/)') 'EXTCDE FILE=', file
        END IF
        IF (PRESENT(line)) THEN
          WRITE(*,'(/A,I8/)') 'EXTCDE LINE=', line
        END IF
        IF (PRESENT(comm)) THEN
          WRITE(*,'(/A,I6/)') 'EXTCDE COMM=', comm
        END IF
        CALL mpi_abort ( mpi_comm_world, iexit, ierr_mpi )
      END IF
    END IF
#endif
    !/
    !/ Handle non-MPI exit
    !/
    CALL exit ( iexit )
    !/
    !/ End of EXTCDE ----------------------------------------------------- /
    !/
  END SUBROUTINE extcde
  !/ ------------------------------------------------------------------- /
  !  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  !  This subroutine turn the wave spectrum by an fixed angle anti-clockwise
  !  so that it may be used in the rotated or stanadard system.
  !  First created:   26 Aug 2005   Jian-Guo Li
  !  Last modified:   21 Feb 2008   Jian-Guo Li
  !
  ! Subroutine Interface:
 
  Subroutine w3spectn( NFreq, NDirc, Alpha, Spectr )
 
    ! Description:
    !   Rotates wave spectrum anticlockwise by angle alpha in degree
    !   This routine is distinct from W3ACTURN since orders spectrum as freq, dirn
    !
    ! Subroutine arguments
    INTEGER, INTENT(IN) :: NFreq, NDirc         ! No. freq and dirn bins
    REAL,    INTENT(IN) :: Alpha                ! Turning angle (degrees)
    REAL, INTENT(INOUT) :: Spectr(NFreq,NDirc)  ! Wave spectrum in/out
 
    ! Local variables
    INTEGER :: ii, jj, kk, nsft
    REAL    :: Ddirc, frac, CNST
    REAL, Dimension(NFreq)      ::  Wrkfrq, Tmpfrq
    REAL, Dimension(NFreq,NDirc)::  Wrkspc
 
    ! Check input bin numbers
    IF( (nfreq .LT. 0) .OR. (ndirc .LT. 0) )  THEN
      print*, " Invalid bin number NF or ND", nfreq, ndirc
      RETURN
    ELSE
      ddirc=360.0/float(ndirc)
    ENDIF
 
    ! Work out shift bin number and fraction
 
    cnst=alpha/ddirc
    nsft=int( cnst )
    frac= cnst - float( nsft )
    !     PRINT*, ' nsft and frac =', nsft, frac
 
    ! Shift nsft bins if >=1
    IF( abs(nsft) .GE. 1 )  THEN
      DO ii=1, ndirc
 
        ! Wave spectral direction bin number is assumed to increase Anti-clockwise from EAST
        ! So shift nsft bins anticlockwise results in local bin number decreasing by nsft
        jj=ii - nsft
 
        ! As nsft may be either positive or negative depends on alpha, wrapping may
        ! happen in either ends of the bin number train
        IF( jj > ndirc )  jj=jj - ndirc
        IF( jj < 1     )  jj=jj + ndirc
 
        ! Copy the selected bin to the loop bin number
        wrkspc(:,ii)=spectr(:,jj)
 
      ENDDO
 
      ! If nsft=0, no need to shift, simply copy
    ELSE
      wrkspc = spectr
    ENDIF
 
    ! Pass fraction of wave energy in frac direction
    ! Wave spectral direction bin number is assumed to increase Anti-clockwise from EAST
    ! So Positive frac or anticlock case, smaller bin upstream
    IF( frac > 0.0 ) THEN
      tmpfrq=wrkspc(:,ndirc)*frac
      DO kk=1, ndirc
        wrkfrq=wrkspc(:,kk)*frac
        spectr(:,kk)=wrkspc(:,kk) - wrkfrq + tmpfrq
        tmpfrq=wrkfrq
      ENDDO
    ELSE
      ! Negative or clockwise case, larger bin upstream
      tmpfrq=wrkspc(:,1)*frac
      DO kk=ndirc, 1, -1
        wrkfrq=wrkspc(:,kk)*frac
        spectr(:,kk)=wrkspc(:,kk) + wrkfrq - tmpfrq
        tmpfrq=wrkfrq
      ENDDO
    ENDIF
 
    ! Spectral turning completed
 
    RETURN
  END SUBROUTINE w3spectn
  !
  !  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  !  This subroutine turn the wave action by an angle (deg) anti-clockwise
  !  so that it may be used in the rotated or stanadard system.
  !  First created:   26 Aug 2005   Jian-Guo Li
  !  Last modified:    9 Oct 2008   Jian-Guo Li
  !
  ! Subroutine Interface:
 
  Subroutine w3acturn( NDirc, NFreq, Alpha, Spectr )
 
    ! Description:
    !   Rotates wave spectrum anticlockwise by angle alpha
    !   Routine is distinct from W3SPECTN since orders spectrum as dirn, freq
    !
    ! Subroutine arguments
    INTEGER, INTENT(IN) :: NFreq, NDirc          ! No. freq and dirn bins
    REAL,    INTENT(IN) :: Alpha                 ! Turning angle (degrees)
    REAL, INTENT(INOUT) :: Spectr(NDirc, NFreq)  ! Wave action in/out
 
    ! Local variables
    INTEGER :: ii, jj, kk, nsft
    REAL    :: Ddirc, frac, CNST
    REAL, Dimension(NFreq)      ::  Wrkfrq, Tmpfrq
    REAL, Dimension(NDirc,NFreq)::  Wrkspc
 
    ! Check input bin numbers
    IF( (nfreq .LT. 0) .OR. (ndirc .LT. 0) )  THEN
      print*, " Invalid bin number NF or ND", nfreq, ndirc
      RETURN
    ELSE
      ddirc=360.0/float(ndirc)
    ENDIF
 
    ! Work out shift bin number and fraction
 
    cnst=alpha/ddirc
    nsft=int( cnst )
    frac= cnst - float( nsft )
    !     PRINT*, ' nsft and frac =', nsft, frac
 
    ! Shift nsft bins if >=1
    IF( abs(nsft) .GE. 1 )  THEN
      DO ii=1, ndirc
 
        ! Wave spectral direction bin number is assumed to increase Anti-clockwise from EAST
        ! So shift nsft bins anticlockwise results in local bin number decreasing by nsft
        jj=ii - nsft
 
        ! As nsft may be either positive or negative depends on alpha, wrapping may
        ! happen in either ends of the bin number train
        IF( jj > ndirc )  jj=jj - ndirc
        IF( jj < 1     )  jj=jj + ndirc
 
        ! Copy the selected bin to the loop bin number
        wrkspc(ii,:)=spectr(jj,:)
 
      ENDDO
 
      ! If nsft=0, no need to shift, simply copy
    ELSE
      wrkspc = spectr
    ENDIF
 
    ! Pass fraction of wave energy in frac direction
    ! Wave spectral direction bin number is assumed to increase anti-clockwise from EAST
    ! So positive frac or anticlock case, smaller bin upstream
    IF( frac > 0.0 ) THEN
      tmpfrq=wrkspc(ndirc,:)*frac
      DO kk=1, ndirc
        wrkfrq=wrkspc(kk,:)*frac
        spectr(kk,:)=wrkspc(kk,:) - wrkfrq + tmpfrq
        tmpfrq=wrkfrq
      ENDDO
    ELSE
      ! Negative or clockwise case, larger bin upstream
      tmpfrq=wrkspc(1,:)*frac
      DO kk=ndirc, 1, -1
        wrkfrq=wrkspc(kk,:)*frac
        spectr(kk,:)=wrkspc(kk,:) + wrkfrq - tmpfrq
        tmpfrq=wrkfrq
      ENDDO
    ENDIF
 
    ! Spectral turning completed
 
    RETURN
  END SUBROUTINE w3acturn
  !
  !Li  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  !Li
  !Li  Merged UM source code for rotated grid, consisting the following
  !Li  original subroutines in UM 6.1
  !Li    LLTOEQ1A  WCOEFF1A  and  LBCROTWINDS1
  !Li  The last subroutine is modified to process only one level winds
  !Li  cpp directives are removed and required header C_Pi.h inserted.
  !Li         Jian-Guo Li     26 May 2005
  !Li
  !Li  The WCOEFF1A subroutine is merged into LLTOEQ to reduce repetition
  !Li  of the same calculations. Subroutine interface changed to
  !Li  LLTOEQANGLE
  !Li         Jian-GUo Li     23 Aug 2005
  !Li
  !Li  Subroutine W3LLTOEQ   --------------------------------------------
  !Li
  !Li  Purpose:  Calculates latitude and longitude on equatorial
  !Li            latitude-longitude (eq) grid used in regional
  !Li            models from input arrays of latitude and
  !Li            longitude on standard grid. Both input and output
  !Li            latitudes and longitudes are in degrees.
  !Li            Also calculate rotation angle in degree to tranform
  !Li            standard wind velocity into equatorial wind.
  !Li            Valid for 0<PHI_POLE<90 or new pole in N. hemisphere.
  !Li
  !* Arguments:--------------------------------------------------------
  SUBROUTINE w3lltoeq ( PHI, LAMBDA, PHI_EQ, LAMBDA_EQ, ANGLED, PHI_POLE, &
       LAMBDA_POLE, POINTS )
 
    INTEGER:: POINTS    !IN  Number of points to be processed
 
    REAL :: PHI_POLE,  & !IN  Latitude of equatorial lat-lon pole
         &        LAMBDA_POLE  !INOUT  Longitude of equatorial lat-lon pole
 
    REAL, DIMENSION(POINTS) ::         &
         &        PHI,       & !IN  Latitude
         &        LAMBDA,    & !IN  Longitude
         &        ANGLED,    & !OUT turning angle in deg for standard wind
         &        LAMBDA_EQ, & !OUT Longitude in equatorial lat-lon coords
         &        PHI_EQ       !OUT Latitude in equatorial lat-lon coords
 
    ! Define local varables:-----------------------------------------------
    REAL(KIND=8) :: a_lambda, a_phi, e_lambda, e_phi, sin_phi_pole, cos_phi_pole, &
         term1, term2, arg, lambda_zero, lambda_pole_keep
    INTEGER      :: I
 
    REAL(KIND=8), parameter :: small=1.0e-6
 
    ! Double precision versions of values in constants.ftn:
    REAL(KIND=8), parameter         :: pi = 3.141592653589793
    REAL(KIND=8), parameter         :: recip_pi_over_180 = 180. / pi
    REAL(KIND=8), parameter         :: pi_over_180   = pi / 180.
 
    !*----------------------------------------------------------------------
 
    ! 1. Initialise local constants
    ! Scale lambda pole to range -180 to 180 degs
    lambda_pole_keep=lambda_pole
    IF (lambda_pole.LE.-180.0) lambda_pole=lambda_pole+360.d0
    IF (lambda_pole.GT. 180.0) lambda_pole=lambda_pole-360.d0
 
    ! Latitude of zeroth meridian
    lambda_zero=lambda_pole+180.d0
    ! Sine and cosine of latitude of eq pole
    IF (phi_pole >= 0.0) THEN
      sin_phi_pole =  sin(pi_over_180*phi_pole)
      cos_phi_pole =  cos(pi_over_180*phi_pole)
    ELSE
      sin_phi_pole = -sin(pi_over_180*phi_pole)
      cos_phi_pole = -cos(pi_over_180*phi_pole)
    ENDIF
 
    ! 2. Transform from standard to equatorial latitude-longitude
 
    DO i= 1, points
 
      ! Scale longitude to range -180 to +180 degs
 
      a_lambda=lambda(i)-lambda_zero
      IF(a_lambda.GT. 180.0) a_lambda=a_lambda-360.d0
      IF(a_lambda.LE.-180.0) a_lambda=a_lambda+360.d0
 
      ! Convert latitude & longitude to radians
 
      a_lambda=pi_over_180*a_lambda
      a_phi=pi_over_180*phi(i)
 
      ! Compute eq latitude using equation (4.4)
 
      arg=-cos_phi_pole*cos(a_phi)*cos(a_lambda) + sin_phi_pole*sin(a_phi)
      arg=min(arg, 1.d0)
      arg=max(arg,-1.d0)
      e_phi=asin(arg)
      phi_eq(i)=recip_pi_over_180*e_phi
 
      ! Compute eq longitude using equation (4.6)
 
      term1 = sin_phi_pole*cos(a_phi)*cos(a_lambda) + cos_phi_pole*sin(a_phi)
      term2 = cos(e_phi)
      IF(term2 .LT. small) THEN
        e_lambda=0.d0
      ELSE
        arg=term1/term2
        arg=min(arg, 1.d0)
        arg=max(arg,-1.d0)
        e_lambda=recip_pi_over_180*acos(arg)
        e_lambda=sign(e_lambda,a_lambda)
      ENDIF
 
      ! Scale longitude to range 0 to 360 degs
 
      IF(e_lambda.GE.360.0) e_lambda=e_lambda-360.d0
      IF(e_lambda.LT.  0.0) e_lambda=e_lambda+360.d0
      lambda_eq(i)=e_lambda
 
      !Li  Calculate turning angle for standard wind velocity
 
      e_lambda=pi_over_180*lambda_eq(i)
 
      ! Formulae used are from eqs (4.19) and (4.21)
 
      term2=sin(e_lambda)
      arg= sin(a_lambda)*term2*sin_phi_pole      &
           &    +cos(a_lambda)*cos(e_lambda)
      arg=min(arg, 1.d0)
      arg=max(arg,-1.d0)
      term1=recip_pi_over_180*acos(arg)
      angled(i)=sign(term1,term2)
      !Li
 
    ENDDO
 
    ! Reset Lambda pole to the setting on entry to subroutine
    lambda_pole=lambda_pole_keep
 
    RETURN
  END SUBROUTINE w3lltoeq
  !
  !Li  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  !Li
  !Li  Merged UM source code for rotated grid, consiting the following
  !Li  original subroutines in UM 6.1
  !Li    EQTOLL1A  WCOEFF1A  and  LBCROTWINDS1
  !Li  The last subroutine is modified to process only one level winds
  !Li  cpp directives are removed and required header C_Pi.h inserted.
  !Li         Jian-Guo Li     26 May 2005
  !Li
  !Li  The WCOEFF1A subroutine is merged into EQTOLL to reduce repetition
  !Li  of the same calculations. Subroutine interface changed to
  !Li  EQTOLLANGLE
  !Li  First created:   Jian-GUo Li     23 Aug 2005
  !Li  Last modified:   Jian-GUo Li     25 Feb 2008
  !Li
  !Li  Subroutine W3EQTOLL  --------------------------------------------
  !Li
  !Li  Purpose:  Calculates latitude and longitude on standard grid
  !Li            from input arrays of latitude and longitude on
  !Li            equatorial latitude-longitude (eq) grid used
  !Li            in regional models. Both input and output latitudes
  !Li            and longitudes are in degrees.
  !Li            Also calculate rotation angle in degree to tranform
  !Li            standard wind velocity into equatorial wind.
  !Li            Valid for 0<PHI_POLE<90 or new pole in N. hemisphere.
  !Li
  !Li  Arguments:--------------------------------------------------------
 
  SUBROUTINE w3eqtoll( PHI_EQ, LAMBDA_EQ, PHI, LAMBDA,   &
       &                 ANGLED, PHI_POLE, LAMBDA_POLE, POINTS )
 
    INTEGER:: POINTS      !IN  Number of points to be processed
 
    REAL :: PHI_POLE,   & !IN  Latitude of equatorial lat-lon pole
         &        LAMBDA_POLE   !IN  Longitude of equatorial lat-lon pole
 
    REAL, DIMENSION(POINTS) ::         &
         &        PHI,       & !OUT Latitude
         &        LAMBDA,    & !OUT Longitude (0 =< LON < 360)
         &        ANGLED,    & !OUT turning angle in deg for standard wind
         &        LAMBDA_EQ, & !IN  Longitude in equatorial lat-lon coords
         &        PHI_EQ       !IN  Latitude in equatorial lat-lon coords
 
    ! Local varables:------------------------------------------------------
    REAL(KIND=8) :: e_lambda, e_phi, a_lambda, a_phi,                 &
         sin_phi_pole, cos_phi_pole,                       &
         term1, term2, arg, lambda_zero
    INTEGER :: I
 
    REAL(KIND=8), parameter :: small=1.0e-6
 
    ! Double precision versions of values in constants.ftn:
    REAL(KIND=8), parameter         :: pi = 3.141592653589793
    REAL(KIND=8), parameter         :: recip_pi_over_180 = 180. / pi
    REAL(KIND=8), parameter         :: pi_over_180   = pi / 180.
 
    ! ----------------------------------------------------------------------
 
    ! 1. Initialise local constants
 
    ! Latitude of zeroth meridian
    lambda_zero=lambda_pole+180.d0
    ! Sine and cosine of latitude of eq pole
    IF (phi_pole >= 0.0) THEN
      sin_phi_pole =  sin(pi_over_180*phi_pole)
      cos_phi_pole =  cos(pi_over_180*phi_pole)
    ELSE
      sin_phi_pole = -sin(pi_over_180*phi_pole)
      cos_phi_pole = -cos(pi_over_180*phi_pole)
    ENDIF
 
    ! 2. Transform from equatorial to standard latitude-longitude
 
    DO i= 1, points
 
      ! Scale eq longitude to range -180 to +180 degs
 
      e_lambda=lambda_eq(i)
      IF(e_lambda.GT. 180.0) e_lambda=e_lambda-360.d0
      IF(e_lambda.LT.-180.0) e_lambda=e_lambda+360.d0
 
      ! Convert eq latitude & longitude to radians
 
      e_lambda=pi_over_180*e_lambda
      e_phi=pi_over_180*phi_eq(i)
 
      ! Compute latitude using equation (4.7)
 
      arg=cos_phi_pole*cos(e_phi)*cos(e_lambda) + sin_phi_pole*sin(e_phi)
      arg=min(arg, 1.d0)
      arg=max(arg,-1.d0)
      a_phi=asin(arg)
      phi(i)=recip_pi_over_180*a_phi
 
      ! Compute longitude using equation (4.8)
 
      term1 = cos(e_phi)*sin_phi_pole*cos(e_lambda) - sin(e_phi)*cos_phi_pole
      term2 = cos(a_phi)
      IF(term2.LT.small) THEN
        a_lambda=0.d0
      ELSE
        arg=term1/term2
        arg=min(arg, 1.d0)
        arg=max(arg,-1.d0)
        a_lambda=recip_pi_over_180*acos(arg)
        a_lambda=sign(a_lambda,e_lambda)
        a_lambda=a_lambda+lambda_zero
      END IF
 
      ! Scale longitude to range 0 to 360 degs
 
      IF(a_lambda.GE.360.0) a_lambda=a_lambda-360.d0
      IF(a_lambda.LT.  0.0) a_lambda=a_lambda+360.d0
      lambda(i)=a_lambda
 
      !Li  Calculate turning angle for standard wind velocity
 
      a_lambda=pi_over_180*(lambda(i)-lambda_zero)
 
      ! Formulae used are from eqs (4.19) and (4.21)
 
      term2=sin(e_lambda)
      arg=sin(a_lambda)*term2*sin_phi_pole     &
           &           +cos(a_lambda)*cos(e_lambda)
      arg=min(arg, 1.d0)
      arg=max(arg,-1.d0)
      term1=recip_pi_over_180*acos(arg)
      angled(i)=sign(term1,term2)
      !Li
 
    ENDDO
 
    RETURN
  END SUBROUTINE w3eqtoll
 
  !Li
  !/ ------------------------------------------------------------------- /
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3thrtn ( NSEA, THETA, AnglD, Degrees )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III            NOAA/NMC |
    !/                  |             A. Saulter            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         01-Mar-2018 |
    !/                  +-----------------------------------+
    !/
    !/    01-Mar-2018 : Added subroutine                   ( version 6.02 )
    !
    !  1. Purpose :
    !     Subroutine to de-rotate directions from rotated to standard pole
    !     reference system
    !
    !  2. Method:
    !   Rotates x,y vectors anticlockwise by angle alpha in radians
    !
    !/ ------------------------------------------------------------------- /
    USE constants, ONLY : dera, tpi, undef
    !
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN) :: NSEA        ! Number of sea points
    REAL,    INTENT(IN) :: AnglD(NSEA) ! Turning angle (degrees)
    LOGICAL, INTENT(IN) :: Degrees     ! Use degrees or radians
    REAL, INTENT(INOUT) :: THETA(NSEA) ! Direction seapoint array
    !
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER :: ISEA
    !
    !/ ------------------------------------------------------------------- /
    ! Apply the rotation
    !
    DO isea=1, nsea
      IF ( theta(isea) .NE. undef ) THEN
        IF ( degrees ) THEN
          theta(isea) = theta(isea) - angld(isea)
          IF ( theta(isea) .LT. 0 ) theta(isea) = theta(isea) + 360.0
        ELSE
          theta(isea) = theta(isea) - angld(isea)*dera
          IF ( theta(isea) .LT. 0 ) theta(isea) = theta(isea) + tpi
        END IF
      ENDIF
    END DO
 
    RETURN
  END SUBROUTINE w3thrtn
  !
  !/ ------------------------------------------------------------------- /
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3xyrtn ( NSEA, XVEC, YVEC, AnglD )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III            NOAA/NMC |
    !/                  |             A. Saulter            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         01-Mar-2018 |
    !/                  +-----------------------------------+
    !/
    !/    01-Mar-2018 : Added subroutine                   ( version 6.02 )
    !
    !  1. Purpose :
    !     Subroutine to de-rotate x,y vectors from rotated to standard pole
    !     reference system
    !
    !  2. Method:
    !   Rotates x,y vectors anticlockwise by angle alpha in radians
    !
    !/ ------------------------------------------------------------------- /
    USE constants, ONLY : dera, tpi, undef
    !
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN) :: NSEA        ! Number of sea points
    REAL,    INTENT(IN) :: AnglD(NSEA) ! Turning angle (degrees)
    REAL, INTENT(INOUT) :: XVEC(NSEA), YVEC(NSEA)
    !
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER :: ISEA
    REAL    :: XVTMP, YVTMP
    !
    !/ ------------------------------------------------------------------- /
    ! Apply the rotation
    !
    DO isea=1, nsea
      IF (( xvec(isea) .NE. undef ) .AND. &
           ( yvec(isea) .NE. undef )) THEN
        xvtmp = xvec(isea)*cos(angld(isea)*dera) + yvec(isea)*sin(angld(isea)*dera)
        yvtmp = yvec(isea)*cos(angld(isea)*dera) - xvec(isea)*sin(angld(isea)*dera)
        xvec(isea) = xvtmp
        yvec(isea) = yvtmp
      END IF
    END DO
 
    RETURN
  END SUBROUTINE w3xyrtn
  !
  !/ ------------------------------------------------------------------- /
  !/ ------------------------------------------------------------------- /
  !/
  SUBROUTINE strsplit(STRING,TAB)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |          M. Accensi               |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         29-Apr-2013 !
    !/                  +-----------------------------------+
    !/
    !/    29-Mar-2013 : Origination.                        ( version 4.10 )
    !/
    !  1. Purpose :
    !
    !     Splits string into words
    !
    !  2. Method :
    !
    !     finds spaces and loops
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       STRING   Str   O   String to be splitted
    !       TAB      Str   O   Array of strings
    !     ----------------------------------------------------------------
    !
    CHARACTER(LEN=*), intent(IN)         :: STRING
    CHARACTER(LEN=100), intent(INOUT)    :: TAB(*)
    INTEGER                              :: cnt, I
    CHARACTER(LEN=1024)                  :: tmp_str, ori_str
 
    ! initializes arrays
    ori_str=adjustl(trim(string))
    tmp_str=ori_str
    cnt=0
 
    ! counts the number of substrings
    DO WHILE ((index(tmp_str,' ').NE.0) .AND. (len_trim(tmp_str).NE.0))
      tmp_str=adjustl(tmp_str(index(tmp_str,' ')+1:))
      cnt=cnt+1
    ENDDO
    !
    ! reinitializes arrays
    !
    tmp_str=ori_str
    ! loops on each substring
    DO i=1,cnt
      tab(i)=tmp_str(:index(tmp_str,' '))
      tmp_str=adjustl(tmp_str(index(tmp_str,' ')+1:))
    END DO
 
    RETURN
    !/
    !/ End of STRSPLIT ----------------------------------------------------- /
    !/
  END SUBROUTINE strsplit
  !/
 
  !/ ------------------------------------------------------------------- /
  SUBROUTINE str_to_upper(STR)
    character(*), intent(inout) :: str
    integer :: i
 
    DO i = 1, len(str)
      select case(str(i:i))
      case("a":"z")
        str(i:i) = achar(iachar(str(i:i))-32)
      end select
    END DO
    !/ End of STR_TO_UPPER
    !/ ------------------------------------------------------------------- /
  END SUBROUTINE str_to_upper
 
  !**********************************************************************
  !*                                                                    *
#ifdef W3_T
  !**********************************************************************
  SUBROUTINE ssort1 (X, Y, N, KFLAG)
    !***BEGIN PROLOGUE  SSORT
    !***PURPOSE  Sort an array and optionally make the same interchanges in
    !            an auxiliary array.  The array may be sorted in increasing
    !            or decreasing order.  A slightly modified QUICKSORT
    !            algorithm is used.
    !***LIBRARY   SLATEC
    !***CATEGORY  N6A2B
    !***TYPE      SINGLE PRECISION (SSORT-S, DSORT-D, ISORT-I)
    !***KEYWORDS  SINGLETON QUICKSORT, SORT, SORTING
    !***AUTHOR  Jones, R. E., (SNLA)
    !           Wisniewski, J. A., (SNLA)
    !***DESCRIPTION
    !
    !   SSORT sorts array X and optionally makes the same interchanges in
    !   array Y.  The array X may be sorted in increasing order or
    !   decreasing order.  A slightly modified quicksort algorithm is used.
    !
    !   Description of Parameters
    !      X - array of values to be sorted   (usually abscissas)
    !      Y - array to be (optionally) carried along
    !      N - number of values in array X to be sorted
    !      KFLAG - control parameter
    !            =  2  means sort X in increasing order and carry Y along.
    !            =  1  means sort X in increasing order (ignoring Y)
    !            = -1  means sort X in decreasing order (ignoring Y)
    !            = -2  means sort X in decreasing order and carry Y along.
    !
    !***REFERENCES  R. C. Singleton, Algorithm 347, An efficient algorithm
    !                 for sorting with minimal storage, Communications of
    !                 the ACM, 12, 3 (1969), pp. 185-187.
    !***REVISION HISTORY  (YYMMDD)
    !   761101  DATE WRITTEN
    !   761118  Modified to use the Singleton quicksort algorithm.  (JAW)
    !   890531  Changed all specific intrinsics to generic.  (WRB)
    !   890831  Modified array declarations.  (WRB)
    !   891009  Removed unreferenced statement labels.  (WRB)
    !   891024  Changed category.  (WRB)
    !   891024  REVISION DATE from Version 3.2
    !   891214  Prologue converted to Version 4.0 format.  (BAB)
    !   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
    !   901012  Declared all variables; changed X,Y to SX,SY. (M. McClain)
    !   920501  Reformatted the REFERENCES section.  (DWL, WRB)
    !   920519  Clarified error messages.  (DWL)
    !   920801  Declarations section rebuilt and code restructured to use
    !           IF-THEN-ELSE-ENDIF.  (RWC, WRB)
    !***END PROLOGUE  SSORT
    !     .. Scalar Arguments ..
    INTEGER KFLAG, N
    !     .. Array Arguments ..
    real*4 x(*), y(*)
    !     .. Local Scalars ..
    real*4 r, t, tt, tty, ty
    INTEGER I, IJ, J, K, KK, L, M, NN
    !     .. Local Arrays ..
    INTEGER IL(21), IU(21)
    !     .. External Subroutines ..
    !     None
    !     .. Intrinsic Functions ..
    INTRINSIC abs, int
    !***FIRST EXECUTABLE STATEMENT  SSORT
    nn = n
    IF (nn .LT. 1) THEN
      WRITE (*,*) 'The number of values to be sorted is not positive.'
      RETURN
    ENDIF
    !
    kk = abs(kflag)
    IF (kk.NE.1 .AND. kk.NE.2) THEN
      WRITE (*,*) 'The sort control parameter, K, is not 2, 1, -1, or -2.'
      RETURN
    ENDIF
    !
    !     Alter array X to get decreasing order if needed
    !
    IF (kflag .LE. -1) THEN
      DO  i=1,nn
        x(i) = -x(i)
      end do
    ENDIF
    !
    IF (kk .EQ. 2) GO TO 100
    !
    !     Sort X only
    !
    m = 1
    i = 1
    j = nn
    r = 0.375e0
    !
20  IF (i .EQ. j) GO TO 60
    IF (r .LE. 0.5898437e0) THEN
      r = r+3.90625e-2
    ELSE
      r = r-0.21875e0
    ENDIF
    !
30  k = i
    !
    !     Select a central element of the array and save it in location T
    !
    ij = i + int((j-i)*r)
    t = x(ij)
    !
    !     If first element of array is greater than T, interchange with T
    !
    IF (x(i) .GT. t) THEN
      x(ij) = x(i)
      x(i) = t
      t = x(ij)
    ENDIF
    l = j
    !
    !     If last element of array is less than than T, interchange with T
    !
    IF (x(j) .LT. t) THEN
      x(ij) = x(j)
      x(j) = t
      t = x(ij)
      !
      !        If first element of array is greater than T, interchange with T
      !
      IF (x(i) .GT. t) THEN
        x(ij) = x(i)
        x(i) = t
        t = x(ij)
      ENDIF
    ENDIF
    !
    !     Find an element in the second half of the array which is smaller
    !     than T
    !
40  l = l-1
    IF (x(l) .GT. t) GO TO 40
    !
    !     Find an element in the first half of the array which is greater
    !     than T
    !
50  k = k+1
    IF (x(k) .LT. t) GO TO 50
    !
    !     Interchange these elements
    !
    IF (k .LE. l) THEN
      tt = x(l)
      x(l) = x(k)
      x(k) = tt
      GO TO 40
    ENDIF
    !
    !     Save upper and lower subscripts of the array yet to be sorted
    !
    IF (l-i .GT. j-k) THEN
      il(m) = i
      iu(m) = l
      i = k
      m = m+1
    ELSE
      il(m) = k
      iu(m) = j
      j = l
      m = m+1
    ENDIF
    GO TO 70
    !
    !     Begin again on another portion of the unsorted array
    !
60  m = m-1
    IF (m .EQ. 0) GO TO 190
    i = il(m)
    j = iu(m)
    !
70  IF (j-i .GE. 1) GO TO 30
    IF (i .EQ. 1) GO TO 20
    i = i-1
    !
80  i = i+1
    IF (i .EQ. j) GO TO 60
    t = x(i+1)
    IF (x(i) .LE. t) GO TO 80
    k = i
    !
90  x(k+1) = x(k)
    k = k-1
    IF (t .LT. x(k)) GO TO 90
    x(k+1) = t
    GO TO 80
    !
    !     Sort X and carry Y along
    !
100 m = 1
    i = 1
    j = nn
    r = 0.375e0
    !
110 IF (i .EQ. j) GO TO 150
    IF (r .LE. 0.5898437e0) THEN
      r = r+3.90625e-2
    ELSE
      r = r-0.21875e0
    ENDIF
    !
120 k = i
    !
    !     Select a central element of the array and save it in location T
    !
    ij = i + int((j-i)*r)
    t = x(ij)
    ty = y(ij)
    !
    !     If first element of array is greater than T, interchange with T
    !
    IF (x(i) .GT. t) THEN
      x(ij) = x(i)
      x(i) = t
      t = x(ij)
      y(ij) = y(i)
      y(i) = ty
      ty = y(ij)
    ENDIF
    l = j
    !
    !     If last element of array is less than T, interchange with T
    !
    IF (x(j) .LT. t) THEN
      x(ij) = x(j)
      x(j) = t
      t = x(ij)
      y(ij) = y(j)
      y(j) = ty
      ty = y(ij)
      !
      !        If first element of array is greater than T, interchange with T
      !
      IF (x(i) .GT. t) THEN
        x(ij) = x(i)
        x(i) = t
        t = x(ij)
        y(ij) = y(i)
        y(i) = ty
        ty = y(ij)
      ENDIF
    ENDIF
    !
    !     Find an element in the second half of the array which is smaller
    !     than T
    !
130 l = l-1
    IF (x(l) .GT. t) GO TO 130
    !
    !     Find an element in the first half of the array which is greater
    !     than T
    !
140 k = k+1
    IF (x(k) .LT. t) GO TO 140
    !
    !     Interchange these elements
    !
    IF (k .LE. l) THEN
      tt = x(l)
      x(l) = x(k)
      x(k) = tt
      tty = y(l)
      y(l) = y(k)
      y(k) = tty
      GO TO 130
    ENDIF
    !
    !     Save upper and lower subscripts of the array yet to be sorted
    !
    IF (l-i .GT. j-k) THEN
      il(m) = i
      iu(m) = l
      i = k
      m = m+1
    ELSE
      il(m) = k
      iu(m) = j
      j = l
      m = m+1
    ENDIF
    GO TO 160
    !
    !     Begin again on another portion of the unsorted array
    !
150 m = m-1
    IF (m .EQ. 0) GO TO 190
    i = il(m)
    j = iu(m)
    !
160 IF (j-i .GE. 1) GO TO 120
    IF (i .EQ. 1) GO TO 110
    i = i-1
    !
170 i = i+1
    IF (i .EQ. j) GO TO 150
    t = x(i+1)
    ty = y(i+1)
    IF (x(i) .LE. t) GO TO 170
    k = i
    !
180 x(k+1) = x(k)
    y(k+1) = y(k)
    k = k-1
    IF (t .LT. x(k)) GO TO 180
    x(k+1) = t
    y(k+1) = ty
    GO TO 170
    !
    !     Clean up
    !
190 IF (kflag .LE. -1) THEN
      DO  i=1,nn
        x(i) = -x(i)
      end do
    ENDIF
    RETURN
  END SUBROUTINE ssort1
 
#endif
 
  !*********************************************************************
  SUBROUTINE diagonalize(a1,d,v,nrot)
    !*********************************************************************
    INTEGER,                          INTENT(out)   :: nrot
    DOUBLE PRECISION, DIMENSION(:)  , INTENT(OUT)   ::d
    DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN)    ::a1  ! Modified from INOUT to IN by F.A. on 2018/01/21
    DOUBLE PRECISION, DIMENSION(:,:), INTENT(OUT)   ::v
 
    INTEGER    i,j,ip,iq,n
    DOUBLE PRECISION       c,g,h,s,sm,t,tau,theta,tresh
    DOUBLE PRECISION    , DIMENSION(size(d)) ::b,z
    DOUBLE PRECISION, DIMENSION(size(d),size(d)) :: a
    LOGICAL, DIMENSION(size(d),size(d)) :: upper_triangle
 
    a=a1
    n=size(d)
    v(:,:)=0.
    upper_triangle(:,:)=.false.
    DO i=1,n
      v(i,i)=1.
      b(i)=a(i,i)
      DO j=i+1,n
        upper_triangle(i,j)=.true.
      ENDDO
    ENDDO
    d(:)=b(:)
    z(:)=0.0
    nrot=0
    DO i=1,50
      sm=sum(abs(a),mask=upper_triangle)
      IF (sm.EQ.0.0) RETURN
      tresh=merge(0.2*sm/n**2,0.0d0,i<4)
      DO ip=1,n-1
        do iq=ip+1,n
          g=100.0*abs(a(ip,iq))
          IF((i > 4).AND.(abs(d(ip))+g.EQ.abs(d(ip))) &
               .AND.(abs(d(iq))+g.EQ.abs(d(iq)))) THEN
            a(ip,iq)=0.0
          ELSE IF (abs(a(ip,iq)) > tresh) THEN
            h=d(iq)-d(ip)
            if (abs(h)+g == abs(h)) THEN
              t=a(ip,iq)/h
            ELSE
              theta=0.5*h/a(ip,iq)
              t=1.0/(abs(theta)+sqrt(1.0+theta**2))
              IF ( theta < 0.0) t=-t
            ENDIF
            c=1.0/sqrt(1+t**2)
            s=t*c
            tau=s/(1.0+c)
            h=t*a(ip,iq)
            z(ip)=z(ip)-h
            z(iq)=z(iq)+h
            d(ip)=d(ip)-h
            d(iq)=d(iq)+h
            a(ip,iq)=0.0
            IF (ip.GE.1) CALL rotate(a(1:ip-1,ip),a(1:ip-1,iq))
            !The IF test was added by F.A. (2005/04/04) because of the following error:
            !Subscript out of range. Location: line 593 column 36 of 'cb_botsc.f90'
            !Subscript number 1 has value 0 in array 'A'
            CALL rotate(a(ip,ip+1:iq-1),a(ip+1:iq-1,iq))
            CALL rotate(a(ip,iq+1:n),a(iq,iq+1:n))
            CALL rotate(v(:,ip),v(:,iq))
            nrot=nrot+1
          ENDIF
        ENDDO
      ENDDO
      b(:)=b(:)+z(:)
      d(:)=b(:)
      z(:)=0.0
    ENDDO
    WRITE(6,*) 'Too many iterations in DIAGONALIZE'
  CONTAINS
    SUBROUTINE rotate(X1,X2)
      DOUBLE PRECISION, DIMENSION(:), INTENT(INOUT) :: X1,X2
      DOUBLE PRECISION, DIMENSION(size(X1)) :: MEM
      mem(:)=x1(:)
      x1(:)=x1(:)-s*(x2(:)+x1(:)*tau)
      x2(:)=x2(:)+s*(mem(:)-x2(:)*tau)
    END SUBROUTINE rotate
  END SUBROUTINE diagonalize
 
  !/ ------------------------------------------------------------------- /
  SUBROUTINE uv_to_mag_dir(U, V, NSEA, MAG, DIR, TOLERANCE, CONV)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           C. Bunney               |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         15-Jan-2021 |
    !/                  +-----------------------------------+
    !/
    !/    15-Jan-2021 : Creation                            ( version 7.12 )
    !/
    !  1. Purpose :
    !
    !     Converts seapoint arrays formulated as U/V vectors into magnitude
    !     and direction arrays.
    !
    !     If MAG and DIR input parameters are not specificed then the
    !     conversion is performed in-place (U => MAG, v => DIR).
    !
    !  2. Parameters
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       U/V       R.Arr  I  Array of U/V components
    !       NSEA      Int    I  Number of sea points
    !       MAG       R.Arr  O  Magnitude array            (Optional)
    !       DIR       R.Arr  O  Direction array (degrees)  (Optional)
    !       TOLERANCE Real   I  Minimum allowed magnitude  (Optional)
    !       CONV      Char   I  Ouput direciton convention (Optional)
    !     ----------------------------------------------------------------
    !
    !  3. Remarks
    !
    !     Optional CONV specifies direction convention. Must be one of:
    !       'N'=Nautical     : North=0, clockwise, direction-from (default)
    !       'O'=Oceangraphic : North=0, clockwise, direction-to
    !       'C'=Cartesian    : North=90, counter-clockwise, direction-to
    !
    !/ ------------------------------------------------------------------- /
    USE constants, ONLY: rade, undef
    REAL, INTENT(INOUT)             :: U(NSEA), V(NSEA)
    INTEGER, INTENT(IN)             :: NSEA
    REAL, INTENT(OUT), OPTIONAL     :: MAG(NSEA), DIR(NSEA)
    REAL, INTENT(IN), OPTIONAL      :: TOLERANCE
    CHARACTER, INTENT(IN), OPTIONAL :: CONV
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !
    REAL :: TOL, SGN, OFFSET, TMP
    CHARACTER :: DIRCONV
    INTEGER :: ISEA
    LOGICAL :: INPLACE
 
    dirconv = 'N'
    tol = 1.0
    inplace = .true.
    IF(PRESENT(tolerance)) tol = tolerance
    IF(PRESENT(conv)) dirconv = conv
    IF(PRESENT(mag) .AND. PRESENT(dir)) inplace = .false.
 
    SELECT CASE (conv)
    CASE('N')
      offset = 630.
      sgn = -1.
    CASE('O')
      offset = 450.
      sgn = -1.
    CASE('C')
      offset = 360.
      sgn = 1.
    CASE DEFAULT
      WRITE(*,*) "UV_TO_MAG_DIR: UNKNOWN DIR CONVENTION: ", dirconv
      CALL extcde(1)
    END SELECT
 
    IF(inplace) THEN
      DO isea=1, nsea
        tmp = sqrt(u(isea)**2 + v(isea)**2)
        IF(tmp .GE. tol) THEN
          v(isea) = mod(offset + (sgn * rade * atan2(v(isea), u(isea))), 360.)
          u(isea) = tmp
        ELSE
          u(isea) = undef
          v(isea) = undef
        END IF
      END DO
    ELSE
      DO isea=1, nsea
        mag(isea) = sqrt(u(isea)**2 + v(isea)**2)
        IF(mag(isea) .GE. tol) THEN
          dir(isea) = mod(offset + (sgn * rade * atan2(v(isea), u(isea))), 360.)
        ELSE
          mag(isea) = undef
          dir(isea) = undef
        END IF
      END DO
    ENDIF
 
  END SUBROUTINE uv_to_mag_dir
 
  !========================================================================
 
  subroutine print_memcheck(iun, msg)
#ifdef W3_MEMCHECK
    USE mallocinfo_m
#endif
    integer          , intent(in) :: iun
    character(len=*) , intent(in) :: msg
 
#ifdef W3_MEMCHECK
    ! local variables
    type(mallinfo_t)        :: mallinfos
 
    write(iun,*) trim(msg)
    call getmallocinfo(mallinfos)
    call printmallinfo(iun, mallinfos)
#endif
  end subroutine print_memcheck
  !/
  !/ End of module W3SERVMD -------------------------------------------- /
  !/
END MODULE w3servmd