text:
C****************************************************************************
C
C WAVEPACK:  routines to compute the wavelet transform of a time series,
C            and significance levels.
C
C Written by:   Christopher Torrence and Gilbert P. Compo
C
C Available from:  http://paos.colorado.edu/research/wavelets/
C
C Requires the following packages:   CFFTPACK, CHISQR
C
C Notes:
C
C  (1) All routines are written in single precision (DOUBLE PRECISION),
C      except for CHISQR, which requires double precision input.
C      Single precision should be sufficient for most applications.
C
C  (2) The CFFTPACK and CHISQR routines were not written by us,
C      and no guarentees are made as to their reliability or efficiency.
C
C  (3) No provision is made for output to files or to graphics packages.
C      The user is expected to call these routines from within
C      their own programs. See sample program "wavetest.f".
C
C  (4) Little error checking is done. Check your input carefully.
C      The programs are not completely ANSI compatible, so use caution.
C
C  (5) Time series are currently limited to 65535 points.
C
C
C Reference: Torrence, C. and G. P. Compo, 1998: A Practical Guide to
C            Wavelet Analysis. Bull. Amer. Meteor. Soc., 79, 61-78.
C
C Notice: Please acknowledge the use of this software in any publications:
C    ``Wavelet software was provided by C. Torrence and G. Compo,
C      and is available at URL: http://paos.colorado.edu/research/wavelets/''.
C
C  Copyright (C) 1998, Christopher Torrence
C This software may be used, copied, or redistributed as long as it is not
C sold and this copyright notice is reproduced on each copy made.  This
C routine is provided as is without any express or implied warranties
C whatsoever.
C
C
C Modified: November 1999 by Arjan van Dijk to include IMPLICIT NONE and
C           to convert all routines to DOUBLE precision.
C****************************************************************************

C****************************************************************************
C WAVELET: computes the wavelet transform of a time series,
C          with appropriate parameters.
C
C
C INPUTS:
C
C  n [INT] = the number of points in "y".
C
C  y [DOUBLE PRECISION] = the time series of length "n".
C
C  dt [DOUBLE PRECISION] = amount of time between each Y value, i.e. the sampling time.
C
C  mother [INT] = An integer giving the mother wavelet to use.
C                   0='Morlet'
C                   1='Paul'
C                   2='DOG' (derivative of Gaussian)
C               If (mother<0 or >2) then default is 'Morlet'.
C
C  param [DOUBLE PRECISION] = mother wavelet parameter. If <0 then default is used.
C            For 'Morlet' this is k0 (wavenumber), default is 6.
C            For 'Paul' this is m (order), default is 4.
C            For 'DOG' this is m (m-th derivative), default is 2.
C
C
C  s0 [DOUBLE PRECISION] = the smallest scale of the wavelet.  Typically = 2*dt.
C         Note: for accurate reconstruction and variance computation
C             set s0=dt for Morlet; s0=dt/4 for Paul
C
C  dj [DOUBLE PRECISION] = the spacing between discrete scales. Typically = 0.25.
C         A smaller # will give better scale resolution, but be slower.
C
C  jtot [INT] = the # of scales.
C             Scales range from s0 up to s0*2^[(jtot-1)*dj],
C             Typically jtot=1+(LOG2(n dt/s0))/dj
C
C  npad [INT] = the total number of points (including padding) to
C             use for the wavelet transform. Typically this is some
C             power of 2. It must be greater or equal to "n".
C             If npad>n, then zeroes are padded onto the end
C             of the time series.
C
C
C OUTPUTS:
C
C  wave [DCMPLX(n,jtot)] = 2D array of the real & imaginary parts
C                 of the wavelet transform, versus time & scale.
C                 CABS(wave) gives the WAVELET amplitude,
C                 ATAN2(AIMAG(wave),DBLE(wave)) gives WAVELET phase.
C                 The WAVELET power spectrum is CABS(wave)**2.
C
C  scale [DOUBLE PRECISION(jtot)] = the wavelet scales that were used.
C
C  period [DOUBLE PRECISION(jtot)] = the "Fourier" periods (in time units) corresponding
C            to "scale".
C
C  coi [DOUBLE PRECISION(n)] = the e-folding factor used for the cone of influence.
C
C
C REQUIRES:   WAVE_FUNCTION, CFFTPACK
C
C  Copyright (C) 1998, Christopher Torrence and Gilbert P. Compo
C This software may be used, copied, or redistributed as long as it is not
C sold and this copyright notice is reproduced on each copy made.  This
C routine is provided as is without any express or implied warranties
C whatsoever.

      SUBROUTINE WAVELET (n,y,dt,mother,param,s0,dj,jtot,npad,
     &             wave,scale,period,coi)
      IMPLICIT none

      INTEGER n,mother,jtot,npad
      DOUBLE PRECISION y(n),dt,param,s0,dj,scale(jtot),period(jtot),
     &  coi(n)
      DOUBLE COMPLEX wave(n,jtot)

      INTEGER i,j,k,nk
      DOUBLE PRECISION ymean,freq1,pi,period1,coi1

C** initialize work arrays
      PARAMETER (nk=65535)
      DOUBLE PRECISION wsave(4*nk+15),kwave(nk)
      DOUBLE COMPLEX yfft(nk),daughter(nk)

      pi = 4.D0*ATAN(1.D0)

      IF (npad.LT.n) THEN
        PRINT*,'**WAVELET: "npad" must be greater than or equal to "n"'
        RETURN
      END IF

      IF ((mother.LT.0).OR.(mother.GT.2)) mother = 0

