speps.f

Go to the documentation of this file.

C> @file
C> @brief Compute utility spectral fields.
C> @author Iredell  @date 92-10-31
 
C> Computes constant fields indexed in the spectral domain
C> in "IBM ORDER" (Zonal wavenumber is the slower index).
C>      
C> If L is the zonal wavenumber and N is the total wavenumber
C> and A is the earth radius, then the fields returned are:
C> - (1) normalizing factor EPSILON=SQRT((N**2-L**2)/(4*N**2-1))
C> - (2) Laplacian factor N*(N+1)/A**2
C> - (3) zonal derivative/Laplacian factor L/(N*(N+1))*A
C> - (4) Meridional derivative/Laplacian factor EPSILON/N*A
C>
C> @param I spectral domain shape (0 for triangular, 1 for rhomboidal)
C> @param M spectral truncation
C> @param EPS ((M+1)*((I+1)*M+2)/2) SQRT((N**2-L**2)/(4*N**2-1))
C> @param EPSTOP (M+1) SQRT((N**2-L**2)/(4*N**2-1)) OVER TOP
C> @param ENN1 ((M+1)*((I+1)*M+2)/2) N*(N+1)/A**2
C> @param ELONN1 ((M+1)*((I+1)*M+2)/2) L/(N*(N+1))*A
C> @param EON ((M+1)*((I+1)*M+2)/2) EPSILON/N*A
C> @param EONTOP (M+1) EPSILON/N*A OVER TOP
C>
C> @author Iredell  @date 92-10-31
      SUBROUTINE speps(I,M,EPS,EPSTOP,ENN1,ELONN1,EON,EONTOP)
      REAL EPS((M+1)*((I+1)*M+2)/2),EPSTOP(M+1)
      REAL ENN1((M+1)*((I+1)*M+2)/2),ELONN1((M+1)*((I+1)*M+2)/2)
      REAL EON((M+1)*((I+1)*M+2)/2),EONTOP(M+1)
      parameter(rerth=6.3712e6,ra2=1./rerth**2)
 
      DO l=0,m
        k=l*(2*m+(i-1)*(l-1))/2+l+1
        eps(k)=0.
        enn1(k)=ra2*l*(l+1)
        elonn1(k)=rerth/(l+1)
        eon(k)=0.
      ENDDO
      DO l=0,m
        DO n=l+1,i*l+m
          k=l*(2*m+(i-1)*(l-1))/2+n+1
          eps(k)=sqrt(float(n**2-l**2)/float(4*n**2-1))
          enn1(k)=ra2*n*(n+1)
          elonn1(k)=rerth*l/(n*(n+1))
          eon(k)=rerth/n*eps(k)
        ENDDO
      ENDDO
      DO l=0,m
        n=i*l+m+1
        epstop(l+1)=sqrt(float(n**2-l**2)/float(4*n**2-1))
        eontop(l+1)=rerth/n*epstop(l+1)
      ENDDO
      RETURN
      END