module spharmt.f90

fortran程序

module spharmt
!
! fortran95 spherical harmonic module.
! For gaussian grids and triangular truncation only.
!
! version 1.1  9/30/2003 (second version - now all fortran using fft99)
! Jeff Whitaker <[e-mail]>

!-----------------------------------------------------------------------------
! to use this module, put "use spharmt" at top of main program
! (just after program statement).
! to initialize arrays needed to compute spherical harmonic transforms
! at T(ntrunc) resolution on a nlon x nlat gaussian grid, 
! do something like this:
! 
! type (sphere) :: sphere_dat
! parameter (nlon=192,nlat=nlon/2)
! parameter (ntrunc=62)
! rsphere = 6.3712e6
! call spharmt_init(sphere_dat,nlon,nlat,ntrunc,re)
 
! derived data type "sphere" contains stuff needed for legendre
! transforms and fft.  By creating multiple instances of sphere
! data types, spherical harmonic transforms on multiple grids can
! be done easily in the same code.
!
! Components of sphere derived data type are:
!
! NLON (integer) - number of longitudes
! NLAT (integer) - number of latitudes
! NTRUNC (integer) - triangular truncation limit
! RE (real) - radius of sphere (m).
! PNM (real pointer array dimension ((ntrunc+1)*(ntrunc+2)/2, nlat) - 
! Associated legendre polynomials.
! HNM (real pointer array, same size as pnm) = -(1 - x*x) * d(pnm)/dx  
! at x = sin(latitude).
! GAULATS (real pointer array dimension nlat) - sin(gaussian latitudes).
! WEIGHTS (real pointer array dimension nlat) - gaussian weights.
! GWRC (real pointer array dimension nlat) - gaussian weights divided by 
! rsphere*(1-x**2).
! INDXM (integer pointer array dimension (ntrunc+1)*(ntrunc+2)/2) - value of
! zonal wavenumber m corresponding to spectral array index 
! nm=1,(ntrunc+1)*(ntrunc+2)/2)
! INDXN (integer pointer array same size as indxm) - value of
! spherical harmonic degree n corresponding to spectral array index 
! nm=1,(ntrunc+1)*(ntrunc+2)/2)
! LAP (real pointer array dimension (ntrunc+1)*(ntrunc+2)/2) - 
! lapacian operator in spectral space = 
! -n*(n+1)/rsphere**2, where n is degree of spherical harmonic.
! ILAP (real pointer array same size as lap) - inverse laplacian operator in 
! spectral space.
! TRIGS, IFAX: arrays needed by Temperton FFT.
! ISINITIALIZED (logical) - true if derived data type has been
! initialized by call to spharm_init, false if not initialized.
! 

! public routines:
!
! SPHARMT_INIT(sphere_dat,nlon,nlat,ntrunc,re):  initialize a sphere object
! (sphere_dat - derived data type "sphere"). inputs are nlon (number of unique
! longitudes), nlat (number of gaussian latitudes), and re 
! (radius of sphere in m). Must be called before anything else.
!
! SPHARMT_DESTROY(sphere_dat):  cleans up pointer arrays allocated by 
! spharmt_init.
!
! SPHARM(sphere_dat, ugrid, anm, idir):  spherical harmonic transform 
! (forward, i.e. grid to spectral, for idir=1 and backward for idir=-1). 
! Arguments are sphere derived data type (sphere_dat), gridded data (ugrid),
! complex spectral coefficients (anm), and flag specifying direction of 
! transform (idir).  See "Import Details" below for information 
! about indexing of grid and spectral arrays.

! GAULW(gaulats,weights): compute sin(gaussian latitudes) and gaussian
! weights.  Number of latitudes determined by size of input arrays.

! LEGEND(x,pnm,hnm): compute associated legendre functions (pnm) and
! their derivates hnm = (X**2-1)*D(PNM)/DX at x = sin(latitude). The
! input arrays pnm and hnm should have dimension (ntrunc+1)*(ntrunc+2)/2, 
! where ntrunc is triangular truncation limit (see Import Details below
! for a description of the spectral indexing).