C** find the time-series mean & remove it
      ymean = 0.D0
      DO 10 i=1,n
        ymean = ymean + y(i)
10    CONTINUE
      ymean = ymean/n
      DO 15 i=1,n
        yfft(i) = y(i) - ymean
15    CONTINUE

C** if desired, pad with extra zeroes
      DO 20 i=n+1,npad
        yfft(i) = 0.D0
20    CONTINUE

C** find the FFT of the time series [Eqn(3)]
      CALL CFFTI(npad,wsave)
      CALL CFFTF(npad,yfft,wsave)
      DO 30 k=1,npad
        yfft(k) = yfft(k)/npad
30    CONTINUE

C** construct the wavenumber array [Eqn(5)]
      freq1 = 2.D0*pi/(DBLE(npad)*dt)
      kwave(1) = 0.D0
      DO 40 i=2,npad/2+1
        kwave(i) = (DBLE(i)-1.D0)*freq1
40    CONTINUE
      DO 50 i=npad/2+2,npad
        kwave(i) = -kwave(npad-i+2)
50    CONTINUE

c**----- Main wavelet loop
c      PRINT '(A8,2A12,/,5X,27("-"))','j','scale','period'
      DO 100 j=1,jtot
        scale(j) = s0*(2.D0**(DBLE(j-1)*dj))
        CALL WAVE_FUNCTION(npad,dt,mother,param,scale(j),
     &            kwave,period1,coi1,daughter)
        period(j) = period1
C**    multiply the daughter by the time-series FFT
        DO 60 k=1,npad
          daughter(k) = daughter(k)*yfft(k)
60      CONTINUE
C**    inverse FFT [Eqn(4)]
        CALL CFFTB(npad,daughter,wsave)
C**    store the wavelet transform, discard zero-padding at end
        DO 70 i=1,n
          wave(i,j) = daughter(i)
70      CONTINUE
c        PRINT '(I8,2F12.3)',j,scale(j),period(j)
100   CONTINUE
c**----- end loop

C** construct the cone of influence
      DO 110 i=1,(n+1)/2
        coi(i) = coi1*dt*(DBLE(i)-1.D0)
        coi(n-i+1) = coi(i)
110   CONTINUE

      RETURN
      END

C****************************************************************************
C WAVE_FUNCTION: computes the daughter wavelets for a particular
C              wavelet function, with appropriate parameters.
C
C
C INPUTS:
C
C  nk [INT] = the number of points in "kwave"
C
C  dt [DOUBLE PRECISION] = amount of time between each Y value, i.e. the sampling time.
C
C  mother [INT] = An integer giving the mother wavelet to use.
C                   0='Morlet'
C                   1='Paul'
C                   2='DOG' (derivative of Gaussian)
C
C  param [DOUBLE PRECISION] = mother wavelet parameter. If <0 then default is used.
C            For 'Morlet' this is k0 (wavenumber), default is 6.
C            For 'Paul' this is m (order), default is 4.
C            For 'DOG' this is m (m-th derivative), default is 2.
C
C  scale1 [DOUBLE PRECISION] = the wavelet scale used to construct the daughter.
C
C  kwave [DOUBLE PRECISION(n)] = vector of wavenumbers, used to construct daughter.
C
C
C OUTPUTS:
C
C  period1 [DOUBLE PRECISION] = the "Fourier" period (in time units) that corresponds
C            to "scale1".
C
C  coi1 [DOUBLE PRECISION] = the e-folding factor used for the cone of influence.
C
C  daughter [DCMPLX(nk)] = real & imaginary parts of the wavelet function
C                    at "scale1" and "kwave".
C
C
C REQUIRES:   FACTORIAL, CHISQR
C
C
C Reference: Tables 1 & 2 in
C            Torrence, C. and G. P. Compo, 1998: A Practical Guide to
C            Wavelet Analysis. Bull. Amer. Meteor. Soc., 79, 61-78.
C
C  Copyright (C) 1998, Christopher Torrence and Gilbert P. Compo
C This software may be used, copied, or redistributed as long as it is not
C sold and this copyright notice is reproduced on each copy made.  This
C routine is provided as is without any express or implied warranties
C whatsoever.

      SUBROUTINE WAVE_FUNCTION (nk,dt,mother,param,scale1,
     &                  kwave,period1,coi1,daughter)
      IMPLICIT none

      INTEGER nk,mother
      DOUBLE PRECISION dt,kwave(nk),param,scale1,period1,coi1
      DOUBLE COMPLEX daughter(nk),norm

      DOUBLE PRECISION expnt,sk,pi,fourier_factor
      INTEGER k,m,factorial
      DOUBLE PRECISION gamma

      pi = 4.D0*ATAN(1.D0)

      IF (mother.EQ.0) THEN
C*******************************************   Morlet wavelet
        IF (param.LT.0) param = 6.D0
        norm = SQRT(2.D0*pi*scale1/dt)*(pi**(-0.25D0))
        DO 10 k=1,nk/2+1
          expnt = -0.5D0*(scale1*kwave(k) - param)**2
          daughter(k) = DCMPLX(norm*EXP(expnt))
10      CONTINUE
         DO 20 k=nk/2+2,nk
          daughter(k) = DCMPLX(0.D0)
20      CONTINUE
        fourier_factor = (4.D0*pi)/(param + SQRT(2.D0+param**2))
        period1 = scale1*fourier_factor
        coi1 = fourier_factor/SQRT(2.D0)
      ELSE IF (mother.EQ.1) THEN
