w3cspcmd.F90

Go to the documentation of this file.

 
 
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
MODULE w3cspcmd
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         01-Nov-2012 |
  !/                  +-----------------------------------+
  !/
  !/    19-Sep-2005 : Origination.                        ( version 3.08 )
  !/    29-May-2009 : Preparing distribution version.     ( version 3.14 )
  !/    01-Nov-2012 : Minor code clean-up (tabs & coments)( version 4.08 )
  !/
  !/    Copyright 2009 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 :
  !
  !     Convert spectra to new discrete spectral grid.
  !
  !  2. Variables and types :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      NCASES    Int.  Private  Number of cases for which interpol.
  !                               data is stored.
  !      IDATA     CASE  Private  Interpolation data.
  !     ----------------------------------------------------------------
  !
  !     Elements of the data structure CASE are given below. The middle
  !     block pf parameters has pointer aliasses with the same name in
  !     the subroutine.
  !
  !      Name      Type  Description
  !     ----------------------------------------------------------------
  !      ICASE     Int.  Number of case.
  !      NFR1, NTH1, NFR2, NTH2, XF1, FR1, TH1, XF2, FR2, TH2
  !                      Same as in parameter list of routine.
  !
  !      DTH1      Real  Directional increment.
  !      DTH2      Real  Directional increment.
  !      IDTH      I.A.  Index information for redistribution of
  !                      energy in direction space.
  !      RDTH      R.A.  Factors corresponding to IDTH.
  !      FRQ1      R.A.  Frequencies.
  !      FRQ2      R.A.  Frequencies.
  !      XDF1      Real  Factor for increments.
  !      XDF2      Real  Factor for increments.
  !      NFR2T     Int.  Frequency to start the tail.
  !      IDFR      I.A.  Idem for frequencies.
  !      RDFR      R.A.  Factors corresponding to IDFR.
  !
  !      NEXT      CASE  Pointer to next data set stored.
  !     ----------------------------------------------------------------
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3CSPC    Subr. Public   Perform conversion for vector of
  !                               spectra.
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !     See subroutine W3CSPC.
  !
  !  5. Remarks :
  !
  !     - Conversion data are sored in an endless linked chain, which
  !       is tested at the beginning of the routine.
  !
  !  6. Switches :
  !
  !     See subroutine.
  !
  !  7. Source code :
  !
  !/ ------------------------------------------------------------------- /
  PUBLIC
  !/
  TYPE case
    INTEGER               :: icase, nfr1, nth1, nfr2, nth2, nfr2t
    REAL                  :: xf1, fr1, th1, xf2, fr2, th2,       &
         dth1, dth2, xdf1, xdf2
    INTEGER, POINTER      :: idth(:,:), idfr(:,:)
    REAL, POINTER         :: rdth(:,:), frq1(:), frq2(:), rdfr(:,:)
    TYPE(case), POINTER   :: next
  END TYPE case
  !/
  INTEGER, PRIVATE        :: ncases = 0
  TYPE(case), PRIVATE, POINTER :: idata
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3cspc ( SP1, NFR1, NTH1, XF1, FR1, TH1,             &
       SP2, NFR2, NTH2, XF2, FR2, TH2,             &
       NSP, NDST, NDSE, FTL )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         01-Nov-2012 !
    !/                  +-----------------------------------+
    !/
    !/    19-Sep-2005 : Origination.                        ( version 3.08 )
    !/    01-Nov-2012 : code clean up (tab spaces, comments)( version 4.08 )
    !/
    !  1. Purpose :
    !
    !     Convert a set of spectra to a new spectral grid.
    !
    !  2. Method :
    !
    !     Conservative distribution of input energies over new grid.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       SP1     R.A.   I   Input spectra.
    !       NFR1    Int.   I   Input number of frequencies.
    !       NTH1    Int.   I   Input number of directions.
    !       XFR     Real   I   Input frequency increment factor.
    !       FR1     Real   I   First input frequency.
    !       TH1     Real   I   First input direction.
    !       SP2     R.A.   O   Output spectra.
    !       NFR2, NTH2, XF2, FR2, TH2
    !                      !   Specral description for output spectra.
    !       NSP     Int.   I   Number of spectra.
    !       NDST    int.   I   Unit number for test output.
    !       NDSE    int.   I   Unit number for error output.
    !       FTAIL   Real   I   Factor for tail description = XF2**N
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Sur.  W3SERVMD Subroutine tracing.
    !      EXTCDE    Sur.    Id     program abort.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3IOBC    Subr. W3IOBCMD Updating boundary conditions.
    !                Subr           Multi scale model bound. data input.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !     - Check on input parameters.
    !
    !  7. Remarks :
    !
    !     - The inner loop of the actual redistribution is over the
    !       individual spectra, optimizing this routine for large numbers
    !       of conversions in a single call.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S    Enable subroutine tracing.
    !
    !     !/T    Enable test output.
    !     !/T1   Test output for searching in stored data.
    !     !/T2   Test output for redistribution data.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    !
    USE w3servmd, ONLY: extcde
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: NSP, NFR1, NTH1, NFR2, NTH2, NDST, NDSE
    REAL, INTENT(IN)        :: SP1(NTH1,NFR1,NSP), XF1, FR1, TH1,   &
         XF2, FR2, TH2, FTL
    REAL, INTENT(OUT)       :: SP2(NTH2,NFR2,NSP)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: I, NRMAX, J, I1, L1, J1, I2, L2, J2, &
         ISP
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: LOW, HGH, RLOW, RHGH, BLOW, BHGH,    &
         FRAC, AUX1, AUX2, R1, R2, FACT
    LOGICAL                 :: FOUND
    TYPE(case), POINTER     :: CURRENT
    !/
    !/ ------------------------------------------------------------------- /
    !/ Pointers for aliases
    !/
    INTEGER, POINTER        :: IDTH(:,:), IDFR(:,:), NFR2T
    REAL, POINTER           :: DTH1, DTH2, RDTH(:,:), FRQ1(:),      &
         FRQ2(:), XDF1, XDF2, RDFR(:,:)
    !/
