module spharmt.f90

fortran程序

 
         PMNJM1   = SQRT( DBLE(2 * M + 3) ) * X * PMNJM2
         IF (NM .NE. NMDIM) PMN(NM+1) = PMNJM1
 
         NP1 = M + 1  
         EPSNMP1 = SQRT( DBLE(NP1*NP1 - M*M) / DBLE(4*NP1*NP1 - 1) )
         EPSPMN   = X * PMNJM1 - EPSNMP1 * PMNJM2 
 
         HMN(NM) = DBLE(M) * EPSNMP1 * PMNJM1 
         IF (NM .NE. NMDIM) &
         HMN(NM+1) = DBLE(M+1) * EPSPMN  -  &
         DBLE(M+2) * EPSNMP1 * PMNJM2
 
!****    LOOP FOR THE 'N' INDEX.
!        NOW APPLY THE RECURSION FORMULAE FOR J .GE. 3.
 
         DO J    = 3, NMAX  
            N = M + J - 1
            NM = NMSTRT + J
            EPS = SQRT( DBLE(N*N - M*M) / DBLE(4*N*N - 1) )
            PMNJ     = EPSPMN / EPS 
            PMN(NM) = PMNJ
 
!        COMPUTE EPS * PMN FOR J+1.
 
            EPSPMN   = X * PMNJ - EPS * PMNJM1  
 
!        COMPUTE THE DERIVATIVE.
 
            HMN(NM) = DBLE(N) * EPSPMN -  &
            DBLE(N+1) * EPS * PMNJM1 
 
            PMNJM2   = PMNJM1 
            PMNJM1   = PMNJ
         ENDDO  
         NMSTRT = NMSTRT + NMAX
     ENDDO

 END SUBROUTINE LEGEND

 subroutine rfft(sphere_dat, data, coeff, idir)

! real multiple fft (uses temperton fft991)

 type (sphere), intent(in) :: sphere_dat
 
 real, dimension(sphere_dat%nlons,sphere_dat%nlats), intent(inout) :: data
 real, dimension(sphere_dat%nlats*(sphere_dat%nlons+2)) :: wrk1
 real, dimension(sphere_dat%nlats*(sphere_dat%nlons+1)) :: wrk2
 complex, dimension(sphere_dat%ntrunc+1,sphere_dat%nlats), intent(inout) :: coeff
 integer ::  nlons,nlats,ntrunc,mwaves,i,j,m,n
 integer, intent(in) :: idir

 if (.not. sphere_dat%isinitialized) then
    print *, 'uninitialized sphere object in rfft!'
    stop
 end if


 nlons = sphere_dat%nlons
 nlats = sphere_dat%nlats
 ntrunc = sphere_dat%ntrunc
 if (ntrunc .gt. nlons/2) then
    print *, 'ntrunc must be less than or equal to nlons in rfft'
    stop
 end if

 mwaves = ntrunc+1

!==> forward transform.

 if (idir .eq. +1) then
 
!==> copy the data into the work array.
!    transforms are computed pairwise using a complex fft.
     
   n = 0
   wrk1 = 0.
   do j=1,nlats
   do i=1,nlons+2  
      n = n + 1
      wrk1(n) = 0.0
      if (i .le. nlons) then
         wrk1(n) = data(i,j)
      end if
   enddo
   enddo
 
   call fft991(wrk1,wrk2,sphere_dat%trigs,sphere_dat%ifax,1,nlons+2,nlons,nlats,-1)
 
   n = -1
   do j=1,nlats
   do m=1,(nlons/2)+1
      n = n + 2
      if (m .le. mwaves) then
        coeff(m,j) = cmplx(wrk1(n),wrk1(n+1)) 
      end if
   enddo
   enddo