C*******************************************   Paul wavelet
        IF (param.LT.0) param = 4.D0
        m = INT(param)
        norm = SQRT(2.D0*pi*scale1/dt)*
     &         (2.D0**m/SQRT(DBLE(m*FACTORIAL(2*m-1))))
        DO 30 k=1,nk/2+1
          expnt = -scale1*kwave(k)
          daughter(k) = DCMPLX(norm*(scale1*kwave(k))**m*EXP(expnt))
30      CONTINUE
         DO 40 k=nk/2+2,nk
          daughter(k) = DCMPLX(0.D0)
40      CONTINUE
        fourier_factor = (4.D0*pi)/(2.D0*m + 1.D0)
        period1 = scale1*fourier_factor
        coi1 = fourier_factor*SQRT(2.D0)
      ELSE IF (mother.EQ.2) THEN
C*******************************************   DOG wavelet
        IF (param.LT.0) param = 2.D0
        m = INT(param)
        norm = SQRT(2.D0*pi*scale1/dt)*SQRT(1.D0/GAMMA(m+0.5D0))
        norm = -norm*(DCMPLX(0.D0,1.D0)**m)
        DO 50 k=1,nk
          sk = scale1*kwave(k)
          daughter(k) = norm*(sk**m)*EXP(-0.5D0*sk**2)
50      CONTINUE
        fourier_factor = 2.D0*pi*SQRT(2.D0/(2.D0*m+1.D0))
        period1 = scale1*fourier_factor
        coi1 = fourier_factor/SQRT(2.D0)
      ELSE
        stop
      END IF
      RETURN
      END

C****************************************************************************
C FACTORIAL: compute the factorial (n!) of an integer n
C  Copyright (C) 1998, Christopher Torrence
      FUNCTION FACTORIAL(n)
      IMPLICIT NONE
      INTEGER factorial,n,i

      factorial = 1
      DO 10 i=1,n
        factorial = factorial*i
10    CONTINUE
      END

C****************************************************************************
C WAVE_SIGNIF: computes the significance levels for a wavelet transform.
C
C
C INPUTS:
C
C  isigtest [INT] = 0, 1, or 2.
C
C          If 0, then just do a regular chi-square test,
C              i.e. Eqn (18) from Torrence & Compo.
C          If 1, then do a "time-average" test, i.e. Eqn (23).
C              In this case, DOF(j) should be set to NA, the number
C              of local wavelet spectra that were averaged together
C              at each scale. For the Global Wavelet Spectrum,
C              this would be dof(j)=N-scale(j),
C              where N is the number of points in your time series.
C          If 2, then do a "scale-average" test, i.e. Eqns (25)-(28).
C              In this case, "dof(1)" and "dof(2)" should be set to the
C              smallest (S1) and largest (S2) scales that were averaged
C              together, respectively.
C              e.g. if you scale-averaged scales between 2 and 8,
C                   then dof(1)=2.0 and dof(2)=8.0
C
C
C  n [INT] = the number of points in "y".
C
C  y [DOUBLE PRECISION] = the time series of length "n".
C
C  dt [DOUBLE PRECISION] = amount of time between each Y value, i.e. the sampling time.
C
C  mother [INT] = An integer giving the mother wavelet to use.
C                   0='Morlet'
C                   1='Paul'
C                   2='DOG' (derivative of Gaussian)
C
C  param [DOUBLE PRECISION] = mother wavelet parameter.
C
C  s0 [DOUBLE PRECISION] = the smallest scale of the wavelet.
C
C  dj [DOUBLE PRECISION] = the spacing between discrete scales.
C
C  jtot [INT] = the # of scales.
C
C  scale [DOUBLE PRECISION(jtot)] = the wavelet scales that were used.
C
C  period [DOUBLE PRECISION(jtot)] = the "Fourier" periods corresponding to "scale".
C
C  lag1 [DOUBLE PRECISION] = lag 1 Autocorrelation, used for signif levels.
C              Default is 0.0, which corresponds to white-noise.
C
C  siglvl [DOUBLE PRECISION] = significance level to use. Default is 0.05 (the "5%" level)
C
C  dof [DOUBLE PRECISION(jtot)] = degrees-of-freedom for signif test.
C     IF SIGTEST=0, then the input dof is ignored.
C     IF SIGTEST=1, then dof(j) = NA, the number of times averaged together.
C     IF SIGTEST=2, then dof(1)=S1, dof(2)=S2, the range of scales averaged.
C
C
C OUTPUTS:
C
C  dof [DOUBLE PRECISION(jtot)] = degrees-of-freedom that were actually used.
C     IF SIGTEST=0, then dof(j) = 2 (or 1 for the 'DOG')
C     IF SIGTEST=1, then dof(j) = degrees-of-freedom versus scale.
C     IF SIGTEST=2, then dof(1)=degrees-of-freedom, dof(2...jtot)=0.0
C
C  fft_theor [DOUBLE PRECISION(jtot)] = theoretical red-noise spectrum vs scale.
C     IF SIGTEST=2, then fft_theor(1) = the average spectrum from S1-->S2
C                   fft_theor(2...jtot) = 0.0
C
C  signif [DOUBLE PRECISION(jtot)] = significance levels vs scale.
C     IF SIGTEST=2, then signif(1) = the significance level
C                   signif(2...jtot) = 0.0
C
C  ymean [DOUBLE PRECISION] = the mean of the time series.
C
C  variance [DOUBLE PRECISION] = the variance of the time series.
C
C  Cdelta [DOUBLE PRECISION] = the constant "Cdelta" for the mother wavelet (Table 2).
C
C  psi0[DOUBLE PRECISION] = the constant 'psi(0)' for the mother wavelet (Table 2)
C
C REQUIRES:   CHISQR
C
C  Copyright (C) 1998, Christopher Torrence and Gilbert P. Compo
C This software may be used, copied, or redistributed as long as it is not
C sold and this copyright notice is reproduced on each copy made.  This
C routine is provided as is without any express or implied warranties
C whatsoever.

      SUBROUTINE WAVE_SIGNIF (isigtest,n,y,dt,mother,param,dj,jtot,
     &       scale,period,lag1,siglvl,dof,fft_theor,signif,
     &       ymean,variance,Cdelta,psi0)
      IMPLICIT none

      INTEGER isigtest,n,mother,jtot
      DOUBLE PRECISION y(n),dt,param,dj,scale(jtot),period(jtot)
      DOUBLE PRECISION lag1,siglvl,dof(jtot),fft_theor(jtot),
     &  signif(jtot)
      DOUBLE PRECISION ymean,variance,Cdelta,psi0

      INTEGER i,j,m,status,javg1,javg2,navg
      DOUBLE PRECISION pi,freq1,dofmin,gammafac,dj0,Savg,Smid
      DOUBLE PRECISION fft_theor1
      DOUBLE PRECISION chisqr,p,q,bound

      pi = 4.D0*ATAN(1.D0)

      IF (siglvl.LE.0.) siglvl = 0.05D0
      IF (lag1.LE.0.D0) lag1 = 0.D0

      Cdelta = -1.D0
      gammafac = -1.D0
      dj0 = -1.D0
      psi0 = -1.D0

      IF (mother.EQ.0) THEN
