        program pwelch_multi
c
chelpini
c_______________________________________________________________
c
c Target:       Welch's periodogram calculation 
c
c Method:       
c               Time series are splitted into a number of 
c               subseries of size N, with an overlapping of 
c               M points. These subseries are filtered with a
c               Hamming window and then their FFT is calculated.
c               Welch's periodogram is the result of averaging  
c               the periodogram of every single subserie.
c
c Comments:     Missing data in the input data are allowed.
c                
c               The program allows the computation of Welch's
c               periodogram for many series at the same time 
c               (nsta)
c
c               Data input and output in .ext format
c               NetCDF files can easily be converted to .ext 
c               using CDO: cdo -f ext copy <filenc> <filext>
c
c
c Program created by Félix García-Pereira
c Last version: October 3, 2021
c Contact info: felgar03@ucm.es
c
c If this program is used, please cite:             
c 
c García-Pereira, F., González-Rouco, J. F., Schmid, T., 
c Melo-Aguilar, C, Vegas-Cañas, C., Steinert, N. J.,       
c Roldán-Gómez, P. J., Cuesta-Valero, F. J., García-García, A., 
c Beltrami, H., and de Vrese, H.: "Thermodynamic and 
c hydrological drivers of the subsurface thermal regime in 
c Central Spain". Earth Surf. Dynam., submitted, 2023.
c________________________________________________________________
chelpfin
c
c
c
	parameter(nser=10,ntp=20000,spval=9E10,pi=3.141592654)
c
        integer       uni,nsta,nvar,ntot,nfr,ssize,ssizecase,ovlap,opsalida
        integer       ovlapcs, windowing
        integer       ia(ntp,4),ib(ntp,4)
        real          avg
        real          f(ntp,20000),h(ntp),g(ntp,1000),gg(ntp,1000)
        real          ampl(ntp),amplis(ntp),ampn(ntp),fase(ntp)
        real          ssamp(ntp,1000),ffass(ntp,1000),fases(ntp)
        real, allocatable :: pwd(:,:),ampliss(:,:),fasess(:,:)
        real, allocatable :: frec(:,:),pe(:,:)
c
        character*100 input,opp,opf,opamp,opfas,fdat,fext
        character*100 opampext,opampdat,opfasdat,opfasext
c
c
      write(*,*)' Minimum sample size required: 50 p'
	    write(*,*)' Input series (.ext required)'
	    read(*,100)input
      write(*,*)' Subsample size (1/sample size): '
      write(*,*)' 10%(1), 25%(2), 40%(3), 50%(4), 75%(5), 100%(6)?'
      read(*,*)ssizecase
      write(*,*)' Allowed overlap (1/subsample size): '
      write(*,*)' 0%(1), 25%(2), 40%(3), 50%(4), 75%(5), 90%(6)?'
      read(*,*)ovlapcs
      write(*,*)' Windowing (0:no/1:Hamming)?'
      read(*,*)windowing
      write(*,*)' Output: amplitudes(0) or power(1)?'
      read(*,*)opsalida
      write(*,*)' Output name for amplitudes/power (NOT append format)'
      read(*,*)opamp
      write(*,*)' Output name for phases (NOT append format)'
      read(*,*)opfas
 100	format(a100)
c
c
c Output files name creation
c
        opp='periods.ext'
        opf='frequencies.ext'
        fdat='.dat'
        fext='.ext'

        opampext = trim(opamp) // trim(fext)
        opampdat = trim(opamp) // trim(fdat)
        opfasext = trim(opfas) // trim(fext)
        opfasdat = trim(opfas) // trim(fdat)

c Output files designation

      open(3, file='pwelch.log', form='formatted')
	    open(10,file=input,form='unformatted')
	    open(20,file=opampext,form='unformatted')
      open(21,file=opfasext,form='unformatted')
      open(22,file=opp,form='unformatted',status='unknown')
      open(23,file=opf,form='unformatted',status='unknown')
      open(30,file=opampdat,form='formatted') 
      open(31,file=opfasdat,form='formatted')