!==> inverse transform.

 else if (idir .eq. -1) then

   wrk1 = 0.
   n = -1
   do j=1,nlats
   do m=1,(nlons/2)+1
      n = n + 2
      if (m .le. mwaves) then
         wrk1(n) = real(coeff(m,j))
         wrk1(n+1) = aimag(coeff(m,j))
      end if
   enddo
   enddo
 
   call fft991(wrk1,wrk2,sphere_dat%trigs,sphere_dat%ifax,1,nlons+2,nlons,nlats,1)
 
   n = 0
   do j=1,nlats
   do i=1,nlons+2  
      n = n + 1
      if (i .le. nlons) then
         data(i,j) = wrk1(n)
      end if
   enddo
   enddo

 else
    write(6,*) ' idir must be +1 or -1 in RFFT!'
    write(6,*) ' execution terminated.'
    stop
 end if

 end subroutine rfft

 subroutine spharm(sphere_dat, ugrid, anm, idir)

! spherical harmonic transform

 type (sphere), intent(in) :: sphere_dat
 
 real, dimension(sphere_dat%nlons,sphere_dat%nlats), intent(inout) :: ugrid
 complex, dimension((sphere_dat%ntrunc+1)*(sphere_dat%ntrunc+2)/2), intent(inout) :: anm
 complex, dimension(sphere_dat%ntrunc+1,sphere_dat%nlats) :: am
 integer ::  nlats,ntrunc,mwaves,nmstrt,nm,m,n,j
 integer, intent(in) :: idir

 if (.not. sphere_dat%isinitialized) then
    print *, 'uninitialized sphere object in spharm!'
    stop
 end if

 nlats = sphere_dat%nlats
 ntrunc = sphere_dat%ntrunc
 mwaves = ntrunc+1

 IF (IDIR .eq. +1) then
 
!==>  GRID SPACE TO SPECTRAL SPACE TRANSFORMATION 
!     FIRST, INITIALIZE ARRAY.  
 
   anm = 0.
 
!==> perform ffts on each latitude.

   call rfft(sphere_dat, ugrid, am, 1)
 
!==>  SUM OVER ALL GAUSSIAN LATITUDES FOR EACH MODE AND EACH WAVE TO
!     OBTAIN THE TRANSFORMED VARIABLE IN SPECTRAL SPACE.
 
   do j=1,nlats
     NMSTRT = 0
     DO m = 1, mwaves
       DO n = 1, mwaves-m+1
          NM = NMSTRT + n
          anm(NM)=anm(NM)+sphere_dat%pnm(nm,j)*sphere_dat%weights(j)*am(m,j)
       enddo
       nmstrt = nmstrt + mwaves-m+1
     enddo
   enddo
 
!==>  SPECTRAL SPACE TO GRID SPACE TRANSFORMATION.
 
 else if (idir .eq. -1) then
 
 DO J = 1, NLATS 
 
!==>  INVERSE LEGENDRE TRANSFORM TO GET VALUES OF THE ZONAL FOURIER
!     TRANSFORM AT LATITUDE j.
 
!==>  SUM THE VARIOUS MERIDIONAL MODES TO BUILD THE FOURIER SERIES
!     COEFFICIENT FOR ZONAL WAVENUMBER M=I-1 AT the GIVEN LATITUDE. 
 
    NMSTRT = 0
    do m = 1, MWAVES
      am(m,j) = cmplx(0., 0.)
      DO n = 1, mwaves-m+1
         NM = NMSTRT + n
         am(m,j) = am(m,j)  +  anm(NM) * sphere_dat%pnm(NM,j) 
      enddo
      NMSTRT = NMSTRT + mwaves-m+1
    enddo
 
 ENDDO

!==>  FOURIER TRANSFORM TO COMPUTE THE VALUES OF THE VARIABLE IN GRID 
!     SPACE at THE J-TH LATITUDE. 

   call rfft(sphere_dat, ugrid, am, -1)
 
 else
   print *, 'error in spharm: idir must be -1 or +1!'
   print *, 'execution terminated in subroutine spharm'
 end if

 end subroutine spharm

 subroutine getuv(sphere_dat,vrtnm,divnm,ug,vg)

