PROGRAM Hough_Functions

USE m_Parameters

USE m_Arrays

USE m_Orthogonality

USE m_Vectorial_Hough_Functions

IMPLICIT NONE
 
! Computes Hough Functions for a Given Equivalent Depth

! This Procedure is Hardwired for Odd Rhomboidal Truncation,
! Then Mmax and na Must Be Even and nm Must Be a Multiple of 6.
! jm Must Be Even Also. k0 Must Be At Least 6.

INTEGER (KIND=ki8) :: s, ks, n, m, isyas, kind, j

REAL (KIND=kr8) :: sumg

CHARACTER (LEN=37) :: hghfnt

! Get Parameters and Array

CALL Get_Parameters

CALL Get_Arrays

! Compute Parameters for Associated Legendre Functions

CALL Gamma_Epsilon

! Compute Gaussian Latitudes and Weigth

CALL Get_Gaussian_Grid

hghfnt='Hough_Functions'//exth//'_h   '
WRITE (UNIT=hghfnt(35:37), FMT='(I3.3)') nHn
WRITE (UNIT=0, FMT='(1X,A)') hghfnt

WRITE (UNIT=*, FMT=*) ' '
WRITE (UNIT=*, FMT=5) Hn, na, nb, 0, Mmax-1, jg

IF (test) THEN
  WRITE (UNIT=*, FMT='(1X,A)') hghfnt//'.num'
  OPEN  (UNIT=86, FILE=TRIM(dirhgf)//hghfnt//'.num', STATUS='UNKNOWN')
  WRITE (UNIT=86, FMT=*) ' '
  WRITE (UNIT=86, FMT=5) Hn, na, nb, 0, Mmax-1, jg
  WRITE (UNIT=*, FMT='(1X,A)') hghfnt//'.anl'
  OPEN  (UNIT=87, FILE=TRIM(dirhgf)//hghfnt//'.anl', STATUS='UNKNOWN')
  WRITE (UNIT=87, FMT=*) ' '
  WRITE (UNIT=87, FMT=5) Hn, na, nb, 0, Mmax-1, jg
  WRITE (UNIT=*, FMT='(1X,A)') hghfnt//'.cft'
  OPEN  (UNIT=88, FILE=TRIM(dirhgf)//hghfnt//'.cft', STATUS='UNKNOWN')
  WRITE (UNIT=88, FMT=*)' '
  WRITE (UNIT=88, FMT=5) Hn, na, nb, 0, Mmax-1, jg
END IF

5 FORMAT ('  Hn = ',1PG12.4,0P,'  na = ',I3,'  nb = ',I3, &
          '  s0 = ',I3,'  sm = ',I3,'  jg = ',I3,/)
WRITE (UNIT=*, FMT=10)
10 FORMAT (10X,'Latitude',11X,'sin(Latitude)',16X,'Gaussian Weight',/)
DO j=1,jg
  WRITE (UNIT=*, FMT=20) rad*ASIN(yg(j)), yg(j), wg(j)
END DO
20 FORMAT (10X,F7.4,5X,F25.21,5X,F25.21)
sumg=2.0_kr8*SUM(wg)
WRITE (UNIT=*, FMT='(/,A,F25.21,1PG13.5,/)') ' SUM(wg) : ', sumg, 2.0_kr8-sumg

! Begin to Storing the Hough Function on a Disk File Unit=77

WRITE (UNIT=*, FMT='(1X,A)') hghfnt
OPEN  (UNIT=77, FILE=TRIM(dirhgf)//hghfnt,STATUS='UNKNOWN', FORM='UNFORMATTED')

WRITE (UNIT=77) epsl, nb, jg
WRITE (UNIT=77) yg, wg

ks=1
s=ks-1

CALL Get_Legendre_Functions (ks)

isyas=0

WRITE (UNIT=*, FMT=*) ' '
WRITE (UNIT=*, FMT=30) s, isyas, epsl, nb
WRITE (UNIT=*, FMT=*) ' '
IF (test) THEN
  WRITE (UNIT=88, FMT=*) ' '
  WRITE (UNIT=88, FMT=30) s, isyas, epsl, nb
END IF
WRITE(UNIT=0, FMT=35) s, isyas, epsl, nb, nHn, Hn

30 FORMAT (' s = ',I3,'  sym = ',I3,'  eps = ',F9.1,'  nb = ',I3)
35 FORMAT (' s = ',I3,'  sym = ',I3,'  eps = ',F9.1,'  nb = ',I3,'  nHn = ', I3,'  Hn = ',1PG12.5)

DO n=1,nm
  DO m=1,no
    ans(n,m)=0.0_kr8
    bns(n,m)=0.0_kr8
    cns(n,m)=0.0_kr8
  END DO
END DO

CALL s0f0sy
CALL s0fnsy 
CALL normsy (0_ki8)
CALL uvzsy (0_ki8)
CALL energ
IF (test) THEN
  CALL Orthogonality_Symmetric (0_ki8)
  CALL Orthogonality_Test (0_ki8, ' sym. ')
  CALL Print_Coefficients (ans, 0_ki8, ' sym. ', 'ans')
  CALL Print_Coefficients (bns, 0_ki8, ' sym. ', 'bns')
  CALL Print_Coefficients (cns, 0_ki8, ' sym. ', 'cns')
END IF

! Storing the Hough Functions on Disk File Unit=77

DO n=1,no
  kind=1+n/nb-(MOD(n-1,nb)+1)/nb
  WRITE (UNIT=77) s, isyas, kind, fr(n), &
                  (uh(n,j),j=1,jg), (vh(n,j),j=1,jg), (zh(n,j),j=1,jg), &
                  (duh(n,j),j=1,jg), (dvh(n,j),j=1,jg), (dzh(n,j),j=1,jg)
!                ,(d2uh(n,j),j=1,jg), (d2vh(n,j),j=1,jg), (d2zh(n,h),j=1,jg)

! Printout of Spectral Energy

  WRITE (UNIT=*, FMT=40) fr(n), et(n), akeu(n), akev(n), ape(n)
  40 FORMAT (' sig=',1PG10.3,' et=',G15.8,' keu=',G10.3, &
             ' kev=',G10.3,' ape=',G10.3)
  IF (ABS(1.0-et(n)) .GT. 1.0E-7) WRITE (UNIT=0, FMT=40) fr(n), et(n), akeu(n), akev(n), ape(n)
END DO

isyas=1

WRITE (UNIT=*, FMT=*) ' '
WRITE (UNIT=*, FMT=30) s, isyas, epsl, nb
WRITE (UNIT=*, FMT=*) ' '
IF (test) THEN
  WRITE (UNIT=88, FMT=*) ' '
  WRITE (UNIT=88, FMT=30) s, isyas, epsl, nb
END IF
WRITE (UNIT=0, FMT=35) s, isyas, epsl, nb, nHn, Hn

DO n=1,nm
  DO m=1,no
    ans(n,m)=0.0_kr8
    bns(n,m)=0.0_kr8
    cns(n,m)=0.0_kr8
  END DO
END DO

CALL s0f0as
CALL s0fnas
CALL normas (0_ki8)
CALL uvzas (0_ki8)
CALL energ
IF (test) THEN
  CALL Orthogonality_Asymmetric(0_ki8)
  CALL Orthogonality_Test(0_ki8, ' asym.')
  CALL Print_Coefficients(ans, 0_ki8, ' asym.', 'ans')
  CALL Print_Coefficients(bns, 0_ki8, ' asym.', 'bns')
  CALL Print_Coefficients(cns, 0_ki8, ' asym.', 'cns')
END IF

! Storing the Hough Functions on Disk File unit=77

DO n=1,no
  kind=1+n/nb-(MOD(n-1,nb)+1)/nb
  WRITE (UNIT=77) s, isyas, kind, fr(n), &
                  (uh(n,j),j=1,jg),(vh(n,j),j=1,jg), (zh(n,j),j=1,jg), &
                  (duh(n,j),j=1,jg), (dvh(n,j),j=1,jg), (dzh(n,j),j=1,jg)
!                ,(d2uh(n,j),j=1,jg), (d2vh(n,j),j=1,jg), (d2zh(n,h),j=1,jg)

! Printout of Spectral Energy

  WRITE (UNIT=*, FMT=40) fr(n), et(n), akeu(n), akev(n), ape(n)
END DO

DO ks=2,Mmax
  s=ks-1

  CALL Get_Legendre_Functions (ks)

  isyas=0

  WRITE (UNIT=*, FMT=*) ' '
  WRITE (UNIT=*, FMT=30) s, isyas, epsl, nb
  WRITE (UNIT=*, FMT=*) ' '
  IF (test) THEN
    WRITE (UNIT=88, FMT=*) ' '
    WRITE (UNIT=88, FMT=30) s, isyas, epsl, nb
  END IF
  WRITE (UNIT=0, FMT=35) s, isyas, epsl, nb, nHn, Hn

  DO n=1,nm
    DO m=1,no
      ans(n,m)=0.0_kr8
      bns(n,m)=0.0_kr8
      cns(n,m)=0.0_kr8
    END DO
  END DO

  CALL snfnsy (s)
  CALL normsy (s)
  CALL uvzsy (s)
  CALL energ
  IF (test) THEN
    CALL Orthogonality_Symmetric (s)
    CALL Orthogonality_Test (s, ' sym. ')
    CALL Print_Coefficients (ans, s, ' sym. ', 'ans')
    CALL Print_Coefficients (bns, s, ' sym. ', 'bns')
    CALL Print_Coefficients (cns, s, ' sym. ', 'cns')
  END IF

! Storing the Hough Functions on Disk File unit=77

  DO n=1,no
    kind=1+n/nb-(MOD(n-1,nb)+1)/nb
    WRITE (UNIT=77) s, isyas, kind,fr(n), &
                    (uh(n,j),j=1,jg), (vh(n,j),j=1,jg), (zh(n,j),j=1,jg), &
                    (duh(n,j),j=1,jg), (dvh(n,j),j=1,jg), (dzh(n,j),j=1,jg)
!                  ,(d2uh(n,j),j=1,jg), (d2vh(n,j),j=1,jg), (d2zh(n,h),j=1,jg)

! Printout of Spectral Energy

    WRITE (UNIT=*, FMT=40) fr(n), et(n), akeu(n), akev(n), ape(n)
  END DO

  isyas=1

  WRITE (UNIT=*, FMT=*) ' '
  WRITE (UNIT=*, FMT=30) s, isyas, epsl, nb
  WRITE (UNIT=*, FMT=*) ' '
  IF (test) THEN
    WRITE (UNIT=88, FMT=*) ' '
    WRITE (UNIT=88, FMT=30) s, isyas, epsl, nb
  END IF
  WRITE (UNIT=0, FMT=35) s, isyas, epsl, nb, nHn, Hn

  DO n=1,nm
    DO m=1,no
      ans(n,m)=0.0_kr8
      bns(n,m)=0.0_kr8
      cns(n,m)=0.0_kr8
    END DO
  END DO

  CALL snfnas (s)
  CALL normas (s)
  CALL uvzas (s)
  CALL energ
  IF (test) THEN
    CALL Orthogonality_Asymmetric (s)
    CALL Orthogonality_Test (s, ' asym.')
    CALL Print_Coefficients (ans, s, ' asym.', 'ans')
    CALL Print_Coefficients (bns, s, ' asym.', 'bns')
    CALL Print_Coefficients (cns, s, ' asym.', 'cns')
  END IF

! Storing the Hough Functions on Disk File unit=77

  DO n=1,no
    kind=1+n/nb-(MOD(n-1,nb)+1)/nb
    WRITE(UNIT=77) s, isyas, kind, fr(n), &
                   (uh(n,j),j=1,jg), (vh(n,j),j=1,jg), (zh(n,j),j=1,jg), &
                   (duh(n,j),j=1,jg), (dvh(n,j),j=1,jg), (dzh(n,j),j=1,jg)
!                 ,(d2uh(n,j),j=1,jg), (d2vh(n,j),j=1,jg), (d2zh(n,h),j=1,jg)

!  Printout of Spectral Energy

    WRITE (UNIT=*, FMT=40) fr(n), et(n), akeu(n), akev(n), ape(n)
  END DO

END DO

CLOSE (UNIT=77)
WRITE (UNIT=*, FMT=*) ' '
IF (test) THEN
  WRITE (UNIT=86, FMT=*) ' '
  CLOSE (UNIT=86)
  WRITE (UNIT=87, FMT=*) ' '
  CLOSE (UNIT=87)
  WRITE (UNIT=88, FMT=*) ' '
  CLOSE (UNIT=88)
END IF

CALL Free_Arrays


CONTAINS


SUBROUTINE Gamma_Epsilon
 
  ! Computes Parameters gams and eps Used for the
  ! Associated Legendre Functions Calcultions

  IMPLICIT NONE

  INTEGER (KIND=ki8) :: k, n, s, m

  gams(1)=SQRT(0.5_kr8)
  DO  k=1,Mmax-1
    s=k-1
    gams(k+1)=gams(k)*SQRT(REAL(2*s+3,kr8))/SQRT(REAL(2*s+2,kr8))
  END DO
  DO  k=1,Mmax
    s=k-1
    DO  n=1,ne
      m=s+n-1
      eps(k,n)=SQRT(REAL(m*m-s*s,kr8)/REAL(4*m*m-1,kr8))
    END DO
  END DO

END SUBROUTINE Gamma_Epsilon


SUBROUTINE Get_Legendre_Functions (ks)

  ! Computes the Normalized Associated Legendre Functions
  !          and its Derivatives

  IMPLICIT NONE

  INTEGER (KIND=ki8), INTENT (IN) :: ks

  INTEGER (KIND=ki8) :: k, j, n, s

  REAL (KIND=kr8) :: p11, rn

  DO  j=1,jg
    DO  k=1,ks
      s=k-1
      psm(1,j)=gams(k)
      IF (s /= 0) THEN
        p11=1.0_kr8
        DO  n=1,s
          p11=p11*xg(j)
        END DO
        psm(1,j)=psm(1,j)*p11
      END IF
    END DO
    psm(2,j)=yg(j)*psm(1,j)/eps(ks,2)
    DO  n=3,np
      psm(n,j)=(yg(j)*psm(n-1,j)-eps(ks,n-1)*psm(n-2,j))/eps(ks,n)
    END DO
    n=1
    rn=REAL(s+n-1,kr8)
    qsm(n,j)=-rn*eps(ks,2)*psm(2,j)
    DO  n=2,nq
      rn=REAL(s+n-1,kr8)
      qsm(n,j)=(rn+1.0_kr8)*eps(ks,n)*psm(n-1,j)-rn*eps(ks,n+1)*psm(n+1,j)
    END DO
  END DO

END SUBROUTINE Get_Legendre_Functions


SUBROUTINE Get_Gaussian_Grid
 
  ! Computes the Gaussian Latitudes and Weights

  IMPLICIT NONE

  INTEGER (KIND=ki8) :: Jmaxh, Jmax, k1, k, Jmaxh2, kk

  REAL (KIND=kr8) :: tol, cw, dyi, y0, dy0, p2, p1, yk, xk, wk

  tol=EPSILON(tol)
  Jmaxh=jg
  Jmax=2*Jmaxh
  k1=Jmax-1
  cw=2.0/REAL(Jmax*Jmax,kr8)
  dyi=MAX(1.0_kr8,0.5_kr8*REAL(jg,kr8)-10.0_kr8)
  dyi=((ATAN(1.0_kr8)/90.0_kr8)/dyi)/10.0_kr8
  y0=0.0_kr8
  DO k=1,Jmaxh
    dy0=dyi
    10 CALL Legendre_Polynomials (Jmax, y0, p2)
    20 p1=p2
    y0=y0+dy0
    CALL Legendre_Polynomials (Jmax, y0, p2)
    IF (SIGN(1.0_kr8,p1) == SIGN(1.0_kr8,p2)) GO TO 20
    IF (dy0 > tol) THEN
      y0=y0-dy0
      dy0=0.25_kr8*dy0
      GO TO 10
    END IF
    yg(k)=y0
    xg(k)=SQRT(1.0_kr8-y0*y0)
    CALL Legendre_Polynomials (k1, y0, p1)
    wg(k)=cw*(1.0_kr8-y0*y0)/(p1*p1)
  END DO
  Jmaxh2=Jmaxh/2
  DO k=1,Jmaxh2
    kk=Jmaxh-k+1
    yk=yg(k)
    xk=xg(k)
    wk=wg(k)
    yg(k)=yg(kk)
    xg(k)=xg(kk)
    wg(k)=wg(kk)
    yg(kk)=yk
    xg(kk)=xk
    wg(kk)=wk
  END DO

END SUBROUTINE Get_Gaussian_Grid


SUBROUTINE Legendre_Polynomials (k, y, p)

  ! Computes the Legendre Polynomials

  IMPLICIT NONE

  INTEGER (KIND=ki8), INTENT (IN) :: k

  REAL (KIND=kr8), INTENT (IN) :: y

  REAL (KIND=kr8), INTENT (OUT) :: p

  INTEGER (KIND=ki8) :: i

  REAL (KIND=kr8) :: y1, y2, fi, g2, g1, y3

  y1=1.0_kr8
  y2=y
  DO i=2,k
    fi=1.0_kr8/REAL(i,kr8)
    g2=y*y2*REAL(i+i-1,kr8)
    g1=y1*REAL(i-1,kr8)
    y3=(g2-g1)*fi
    y1=y2
    y2=y3
  END DO
  p=y3

END SUBROUTINE Legendre_Polynomials


END PROGRAM Hough_Functions