c c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c . . c . copyright (c) 1998 by UCAR . c . . c . University Corporation for Atmospheric Research . c . . c . all rights reserved . c . . c . . c . SPHEREPACK . c . . c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c c c c ... file shaes.f c c this file contains code and documentation for subroutines c shaes and shaesi c c ... files which must be loaded with shaes.f c c sphcom.f, hrfft.f c c subroutine shaes(nlat,nlon,isym,nt,g,idg,jdg,a,b,mdab,ndab, c + wshaes,lshaes,work,lwork,ierror) c c subroutine shaes performs the spherical harmonic analysis c on the array g and stores the result in the arrays a and b. c the analysis is performed on an equally spaced grid. the c associated legendre functions are stored rather than recomputed c as they are in subroutine shaec. the analysis is described c below at output parameters a,b. c c sphcom.f, hrfft.f c c c input parameters c c nlat the number of colatitudes on the full sphere including the c poles. for example, nlat = 37 for a five degree grid. c nlat determines the grid increment in colatitude as c pi/(nlat-1). if nlat is odd the equator is located at c grid point i=(nlat+1)/2. if nlat is even the equator is c located half way between points i=nlat/2 and i=nlat/2+1. c nlat must be at least 3. note: on the half sphere, the c number of grid points in the colatitudinal direction is c nlat/2 if nlat is even or (nlat+1)/2 if nlat is odd. c c nlon the number of distinct londitude points. nlon determines c the grid increment in longitude as 2*pi/nlon. for example c nlon = 72 for a five degree grid. nlon must be greater c than or equal to 4. the efficiency of the computation is c improved when nlon is a product of small prime numbers. c c isym = 0 no symmetries exist about the equator. the analysis c is performed on the entire sphere. i.e. on the c array g(i,j) for i=1,...,nlat and j=1,...,nlon. c (see description of g below) c c = 1 g is antisymmetric about the equator. the analysis c is performed on the northern hemisphere only. i.e. c if nlat is odd the analysis is performed on the c array g(i,j) for i=1,...,(nlat+1)/2 and j=1,...,nlon. c if nlat is even the analysis is performed on the c array g(i,j) for i=1,...,nlat/2 and j=1,...,nlon. c c c = 2 g is symmetric about the equator. the analysis is c performed on the northern hemisphere only. i.e. c if nlat is odd the analysis is performed on the c array g(i,j) for i=1,...,(nlat+1)/2 and j=1,...,nlon. c if nlat is even the analysis is performed on the c array g(i,j) for i=1,...,nlat/2 and j=1,...,nlon. c c nt the number of analyses. in the program that calls shaes, c the arrays g,a and b can be three dimensional in which c case multiple analyses will be performed. the third c index is the analysis index which assumes the values c k=1,...,nt. for a single analysis set nt=1. the c discription of the remaining parameters is simplified c by assuming that nt=1 or that the arrays g,a and b c have only two dimensions. c c g a two or three dimensional array (see input parameter c nt) that contains the discrete function to be analyzed. c g(i,j) contains the value of the function at the colatitude c point theta(i) = (i-1)*pi/(nlat-1) and longitude point c phi(j) = (j-1)*2*pi/nlon. the index ranges are defined c above at the input parameter isym. c c c idg the first dimension of the array g as it appears in the c program that calls shaes. if isym equals zero then idg c must be at least nlat. if isym is nonzero then idg c must be at least nlat/2 if nlat is even or at least c (nlat+1)/2 if nlat is odd. c c jdg the second dimension of the array g as it appears in the c program that calls shaes. jdg must be at least nlon. c c mdab the first dimension of the arrays a and b as it appears c in the program that calls shaes. mdab must be at least c min0(nlat,(nlon+2)/2) if nlon is even or at least c min0(nlat,(nlon+1)/2) if nlon is odd. c c ndab the second dimension of the arrays a and b as it appears c in the program that calls shaes. ndab must be at least nlat c c wshaes an array which must be initialized by subroutine shaesi. c once initialized, wshaes can be used repeatedly by shaes c as long as nlon and nlat remain unchanged. wshaes must c not be altered between calls of shaes. c c lshaes the dimension of the array wshaes as it appears in the c program that calls shaes. define c c l1 = min0(nlat,(nlon+2)/2) if nlon is even or c l1 = min0(nlat,(nlon+1)/2) if nlon is odd c c and c c l2 = nlat/2 if nlat is even or c l2 = (nlat+1)/2 if nlat is odd c c then lshaes must be at least c c (l1*l2*(nlat+nlat-l1+1))/2+nlon+15 c c work a work array that does not have to be saved. c c lwork the dimension of the array work as it appears in the c program that calls shaes. define c c l2 = nlat/2 if nlat is even or c l2 = (nlat+1)/2 if nlat is odd c c if isym is zero then lwork must be at least c (nt+1)*nlat*nlon. if isym is not zero then c lwork must be at least (nt+1)*l2*nlon. c c c ************************************************************** c c output parameters c c a,b both a,b are two or three dimensional arrays (see input c parameter nt) that contain the spherical harmonic c coefficients in the representation of g(i,j) given in the c discription of subroutine shses. for isym=0, a(m,n) and c b(m,n) are given by the equations listed below. symmetric c versions are used when isym is greater than zero. c c c c definitions c c 1. the normalized associated legendre functions c c pbar(m,n,theta) = sqrt((2*n+1)*factorial(n-m)/(2*factorial(n+m))) c *sin(theta)**m/(2**n*factorial(n)) times the c (n+m)th derivative of (x**2-1)**n with respect c to x=cos(theta) c c 2. the normalized z functions for m even c c zbar(m,n,theta) = 2/(nlat-1) times the sum from k=0 to k=nlat-1 of c the integral from tau = 0 to tau = pi of c cos(k*theta)*cos(k*tau)*pbar(m,n,tau)*sin(tau) c (first and last terms in this sum are divided c by 2) c c 3. the normalized z functions for m odd c c zbar(m,n,theta) = 2/(nlat-1) times the sum from k=0 to k=nlat-1 of c of the integral from tau = 0 to tau = pi of c sin(k*theta)*sin(k*tau)*pbar(m,n,tau)*sin(tau) c c 4. the fourier transform of g(i,j). c c c(m,i) = 2/nlon times the sum from j=1 to j=nlon c of g(i,j)*cos((m-1)*(j-1)*2*pi/nlon) c (the first and last terms in this sum c are divided by 2) c c s(m,i) = 2/nlon times the sum from j=2 to j=nlon c of g(i,j)*sin((m-1)*(j-1)*2*pi/nlon) c c 5. the maximum (plus one) longitudinal wave number c c mmax = min0(nlat,(nlon+2)/2) if nlon is even or c mmax = min0(nlat,(nlon+1)/2) if nlon is odd. c c then for m=0,...,mmax-1 and n=m,...,nlat-1 the arrays a,b are c given by c c a(m+1,n+1) = the sum from i=1 to i=nlat of c c(m+1,i)*zbar(m,n,theta(i)) c (first and last terms in this sum are c divided by 2) c c b(m+1,n+1) = the sum from i=1 to i=nlat of c s(m+1,i)*zbar(m,n,theta(i)) c c c ierror = 0 no errors c = 1 error in the specification of nlat c = 2 error in the specification of nlon c = 3 error in the specification of isym c = 4 error in the specification of nt c = 5 error in the specification of idg c = 6 error in the specification of jdg c = 7 error in the specification of mdab c = 8 error in the specification of ndab c = 9 error in the specification of lshaes c = 10 error in the specification of lwork c c c **************************************************************** c subroutine shaesi(nlat,nlon,wshaes,lshaes,work,lwork,dwork, c + ldwork,ierror) c c subroutine shaesi initializes the array wshaes which can then c be used repeatedly by subroutine shaes c c input parameters c c nlat the number of colatitudes on the full sphere including the c poles. for example, nlat = 37 for a five degree grid. c nlat determines the grid increment in colatitude as c pi/(nlat-1). if nlat is odd the equator is located at c grid point i=(nlat+1)/2. if nlat is even the equator is c located half way between points i=nlat/2 and i=nlat/2+1. c nlat must be at least 3. note: on the half sphere, the c number of grid points in the colatitudinal direction is c nlat/2 if nlat is even or (nlat+1)/2 if nlat is odd. c c nlon the number of distinct londitude points. nlon determines c the grid increment in longitude as 2*pi/nlon. for example c nlon = 72 for a five degree grid. nlon must be greater c than or equal to 4. the efficiency of the computation is c improved when nlon is a product of small prime numbers. c c lshaes the dimension of the array wshaes as it appears in the c program that calls shaesi. define c c l1 = min0(nlat,(nlon+2)/2) if nlon is even or c l1 = min0(nlat,(nlon+1)/2) if nlon is odd c c and c c l2 = nlat/2 if nlat is even or c l2 = (nlat+1)/2 if nlat is odd c c then lshaes must be at least c c (l1*l2*(nlat+nlat-l1+1))/2+nlon+15 c c work a real work array that does not have to be saved. c c lwork the dimension of the array work as it appears in the c program that calls shaesi. define c c l1 = min0(nlat,(nlon+2)/2) if nlon is even or c l1 = min0(nlat,(nlon+1)/2) if nlon is odd c c and c c l2 = nlat/2 if nlat is even or c l2 = (nlat+1)/2 if nlat is odd c c then lwork must be at least c c 5*nlat*l2+3*((l1-2)*(nlat+nlat-l1-1))/2 c c c dwork a double precision work array that does not have to be saved. c c ldwork the dimension of the array dwork as it appears in the c program that calls shaesi. ldwork must be at least nlat+1 c c c output parameters c c wshaes an array which is initialized for use by subroutine shaes. c once initialized, wshaes can be used repeatedly by shaes c as long as nlon and nlat remain unchanged. wshaes must c not be altered between calls of shaes. c c ierror = 0 no errors c = 1 error in the specification of nlat c = 2 error in the specification of nlon c = 3 error in the specification of lshaes c = 4 error in the specification of lwork c = 5 error in the specification of ldwork c c c ****************************************************************