! compute U,V (winds times cos(lat)) from vrtnm,divnm
! (spectral coeffs of vorticity and divergence).

 type (sphere), intent(in) :: sphere_dat
 real, dimension(sphere_dat%nlons,sphere_dat%nlats), intent(out) ::  ug,vg
 complex, dimension((sphere_dat%ntrunc+1)*(sphere_dat%ntrunc+2)/2), intent(in) :: vrtnm,divnm
 complex, dimension(sphere_dat%ntrunc+1,sphere_dat%nlats)  :: um,vm
 integer :: nlats,ntrunc,mwaves,m,j,n,nm,nmstrt
 real ::  rm

 if (.not. sphere_dat%isinitialized) then
    print *, 'uninitialized sphere object in getuv!'
    stop
 end if


 nlats = sphere_dat%nlats
 ntrunc = sphere_dat%ntrunc
 mwaves = ntrunc+1

 do j=1,nlats
    nmstrt = 0
    do m=1,mwaves
       rm = m-1
       um(m,j) = cmplx(0.,0.)
       vm(m,j) = cmplx(0.,0.)
       do n=1,mwaves-m+1
          nm = nmstrt + n
          um(m,j) = um(m,j) + (sphere_dat%ilap(nm)/sphere_dat%rsphere)*( &
                    cmplx(0.,rm)*divnm(nm)*sphere_dat%pnm(nm,j) + &
                    vrtnm(nm)*sphere_dat%hnm(nm,j) )
          vm(m,j) = vm(m,j) + (sphere_dat%ilap(nm)/sphere_dat%rsphere)*( &
                    cmplx(0.,rm)*vrtnm(nm)*sphere_dat%pnm(nm,j) - &
                    divnm(nm)*sphere_dat%hnm(nm,j) )
       enddo
       nmstrt = nmstrt + mwaves-m+1
    enddo
 enddo

 call rfft(sphere_dat, ug, um, -1)
 call rfft(sphere_dat, vg, vm, -1)
 
 end subroutine getuv

 subroutine getvrtdiv(sphere_dat,vrtnm,divnm,ug,vg)

! compute vrtnm,divnm (spectral coeffs of vorticity and
! divergence) from U,V (winds time cos(lat)).

 type (sphere), intent(in) :: sphere_dat
 
 real, dimension(sphere_dat%nlons,sphere_dat%nlats), intent(inout) ::  ug,vg
 complex, dimension((sphere_dat%ntrunc+1)*(sphere_dat%ntrunc+2)/2), intent(in) :: vrtnm,divnm
 complex, dimension(sphere_dat%ntrunc+1,sphere_dat%nlats) :: um,vm

 if (.not. sphere_dat%isinitialized) then
    print *, 'uninitialized sphere object in getvrtdiv!'
    stop
 end if

 call rfft(sphere_dat, ug, um, 1)
 call rfft(sphere_dat, vg, vm, 1)

 call sumnm(sphere_dat,um,vm,divnm,1,1)
 call sumnm(sphere_dat,vm,um,vrtnm,1,-1)
 
 end subroutine getvrtdiv

 subroutine sumnm(sphere_dat,am,bm,anm,isign1,isign2)
