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 shagc.f c c this file contains code and documentation for subroutines c shagc and shagci c c ... files which must be loaded with shagc.f c c sphcom.f, hrfft.f, gaqd.f c c c subroutine shagc(nlat,nlon,isym,nt,g,idg,jdg,a,b,mdab,ndab, c + wshagc,lshagc,work,lwork,ierror) c c subroutine shagc 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 a gaussian grid in colatitude c and an equally spaced grid in longitude. the associated c legendre functions are recomputed rather than stored as they c are in subroutine shags. the analysis is described below c at output parameters a,b. c c input parameters c c nlat the number of points in the gaussian colatitude grid on the c full sphere. these lie in the interval (0,pi) and are compu c in radians in theta(1),...,theta(nlat) by subroutine gaqd. c if nlat is odd the equator will be included as the grid poi c theta((nlat+1)/2). if nlat is even the equator will be c excluded as a grid point and will lie half way between c theta(nlat/2) and theta(nlat/2+1). nlat must be at least 3. c note: on the half sphere, the number of grid points in the c colatitudinal direction is nlat/2 if nlat is even or c (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 shagc, 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 gaussian c point theta(i) and longitude point phi(j) = (j-1)*2*pi/nlon c the index ranges are defined above at the input parameter c isym. c c idg the first dimension of the array g as it appears in the c program that calls shagc. if isym equals zero then idg c must be at least nlat. if isym is nonzero then idg must c be at least nlat/2 if nlat is even or at least (nlat+1)/2 c if nlat is odd. c c jdg the second dimension of the array g as it appears in the c program that calls shagc. 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 shagc. mdab must be at least c min0((nlon+2)/2,nlat) if nlon is even or at least c min0((nlon+1)/2,nlat) 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 shaec. ndab must be at least nlat c c wshagc an array which must be initialized by subroutine shagci. c once initialized, wshagc can be used repeatedly by shagc. c as long as nlat and nlon remain unchanged. wshagc must c not be altered between calls of shagc. c c lshagc the dimension of the array wshagc as it appears in the c program that calls shagc. 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 lshagc must be at least c c nlat*(2*l2+3*l1-2)+3*l1*(1-l1)/2+nlon+15 c 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 shagc. 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 if isym is zero then lwork must be at least c c nlat*(nlon*nt+max0(3*l2,nlon)) c c if isym is not zero then lwork must be at least c c l2*(nlon*nt+max0(3*nlat,nlon)) 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 shagc. 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 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 fourier transform of g(i,j). c c c(m,i) = 2/nlon times the sum from j=1 to j=nlon of c 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 of c g(i,j)*sin((m-1)*(j-1)*2*pi/nlon) c c c 3. the gaussian points and weights on the sphere c (computed by subroutine gaqd). c c theta(1),...,theta(nlat) (gaussian pts in radians) c wts(1),...,wts(nlat) (corresponding gaussian weights) c c 4. 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 c then for m=0,...,mmax-1 and n=m,...,nlat-1 the arrays a,b c are given by c c a(m+1,n+1) = the sum from i=1 to i=nlat of c c(m+1,i)*wts(i)*pbar(m,n,theta(i)) c c b(m+1,n+1) = the sum from i=1 to nlat of c s(m+1,i)*wts(i)*pbar(m,n,theta(i)) 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 lshagc c = 10 error in the specification of lwork c c c **************************************************************** c c subroutine shagci(nlat,nlon,wshagc,lshagc,dwork,ldwork,ierror) c c subroutine shagci initializes the array wshagc which can then c be used repeatedly by subroutines shagc. it precomputes c and stores in wshagc quantities such as gaussian weights, c legendre polynomial coefficients, and fft trigonometric tables. c c input parameters c c nlat the number of points in the gaussian colatitude grid on the c full sphere. these lie in the interval (0,pi) and are compu c in radians in theta(1),...,theta(nlat) by subroutine gaqd. c if nlat is odd the equator will be included as the grid poi c theta((nlat+1)/2). if nlat is even the equator will be c excluded as a grid point and will lie half way between c theta(nlat/2) and theta(nlat/2+1). nlat must be at least 3. c note: on the half sphere, the number of grid points in the c colatitudinal direction is nlat/2 if nlat is even or c (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 wshagc an array which must be initialized by subroutine shagci. c once initialized, wshagc can be used repeatedly by shagc c as long as nlat and nlon remain unchanged. wshagc must c not be altered between calls of shagc. c c lshagc the dimension of the array wshagc as it appears in the c program that calls shagc. 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 lshagc must be at least c c nlat*(2*l2+3*l1-2)+3*l1*(1-l1)/2+nlon+15 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 shagci. ldwork must be at least c c nlat*(nlat+4) c c output parameter c c wshagc an array which must be initialized before calling shagc or c once initialized, wshagc can be used repeatedly by shagc or c as long as nlat and nlon remain unchanged. wshagc must not c altered between calls of shagc. 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 lshagc c = 4 error in the specification of ldwork c = 5 failure in gaqd to compute gaussian points c (due to failure in eigenvalue routine) c c c ****************************************************************