#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 ! 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 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(:,:) !/ ! ! -------------------------------------------------------------------- / ! 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 ! ! -------------------------------------------------------------------- / ! 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 ! 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 ! 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 ! ! 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 ! ! -------------------------------------------------------------------- / ! 3. Convert ! 3.a Discrete energies ! 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 ! 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 ! 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 ***'/) ! !/ !/ End of W3CSPC ----------------------------------------------------- / !/ END SUBROUTINE W3CSPC !/ !/ End of module W3CSPCMD -------------------------------------------- / !/ END MODULE W3CSPCMD