c
c.....  Input data reading 
c
        i=0
 10     i=i+1
c
        uni=10
        read(uni,end=999) (ia(i,j),j=1,4)
        nvar=ia(i,3)
        nsta=ia(i,4)
        read(uni)(f(i,j),j=1,nsta)
c        write(3,*)ia(i,1),nsta,f(i,1)
        goto 10
 999    ntot=i-1
 20     continue
c        write(3,*)'ntot,nsta ',ntot,nsta
c

c Output data size allocation

        allocate(ampliss(ntp,nsta))
        allocate(fasess(ntp,nsta))
        allocate(pwd(ntp,nsta))
        allocate(frec(ntp,nsta))
        allocate(pe(ntp,nsta))
c
c
        
        do st=1,nsta         !Loop ending in line 257

        write(*,*) '                Station/gridpoint number'
        write(*,*)st

        do 40 i=1,ntot  
        h(i)=f(i,st)                   
 40     continue

c.....  Serie splitting according to ssize and ovlap parameters
        

        SELECT CASE (ssizecase)
          CASE (1)
            ssize=int(0.1*ntot)    
          CASE (2)
            ssize=int(0.25*ntot)
          CASE (3)
            ssize=int(0.4*ntot)
          CASE (4)
            ssize=int(0.5*ntot)
          CASE (5)
            ssize=int(0.75*ntot)
          CASE (6)
            ssize=int(1.00*ntot)
          CASE DEFAULT
            WRITE(*,*)  "Unknown subsample size case"
        END SELECT


        SELECT CASE (ovlapcs)
          CASE (1)
            ovlap=int(0.0*ssize) 
          CASE (2)
            ovlap=int(0.25*ssize)
          CASE (3)
            ovlap=int(0.4*ssize)    
          CASE (4)
            ovlap=int(0.5*ssize)
          CASE (5)
            ovlap=int(0.75*ssize)
          CASE (6)
            ovlap=int(0.9*ssize)
          CASE DEFAULT
            WRITE(*,*)  "Unknown allowed overlapping"
        END SELECT

        k=0
        j=1

        do while (k.le.ntot-ssize+1)
          do i=k+1,k+ssize
            g(i-(ssize-ovlap)*(j-1),j)=h(i)
          enddo
          k=k+ssize-ovlap
          j=j+1
        enddo

c        write(*,*)j,k,ssize,ovlap,ntot
        write(*,*) 'Number of subseries', j-1


c..... Subseries with NaN data >25 % are discarded

        k=0

        do l=1,j-1
          nnnans=0

          do i=1,ssize
            if (g(i,l).ne.spval) then
              nnnans=nnnans+1
            else
              nnnans=nnnans
            endif
          enddo

          if (nnnans.ge.int(0.75*ssize)) then
            k=k+1
            do i=1,ssize
              g(i,k)=g(i,l)
            enddo
          else
            k=k
          endif
        enddo

        je=k+1

        if (je == 1) then
          do i=1,ssize/2+1
            pwd(i,st)=spval
            ampliss(i,st)=spval
          end do
          cycle        
        end if


c.....  Computing anomalies with respect the mean for every subserie

        do l=1,j-1

          call med(g(:,l),ssize,avg)   

          do i=1,ssize
            if (abs(gg(i,l)-spval).gt.10000) then
              g(i,l)=g(i,l)-avg
            else
              g(i,l)=spval
            endif
          enddo

        enddo