C*******************************************   Morlet wavelet
        dofmin = 2.D0
        IF (param.EQ.6.D0) THEN
          Cdelta = 0.776D0
          gammafac = 2.32D0
          dj0 = 0.60D0
          psi0 = pi**(-0.25D0)
        END IF
      ELSE IF (mother.EQ.1) THEN
C*******************************************   Paul wavelet
        m = INT(param)
        dofmin = 2.D0
        IF (m.EQ.4) THEN
          Cdelta = 1.132D0
          gammafac = 1.17D0
          dj0 = 1.5D0
          psi0 = 1.079D0
        END IF
      ELSE IF (mother.EQ.2) THEN
C*******************************************   DOG wavelet
        m = INT(param)
        dofmin = 1.D0
        IF (m.EQ.2) THEN
          Cdelta = 3.541D0
          gammafac = 1.43D0
          dj0 = 1.4D0
          psi0 = 0.867D0
        ELSE IF (m.EQ.6) THEN
          Cdelta = 1.966D0
          gammafac = 1.37D0
          dj0 = 0.97D0
          psi0 = 0.884D0
        END IF
      ELSE
        stop
      END IF

C** find the time-series variance
      ymean = 0.D0
      DO 10 i=1,n
        ymean = ymean + y(i)
10    CONTINUE
      ymean = ymean/n
      variance = 0.D0
      DO 15 i=1,n
        variance = variance + (y(i) - ymean)**2
15    CONTINUE
      variance = variance/(DBLE(n)) ! - 1.D0)

C** construct theoretical red(white)-noise power spectrum [Eqn(16)]
      DO 20 j=1,jtot
        freq1 = dt/period(j)
        fft_theor(j) = variance*(1.D0-lag1**2)/
     &          (1.D0 - 2.D0*lag1*COS(freq1*2.D0*pi) + lag1**2)
20    CONTINUE

      q = DBLE(siglvl)
      p = 1d0 - q

      IF (isigtest.EQ.0) THEN
C*******************************************   no smoothing, dof=dofmin
C   see Eqn(18)
        DO 30 j=1,jtot
          dof(j) = dofmin
          CALL CDFCHI(2,p,q,chisqr,DBLE(dofmin),status,bound)
          signif(j) = fft_theor(j)*chisqr/dofmin
30      CONTINUE
      ELSE IF (isigtest.EQ.1) THEN
C***********************************   time-averaged, dof depend on scale
        IF (gammafac.LE.0.D0) THEN
        PRINT*,'**WAVE_SIGNIF: "gammafac" undefined for this wavelet'
          RETURN
        END IF
        DO 40 j=1,jtot
            IF (dof(j).LT.1.) dof(j) = 1.D0
C   see Eqn(23)
            dof(j) = dofmin*SQRT(1.D0+(dof(j)*dt/gammafac/scale(j))**2)
            IF (dof(j).LT.dofmin) dof(j) = dofmin
            CALL CDFCHI(2,p,q,chisqr,DBLE(dof(j)),status,bound)
            signif(j) = fft_theor(j)*chisqr/dof(j)
40      CONTINUE
       ELSE IF (isigtest.EQ.2) THEN
C***********************************   scale-averaged, dof depend on scale
        IF (Cdelta.LE.0.) THEN
          PRINT*,'**WAVE_SIGNIF: "Cdelta" and "dj0" '//
     &     'undefined for this wavelet'
          RETURN
        END IF
        javg1 = 0
        javg2 = 0
        DO 50 j=1,jtot
          IF ((scale(j).GE.dof(1)).AND.(javg1.EQ.0)) javg1 = j
          IF (scale(j).LE.dof(2)) javg2 = j
