program timeanal
      PARAMETER (NF=7,ng=(NF-3)/2,nout=3)
      parameter (ntime=248,ilat=36,kgit=10)
      implicit real*4(a-h,o-z)
      real amph(ntime,ilat,kgit),amp(ilat,kgit)
c      real phas(ntime,ilat,kgit)
      real ghr(ntime),ghrinp(ntime),ghi(ntime),ghiinp(ntime),
     $         tHour(ntime) !, tDays(ntime)
      real ghrDiff(ntime), ghiDiff(ntime)
      real FUN(NF,ntime),SQ(NF),A(NF,NF),AM1(NF,NF),X(NF),Y(NF)
      real Xstw(ilat,kgit), Ystw(ilat,kgit)!,stwx(),stwy
      real SIGMA(NF), WORK(2*NF), errx(nf),erry(nf),ierr
      real period(nf), amplr(nf), ampli(nf)
c      real Aw(ng,ilat,kgit),Ae(ng,ilat,kgit)
c      real Pw(ng,ilat,kgit),Pe(ng,ilat,kgit)
      real  apwe(ng,ilat,kgit,2)
c      character*13 wFile, OutFile
c      character*4 filename
      pi = 4. * atan(1.0)
      dolg=0.

c--------------------------m=0
      Period1 = 24.       ! 24h
      Period2 = 12.       ! 12h

      om1 = 2. * pi  / Period1 
      om2 = 2. * pi  / Period2 

      Period(1)=period1
      Period(2)=period2 

      print*,'period1=',period1,' period2=',period2

c---- read  data
      open (11, FILE='LH_a0_La.dx' ,
     $      form='unformatted',
     *      access='direct', status='old',recl=4*ILAT*kgit)
c      open (12, FILE='mw_a0.dx' ,
c     $      form='unformatted',
c     *      access='direct', status='old',recl=4*ntime*ILAT*kgit)
c      open (13, FILE='tp_a0.dx' ,
c     $      form='unformatted',
c     *      access='direct', status='old',recl=4*ntime*ILAT*kgit)

C_________Output
      open (15, FILE='LH_La_m0.dx' ,
     $      form='unformatted',
     *      access='direct', status='unknown',recl=4*(ng)*ILAT*kgit*2)
c      open (16, FILE='freq_ave_v.dx' ,
c     $      form='unformatted',
c     *      access='direct', status='unknown',recl=4*(ng)*ILAT*kgit*4)
c      open (17, FILE='freq_ave_T.dx' ,
c     $      form='unformatted',
c     *      access='direct', status='unknown',recl=4*(ng)*ILAT*kgit*4)

      do il=1,1   !files
      do j=1,ntime
       read(10+il, rec=j) amp
       do io=1,kgit
       do ii=1,ilat   
       amph(j,ii,io)=amp(ii,io)
       end do
       end do
      end do 
       shour=0.
      print*,'NF=',NF,', ng=',ng, ', m=',m, ', nout=',nout
      do j=1,ntime
       tHour(j) = float(j*nout)   
      shour=shour+tHour(j)
       enddo   
      hourav=shour/float(ntime)

      do io=1,kgit
      do ii=1,ilat

c...calculating Cos and Sinus components

      do j=1,ntime
      ghrinp(j)=amph(j,ii,io)
      ghiinp(j)=amph(j,ii,io)
      end do

      j0=1
      j1=1

        DO J=1,ntime
      ghr(j)=ghrinp(j+j1-1) 
      ghi(j)=ghiinp(j+j1-1)
      FUN(1,J)=1.
      FUN(2,J)=(tHour(j+j1-1)-hourav)/hourav
      FUN(3,J)=FUN(2,J)*FUN(2,J)
      FUN(4,J)=sin(om1 * (tHour(j+j1-1)-tHour(j1))) 
      FUN(5,J)=COS(om1 * (tHour(j+j1-1)-tHour(j1)))  
      FUN(6,J)=sin(om2 * (tHour(j+j1-1)-tHour(j1)))  
      FUN(7,J)=COS(om2 * (tHour(j+j1-1)-tHour(j1)))  
