splaplac.f

Go to the documentation of this file.

C> @file
C> @brief Compute laplacian in spectral space.
C> @author Iredell @date 92-10-31
 
C> Computes the laplacian or the inverse laplacian
C> of a scalar field in spectral space.
C>      
C> Subprogram speps() should be called already.
C>      
C> The Laplacian of Q(L,N) is simply -N*(N+1)/A**2*Q(L,N)
C>
C> @param I spectral domain shape
C> (0 for triangular, 1 for rhomboidal)
C> @param M spectral truncation
C> @param ENN1 N*(N+1)/A**2
C> @param[out] Q if IDIR > 0, scalar field
C> (Q(0,0) is not computed)
C> @param[out] QD2 if IDIR < 0, Laplacian
C> @param IDIR flag
C> - IDIR > 0 to take Laplacian
C> - IDIR < 0 to take inverse Laplacian
C>
C> @author Iredell @date 92-10-31
      SUBROUTINE splaplac(I,M,ENN1,Q,QD2,IDIR)
 
      REAL ENN1((M+1)*((I+1)*M+2)/2)
      REAL Q((M+1)*((I+1)*M+2))
      REAL QD2((M+1)*((I+1)*M+2))
 
C  TAKE LAPLACIAN
      IF(idir.GT.0) THEN
        k=1
        qd2(2*k-1)=0.
        qd2(2*k)=0.
        DO k=2,(m+1)*((i+1)*m+2)/2
          qd2(2*k-1)=q(2*k-1)*(-enn1(k))
          qd2(2*k)=q(2*k)*(-enn1(k))
        ENDDO
 
C  TAKE INVERSE LAPLACIAN
      ELSE
        DO k=2,(m+1)*((i+1)*m+2)/2
          q(2*k-1)=qd2(2*k-1)/(-enn1(k))
          q(2*k)=qd2(2*k)/(-enn1(k))
        ENDDO
      ENDIF
 
      RETURN
      END