50      CONTINUE
        IF ((javg1.EQ.0).OR.(javg2.EQ.0).OR.(javg1.GT.javg2)) THEN
          PRINT*,'**WAVE_SIGNIF: Scales in "dof(1)" & "dof(2)" '//
     &     'are out of range.'
          RETURN
        END IF
        navg = javg2 - javg1 + 1
C   see Eqn(25)
        Savg = 0.D0
        DO 60 j=javg1,javg2
          Savg = Savg + 1.D0/scale(j)
60      CONTINUE
        Savg = 1.D0/Savg
C   see Eqn(27)
        fft_theor1 = 0.D0
        DO 70 j=javg1,javg2
          fft_theor1 = fft_theor1 + fft_theor(j)/scale(j)
70      CONTINUE
        fft_theor(1) = Savg*fft_theor1
C   see Eqn(28)
        Smid = EXP(0.5D0*(LOG(scale(javg1)) + LOG(scale(javg2))))
        dof(1) = (dofmin*navg*Savg/Smid)*SQRT(1 + (navg*dj/dj0)**2)
C   see Eqn(26)
        CALL CDFCHI(2,p,q,chisqr,DBLE(dof(1)),status,bound)
        signif(1)=(dj*dt/Cdelta/Savg)*fft_theor(1)*chisqr/dof(1)
        DO 80 j=2,jtot
          dof(j) = 0.D0
          fft_theor(j) = 0.D0
          signif(j) = 0.D0
80      CONTINUE
      ELSE
        stop
      END IF

      END


text:
	SUBROUTINE AF_IFPT  ( rimnem, cimnem, iret )
C************************************************************************
C* AF_IFPT								*
C*									*
C* This subroutine prints the interface values arrays to the decoder	*
C* log if the verbosity level is 3 or higher.				*
C*									*
C* AF_IFPT  ( RIMNEM, CIMNEM, IRET )					*
C*									*
C* Input parameters:							*
C*	RIMNEM (*)	CHAR*		Interface mnemonics for reals   *
C*	CIMNEM (*)	CHAR*		Interface mnemonics for chars   *
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code 			*
C*					  0 = normal return 		*
C*									*
C**									*
C* Log:									*
C* J. Ator/NP12		08/97						*
C* J. Ator/NCEP		03/98	Use DC_IFPC and DC_IFPR			*
C* D. Kidwell/NCEP	10/98	Added intf mnemonics to calling sequence*
C* D. Kidwell/NCEP	 6/99	Added more icing, turb, sky cover parms *
C* D. Kidwell/NCEP	 8/99	(irhocb,irhoct) -> (irhcbf,irhctf)      *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'afcmn.cmn'
C*
	CHARACTER*(*)	rimnem (*), cimnem (*)
C-----------------------------------------------------------------------
	iret = 0
C
C*	Print out all of the single-level character interface values
C*	that are not "missing".
C
	ncsimn = ( icwcod (1) ) - 1
	CALL DC_IFPC  ( 3, cimnem, civals, ncsimn, ierfpc )
C
C*	Print out all of the single-level real interface values
C*	that are not "missing".
C
	CALL DC_IFPR  ( 3, rimnem, rivals, NRSIMN, ierfpr )
C
C*	Print out all of the multi-level interface values.
C
C*	Turbulence.
C
	IF  ( rivals ( irntrb ) .gt. 0 )  THEN
C
	    WRITE ( UNIT = logmsg, FMT = '( 6A10 )' )
     +		rimnem ( irntrb ),
     +		rimnem ( irdgot (1) ),
     +		rimnem ( irhbot (1) ), rimnem ( irhtot (1) ),
     +		rimnem ( irfqot (1) ), rimnem ( irtpot (1) )
	    CALL DC_WLOG  ( 3, 'AF', 1, logmsg, ierwlg )
	    WRITE ( UNIT = logmsg, FMT = '( 6F10.2 )' )
     +		rivals ( irntrb ),
     +		rivals ( irdgot (1) ),
     +		rivals ( irhbot (1) ), rivals ( irhtot (1) ),
     +		rivals ( irfqot (1) ), rivals ( irtpot (1) )
	    CALL DC_WLOG  ( 3, 'AF', 1, logmsg, ierwlg )
C
	    IF  ( rivals ( irntrb ) .gt. 1 )  THEN
		DO jj = 2, rivals ( irntrb )
		    WRITE ( UNIT = logmsg,
     +			    FMT = '( 10X, 5F10.2 )' )
     +		     rivals ( irdgot ( jj ) ),
     +		     rivals ( irhbot ( jj ) ), rivals ( irhtot ( jj ) ),
     +		     rivals ( irfqot ( jj ) ), rivals ( irtpot ( jj ) )
		    CALL DC_WLOG  ( 3, 'AF', 1, logmsg, ierwlg )
		END DO
	    END IF
C
	END IF
C
C*	Airframe icing.
C
	IF  ( rivals ( irnicg ) .gt. 0 )  THEN
C
	    WRITE ( UNIT = logmsg, FMT = '( 5A10 )' )
     +		rimnem ( irnicg ),
     +		rimnem ( irafic (1) ),
     +		rimnem ( irhboi (1) ), rimnem ( irhtoi (1) ),
     +		rimnem ( irtpoi (1) )
	    CALL DC_WLOG  ( 3, 'AF', 1, logmsg, ierwlg )
	    WRITE ( UNIT = logmsg, FMT = '( 5F10.2 )' )
     +		rivals ( irnicg ),
     +		rivals ( irafic (1) ),
     +		rivals ( irhboi (1) ), rivals ( irhtoi (1) ),
     +		rivals ( irtpoi (1) )
	    CALL DC_WLOG  ( 3, 'AF', 1, logmsg, ierwlg )