c      FUN(8,J)=sin(om7 * (tHour(j+j1-1)-tHour(j1)))   !16dh
c      FUN(9,J)=COS(om7 * (tHour(j+j1-1)-tHour(j1)))   !16dh
       ENDDO

      DO I=1,NF
          S=0.
          DO J=1,ntime
          S=S+FUN(I,J)*FUN(I,J)
          ENDDO
      SQ(I)=SQRT(S)
            DO J=1,ntime
            FUN(I,J)=FUN(I,J)/SQ(I)
            ENDDO
      ENDDO
      DO I=1,NF
          DO K=1,NF
          S=0.
             DO J=1,ntime
             S=S+FUN(I,J)*FUN(K,J)
             ENDDO
          AM1(I,K)=S
          ENDDO
      ENDDO

      tol=1.e-4 
      call LGINF(AM1,nf,nf,nf,tol,A,nf,SIGMA,WORK,ierr)
           k=1

           DO I=1,NF
              S=0.
              DO J=1,ntime
              S=S+ghr(J)*FUN(I,J)
              ENDDO
           WORK(I)=S
           ENDDO
      call PRMV(nf,nf,nf,A,WORK,X)
         DO I=1,NF
         X(I)=X(I)/SQ(I)
         if(I.eq.1) Xstw(ii,io)=X(I)
         ENDDO
           DO I=1,NF
              S=0.
              DO J=1,ntime
              S=S+ghi(J)*FUN(I,J)
              ENDDO
           WORK(I)=S
           ENDDO
      call PRMV(nf,nf,nf,A,WORK,Y)
         DO I=1,NF
         Y(I)=Y(I)/SQ(I)
         if(I.eq.1) Ystw(ii,io)=Y(I)
         ENDDO
C--------------- CALCULATIONS OF ERRORS IN COEFFICIENTS
           DO I=1,NF
                     SSu=0.          
                     SSv=0.          
                     DO J=1,ntime
                                 Su=0.
                                 Sv=0.
                                 DO L=1,NF
                                 Su=Su+A(I,L)*FUN(L,J)*0.1
                                 Sv=Sv+A(I,L)*FUN(L,J)*0.1
                                 END DO
                     SSu=SSu+Su*Su
                     SSv=SSv+Sv*Sv
                     END DO 
           ERRX(I)= SQRT(SSu)/SQ(I)                           
           ERRY(I)= SQRT(SSv)/SQ(I)                           
          END  DO 

c        ng=(nf-2)/2
        do i=1,ng 
        is=4+2*(i-1)
        ic=5+2*(i-1)
        ys=y(is)
        yc=y(ic)
        xs=x(is)
        xc=x(ic)
        amplI(i)=sqrt(ys*ys+yc*yc)
        amplR(i)=sqrt(xs*xs+xc*xc)
        ERRrS=ERRX(IS)
        ERRrC=ERRX(IC)
        ERRiS=ERRY(IS)
        ERRiC=ERRY(IC)
        ERRi=SQRT(ERRiS*ERRiS*YS*YS+ERRiC*ERRiC*YC*YC)/AMPLI(i)
        ERRr=SQRT(ERRrS*ERRrS*XS*XS+ERRrC*ERRrC*XC*XC)/AMPLR(i)
        Frequ= 1./period(i) 
        phaseI=180.*ATAN2(ys,yc)/pi
        if(phaseI.lt.0.)phaseI=phaseI+360.
        phaseR=180.*ATAN2(xs,xc)/pi        
        if(phaseR.lt.0.)phaseR=phaseR+360.
c----------------------------------------------------------------
c        phaseR=phaseR+(thour(j1)-thour(j0))*360./period(i)
c        phaseI=phaseI+(thour(j1)-thour(j0))*360./period(i)

  601 continue
       if(phaser.gt.360.)then
       phaser=phaser-360.
       Go to 601
       End if    
  602 continue
       if(phasei.gt.360.)then
       phasei=phasei-360.
       Go to 602
       End if    
        errphR=180./(xs*xs+xc*xc)*SQRT(xc*xc*errx(is)*errx(is)+
     $                                 xs*xs*errx(ic)*errx(ic))/pi        
        errphI=180./(ys*ys+yc*yc)*SQRT(yc*yc*erry(is)*erry(is)+
     $                                 ys*ys*erry(ic)*erry(ic))/pi        


        ampl=sqrt(amplr(i)*amplr(i)+ampli(i)*ampli(i))
        phas=phaseR
       if( phas.gt.360) phas = phas -360.

       apwe(i,ii,io,1)=ampl
       apwe(i,ii,io,2)=phas

         end do        !von ng

       jout=ntime
       DO J=1,Jout
       tDout=thour(j0+j-1)
            ghrAPPROX=0.
            ghiAPPROX=0.
            DO I=1,nf
            ghrApprox = ghrAPPROX+X(I)*SQ(I)*FUN(I,J)
            ghiApprox = ghiAPPROX+Y(I)*SQ(I)*FUN(I,J)
            ENDDO
            ghrDiff(j) = ghr(j) - ghrApprox
            ghiDiff(j) = ghi(j) - ghiApprox

       enddo

      enddo !breiten
      enddo !Höhen


      write(14+il,rec=1) apwe
      enddo !files
      end