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=20)
c------------  NF has to be < or = NDIM
      PARAMETER (NF=1+2*NG)
c *****************************
c------------------------------ 4-hourly output for 30 days      
      PARAMETER (IDAY=1201)  
      PARAMETER (ILON=64)
      PARAMETER (ILAT=36)
      parameter (kgit=56)
      parameter (nvar=3)
c
      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  aus(ILON,ILAT,kgit,nvar)
      REAL  AM1(NDIM,NF),A(NDIM,NF),U(NDIM,NF),V(NDIM,NF),WORK(NF)
      REAL  FUN(NF,ILON),SIGMA(NF),ALFA(NF),SQ(NF),AG(NG),FG(NG)

      PI=2.*ASIN(1.)

C....OUTPUTFILES
          OPEN (100,FILE='tp_a0_grads.dx',   
     $      form='unformatted',
     *      access='direct', status='unknown',recl=4*ILAT*kgit) 
          OPEN (1,  FILE='tp_m1_grads.dx' ,
     $      form='unformatted',
     *      access='direct', status='unknown',recl=4*ILAT*kgit*2)
          OPEN (2,  FILE='tp_m2_grads.dx',
     $      form='unformatted',
     *      access='direct', status='unknown',recl=4*ILAT*kgit*2) 
          OPEN (3,  FILE='tp_m3_grads.dx',
     $      form='unformatted',
     *      access='direct', status='unknown',recl=4*ILAT*kgit*2)
          OPEN (4,  FILE='tp_m4_grads.dx'  ,
     $      form='unformatted',
     *      access='direct', status='unknown',recl=4*ILAT*kgit*2)  
          open(99,file='err') 

c....INPUTFILE      
          OPEN (101, FILE='uvt400_JF_NF_wQBO_3.dx',
     $      form='unformatted',
     *      access='direct', status='old',recl=4*ILON*ILAT*kgit*nvar)    
c
      KLEF=0
         mm=3
         print*,'mm=',mm
      DO  l=1,IDAY       !alle Zeiten
c     Einlesen der Zeitreihe
      read(101,rec=l) aus
      DO  io=1,kgit   !alle H÷hen
      do 2001 j=1,ILAT
      do 2002 i=1,ILON   
      Geopot(i)=aus(i,j,io,mm)
      dolg(i)=(float(i)-0.9999)*5.625    
2002  continue
c ====================================================
      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     