C
	    IF  ( rivals ( irnicg ) .gt. 1 )  THEN
		DO jj = 2, rivals ( irnicg )
		    WRITE ( UNIT = logmsg,
     +			    FMT = '( 10X, 4F10.2 )' )
     +		     rivals ( irafic ( jj ) ),
     +		     rivals ( irhboi ( jj ) ), rivals ( irhtoi ( jj ) ),
     +		     rivals ( irtpoi ( jj ) )
		    CALL DC_WLOG  ( 3, 'AF', 1, logmsg, ierwlg )
		END DO
	    END IF
C
	END IF
C
C*	Present weather.
C
	IF  ( rivals ( irnpwx ) .gt. 0 )  THEN
C
	    WRITE ( UNIT = logmsg, FMT = '( 4A10 )' )
     +		rimnem ( irnpwx ),
     +		cimnem ( icwcod (1) ),
     +		rimnem ( irhbwx (1) ), rimnem ( irhtwx (1) )
	    CALL DC_WLOG  ( 3, 'AF', 1, logmsg, ierwlg )
	    WRITE ( UNIT = logmsg, FMT = '( F10.2, A10, 2F10.2 )' )
     +		rivals ( irnpwx ),
     +		civals ( icwcod (1) ),
     +		rivals ( irhbwx (1) ), rivals ( irhtwx (1) )
	    CALL DC_WLOG  ( 3, 'AF', 1, logmsg, ierwlg )
C
	    IF  ( rivals ( irnpwx ) .gt. 1 )  THEN
		DO jj = 2, rivals ( irnpwx )
		    WRITE ( UNIT = logmsg,
     +			    FMT = '( 10X, A10, 2F10.2 )' )
     +		     civals ( icwcod ( jj ) ),
     +		     rivals ( irhbwx ( jj ) ), rivals ( irhtwx ( jj ) )
		    CALL DC_WLOG  ( 3, 'AF', 1, logmsg, ierwlg )
		END DO
	    END IF
C
	END IF
C
C*	Cloud data.
C
	IF  ( rivals ( irncld ) .gt. 0 )  THEN
C
	    WRITE ( UNIT = logmsg, FMT = '( 5A10 )' )
     +		rimnem ( irncld ),
     +		rimnem ( irclam (1) ), rimnem ( ircltp (1) ),
     +		rimnem ( irhcbf (1) ), rimnem ( irhctf (1) )
	    CALL DC_WLOG  ( 3, 'AF', 1, logmsg, ierwlg )
	    WRITE ( UNIT = logmsg, FMT = '( 5F10.2 )' )
     +		rivals ( irncld ),
     +		rivals ( irclam (1) ), rivals ( ircltp (1) ),
     +		rivals ( irhcbf (1) ), rivals ( irhctf (1) )
	    CALL DC_WLOG  ( 3, 'AF', 1, logmsg, ierwlg )
C
	    IF  ( rivals ( irncld ) .gt. 1 )  THEN
		DO jj = 2, rivals ( irncld )
		    WRITE ( UNIT = logmsg,
     +			    FMT = '( 10X, 4F10.2 )' )
     +		     rivals ( irclam ( jj ) ), rivals ( ircltp ( jj ) ),
     +		     rivals ( irhcbf ( jj ) ), rivals ( irhctf ( jj ) )
		    CALL DC_WLOG  ( 3, 'AF', 1, logmsg, ierwlg )
		END DO
	    END IF
C
	END IF
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
c  de-redshift a spectrum by z, in log space.  this means
c  shifting wl scale by -log(1+z).

c do the shift in place, or not
c      subroutine deshiftz(spec,nsp,dw,z,specout)
      subroutine deshiftz(spec,nsp,dwlog,z)

      real spec(nsp)
      real specout(nsp)

      include 'pcredshift.h'

      dzlog = log10(1.+z)
      dindx = dzlog / dwlog
      ishift = int(dindx)
      fshift = dindx - ishift
c now ishift and fshift are the positive integer and fractional
c shifts (if z is positive).  we want to shift to lower w

      do i=1,nsp
         j1 = i + ishift
         if (spec(j1) .gt. bad .and. spec(j1+1) .gt. bad) then
            specout(i) = (1.-fshift)*spec(j1) + fshift*spec(j1+1)
         else if (spec(j1) .gt. bad) then
            specout(i) = spec(j1)
         else if (spec(j1+1) .gt. bad) then
            specout(i) = spec(j1+1)
         else
            specout(i) = bad
         end if
      end do
c to deshift in place
      do i=1,nsp
         spec(i) = specout(i)
      end do

      return
      end