#ifdef W3_S
    CALL strace (ient, 'W3CSPC')
#endif
    !
    ! -------------------------------------------------------------------- /
    ! 0.  Initializations
    ! 0.a Check input
    !
    IF ( nfr1.LT.3 .OR. nth1.LT.4 .OR. xf1.LE.1. .OR. fr1.LE.0. .OR.&
         nfr2.LT.3 .OR. nth2.LT.4 .OR. xf2.LE.1. .OR. fr2.LE.0. ) THEN
      WRITE (ndse,900) nfr1, nth1, xf1, fr1, nfr2, nth2, xf2, fr2
      CALL extcde ( 1 )
    END IF
    !
    IF ( nsp .LT. 0 ) THEN
      WRITE (ndse,901)
      CALL extcde ( 2 )
    END IF
    !
    IF ( nsp .EQ. 0 ) THEN
      WRITE (ndse,902)
      RETURN
    END IF
    !
    ! 0.b Test output
    !
#ifdef W3_T
    WRITE (ndst,9000) nsp, nfr1, nth1, xf1, fr1, th1*rade,       &
         nfr2, nth2, xf2, fr2, th2*rade, ftl
#endif
    !
    ! -------------------------------------------------------------------- /
    ! 1.  Search stored interpolation data for match
    !
    found  = .false.
    !
    DO i=1, ncases
      !
      IF ( i .EQ. 1 ) THEN
        current => idata
      ELSE
        current => current%NEXT
      END IF
      !
#ifdef W3_T1
      WRITE (ndst,9010) i, current%NFR1, current%NTH1,          &
           current%XF1, current%FR1, current%TH1*rade,  &
           current%NFR2, current%NTH2,          &
           current%XF2, current%FR2, current%TH2*rade
#endif
      !
      found = current%NFR1.EQ.nfr1 .AND. current%NFR2.EQ.nfr2 .AND. &
           current%NTH1.EQ.nth1 .AND. current%NTH2.EQ.nth2 .AND. &
           current%XF1 .EQ.xf1  .AND. current%XF2 .EQ.xf2  .AND. &
           current%FR1 .EQ.fr1  .AND. current%FR2 .EQ.fr2  .AND. &
           current%TH1 .EQ.th1  .AND. current%TH2 .EQ.th2
      IF ( found ) EXIT
      !
    END DO
    !
    ! -------------------------------------------------------------------- /
    ! 2.  Link or compute interpolation data
    ! 2.a Link
    !
    IF ( found ) THEN
      !
#ifdef W3_T
      WRITE (ndst,9020) i
#endif
      !
      dth1   => current%DTH1
      dth2   => current%DTH2
      idth   => current%IDTH
      rdth   => current%RDTH
      !
      frq1   => current%FRQ1
      frq2   => current%FRQ2
      xdf1   => current%XDF1
      xdf2   => current%XDF2
      nfr2t  => current%NFR2T
      idfr   => current%IDFR
      rdfr   => current%RDFR
      !
      ! 2.b Compute
      !
    ELSE
      !
      ncases = ncases + 1
#ifdef W3_T
      WRITE (ndst,9021) ncases