c.....  Hamming filtering of the subserie
  
        do l=1,j-1
          nnnans=0

          do i=1,ssize
            if (abs(g(i,l)-spval).gt.10000) then
              nnnans=nnnans+1
            else
              nnnans=nnnans
            endif
          enddo

          do i=1,ssize
            if (windowing.eq.0) then
              if (abs(g(i,l)-spval).gt.10000) then
                gg(i,l)=g(i,l)
              else
                gg(i,l)=-9999.9004
              endif
            elseif (windowing.eq.1) then
              if (abs(g(i,l)-spval).gt.10000) then
                gg(i,l)=g(i,l)*(0.538-0.462*cos(2*pi*i/(nnnans-1)))
              else
                gg(i,l)=-9999.9004
              endif
            endif
          enddo

        enddo

c.....  Fractioning and filtering of the original input serie completed in line 257

c.....  Welch's periodogram computing      
c.....  FFT of every subserie

        do l=1,j-1

          call armo(ssize,nfr,gg(:,l),ampl,amplm,ampn,ampnm,nonans,fase) 

          ssamp(:,l)=ampl
          ffass(:,l)=fase

        enddo

c.....  Subserie averaging (Welch)
c.....  NOTE: ssize/2 is taken since for a sample size of n the FFT size is n/2 (F_Nyquist=1/2)

        do i=1,int(ssize/2)
          sum=0
          do l=1,j-1
            sum=sum+ssamp(i,l)
          enddo
          amplis(i)=sum/(j-1)
        enddo

        do i=1,int(ssize/2)
          sum=0
          do l=1,j-1
            sum=sum+ffass(i,l)
          enddo
          fases(i)=sum/(j-1)
        enddo
c
        do 60 i=1,nfr
        frec(i,st)=i/float(ssize)
        pe(i,st)=float(ssize)/i
        pwd(i,st)=2*amplis(i)/ssize
        ampliss(i,st)=sqrt(2*amplis(i)/ssize)
        fasess(i,st)=fases(i)

 60     continue

        end do          !LOOP end for series in nsta


c 

        do i=1,nfr
          ib(i,1)=i
          ib(i,2)=nfr+1
          ib(i,3)=nvar
          ib(i,4)=nsta
        enddo
       

        if (opsalida.eq.0) then
        do i=1,nfr
          write(20) (ib(i,j),j=1,4)
          write(20)ampliss(i,:)
        enddo
        elseif (opsalida.eq.1) then
        do i=1,nfr
          write(20) (ib(i,j),j=1,4)
          write(20)pwd(i,:)
        enddo
        endif

        do i=1,nfr
          write(21) (ib(i,j),j=1,4)
          write(21)fasess(i,:)
        enddo        

        do i=1,nfr
          write(22) (ib(i,j),j=1,4)
          write(22)pe(i,:)
        enddo

        do i=1,nfr
          write(23) (ib(i,j),j=1,4)
          write(23)frec(i,:)
        enddo


        if (opsalida.eq.0) then
        do j=1,nfr
        write(30,*)pe(j,1),frec(j,1),ampliss(j,:)
        end do
        elseif (opsalida.eq.1) then
        do j=1,nfr
        write(30,*)pe(j,1),frec(j,1),pwd(j,:)
        enddo
        endif

        do j=1,nfr
        write(31,*)pe(j,1),frec(j,1),fasess(j,:)
        end do

c--
c  30     continue
c
c        stop
        end
c
c===========================================
c
         SUBROUTINE MED(vec,l,mean)
c
c        vec: input vector
c        l: length vec
c        mean: =
c===========================================
c
         integer l,c,nnans
         real mean,sum
         real vec(200000)
        parameter(nser=10,ntp=200000,spval=9E10)
c.....
         mean=0.0
         sum=0.
         nnans=0
c.....
         do 300 c=1,l
         if (abs(vec(c)-spval).gt.10000) then
         sum=sum+vec(c)
         nnans=nnans+1
         else
         sum=sum
         nnans=nnans
         endif
c
 300     continue
c
         mean=sum/nnans
         sum=0.
c.....
        return
        end
c

C========================================
c Harmonic analysis (FFT) subroutine
       SUBROUTINE ARMO(nmaf,m,vec,ampl,amplm,ampn,ampnm,nonans,fas)
