!C==============================================================
!C
!c    calculation of longitudinal harmonics    
!c
!c    in NCAR/NCEP Geopotential height data
!C
!c==============================================================
!c
!C------------------------- the number of zonal harmonics
      PARAMETER (NG=4)
      Parameter (NDIM=9)
!c------------  NF has to be < or = NDIM
      PARAMETER (NF=1+2*NG)
!c *****************************
!c------------------------------ 3-hourly data for 31 days      
      PARAMETER (IDAY=8)  
      PARAMETER (ILON=64)
      PARAMETER (ILAT=36)
      parameter (kgit=10)

      REAL  dolg(ILON),Geopot(ILON),gh_A0(ILAT,kgit)
      REAL  gh_A1(ILAT,kgit),  gh_A2(ILAT,kgit)
      REAL  gh_A3(ILAT,kgit)
      REAL  gh_F1(ILAT,kgit),  gh_F2(ILAT,kgit)
      REAL  gh_F3(ILAT,kgit)
      REAL  gh_A4(ILAT,kgit),  gh_F4(ILAT,kgit)
      REAL  aus1(ILON,ILAT,kgit)

      REAl  AM1(NDIM,NF),A(NDIM,NF),U(NDIM,NF),V(NDIM,NF),WORK(2*NF)
      REAl  FUN(NF,ilon),SIGMA(NF),AlFA(NF),SQ(NF),AG(NG),FG(NG)


      PI=2.*ASIN(1.)

!....OUTPUTFILES
          OPEN (100,FILE='lh_a0_oct.dx', &  
     &      form='unformatted',access='direct', status='unknown',recl=4*ILAT*kgit) 
          OPEN (1,  FILE='lh_m1_oct.dx', &
     &      form='unformatted',access='direct', status='unknown',recl=4*ILAT*kgit*2)
          OPEN (2,  FILE='lh_m2_oct.dx', &
     &      form='unformatted',access='direct', status='unknown',recl=4*ILAT*kgit*2) 
          OPEN (3,  FILE='lh_m3_oct.dx', &
     &      form='unformatted',access='direct', status='unknown',recl=4*ILAT*kgit*2)
          OPEN (4,  FILE='lh_m4_oct.dx', &
     &      form='unformatted',access='direct', status='unknown',recl=4*ILAT*kgit*2)  
          open(99,file='err') 

!....INPUTFILE      
          OPEN (101, FILE='oct_DT.dat', &
     &      form='unformatted',access='direct', status='old',recl=4*ILON*ILAT*kgit)    

      KLEF=0

      DO  l=1,IDAY       !alle Zeiten
!     Einlesen der Zeitreihe
      read(101,rec=l) aus1

      DO  io=1,kgit   !alle H÷hen
      do 2001 j=1,ILAT
      do 2002 i=1,ILON   
      Geopot(i)=aus1(i,j,io)
      if(Geopot(i).eq.0.) Geopot(i)= 1.e-4
      dolg(i)=(float(i)-0.9999)*5.625    
2002  continue
! ====================================================
      CALL HARMON(NDIM,NF,AM1,A,U,V,FUN,SIGMA,ALFA,SQ,ILON,DOLG,Geopot,&
     &               NG,ACONST,AG,FG,WORK,KLEF)
      KLEF=1
      gh_A0(j,io)=ACONST
      gh_A1(j,io)=AG(1)
      gh_F1(j,io)=FG(1)
      gh_A2(j,io)=AG(2)
      gh_F2(j,io)=FG(2)
      gh_A3(j,io)=AG(3)
      gh_F3(j,io)=FG(3)
      gh_A4(j,io)=AG(4)
      gh_F4(j,io)=FG(4)
2001  continue
      ENDDO   !zeiten
      WRITE(100,rec=l) gh_A0
      WRITE  (1,rec=l) gh_A1,gh_F1
      WRITE  (2,rec=l) gh_A2,gh_F2
      WRITE  (3,rec=l) gh_A3,gh_F3
      WRITE  (4,rec=l) gh_A4,gh_F4
      Enddo   !h÷hen
      CLOSE(100)
      CLOSE(1)
      CLOSE(2)
      CLOSE(3)
      CLOSE(4)
      CLOSE(99)
       STOP
       END     
!============================= SUBROUTINE PART =====================     
      SUBROUTINE HARMON(NDIM,NF,AM1,A,U,V,FUN,SIGMA,ALFA,SQ, &
     &                  NDOLG,DOLG,YDOLG,NG,ACONST,AG,FG,WORK,KLEF)
      DIMENSION AM1(NDIM,NF),A(NDIM,NF),U(NDIM,NF),V(NDIM,NF), &
     &          FUN(NF,NDOLG),SIGMA(NF),ALFA(NF),SQ(NF),DOLG(NDOLG), &
     &          YDOLG(NDOLG),AG(NG),FG(NG),WORK(2*NF)
      PI=2.*ASIN(1.) 
      IF (KLEF.GT.1) GO TO 35
      RELERR=0.0001
      DO 10 K=1,NDOLG
      FUN(1,K)=1.
      DO 20 J=1,NG
      S=float(J)*DOLG(K)*PI/180.
      FUN(J+1,K)=COS(S)
   20 FUN(J+1+NG,K)=SIN(S)
   10 CONTINUE
!C------------ normalization of basic functions
      DO 200 J=1,NF
      S=0.
      DO 400 K=1,NDOLG
  400 S=S+FUN(J,K)*FUN(J,K)
      SQ(J)=SQRT(S)
      DO 500 K=1,NDOLG
  500 FUN(J,K)=FUN(J,K)/SQ(J)
  200 CONTINUE
!C------ least square approximation
      DO 30 I=1,NF
      DO 30 J=1,NF
      S=0.
      DO 40 K=1,NDOLG
   40 S=S+FUN(I,K)*FUN(J,K)
      AM1(I,J)=S
   30 CONTINUE

      tol=1.e-4 
      call LGINF(AM1,NDIM,nf,nf,tol,A,nf,SIGMA,WORK,ierr)
  
   35 CONTINUE
      DO 50 I=1,NF
      S=0.
      DO 60 K=1,NDOLG
   60 S=S+YDOLG(K)*FUN(I,K)
      WORK(I)=S
   50 CONTINUE  
      CALL PRMV(NDIM,NF,NF,A,WORK,ALFA)
      DO 70 I=1,NF
   70 ALFA(I)=ALFA(I)/SQ(I)
      ACONST=ALFA(1)
!C----- amplitudes ad phases calculation 
      DO 80 J=1,NG
      AR=ALFA(1+J)
      AI=ALFA(1+J+NG)
      AG(J)=SQRT(AR*AR+AI*AI)
!c------ the phase - position of ridge
      FG(J)=360.+ATAN2(AI,AR)*180./PI
      IF(FG(J).GT.360.)FG(J)=FG(J)-360.
      FG(J)=FG(J)/FLOAT(J)
   80 CONTINUE
      RETURN
      END