#endif
      !
      ! 2.b.1 Point and allocate as necessary
      !
      IF ( ncases .EQ. 1 ) THEN
        ALLOCATE ( idata )
        current => idata
      ELSE
        ALLOCATE ( current%NEXT )
        current => current%NEXT
      END IF
      !
      ! 2.b.2 Store test data
      !
      current%ICASE = ncases
      current%NFR1  = nfr1
      current%NTH1  = nth1
      current%XF1   = xf1
      current%FR1   = fr1
      current%TH1   = th1
      current%NFR2  = nfr2
      current%NTH2  = nth2
      current%XF2   = xf2
      current%FR2   = fr2
      current%TH2   = th2
      !
      ! 2.b.3 Directional redistribution data
      !
      dth1   => current%DTH1
      dth1   = tpi / real(nth1)
      dth2   => current%DTH2
      dth2   = tpi / real(nth2)
      !
      IF ( dth1 .LE. dth2 ) THEN
        nrmax  = 2
      ELSE
        nrmax  = 2 + int(dth1/dth2)
      END IF
      !
      ALLOCATE (current%IDTH(0:nrmax,nth1),current%RDTH(nrmax,nth1))
      idth   => current%IDTH
      rdth   => current%RDTH
      idth   = 0
      rdth   = 0.
      !
      DO i=1, nth1
        low    = th1 + real(i-1)*dth1 - 0.5*dth1
        hgh    = low + dth1
        rlow   = 1. + (low-th2)/dth2
        rhgh   = 1. + (hgh-th2)/dth2
        DO j=nint(rlow), nint(rlow)+nrmax-1
          blow   = th2 + real(j-1)*dth2 - 0.5*dth2
          bhgh   = blow + dth2
          frac   = (min(bhgh,hgh)-max(blow,low)) / (hgh-low)
          IF ( frac .GT. 1.e-5 ) THEN
            idth(0,i) = idth(0,i) + 1
            idth(idth(0,i),i) = 1 + mod(j-1+nth2,nth2)
            rdth(idth(0,i),i) = frac
          END IF
        END DO
      END DO
      !
      ! 2.b.4 Frequency redistribution data
      !
      ALLOCATE ( current%FRQ1(nfr1), current%FRQ2(nfr2) )
      frq1   => current%FRQ1
      frq2   => current%FRQ2
      !
      frq1(1) = fr1
      DO i=2, nfr1
        frq1(i) = xf1 * frq1(i-1)
      END DO
      !
      frq2(1) = fr2
      DO i=2, nfr2
        frq2(i) = xf2 * frq2(i-1)
      END DO
      !
      xdf1   => current%XDF1
      xdf1   = 0.5 * ( xf1 - 1./xf1 )
      xdf2   => current%XDF2
      xdf2   = 0.5 * ( xf2 - 1./xf2 )
      !
      IF ( xdf1 .LE. xdf2 ) THEN
        nrmax  = 2
      ELSE
        nrmax  = 1
        aux1   = xdf1
        aux2   = xdf2
        DO
          nrmax  = nrmax + 1
          aux1   = aux1 - aux2
          aux2   = aux2 / xf2
          IF ( aux1 .LT. 0. ) EXIT
        END DO
      END IF
      !
      ALLOCATE (current%IDFR(0:nrmax,nfr1),current%RDFR(nrmax,nfr1))
      idfr   => current%IDFR
      rdfr   => current%RDFR
      idfr   = 0
      rdfr   = 0.
      !
      DO i=1, nfr1
        IF ( i .EQ. 1 ) THEN
          hgh    = 0.5 * ( frq1(i) + frq1(i+1) )
          low    = hgh - xdf1*frq1(i)
        ELSE
          low    = 0.5 * ( frq1(i) + frq1(i-1) )
          hgh    = low + xdf1*frq1(i)
        END IF
        DO j=1, nfr2
          IF ( j .EQ. 1 ) THEN
            bhgh   = 0.5 * ( frq2(j) + frq2(j+1) )
            blow   = bhgh - xdf2*frq2(j)
          ELSE
            blow   = 0.5 * ( frq2(j) + frq2(j-1) )
            bhgh   = blow + xdf2*frq2(j)
          END IF
          IF ( bhgh .LE. low ) cycle
          IF ( blow .GE. hgh ) EXIT
          frac   = (min(bhgh,hgh)-max(blow,low)) / (hgh-low)
          IF ( frac .LT. 1.e-5 ) cycle
          idfr(0,i) = idfr(0,i) + 1
          idfr(idfr(0,i),i) = j
          rdfr(idfr(0,i),i) = frac
        END DO
      END DO
      !
      nfr2t  => current%NFR2T
      nfr2t  = nfr2 + 1
      DO j=nfr2, 1, -1
        IF ( j .EQ. 1 ) THEN
          bhgh   = 0.5 * ( frq2(j) + frq2(j+1) )
        ELSE
          blow   = 0.5 * ( frq2(j) + frq2(j-1) )
          bhgh   = blow + xdf2*frq2(j)
        END IF
        IF ( bhgh .GT. hgh ) THEN
          nfr2t  = j
        ELSE
          EXIT
        END IF
      END DO
      !
    END IF
    !
    ! 2.c Test output
    !