c_________
c  xvec.-  data vector
c  nmaf.-  xvec size
c  m.-     number of harmonic frequencies
c  ampl .-  amplitude vector
c  amplm.- if n is an even number, amplitude of freq 0.5
c  ampn.-   normalized amp
c  ampnm.- normalized ampn2
c  fas.- phase
C========================================
C
c
      integer c,cc,nn,nmaf,m,nnn,sn
c
      real am,amplm,ampnm
      real A(150000),B(150000),AMPl(150000),ampn(150000),fas(150000)
      real vec(300000),frec(150000),pe(150000)
      parameter(nser=10,ntp=200000,spval=-9999.9004)
c
c..... Number of spectral estimates
c
c
      sn=nmaf/2.
c
      nn=int(sn)        ! Half of the data size
c
      if (sn.eq.nn) then           !n even
c
              m=(nmaf-1)/2
              nnn=1
c
      else                         !n odd
c
              m=(nmaf-1)/2
              nnn=2
      endif
c      write(3,*)nmaf,m
c
c..... Harmonic analysis
c
       CALL SUMAI(vec,nmaf,nnn,m,a,am)    !calculates ai
       CALL SUMBI(vec,nmaf,m,b,nonans)        !calculates bi
       CALL AMPLI(a,am,b,m,nonans,ampl,amplm)  !calculates amplitudes
       CALL NORM(amplm,ampl,nonans,m,ampn,ampnm)  !normalizes ampli
       CALL FASE(a,b,m,nonans,fas)
        do 800 cc=1,m
c         write(3,*) a(cc),am,b(cc),a(cc),ampn(cc)
 800    continue
c
        return
        end
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c

c
c----------------------------------------------------------
      SUBROUTINE SUMAI (vec,nmaf,nnn,m,a,am)
c
C========================================
c Subroutine for AI harmonic coefficients calculation
c_________
c  xvec.-  data vector
c  nmaf.-  xvec size
c  m   .-  number of harmonic frequencies
c  a   .-  ai coefficients
c  anm .-  sum vector
C========================================
c
c
      integer nmaf,m,c,cc,nonans
c
      real pi,suma0,sumam,suma,arg,am
      real vec(300000),xxvec(300000),A(150000)
      parameter(nser=10,ntp=200000,spval=-9999.9004)
c
c..... Initial values for the loop
c
         suma0=0.
         sumam=0.
         pi=3.141592654
c
c..... The loop starts
c

      do 400 c=1,m
c
         suma=0.
         nonans=0
c
      do 410 cc=1,nmaf
c

         arg=(2*pi*c*cc/(float(nmaf))-(2*pi*int(c*cc/float(nmaf))))
         if (abs(vec(cc)-spval).gt.0.01) then
         xxvec(cc)=vec(cc)*COS(arg)
         suma=suma+xxvec(cc)
         nonans=nonans+1
         else
         xxvec(cc)=spval
         suma=suma
         endif

c
 410    continue
c
         a(c)=suma*2/float(nonans)
c
 400    continue
c

       do 430 cc=1,nmaf
c
         if (abs(vec(cc)-spval).gt.0.01) then
         sumam=sumam+(((-1)**(cc))*vec(cc))
         else
         sumam=sumam
         endif
c
 430    continue
c

         if(nnn.eq.1) then
         am=sumam
         else
         am=0
         endif
c
      return
      end
c
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C              END SUBROUTINE SUMAI
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      SUBROUTINE SUMBI(vec,nmaf,m,b,nonans)
c
c
C========================================
c Subroutine for BI harmonic coefficients calculation
c_________
c  xvec.-  data vector
c  nmaf.-  xvec size
c  m   .-  number of harmonic frequencies
c  b   .-  Bi coefficients
C========================================
c
      integer nmaf,m,c,cc,nonans
c
      real pi,sumb,arg
      real vec(300000),xxvec(300000),b(150000)
      parameter(nser=10,ntp=200000,spval=-9999.9004)
