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 ... file shsgc.f c c this file contains code and documentation for subroutines c shsgc and shsgci c c ... files which must be loaded with shsgc.f c c sphcom.f, hrfft.f, gaqd.f c c subroutine shsgc(nlat,nlon,isym,nt,g,idg,jdg,a,b,mdab,ndab, c + wshsgc,lshsgc,work,lwork,ierror) c c subroutine shsgc 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 recomputed rather than stored as they are in subroutine c shsgs. the synthesis is described below at output parameter c g. 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 shsgc, 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 shsgc. 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 shsgc. 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 shsgc. 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 shsgc. ndab must be at least nlat 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 wshsgc an array which must be initialized by subroutine shsgci. c once initialized, wshsgc can be used repeatedly by shsgc c as long as nlat and nlon remain unchanged. wshsgc must c not be altered between calls of shsgc. c c lshsgc the dimension of the array wshsgc as it appears in the c program that calls shsgc. 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 lshsgc must be at least c c nlat*(2*l2+3*l1-2)+3*l1*(1-l1)/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 shsgc. 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 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 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 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 lwshig c = 10 error in the specification of lwork c c c **************************************************************** c c subroutine shsgci(nlat,nlon,wshsgc,lshsgc,dwork,ldwork,ierror) c c subroutine shsgci initializes the array wshsgc which can then c be used repeatedly by subroutines shsgc. it precomputes c and stores in wshsgc 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 wshsgc an array which must be initialized by subroutine shsgci. c once initialized, wshsgc can be used repeatedly by shsgc c as long as nlat and nlon remain unchanged. wshsgc must c not be altered between calls of shsgc. c c lshsgc the dimension of the array wshsgc as it appears in the c program that calls shsgc. 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 lshsgc 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 shsgci. ldwork must be at least c c nlat*(nlat+4) c c output parameter c c wshsgc an array which must be initialized before calling shsgc. c once initialized, wshsgc can be used repeatedly by shsgc c as long as nlat and nlon remain unchanged. wshsgc must not c altered between calls of shsgc. 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 lshsgc 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 ****************************************************************