! GETUV(sphere_dat,vrtnm,divnm,ug,vg): get U,V (u*cos(lat),v*cos(lat) on grid) 
! from spectral coefficients of vorticity and divergence. 
! Input:  sphere_dat, spectral coeffs. of vort and div (vrtnm,divnm).
! Output: gridded U,V (ug,vg).
!
! GETVRTDIV(sphere_dat,vrtnm,divnm,ug,vg): get spectral coefficients of 
! vorticity and divergence (vrtnm,divnm) from gridded U,V (ug,vg).
!
! COSGRAD(sphere_dat,chinm,uchig,vchig): compute cos(lat)*vector gradient.
! Inputs: sphere_dat, spectral coefficient array (chinm).
! Outputs: longitudinal and latitudinal components of gradient on the gaussian
! grid (uchig,vchig).

! SPECSMOOTH(sphere_dat,datagrid,smooth): isotropic spectral smoothing.
! Inputs: sphere_dat, smooth(ntrunc+1) (smoothing factor as a function
! of spherical harmonic degree), datagrid (gridded data to be smoothed).
! Outputs: datagrid (smoothed gridded data).
!
! SUMNM(sphere_dat,am,bm,anm,isign1,isign2):  
! given the arrays of fourier coeffs, am and bm,
! compute the complex spherical harmonic coeffs (anm) of:  
! isign1*( (1./rsphere*(1.-x**2))*d(ag)/d(lon) + (isign2/rsphere)*d(bg)/dx )
! where ag and bg are the grid pt counterparts of am,bm,
! isign1,isign2 are +1 or -1, rsphere is radius of sphere, x=sin(lat))
! am,bm can be computed from gridded data (ag,bg) using RFFT.
!
! RFFT(sphere_dat, data, coeff, idir): computes fourier harmonics in zonal
! direction of a gridded array. idir=+1 for forward (grid
! to fourier) and -1 for backward (fourier to grid) transform.
! data (nlon,nlat) contains gridded data, coeff (ntrunc+1,nlat) contains
! complex zonal fourier harmonics.
! 
! Important Details:
!
! The gridded data is assumed to be oriented such that i=1 is the Greenwich
! meridian and j=1 is the northernmost point. Grid indices increase eastward 
! and southward. The grid increment in longitude is 2*pi/nlon radians.
! For example nlon = 72 for a five degree grid. nlon must be greater than or
! equal to 4. The efficiency of the computation is improved when nlon is a
! product of small prime numbers.

! The spectral data is assumed to be in a complex array of dimension
! (NTRUNC+1)*(NTRUNC+2)/2. NTRUNC is the triangular truncation limit (NTRUNC=42
! for T42). NTRUNC must be <= nlon/2. Coefficients are ordered so that first
! (nm=1) is m=0,n=0, second is m=0,n=1, nm=mtrunc is m=0,n=mtrunc, nm=mtrunc+1 
! is m=1,n=1, etc. In Fortran90 syntax, values of m (degree) and n (order)
! corresponding as a function of the index nm are:

! indxm = (/((m,n=m,mtrunc),m=0,mtrunc)/)
! indxn = (/((n,n=m,mtrunc),m=0,mtrunc)/)

! Conversely, the index nm as a function of m and n is:

! nm = sum((/(i,i=mtrunc+1,mtrunc-m+2,-1)/))+n-m+1

! The associated legendre polynomials are normalized so that the integral  
! (pbar(n,m,theta)**2)*sin(theta) on the interval theta=0 to theta=pi is 1, 
! where: pbar(m,n,theta) = sqrt((2*n+1)*factorial(n-m)/(2*factorial(n+m)))
! *sin(theta)**m/(2**n*factorial(n)) times the (n+m)th derivative of 
! (x**2-1)**n with respect to x=cos(theta).

! note: theta = 0.5*pi - phi, where phi is latitude and theta is colatitude,
! Therefore, cos(theta) = sin(phi) and sin(theta) = cos(phi).

! Note that pbar(0,0,theta)=sqrt(2)/2, 
! and pbar(1,0,theta)=0.5*sqrt(6.)*sin(lat).
!----------------------------------------------------------------------------

! explicit typing

 implicit none

! everything private to module, unless otherwise specified.

 private
 public :: sphere,spharmt_init,spharmt_destroy,spharm,rfft,cosgrad,&
           specsmooth,getuv,getvrtdiv,sumnm,gaulw,legend
 
 type sphere

! rsphere is radius of sphere in m.
   real    :: rsphere = 0.

! nlons is number of longitudes (not including cyclic point).
   integer :: nlons = 0

! nlats is number of gaussian latitudes.
   integer :: nlats = 0

! for fft.
   real, dimension(:), pointer :: trigs
   integer, dimension(:), pointer :: ifax