text:
module ttgt_descriptor_module
    use, intrinsic :: iso_fortran_env, only : int64
    use :: tensor_api, only : tensor

    use :: tc_descriptor_module, only : tc_descriptor

    implicit none
    private

    public :: ttgt_descriptor

    type :: ttgt_descriptor
        logical :: switch_ab
        integer, dimension(:), allocatable :: m, n, k
        integer, dimension(:), allocatable :: corder, aorder, border, cinverseorder
        logical :: permute_c, permute_a, permute_b
        integer, dimension(:), allocatable :: cinverse, ainverse, binverse
    contains
        procedure :: get_matrix_dimensions => get_matrix_dimensions
        procedure :: get_dimensions => get_dimensions
        procedure :: setup => setup
        procedure :: setup_inverse_arrays => setup_inverse_arrays
        procedure :: setup_permute_orders => setup_permute_orders
        procedure :: get_order => get_order
        procedure :: get_inverse_order => get_inverse_order
        procedure :: is_permuted => is_permuted
        procedure :: cleanup => cleanup
        procedure :: clear => clear
    end type ttgt_descriptor

    interface ttgt_descriptor
        module procedure constructor_empty
        module procedure constructor
    end interface ttgt_descriptor

contains
    function constructor_empty() result(this)
        type(ttgt_descriptor) :: this

        call this%clear()
    end function constructor_empty

    function constructor(descr, m, n, k, switch_ab) result(this)
        type(tc_descriptor), intent(in) :: descr
        integer, dimension(:), intent(in) :: m, n, k
        logical, intent(in), optional :: switch_ab
        type(ttgt_descriptor) :: this

        this = ttgt_descriptor()

        this%m = m
        this%n = n
        this%k = k

        if ( present(switch_ab)) this%switch_ab = switch_ab

        call this%setup(descr)

    end function constructor

    function get_matrix_dimensions(this, a, b) result(dimensions)
        class(ttgt_descriptor), intent(in) :: this
        class(tensor), intent(in) :: a, b
        integer(int64), dimension(3) :: dimensions

        integer(int64), dimension(:), allocatable :: all_dims
        integer :: idx

        all_dims = this%get_dimensions(a, b)
        dimensions = 1
        do idx = 1, size(this%m)
            dimensions(1) = dimensions(1)*all_dims(this%m(idx))
        end do
        do idx = 1, size(this%n)
            dimensions(2) = dimensions(2)*all_dims(this%n(idx))
        end do
        do idx = 1, size(this%k)
            dimensions(3) = dimensions(3)*all_dims(this%k(idx))
        end do
    end function get_matrix_dimensions

    function get_dimensions(this, a, b) result(dimensions)
        class(ttgt_descriptor), intent(in) :: this
        class(tensor), intent(in) :: a, b
        integer(int64), dimension(:), allocatable :: dimensions

        integer :: idx

        allocate(dimensions(size(this%ainverse)))
        do idx = 1, size(this%ainverse)
            if ( this%ainverse(idx) > 0 ) dimensions(idx) = a%dims(this%ainverse(idx))
        end do

        do idx = 1, size(this%binverse)
            if ( this%binverse(idx) > 0 ) dimensions(idx) = b%dims(this%binverse(idx))
        end do
    end function get_dimensions

    subroutine setup(this, descr)
        class(ttgt_descriptor), intent(inout) :: this
        type(tc_descriptor), intent(in) :: descr

        call this%setup_inverse_arrays(descr)
        call this%setup_permute_orders()
        this%cinverseorder = this%get_inverse_order(this%corder)
    end subroutine setup

    function get_inverse_order(this, order) result(inverse)
        class(ttgt_descriptor), intent(in) :: this
        integer, dimension(:), intent(in) :: order
        integer, dimension(:), allocatable :: inverse

        integer :: idx

        allocate(inverse(size(order)))
        do idx = 1, size(order)
            inverse(order(idx)) = idx
        end do
    end function get_inverse_order

    subroutine setup_inverse_arrays(this, descr)
        class(ttgt_descriptor), intent(inout) :: this
        type(tc_descriptor), intent(in) :: descr

        this%cinverse = descr%get_inverse(descr%c_indices)
        this%ainverse = descr%get_inverse(descr%a_indices)
        this%binverse = descr%get_inverse(descr%b_indices)
    end subroutine setup_inverse_arrays

    subroutine setup_permute_orders(this)
        class(ttgt_descriptor), intent(inout) :: this

        this%corder = this%get_order(this%cinverse, this%m, this%n)
        if ( this%switch_ab) then
            this%aorder = this%get_order(this%ainverse, this%k, this%n)
            this%border = this%get_order(this%binverse, this%m, this%k)
        else
            this%aorder = this%get_order(this%ainverse, this%m, this%k)
            this%border = this%get_order(this%binverse, this%k, this%n)
        end if

        this%permute_c = this%is_permuted(this%corder)
        this%permute_a = this%is_permuted(this%aorder)
        this%permute_b = this%is_permuted(this%border)
    end subroutine setup_permute_orders

    logical function is_permuted(this, order)
        class(ttgt_descriptor), intent(in) :: this
        integer, dimension(:), intent(in) :: order

        integer :: idx

        is_permuted = .not. all(order == [(idx, idx=1, size(order))])
    end function is_permuted

    function get_order(this, inverse, m, n) result(order_array)
        class(ttgt_descriptor), intent(in) :: this
        integer, dimension(:), intent(in) :: inverse, m, n
        integer, dimension(:), allocatable :: order_array

        integer :: idx, offset

        allocate(order_array(size(m) + size(n)))

        offset = 0
        do idx = 1, size(m)
            order_array(offset + idx) = inverse(m(idx))
        end do

        offset = size(m)
        do idx = 1, size(n)
            order_array(offset + idx) = inverse(n(idx))
        end do
    end function get_order

    subroutine cleanup(this)
        class(ttgt_descriptor), intent(inout) :: this

        if ( allocated(this%m) ) deallocate(this%m)
        if ( allocated(this%n) ) deallocate(this%n)
        if ( allocated(this%k) ) deallocate(this%k)
        call this%clear()
    end subroutine cleanup

    subroutine clear(this)
        class(ttgt_descriptor), intent(inout) :: this

        this%switch_ab = .false.
        this%permute_c = .false.
        this%permute_a = .false.
        this%permute_b = .false.
    end subroutine clear