!
!  given the arrays of fourier coeffs, am and bm,
!  compute the complex spherical harmonic coeffs of:  
!
!  isign1*( (1./re*(1.-x**2))*d(ag)/d(lon) + (isign2/re)*d(bg)/dx )
!
!  where x = sin(lat), isign1 and isign2 are either +1 or -1,
!  ag and bg are the physical space values of am,bm and re
!  is the radius of the sphere.
!
!  the result is returned in anm.
!
!  for example on how to use this routine, see subroutine getvrtdiv.
!
!
 type (sphere), intent(in) :: sphere_dat
 
 integer, intent(in) :: isign1,isign2
 complex, dimension((sphere_dat%ntrunc+1)*(sphere_dat%ntrunc+2)/2) :: anm
 complex, dimension(sphere_dat%ntrunc+1,sphere_dat%nlats), intent(in) :: am,bm
 integer :: nlats,ntrunc,mwaves,j,m,n,nm,nmstrt
 real ::  sign1,sign2,rm

 if (.not. sphere_dat%isinitialized) then
    print *, 'uninitialized sphere object in sumnm!'
    stop
 end if


 sign1 = float(isign1)
 sign2 = float(isign2)
 if (isign2 .ne. 1 .and. isign2 .ne. -1) then
    print *, ' isign2 must be +1 or -1 in sumnm!'
    print *, ' execution terminated in sumnm'
 end if
 if (isign1 .ne. 1 .and. isign1 .ne. -1) then
    print *, ' isign1 must be +1 or -1 in sumnm!'
    print *, ' execution terminated in sumnm'
    stop
 end if
 nlats = sphere_dat%nlats
 ntrunc = sphere_dat%ntrunc
 mwaves = ntrunc+1
 
 anm = 0.
 DO J=1,NLATS
    NMSTRT = 0 
    DO m = 1, MWAVES 
       RM = M-1 
       DO n   = 1, mwaves-m+1
          NM = NMSTRT + n
          aNM(NM) = aNM(NM) + sign1*sphere_dat%GWrC(j)*(CMPLX(0.,RM) &
                    * sphere_dat%PNM(NM,J) * am(m,j) &
                    + sign2 * sphere_dat%HNM(NM,J) * bm(m,j))
       enddo
       nmstrt = nmstrt + mwaves - m +1
    enddo
 enddo

 end subroutine sumnm

 subroutine cosgrad(sphere_dat,divnm,ug,vg)

! compute coslat * gradient of spectral coefficients (divnm)
! vector gradient returned on grid as (ug,vg)

 type (sphere), intent(in) :: sphere_dat
 
 real, dimension(sphere_dat%nlons,sphere_dat%nlats), intent(out) ::  ug,vg
 complex, dimension((sphere_dat%ntrunc+1)*(sphere_dat%ntrunc+2)/2), intent(in) :: divnm
 complex, dimension(sphere_dat%ntrunc+1,sphere_dat%nlats) :: um,vm
 integer :: nlats,ntrunc,mwaves,j,m,n,nm,nmstrt
 real :: rm

 if (.not. sphere_dat%isinitialized) then
    print *, 'uninitialized sphere object in cosgrad!'
    stop
 end if


 nlats = sphere_dat%nlats
 ntrunc = sphere_dat%ntrunc
 mwaves = ntrunc+1
 
 do j=1,nlats
    nmstrt = 0
    do m=1,mwaves
       rm = (m-1)
       um(m,j) = cmplx(0.,0.)
       vm(m,j) = cmplx(0.,0.)
       do n=1,mwaves-m+1
          nm = nmstrt + n
          um(m,j) = um(m,j) + (1./sphere_dat%rsphere)* &
                    cmplx(0.,rm)*divnm(nm)*sphere_dat%pnm(nm,j) 
          vm(m,j) = vm(m,j) - (1./sphere_dat%rsphere)* &
                    divnm(nm)*sphere_dat%hnm(nm,j) 
       enddo
       nmstrt = nmstrt + mwaves - m +1
    enddo
 enddo

 call rfft(sphere_dat, ug, um, -1)
 call rfft(sphere_dat, vg, vm, -1)

 end subroutine cosgrad

 subroutine specsmooth(sphere_dat,datagrid,smooth)
 
! isotropic spectral smoothing of datagrid.
! input: smooth(sphere_dat%ntrunc+1) - smoothing factor as a
! function of degree (sphere_dat%indxn).

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