C==============================================================
C
c    calculatkn 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 *****************************      
      PARAMETER (NTIME=1201)  
      PARAMETER (NlON=64)
      PARAMETER (NlAT=36)
      parameter (kgit=56)
      parameter (nvar=1)
c
      REAl  dolg(nlon),Geopot(nlon),gh_A0(nlat,kgit)
      REAl  gh_A1(nlat,kgit),  gh_A2(nlat,kgit)
      REAl  gh_A3(nlat,kgit)
      REAl  gh_F1(nlat,kgit),  gh_F2(nlat,kgit)
      REAl  gh_F3(nlat,kgit)
      REAl  gh_A4(nlat,kgit),  gh_F4(nlat,kgit)
      REAl  aus(nlon,nlat,kgit,nvar),out1(nlat,kgit,2)
      REAl  AM1(NDIM,NF),A(NDIM,NF),U(NDIM,NF),V(NDIM,NF),WORK(NF)
      REAl  FUN(NF,nlon),SIGMA(NF),AlFA(NF),SQ(NF),AG(NG),FG(NG)

c      PI=2.*ASIN(1.)

C....OUTPUTFIlES
          OPEN (100,FIlE='gh_a0_grads.dx',   
     $      form='unformatted',
     *      access='direct', status='unknown',recl=4*nlat*kgit) 
          OPEN (1,  FIlE='gh_m1_grads.dx' ,
     $      form='unformatted',
     *      access='direct', status='unknown',recl=4*nlat*kgit*2)
c      open (11, FIlE='gh_m1a_Jan1_27km.dat')
c      open (12, FIlE='gh_m1f_Jan1_27km.dat')
c      open (13, FIlE='gh_m0_Jan1_27km.dat')
c      open (21, FIlE='gh_m1a_Jan1_50km.dat')
c      open (22, FIlE='gh_m1f_Jan1_50km.dat')
c      open (23, FIlE='gh_m0_Jan1_50km.dat')
          OPEN (2,  FIlE='gh_m2_grads.dx',
     $      form='unformatted',
     *      access='direct', status='unknown',recl=4*nlat*kgit*2) 
          OPEN (3,  FIlE='gh_m3_grads.dx',
     $      form='unformatted',
     *      access='direct', status='unknown',recl=4*nlat*kgit*2)
          OPEN (4,  FIlE='gh_m4_grads.dx'  ,
     $      form='unformatted',
     *      access='direct', status='unknown',recl=4*nlat*kgit*2)  
          open(99,file='err') 

c....INPUTFIlE      
          OPEN (101, FIlE='phi.dx',
     $      form='unformatted',
     *      access='direct', status='old',recl=4*nlon*nlat*kgit*nvar)    
c
      KlEF=0
      mm=1
         print*,'mm=',mm
      DO  l=1,ntime     !time
      read(101,rec=l) aus
         DO  k=1,kgit   !altitudes
c      read(101,rec=l) aus
       do 2001 j=1,nlat
       do 2002 i=1,nlon   
       Geopot(i)=aus(i,j,k,mm)
       dolg(i)=(float(i)-0.9999)*5.625    
2002  continue
c ====================================================
      CAll HARMON(NDIM,NF,AM1,A,U,V,FUN,SIGMA,AlFA,SQ,nlon,DOlG,Geopot,
     $               NG,ACONST,AG,FG,WORK,KlEF)
      KlEF=1
      gh_A0(j,k)=ACONST
      gh_A1(j,k)=AG(1)
      gh_F1(j,k)=FG(1)
c      out1(j,k,1)=AG(1)
c      out1(j,k,2)=FG(1)
      gh_A2(j,k)=AG(2)
      gh_F2(j,k)=FG(2)
      gh_A3(j,k)=AG(3)
      gh_F3(j,k)=FG(3)
      gh_A4(j,k)=AG(4)
      gh_F4(j,k)=FG(4)
2001  continue

      Enddo   !altitudes
c      ENDDO   !time

 
      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   !time
c----------------------------  write the output amplitudes and phases
c      DO l=1,ntime 
c      thour=4.*Float(l-1)
c      td=-31.+thour/24. 
c       write(11,1011) td,(gh_A1(nlat-j+1,10,l),j=1,nlat)
c       write(12,1011) td,(gh_F1(nlat-j+1,10,l),j=1,nlat)
c       write(13,1011) td,(gh_A0(nlat-j+1,10,l),j=1,nlat)
c       write(21,1011) td,(gh_A1(nlat-j+1,18,l),j=1,nlat)
c       write(22,1011) td,(gh_F1(nlat-j+1,18,l),j=1,nlat)
c       write(23,1011) td,(gh_A0(nlat-j+1,18,l),j=1,nlat)
 1011 format(1x,37f9.2)
c      end do
      ClOSE(100)
      ClOSE(1)
      ClOSE(2)
      ClOSE(3)
      ClOSE(4)
c      close (11)
c      close(12)
      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------------ normalizatkn of basic functkns
      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 approximatkn
      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 calculatkn 
      DO 80 J=1,NG
      AR=AlFA(1+J)
      AI=AlFA(1+J+NG)
      AG(J)=SQRT(AR*AR+AI*AI)
c------ the phase - positkn 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.30) 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