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 vhags.f c c this file contains code and documentation for subroutines c vhags and vhagsi c c ... files which must be loaded with vhags.f c c sphcom.f, hrfft.f, gaqd.f c c subroutine vhags(nlat,nlon,ityp,nt,v,w,idvw,jdvw,br,bi,cr,ci, c + mdab,ndab,wvhags,lvhags,work,lwork,ierror) c c subroutine vhags performs the vector spherical harmonic analysis c on the vector field (v,w) and stores the result in the arrays c br, bi, cr, and ci. v(i,j) and w(i,j) are the colatitudinal c (measured from the north pole) and east longitudinal components c respectively, located at the gaussian colatitude point theta(i) c and longitude phi(j) = (j-1)*2*pi/nlon. the spectral c representation of (v,w) is given at output parameters v,w in c subroutine vhses. 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 computed c in radians in theta(1) <...< theta(nlat) by subroutine gaqd. c if nlat is odd the equator will be included as the grid point 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 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 zero. the axisymmetric case corresponds to nlon=1. c the efficiency of the computation is improved when nlon c is a product of small prime numbers. c c ityp = 0 no symmetries exist about the equator. the analysis c is performed on the entire sphere. i.e. on the c arrays v(i,j),w(i,j) for i=1,...,nlat and c j=1,...,nlon. c c = 1 no symmetries exist about the equator. the analysis c is performed on the entire sphere. i.e. on the c arrays v(i,j),w(i,j) for i=1,...,nlat and c j=1,...,nlon. the curl of (v,w) is zero. that is, c (d/dtheta (sin(theta) w) - dv/dphi)/sin(theta) = 0. c the coefficients cr and ci are zero. c c = 2 no symmetries exist about the equator. the analysis c is performed on the entire sphere. i.e. on the c arrays v(i,j),w(i,j) for i=1,...,nlat and c j=1,...,nlon. the divergence of (v,w) is zero. i.e., c (d/dtheta (sin(theta) v) + dw/dphi)/sin(theta) = 0. c the coefficients br and bi are zero. c c = 3 v is symmetric and w is antisymmetric about the c equator. the analysis is performed on the northern c hemisphere only. i.e., if nlat is odd the analysis c is performed on the arrays v(i,j),w(i,j) for c i=1,...,(nlat+1)/2 and j=1,...,nlon. if nlat is c even the analysis is performed on the the arrays c v(i,j),w(i,j) for i=1,...,nlat/2 and j=1,...,nlon. c c = 4 v is symmetric and w is antisymmetric about the c equator. the analysis is performed on the northern c hemisphere only. i.e., if nlat is odd the analysis c is performed on the arrays v(i,j),w(i,j) for c i=1,...,(nlat+1)/2 and j=1,...,nlon. if nlat is c even the analysis is performed on the the arrays c v(i,j),w(i,j) for i=1,...,nlat/2 and j=1,...,nlon. c the curl of (v,w) is zero. that is, c (d/dtheta (sin(theta) w) - dv/dphi)/sin(theta) = 0. c the coefficients cr and ci are zero. c c = 5 v is symmetric and w is antisymmetric about the c equator. the analysis is performed on the northern c hemisphere only. i.e., if nlat is odd the analysis c is performed on the arrays v(i,j),w(i,j) for c i=1,...,(nlat+1)/2 and j=1,...,nlon. if nlat is c even the analysis is performed on the the arrays c v(i,j),w(i,j) for i=1,...,nlat/2 and j=1,...,nlon. c the divergence of (v,w) is zero. i.e., c (d/dtheta (sin(theta) v) + dw/dphi)/sin(theta) = 0. c the coefficients br and bi are zero. c c = 6 v is antisymmetric and w is symmetric about the c equator. the analysis is performed on the northern c hemisphere only. i.e., if nlat is odd the analysis c is performed on the arrays v(i,j),w(i,j) for c i=1,...,(nlat+1)/2 and j=1,...,nlon. if nlat is c even the analysis is performed on the the arrays c v(i,j),w(i,j) for i=1,...,nlat/2 and j=1,...,nlon. c c = 7 v is antisymmetric and w is symmetric about the c equator. the analysis is performed on the northern c hemisphere only. i.e., if nlat is odd the analysis c is performed on the arrays v(i,j),w(i,j) for c i=1,...,(nlat+1)/2 and j=1,...,nlon. if nlat is c even the analysis is performed on the the arrays c v(i,j),w(i,j) for i=1,...,nlat/2 and j=1,...,nlon. c the curl of (v,w) is zero. that is, c (d/dtheta (sin(theta) w) - dv/dphi)/sin(theta) = 0. c the coefficients cr and ci are zero. c c = 8 v is antisymmetric and w is symmetric about the c equator. the analysis is performed on the northern c hemisphere only. i.e., if nlat is odd the analysis c is performed on the arrays v(i,j),w(i,j) for c i=1,...,(nlat+1)/2 and j=1,...,nlon. if nlat is c even the analysis is performed on the the arrays c v(i,j),w(i,j) for i=1,...,nlat/2 and j=1,...,nlon. c the divergence of (v,w) is zero. i.e., c (d/dtheta (sin(theta) v) + dw/dphi)/sin(theta) = 0. c the coefficients br and bi are zero. c c c nt the number of analyses. in the program that calls vhags, c the arrays v,w,br,bi,cr, and ci can be three dimensional c in which case multiple analyses will be performed. c the third index is the analysis index which assumes the c values 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 all the arrays are two c dimensional. c c v,w two or three dimensional arrays (see input parameter nt) c that contain the vector function to be analyzed. c v is the colatitudnal component and w is the east c longitudinal component. v(i,j),w(i,j) contain the c components at the gaussian colatitude point theta(i) c and longitude phi(j) = (j-1)*2*pi/nlon. the index ranges c are defined above at the input parameter ityp. c c idvw the first dimension of the arrays v,w as it appears in c the program that calls vhags. if ityp .le. 2 then idvw c must be at least nlat. if ityp .gt. 2 and nlat is c even then idvw must be at least nlat/2. if ityp .gt. 2 c and nlat is odd then idvw must be at least (nlat+1)/2. c c jdvw the second dimension of the arrays v,w as it appears in c the program that calls vhags. jdvw must be at least nlon. c c mdab the first dimension of the arrays br,bi,cr, and ci as it c appears in the program that calls vhags. mdab must be at c least min0(nlat,nlon/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 br,bi,cr, and ci as it c appears in the program that calls vhags. ndab must be at c least nlat. c c wvhags an array which must be initialized by subroutine vhgsi. c once initialized, wvhags can be used repeatedly by vhags c as long as nlon and nlat remain unchanged. wvhags must c not be altered between calls of vhags. c c lvhags the dimension of the array wvhags as it appears in the c program that calls vhags. define c c l1 = min0(nlat,nlon/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 lvhags must be at least c c (nlat+1)*(nlat+1)*nlat/2 + nlon + 15 c 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 vhags. define c c l2 = nlat/2 if nlat is even or c l2 = (nlat+1)/2 if nlat is odd c c if ityp .le. 2 then lwork must be at least c the larger of the two quantities c c 3*nlat*(nlat+1)+2 (required by vhagsi) c c and c c (2*nt+1)*nlat*nlon c c if ityp .gt. 2 then lwork must be at least c the larger of the two quantities c c 3*nlat*(nlat+1)+2 (required by vhagsi) c c and c c (2*nt+1)*l2*nlon c c c ************************************************************** c c output parameters c c br,bi two or three dimensional arrays (see input parameter nt) c cr,ci that contain the vector spherical harmonic coefficients c in the spectral representation of v(i,j) and w(i,j) given c in the discription of subroutine vhses. br(mp1,np1), c bi(mp1,np1),cr(mp1,np1), and ci(mp1,np1) are computed c for mp1=1,...,mmax and np1=mp1,...,nlat except for np1=nlat c and odd mp1. mmax=min0(nlat,nlon/2) if nlon is even or c mmax=min0(nlat,(nlon+1)/2) if nlon is odd. 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 ityp c = 4 error in the specification of nt c = 5 error in the specification of idvw c = 6 error in the specification of jdvw c = 7 error in the specification of mdab c = 8 error in the specification of ndab c = 9 error in the specification of lvhags c = 10 error in the specification of lwork c c c subroutine vhagsi(nlat,nlon,wvhags,lvhags,work,lwork,ierror) c c subroutine vhagsi initializes the array wvhags which can then be c used repeatedly by subroutine vhags until nlat or nlon is changed. 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 computed c in radians in theta(1) <...< theta(nlat) by subroutine gaqd. c if nlat is odd the equator will be included as the grid point 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 zero. the axisymmetric case corresponds to nlon=1. c the efficiency of the computation is improved when nlon c is a product of small prime numbers. c c lvhags the dimension of the array wvhags as it appears in the c program that calls vhagsi. define c c l1 = min0(nlat,nlon/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 lvhags must be at least c c 3*nlat*(nlat+1)+2 (required by vhagsi) c c dwork a double precision work space that does not need to be saved c c ldwork the dimension of the array dwork as it appears in the c program that calls vhagsi. ldwork must be at least c c (3*nlat*(nlat+3)+2)/2 c c ************************************************************** c c output parameters c c wvhags an array which is initialized for use by subroutine vhags. c once initialized, wvhags can be used repeatedly by vhags c as long as nlat and nlon remain unchanged. wvhags must not c be altered between calls of vhags. 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 lvhags c = 4 error in the specification of ldwork c