c
c..... Initial values for the loop
c

      pi=3.141592654
c
c..... The loop starts
c
      DO 500 c=1,m

      SUMB=0
      nonans=0
c
        do 510 cc=1,nmaf
c
         ARG=(2*PI*c*cc/(FLOAT(Nmaf))-(2*PI*INT(c*cc/FLOAT(Nmaf))))
         if (abs(vec(cc)-spval).gt.0.01) then
         xxvec(cc)=vec(cc)*SIN(ARG)
         SUMB=xxvec(cc)+SUMB
         nonans=nonans+1
         else
         xxvec(cc)=spval
         SUMB=SUMB
         endif
c
 510     continue
c
         B(c)=SUMB*2/FLOAT(nonans)
c
 500   CONTINUE
c
      return
      end
c
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C              END SUBROUTINE SUMBI
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C

      SUBROUTINE AMPLI(a,am,b,m,nmaf,ampl,amplm)
c       
C========================================
c Harmonic amplitudes subroutine
c_________
c  A.-  ai coefficients
c  b.-  bi coefficients
c  amp.- amplitudes
c  anm.- sumai cummulative addition
c  ampn2.- anm squared
c  m   .-  number of harmonic frequencies
C========================================
c
c
      integer m,c,cc
c
      real am,amplm
      real a(150000),b(150000),ampl(150000)
      parameter(nser=10,ntp=200000,spval=-9999.9004)
c
c..... Amplitudes --> explained variance of every single harmonic
c
      do 600 cc=1,m
c
        ampl(cc)=(((a(cc)**2.0)+(b(cc)**2.0))*nmaf)/2.0
c
 600  continue
c
        amplm=(am**2)/float(nmaf)
c
      return
      end
c
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C              END SUBROUTINE AMPLI
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C
c
      SUBROUTINE NORM(amplm,ampl,nmaf,m,ampn,ampnm)
c
C========================================
c Normalized harmonic amplitudes subroutine
c_________
c  amp.- amplitudes
c  anm.- sumai cummulative addition
c  ampn2.- anm squared
c  n   .-  data size
c  m   .-  number of harmonic frequencies
c  g   .-  normalized amplitudes
c  gn2 .-  normalized anm
C========================================
c
c
      integer M,N,c,cc
      real AMPl(300000),ampn(300000)
      real amplm,ampnm
      real zma,szma,num
      parameter(nser=10,ntp=200000,spval=-9999.9004)
c
c..... Initial values for the loop
c
      SZMA=0
      ZMA=0
c
c.....
c
       DO 650 c=1,m
c
          ZMA=zma+ampl(c)
c
c       write(3,*)c,ampl(c),zma
c
 650  CONTINUE
C
         szma=zma+amplm
c       write(3,*)'szma',szma
c
       DO 660 c=1,m
c
        num=ampl(c)
        ampn(c)=num/szma

c       write(3,'(i4,2f15.4)')c,ampl(c),ampn(c)
c
 660  CONTINUE
c
       ampnm=amplm/szma
c       write(3,*)'ampnm',ampnm
C
      return
      end
c
C

      SUBROUTINE FASE(a,b,m,nmaf,fas)
c       
C========================================
c Harmonic phase calculation subroutine
c_________
c  A.-  ai coefficients
c  b.-  bi coefficients
c  fas.- phase values
c  m   .-  number of harmonic frequencies
C========================================
c
c
      integer m,c,cc
c
      real a(150000),b(150000),fas(150000)
      parameter(nser=10,ntp=200000,spval=-9999.9004)
c
c
      do 600 cc=1,m
c
        fas(cc)=atan2(b(cc),a(cc))

        if (fas(cc).lt.0) then
          fas(cc)=fas(cc)+6.2832
        endif
c
 600  continue
c
      return
      end
c
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C              END SUBROUTINE PHASE
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C