! ntrunc is triangular truncation limit (e.g. 42 for T42 truncation)
   integer :: ntrunc = 0

! gaulats is sin(gaussian lats), weights are gaussian weights,
! gwrc is gaussian weights divided by re*cos(lat)**2, lap is 
! the laplacian operatore -n*(n+1)/re**2 (where n is the degree 
! of the spherical harmonic),
! ilap is the inverse laplacian (1/lap, set to zero for n=0).
   real, dimension(:), pointer :: gaulats, weights, gwrc, lap, ilap

! pnm is associated legendre polynomials, hnm is -(1-x**2)d(pnm)/dx,
! where x = gaulats.
   real, dimension(:,:), pointer :: pnm, hnm

! indxm is zonal wavenumber index, indxn is spherical harmonic degree index.
   integer, dimension(:), pointer :: indxm, indxn
   logical :: isinitialized=.FALSE.

 end type sphere

 contains

 subroutine spharmt_init(sphere_dat,nlon,nlat,ntrunc,re)

! initialize a sphere object.

   integer, intent(in) :: nlon,nlat,ntrunc
   real, intent(in) :: re
   type (sphere), intent(inout) :: sphere_dat
   real, dimension(:), allocatable :: pnm_tmp,hnm_tmp
   double precision, dimension(:), allocatable :: gaulats_tmp,weights_tmp
   integer :: nmdim,m,n,j

   sphere_dat%nlons = nlon
   sphere_dat%nlats = nlat
   sphere_dat%ntrunc = ntrunc
   sphere_dat%rsphere = re

   nmdim = (ntrunc+1)*(ntrunc+2)/2

   allocate(gaulats_tmp(nlat))
   allocate(weights_tmp(nlat))
   call gaulw(gaulats_tmp,weights_tmp)
   allocate(sphere_dat%gwrc(nlat))
   sphere_dat%gwrc = weights_tmp/(dble(re)*(1.d0-gaulats_tmp**2))
   allocate(sphere_dat%gaulats(nlat))
   allocate(sphere_dat%weights(nlat))
   sphere_dat%gaulats = gaulats_tmp
   sphere_dat%weights = weights_tmp
   deallocate(weights_tmp)

   allocate(sphere_dat%indxm(nmdim))
   allocate(sphere_dat%indxn(nmdim))
   allocate(sphere_dat%lap(nmdim))
   allocate(sphere_dat%ilap(nmdim))

   sphere_dat%indxm = (/((m,n=m,ntrunc),m=0,ntrunc)/)
   sphere_dat%indxn = (/((n,n=m,ntrunc),m=0,ntrunc)/)
   sphere_dat%lap(:)=-real(sphere_dat%indxn(:))*real(sphere_dat%indxn(:)+1)/re**2
   sphere_dat%ilap(1) = 0.
   sphere_dat%ilap(2:nmdim) = 1./sphere_dat%lap(2:nmdim)

   allocate(sphere_dat%pnm(nmdim,nlat))
   allocate(sphere_dat%hnm(nmdim,nlat))
   allocate(pnm_tmp(nmdim))
   allocate(hnm_tmp(nmdim))
   do j=1,nlat
      call legend(gaulats_tmp(j),pnm_tmp,hnm_tmp)
      sphere_dat%pnm(:,j) = pnm_tmp(:)
      sphere_dat%hnm(:,j) = hnm_tmp(:)
   enddo
   deallocate(gaulats_tmp)
   deallocate(pnm_tmp)
   deallocate(hnm_tmp)

   allocate(sphere_dat%trigs((3*nlon/2)+1))
   allocate(sphere_dat%ifax(13))
   call set99(sphere_dat%trigs,sphere_dat%ifax,nlon)

   sphere_dat%isinitialized = .TRUE.

 end subroutine spharmt_init

 subroutine spharmt_destroy(sphere_dat)

! deallocate pointers in sphere object.

   type (sphere), intent(inout) :: sphere_dat

   if (.not. sphere_dat%isinitialized) return

   deallocate(sphere_dat%gaulats)
   deallocate(sphere_dat%weights)
   deallocate(sphere_dat%gwrc)
   deallocate(sphere_dat%pnm)
   deallocate(sphere_dat%hnm)
   deallocate(sphere_dat%indxm)
   deallocate(sphere_dat%indxn)
   deallocate(sphere_dat%lap)
   deallocate(sphere_dat%ilap)
   deallocate(sphere_dat%trigs)
   deallocate(sphere_dat%ifax)

 end subroutine spharmt_destroy

 subroutine gaulw(sinlats, wts)