C ============================= 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(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
      CALL INVSVD(NDIM,NF,NF,AM1,A,U,V,SIGMA,
     *                           WORK,RELERR)  
   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
c*********************************************************************
C
C     " I N V S V D "
C
C
      SUBROUTINE INVSVD(NDIM,M,N,A,B,U,V,
     *                  SIGMA,WORK,RELERR)
C
      REAL A(NDIM,N),U(NDIM,N),V(NDIM,N),
     *     SIGMA(N),WORK(N),B(NDIM,M)
C

      CALL SVD(NDIM,M,N,A,SIGMA,.TRUE.,U,.TRUE.,V,
     *                                   IERR,WORK)
C
      IF(IERR.NE.0) WRITE(99,1) IERR
   1  FORMAT(/'   /INVSVD/    TROUBLE.IERR = ',I4/)
C
      SIGMA1=0.
      DO 2 I=1,N
      IF(SIGMA(I).GT.SIGMA1) SIGMA1=SIGMA(I)
   2  CONTINUE
C
      TAU=RELERR*SIGMA1
C
      DO 4 I=1,N
      DO 4 K=1,M
C
      S=0.
      DO 3 J=1,N
      IF(SIGMA(J).LE.TAU) GO TO 3
      S=S+V(I,J)*U(K,J)/SIGMA(J)
   3  CONTINUE
C
      B(I,K)=S
C
   4  CONTINUE
C
      RETURN
      END
c *********************************************************************
      SUBROUTINE PRMV(NDIM,NI,NJ,A,B,C)
      REAL A(NDIM,NJ),B(NJ),C(NI),S
      DO 1 I=1,NI
      S=0.
      DO 2 J=1,NJ
    2 S=S+A(I,J)*B(J)
      C(I)=S
    1 CONTINUE
      RETURN
      END
c *********************************************************************
C
C     " S V D "
C
C
      SUBROUTINE SVD(NM,M,N,A,W,MATU,U,MATV,V,IERR,RV1)
C
      REAL A(NM,N),W(N),U(NM,N),V(NM,N),RV1(N)
C
      LOGICAL MATU,MATV
C
                       IERR=0
      DO 100 I=1,M
      DO 100 J=1,N
  100 U(I,J)=A(I,J)
      G=0.
      SCALE=0.
      ANORM=0.
      DO 300 I=1,N
      L=I+1
      RV1(I)=SCALE*G
      G=0.
      S=0.
      SCALE=0.
      IF(I.GT.M) GO TO 210
      DO 120 K=I,M
      SCALE=SCALE+ABS(U(K,I))
  120 CONTINUE
      IF(SCALE.EQ.0.) GO TO 210
      DO 130 K=I,M
      U(K,I)=U(K,I)/SCALE
      S=S+U(K,I)**2
  130 CONTINUE
      F=U(I,I)
      G=-SIGN(SQRT(S),F)
      H=F*G-S
      U(I,I)=F-G
      IF(I.EQ.N) GO TO 190
      DO 150 J=L,N
      S=0.
      DO 140 K=I,M
      S=S+U(K,I)*U(K,J)
  140 CONTINUE
      F=S/H
      DO 150 K=I,M
      U(K,J)=U(K,J)+F*U(K,I)
  150 CONTINUE
  190 DO 200 K=I,M
  200 U(K,I)=SCALE*U(K,I)
  210 W(I)=SCALE*G
      G=0.
      S=0.
      IF(I.GT.M.OR.I.EQ.N) GO TO 290
      DO 220 K=L,N
  220 SCALE=SCALE+ABS(U(I,K))
C
      IF(SCALE.EQ.0.) GO TO 290
      DO 230 K=L,N
      U(I,K)=U(I,K)/SCALE
      S=S+U(I,K)**2
  230 CONTINUE
      F=U(I,L)
      G=-SIGN(SQRT(S),F)
      H=F*G-S
      U(I,L)=F-G
      DO 240 K=L,N
  240 RV1(K)=U(I,K)/H
      IF(I.EQ.M) GO TO 270
      DO 260 J=L,M
      S=0.
      DO 250 K=L,N
  250 S=S+U(J,K)*U(I,K)
      DO 260 K=L,N
  260 U(J,K)=U(J,K)+S*RV1(K)
  270 DO 280 K=L,N
  280 U(I,K)=SCALE*U(I,K)
  290 ANORM=AMAX1(ANORM,ABS(W(I))+ABS(RV1(I)))
  300 CONTINUE
      IF(.NOT.MATV) GO TO 410
      DO 400 II=1,N
      I=N+1-II
      IF(I.EQ.N ) GO TO 390
      IF(G.EQ.0.) GO TO 360
      DO 320 J=L,N
  320 V(J,I)=(U(I,J)/U(I,L))/G
      DO 350 J=L,N
      S=0.
      DO 340 K=L,N
  340 S=S+U(I,K)*V(K,J)
      DO 350 K=L,N
  350 V(K,J)=V(K,J)+S*V(K,I)
  360 DO 380 J=L,N
      V(I,J)=0.
      V(J,I)=0.
  380 CONTINUE
  390 V(I,I)=1.
      G=RV1(I)
      L=I
  400 CONTINUE
  410 IF(.NOT.MATU) GO TO 510
      MN=N
      IF(M.LT.N)MN=M
      DO 500 II=1,MN
      I=MN+1-II
      L=I+1
      G=W(I)
      IF(I.EQ.N)  GO TO 430
      DO 420 J=L,N
  420 U(I,J)=0.
  430 IF(G.EQ.0.) GO TO 475
      IF(I.EQ.MN) GO TO 460
      DO 450 J=L,N
      S=0.
      DO 440 K=L,M
  440 S=S+U(K,I)*U(K,J)
      F=(S/U(I,I))/G
      DO 450 K=I,M
      U(K,J)=U(K,J)+F*U(K,I)
  450 CONTINUE
  460 DO 470 J=I,M
  470 U(J,I)=U(J,I)/G
      GO TO 490
  475 DO 480 J=I,M
  480 U(J,I)=0.
  490 U(I,I)=U(I,I)+1.
  500 CONTINUE
  510 DO 700 KK=1,N
      K1=N-KK
      K=K1+1
      ITS=0
  520 DO 530 LL=1,K
      L1=K-LL
      L=L1+1
      IF(ABS(RV1(L))+ANORM.EQ.ANORM)  GO TO 565
      IF(ABS(W(L1)) +ANORM.EQ.ANORM)  GO TO 540
  530 CONTINUE
  540 C=0.
      S=1.
      DO 560 I=L,K
      F=S*RV1(I)
      RV1(I)=C*RV1(I)
      IF(ABS(F)+ANORM.EQ.ANORM) GO TO 565
      G=W(I)
      H=SQRT(F*F+G*G)
      W(I)=H
      C= G/H
      S=-F/H
      IF(.NOT.MATU) GO TO 560
      DO 550 J=1,M
      Y=U(J,L1)
      Z=U(J,I)
      U(J,L1)=Y*C+Z*S
      U(J,I)=-Y*S+Z*C
  550 CONTINUE
  560 CONTINUE
  565 Z=W(K)
      IF( L .EQ. K) GO TO  650
      IF(ITS.EQ.300) GO TO 1000
      ITS=ITS+1
      X=W(L)
      Y=W(K1)
      G=RV1(K1)
      H=RV1(K)
      F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.*H*Y)
      G=SQRT(F*F+1.)
      F=((X-Z)*(X+Z)+H*(Y/(F+SIGN(G,F))-H))/X
C
      C=1.
      S=1.
      DO 600 I1=L,K1
      I=I1+1
      G=RV1(I)
      Y=W(I)
      H=S*G
      G=C*G
      Z=SQRT(F*F+H*H)
      RV1(I1)=Z
      C=F/Z
      S=H/Z
      F= X*C+G*S
      G=-X*S+G*C
      H=Y*S
      Y=Y*C
      IF(.NOT.MATV) GO TO 575
      DO 570 J=1,N
      X=V(J,I1)
      Z=V(J,I)
      V(J,I1)=X*C+Z*S
      V(J,I)=-X*S+Z*C
  570 CONTINUE
  575 Z=SQRT(F*F+H*H)
      W(I1)=Z
      IF(Z.EQ.0.)   GO TO 580
      C=F/Z
      S=H/Z
  580 F=C*G+S*Y
      X=-S*G+C*Y
      IF(.NOT.MATU) GO TO 600
      DO 590 J=1,M
      Y=U(J,I1)
      Z=U(J,I)
      U(J,I1)=Y*C+Z*S
      U(J,I)=-Y*S+Z*C
  590 CONTINUE
  600 CONTINUE
      RV1(L)=0.
      RV1(K)=F
      W(K)  =X
      GO TO 520
  650 IF(Z.GE.0.)   GO TO 700
      W(K)=-Z
      IF(.NOT.MATV) GO TO 700
      DO 690 J=1,N
  690 V(J,K)=-V(J,K)
  700 CONTINUE
      GO TO 1001
 1000 IERR=K
C
 1001 RETURN
      END