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 shsgs.f c c this file contains code and documentation for subroutines c shsgs and shsgsi c c ... files which must be loaded with shsgs.f c c sphcom.f, hrfft.f, gaqd.f c c subroutine shsgs(nlat,nlon,isym,nt,g,idg,jdg,a,b,mdab,ndab, c 1 wshsgs,lshsgs,work,lwork,ierror) c c subroutine shsgs performs the spherical harmonic synthesis c on the arrays a and b and stores the result in the array g. c the synthesis is performed on an equally spaced longitude grid c and a gaussian colatitude grid. the associated legendre functions c are stored rather than recomputed as they are in subroutine c shsgc. the synthesis is described below at output parameter c g. c 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 synthesis 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 synthesis c is performed on the northern hemisphere only. i.e. c if nlat is odd the synthesis 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 synthesis 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 synthesis is c performed on the northern hemisphere only. i.e. c if nlat is odd the synthesis 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 synthesis is performed on the c array g(i,j) for i=1,...,nlat/2 and j=1,...,nlon. c c nt the number of syntheses. in the program that calls shsgs, c the arrays g,a and b can be three dimensional in which c case multiple synthesis will be performed. the third c index is the synthesis index which assumes the values c k=1,...,nt. for a single synthesis 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 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 a,b two or three dimensional arrays (see the input parameter c nt) that contain the coefficients in the spherical harmonic c expansion of g(i,j) given below at the definition of the c output parameter g. a(m,n) and b(m,n) are defined for c indices m=1,...,mmax and n=m,...,nlat where mmax is the c maximum (plus one) longitudinal wave number given by 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 mdab the first dimension of the arrays a and b as it appears c in the program that calls shsgs. 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 shsgs. ndab must be at least nlat c c wshsgs an array which must be initialized by subroutine shsgsi. c once initialized, wshsgs can be used repeatedly by shsgs c as long as nlat and nlon remain unchanged. wshsgs must c not be altered between calls of shsgs. c c lshsgs the dimension of the array wshsgs as it appears in the c program that calls shsgs. 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 lshsgs must be at least c c nlat*(3*(l1+l2)-2)+(l1-1)*(l2*(2*nlat-l1)-3*l1)/2+nlon+15 c c c lwork the dimension of the array work as it appears in the c program that calls shsgs. define c c l2 = nlat/2 if nlat is even or c l2 = (nlat+1)/2 if nlat is odd c c c if isym is zero then lwork must be at least c c nlat*nlon*(nt+1) c c if isym is nonzero then lwork must be at least c c l2*nlon*(nt+1) c c c ************************************************************** c c output parameters c c g a two or three dimensional array (see input parameter nt) c that contains the discrete function which is synthesized. c g(i,j) contains the value of the function at the gaussian c colatitude point theta(i) and longitude point c phi(j) = (j-1)*2*pi/nlon. the index ranges are defined c above at the input parameter isym. for isym=0, g(i,j) c is given by the the equations listed below. symmetric c versions are used when isym is greater than zero. c c the normalized associated legendre functions are given by 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 define the maximum (plus one) longitudinal wave number c as 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 g(i,j) = the sum from n=0 to n=nlat-1 of c c .5*pbar(0,n,theta(i))*a(1,n+1) c c plus the sum from m=1 to m=mmax-1 of c c the sum from n=m to n=nlat-1 of c c pbar(m,n,theta(i))*(a(m+1,n+1)*cos(m*phi(j)) c -b(m+1,n+1)*sin(m*phi(j))) 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 lshsgs c = 10 error in the specification of lwork c c c **************************************************************** c c subroutine shsgsi(nlat,nlon,wshsgs,lshsgs,work,lwork,dwork,ldwork, c + ierror) c c subroutine shsgsi initializes the array wshsgs which can then c be used repeatedly by subroutines shsgs. it precomputes c and stores in wshsgs 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 wshsgs an array which must be initialized by subroutine shsgsi. c once initialized, wshsgs can be used repeatedly by shsgs c as long as nlat and nlon remain unchanged. wshsgs must c not be altered between calls of shsgs. c c lshsgs the dimension of the array wshsgs as it appears in the c program that calls shsgs. 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 lshsgs must be at least c c nlat*(3*(l1+l2)-2)+(l1-1)*(l2*(2*nlat-l1)-3*l1)/2+nlon+15 c c work a real work space which need not be saved c c lwork the dimension of the array work as it appears in the c program that calls shsgsi. lwork must be at least c 4*nlat*(nlat+2)+2 in the routine calling shsgsi c c dwork a double precision work array that does not have to be saved. c c ldwork the length of dwork in the calling routine. ldwork must c be at least nlat*(nlat+4) c c output parameter c c wshsgs an array which must be initialized before calling shsgs or c once initialized, wshsgs can be used repeatedly by shsgs or c as long as nlat and nlon remain unchanged. wshsgs must not c altered between calls of shsgs. 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 lshsgs c = 4 error in the specification of lwork c = 5 error in the specification of ldwork c = 5 failure in gaqd to compute gaussian points c (due to failure in eigenvalue routine) c c