! compute sin of gaussian latitudes and gaussian weights.
! uses the iterative method presented in Numerical Recipes.

 double precision, intent (inout), dimension(:) :: sinlats, wts

 integer :: itermax
 integer :: i, iter, j, nlat, nprec
 double precision :: pi, pp, p1, p2, p3, z, z1, converg, ten


 ten = 10.d0
 converg = ten * EPSILON(ten)

 nprec = precision(converg)
 converg = .1**nprec

 itermax = 10

 pi = 4.0*ATAN(1.0)

 nlat = size(sinlats)
 if (size(sinlats) .ne. size(wts)) then
   print *, 'sinlats and wts must be same size in gaulw!'
   stop
 end if

 do i=1,nlat
   z = cos(pi*(i - 0.25)/(nlat + 0.5))
   do iter=1,itermax
     p1 = 1.0
     p2 = 0.0

     do j=1,nlat
        p3 = p2
        p2 = p1
        p1 = ((2.0*j - 1.0)*z*p2 - (j - 1.0)*p3)/j
     end do

     pp = nlat*(z*p1 - p2)/(z*z - 1.0E+00)
     z1 = z
     z  = z1 - p1/pp
     if(ABS(z - z1) .LT. converg) go to 10
   end do
   print *, 'abscissas failed to converge in itermax iterations'
   print *, 'stopping in gaulw!'
   stop

10 continue

   sinlats (i)     = z
   wts (i)     = 2.0/((1.0 - z*z)*pp*pp)

 end do

 return
 end subroutine gaulw

 SUBROUTINE LEGEND(X,PMN,HMN)
!
!     THIS SUBROUTINE COMPUTES ASSOCIATED LEGENDRE  
!     FUNCTIONS, PMN, AND THE DERIVATIVE QUANTITY 
!     HMN = -(1 - X*X) * D(PMN)/DX  AT X = COS( COLATITUDE )  
!     GIVEN SOME COLATITUDE IN THE DOMAIN.
!
!     ACCURACY: 
!
!     THE RECURSION RELATION EMPLOYED IS LESS ACCURATE
!     NEAR THE POLES FOR M + N .GE. ROUGHLY 75 BECAUSE THE RECURSION
!     INVOLVES THE DIFFERENCE OF NEARLY EQUAL NUMBERS, SO
!     THE ASSOCIATED LEGENDRE FUNCTIONS ARE COMPUTED USING DOUBLE
!     PRECISION ACCUMULATORS. 
!     SEE BELOUSOV, TABLES OF NORMALIZED ASSOCIATED LEGENDRE  
!     POLYNOMIALS, C. 1962, FOR INFORMATION ON THE ACCURATE 
!     COMPUTATION OF ASSOCIATED LEGENDRE FUNCTIONS.
 
      real, dimension(:), intent(inout) ::  PMN,  HMN
      double precision, intent(in) :: x
      integer :: m,n,nm,i,nmax,np1,nmstrt,j,nmdim,ntrunc
      double precision :: A, B, PROD, SINSQ, &
      eps,epsnmp1,EPSPMN, PMNJ, PMNJM1, PMNJM2
 
 
!**** SET PARAMETERS FOR ENTRY INTO THE RECURSIVE FORMULAE. 
 
 
      SINSQ = 1.D0 - X * X
 
      A     = 1.D0
      B     = 0.D0
      PROD  = 1.D0
 
      nmdim = size(pmn)
      ntrunc = nint((-1.+sqrt(1+8*float(nmdim)))/2.)-1
      if (size(pmn) .ne. size(hmn)) then
         print *, 'pnm and hnm must be same size in subroutine legend!'
         stop
      end if
 
!**** LOOP FOR THE 'M' INDEX. 
 
      NMSTRT = 0
      DO I = 1, NTRUNC+1
 
         M    = (I - 1)
         NMAX = NTRUNC+1-M
 
!        WHEN M=0 (I=1), STANDARD LEGENDRE POLYNOMIALS ARE
!        GENERATED. 
 
         IF (M .NE. 0) THEN
               A    = A + 2.D0
               B    = B + 2.D0
               PROD = PROD * SINSQ * A / B 
         END IF
 
!****    GENERATE PMN AND HMN FOR J = 1 AND 2.
 
         PMNJM2   = SQRT(0.5D0 * PROD)
         NM = NMSTRT + 1
         PMN(NM) = PMNJM2