end module ttgt_descriptor_module

max_stars_repo_name: cheshyre/ntcl-algorithms

max_issues_repo_name: cheshyre/ntcl-algorithms

max_forks_repo_name: cheshyre/ntcl-algorithms


text:
module mod_preconditioner
#include <messenger.h>
    use mod_kinds,              only: rk, ik
    use type_preconditioner,    only: preconditioner_t

    ! Import preconditioner types
    use precon_identity,            only: precon_identity_t
    use precon_jacobi,              only: precon_jacobi_t
    use precon_ILU0,                only: precon_ILU0_t
    use precon_RASILU0,             only: precon_RASILU0_t
    implicit none

    ! Instantiate preconditioner types for sourcing
    type(precon_identity_t)             :: IDENTITY
    type(precon_jacobi_t)               :: BLOCKJACOBI
    type(precon_ILU0_t)                 :: ILU0
    type(precon_RASILU0_t)              :: RASILU0

contains

    !>
    !!
    !!  @author Nathan A. Wukie
    !!  @date   11/22/2016
    !!
    !!
    !-----------------------------------------------------------------------------------
    subroutine create_preconditioner(precon_string,instance)
        character(*),                           intent(in)      :: precon_string
        class(preconditioner_t), allocatable,   intent(inout)   :: instance

        character(:),   allocatable :: user_msg

        select case (trim(precon_string))
            case ('identity','Identity','IDENTITY')
                allocate(instance, source=IDENTITY)

            case ('jacobi','Jacobi','blockjacobi','BlockJacobi')
                allocate(instance, source=BLOCKJACOBI)

            case('ilu0','ILU0')
                allocate(instance, source=ILU0)

            case('rasilu0','RASILU0', 'ras-ilu0', 'RAS-ILU0', 'RAS+ILU0', 'ras+ilu0', 'RAS ILU0', 'ras ilu0')
                allocate(instance, source=RASILU0)

            case default
                user_msg = "create_preconditioner: It seems like we can't find the preconditioner &
                            that was specified in the input string. You might check the spelling &
                            and also make sure the preconditioner is registered in the file &
                            mod_preconditioner.f90."
                call chidg_signal_one(FATAL,user_msg,trim(precon_string))

        end select

        !
        ! Make sure the preconditioner was allocated
        !
        user_msg = "create_preconditioner: preconditioner was not allocated. &
                    Check that the desired solver was registered and instantiated &
                    in the mod_preconditioner module."
        if (.not. allocated(instance)) call chidg_signal_one(FATAL,user_msg,trim(precon_string))

    end subroutine create_preconditioner
    !***********************************************************************************

end module mod_preconditioner

max_stars_repo_name: wanglican/ChiDG

max_issues_repo_name: haohb/ChiDG

max_forks_repo_name: haohb/ChiDG


text:
!> \file chk_tools.f90
!! \brief checking module for the tools module
!! @author Innocent Souopgui
!! \details
!<
MODULE chk_tools
  USE general_constant
IMPLICIT NONE

CONTAINS

  !> \brief check argument for the function/subroutine regular_coord_count
  !! \see regular_coord_count
  !<
  SUBROUTINE chk_regular_coord_count(ida_step, ida_max, ida_min, ida_ncoord)
    INTEGER, DIMENSION(:), INTENT(IN) :: ida_step, ida_max
    INTEGER, DIMENSION(:), OPTIONAL, INTENT(IN)    :: ida_min
      INTEGER, DIMENSION(:), OPTIONAL, INTENT(OUT) :: ida_ncoord
    INTEGER, DIMENSION(SIZE(ida_max)) :: ila_min, ila_count
    LOGICAL :: ll_conform_arg !> used to check the conformity of arguments

    IF( PRESENT(ida_min) ) THEN
      ll_conform_arg = ( (SIZE(ida_step)==SIZE(ida_max)).AND.(SIZE(ida_step)==SIZE(ida_min)) )
    ELSE
      ll_conform_arg = ( SIZE(ida_step)==SIZE(ida_max) )
    END IF
    IF( PRESENT(ida_ncoord) ) THEN
      ll_conform_arg = ( ll_conform_arg.AND.(SIZE(ida_step)==SIZE(ida_ncoord)) )
    END IF
    !Check size
    IF(.NOT.ll_conform_arg)THEN
      CALL stop_progam('', 'In chk_tools::regular_coord_count: different sizes for arrays ida_step, ida_max [, ida_min] [and ida_ncoord] ')
    END IF
  END SUBROUTINE chk_regular_coord_count
END MODULE chk_tools

max_stars_repo_name: souopgui/da-toolbox

max_issues_repo_name: souopgui/da-toolbox

max_forks_repo_name: souopgui/da-toolbox


text:
program prog
  use mod2
  implicit none

  if (modval1 + modval2 /= 3) then
    stop 1
  end if

end program prog

max_stars_repo_path: test cases/fortran/7 generated/prog.f90

max_issues_repo_path: test cases/fortran/7 generated/prog.f90

max_forks_repo_path: test cases/fortran/7 generated/prog.f90