#ifdef W3_T2
    WRITE (ndst,9022)
    DO i=1, nth1
      WRITE (ndst,9024) i, idth(0,i),                           &
           (idth(j,i),rdth(j,i),j=1,idth(0,i))
    END DO
    WRITE (ndst,9023) nfr2t
    DO i=1, nfr1
      WRITE (ndst,9024) i, idfr(0,i),                           &
           (idfr(j,i),rdfr(j,i),j=1,idfr(0,i))
    END DO
#endif
    !
    ! -------------------------------------------------------------------- /
    ! 3.  Convert
    ! 3.a Discrete energies
    !
#ifdef W3_T
    WRITE (ndst,9030)
#endif
    !
    sp2    = 0.
    !
    DO i2=1, nfr1
      DO l2=1, idfr(0,i2)
        j2   = idfr(l2,i2)
        r2   = rdfr(l2,i2)
        DO i1=1,nth1
          DO l1=1, idth( 0,i1)
            j1   = idth(l1,i1)
            r1   = rdth(l1,i1)
            frac   = r2 * frq1(i2) * xdf1 * r1 * dth1
            sp2(j1,j2,:) = sp2(j1,j2,:) + frac * sp1(i1,i2,:)
          END DO
        END DO
      END DO
    END DO
    !
    ! 3.b Energy densities
    !
#ifdef W3_T
    WRITE (ndst,9031)
#endif
    !
    DO j2=1, nfr2
      DO j1=1, nth2
        fact   = 1. / ( frq2(j2) * xdf2 * dth2 )
        sp2(j1,j2,:) = fact * sp2(j1,j2,:)
      END DO
    END DO
    !
    ! 3.c Add the tail
    !
#ifdef W3_T
    WRITE (ndst,9032)
#endif
    !
    DO j2=nfr2t, nfr2
      sp2(:,j2,:) = ftl * sp2(:,j2-1,:)
    END DO
    !
    RETURN
    !
    ! Formats
    !
900 FORMAT (/' *** ERROR W3CSPC: ILLEGAL INPUT PARAMETERS ***'/     &
         '                   INPUT  : ',2i8,2f10.4/             &
         '                   OUTPUT : ',2i8,2f10.4)
901 FORMAT (/' *** ERROR W3CSPC: NEGATIVE NUMBER OF SPECTRA ***'/)
902 FORMAT (/' *** WARNING W3CSPC: NO SPECTRA ***'/)
    !
#ifdef W3_T
9000 FORMAT ( ' TEST W3CSPC : NR. OF SPECTRA : ',i8/             &
         '               INPUT SPECTRA  : ',2i4,2f8.4,f6.1/ &
         '               OUTPUT SPECTRA : ',2i4,2f8.4,f6.1/ &
         '               TAIL FACTOR    : ',f8.5)
#endif
    !
#ifdef W3_T1
9010 FORMAT ( ' TEST W3CSPC : TEST INFO CASE : ',i8/             &
         '               INPUT SPECTRA  : ',2i4,2f8.4,f6.1/ &
         '               OUTPUT SPECTRA : ',2i4,2f8.4,f6.1)
#endif
    !
#ifdef W3_T
9020 FORMAT ( ' TEST W3CSPC : USING STORED DATA FOR CASE',i4)
9021 FORMAT ( ' TEST W3CSPC : COMPUTING DATA FOR CASE',i4)
#endif
#ifdef W3_T2
9022 FORMAT ( ' TEST W3CSPC : DIRECTIONAL DISTRIBUTION DATA')
9023 FORMAT ( ' TEST W3CSPC : FREQUENCY DISTRIBUTION DATA, ',    &
         'TAIL AT',i4)
9024 FORMAT ( '           ',i4,i4,' :',10(i4,f5.2) )
#endif
    !
#ifdef W3_T
9030 FORMAT ( ' TEST W3CSPC : STARTING CONVERSION')
9031 FORMAT ( ' TEST W3CSPC : ENERGIES TO DENSITIES')
9032 FORMAT ( ' TEST W3CSPC : ADD TAIL')
#endif
    !/
    !/ End of W3CSPC ----------------------------------------------------- /
    !/
  END SUBROUTINE w3cspc
  !/
  !/ End of module W3CSPCMD -------------------------------------------- /
  !/
END MODULE w3cspcmd