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 vtsec.f c c this file includes documentation and code for c subroutines vtsec and vtseci c c ... files which must be loaded with vtsec.f c c sphcom.f, hrfft.f, vhaec.f, vhsec.f c c subroutine vtsec(nlat,nlon,ityp,nt,vt,wt,idvw,jdvw,br,bi,cr,ci, c + mdab,ndab,wvts,lwvts,work,lwork,ierror) c c given the vector harmonic analysis br,bi,cr, and ci (computed c by subroutine vhaec) of some vector function (v,w), this c subroutine computes the vector function (vt,wt) which is c the derivative of (v,w) with respect to colatitude theta. vtsec c is similar to vhsec except the vector harmonics are replaced by c their derivative with respect to colatitude with the result that c (vt,wt) is computed instead of (v,w). vt(i,j) is the derivative c of the colatitudinal component v(i,j) at the point theta(i) = c (i-1)*pi/(nlat-1) and longitude phi(j) = (j-1)*2*pi/nlon. the c spectral representation of (vt,wt) is given below at output c parameters vt,wt. 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 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 synthesis c is performed on the entire sphere. i.e. the arrays c vt(i,j),wt(i,j) are computed for i=1,...,nlat and c j=1,...,nlon. c c = 1 no symmetries exist about the equator however the c the coefficients cr and ci are zero. the synthesis c is performed on the entire sphere. i.e. the arrays c vt(i,j),wt(i,j) are computed for i=1,...,nlat and c j=1,...,nlon. c c = 2 no symmetries exist about the equator however the c the coefficients br and bi are zero. the synthesis c is performed on the entire sphere. i.e. the arrays c vt(i,j),wt(i,j) are computed for i=1,...,nlat and c j=1,...,nlon. c c = 3 vt is odd and wt is even about the equator. the c synthesis is performed on the northern hemisphere c only. i.e., if nlat is odd the arrays vt(i,j),wt(i,j) c are computed for i=1,...,(nlat+1)/2 and j=1,...,nlon. c if nlat is even the arrays vt(i,j),wt(i,j) are computed c for i=1,...,nlat/2 and j=1,...,nlon. c c = 4 vt is odd and wt is even about the equator and the c coefficients cr and ci are zero. the synthesis is c performed on the northern hemisphere only. i.e. if c nlat is odd the arrays vt(i,j),wt(i,j) are computed c for i=1,...,(nlat+1)/2 and j=1,...,nlon. if nlat is c even the arrays vt(i,j),wt(i,j) are computed for c i=1,...,nlat/2 and j=1,...,nlon. c c = 5 vt is odd and wt is even about the equator and the c coefficients br and bi are zero. the synthesis is c performed on the northern hemisphere only. i.e. if c nlat is odd the arrays vt(i,j),wt(i,j) are computed c for i=1,...,(nlat+1)/2 and j=1,...,nlon. if nlat is c even the arrays vt(i,j),wt(i,j) are computed for c i=1,...,nlat/2 and j=1,...,nlon. c c = 6 vt is even and wt is odd about the equator. the c synthesis is performed on the northern hemisphere c only. i.e., if nlat is odd the arrays vt(i,j),wt(i,j) c are computed for i=1,...,(nlat+1)/2 and j=1,...,nlon. c if nlat is even the arrays vt(i,j),wt(i,j) are computed c for i=1,...,nlat/2 and j=1,...,nlon. c c = 7 vt is even and wt is odd about the equator and the c coefficients cr and ci are zero. the synthesis is c performed on the northern hemisphere only. i.e. if c nlat is odd the arrays vt(i,j),wt(i,j) are computed c for i=1,...,(nlat+1)/2 and j=1,...,nlon. if nlat is c even the arrays vt(i,j),wt(i,j) are computed for c i=1,...,nlat/2 and j=1,...,nlon. c c = 8 vt is even and wt is odd about the equator and the c coefficients br and bi are zero. the synthesis is c performed on the northern hemisphere only. i.e. if c nlat is odd the arrays vt(i,j),wt(i,j) are computed c for i=1,...,(nlat+1)/2 and j=1,...,nlon. if nlat is c even the arrays vt(i,j),wt(i,j) are computed for c i=1,...,nlat/2 and j=1,...,nlon. c c nt the number of syntheses. in the program that calls vtsec, c the arrays vt,wt,br,bi,cr, and ci can be three dimensional c in which case multiple syntheses will be performed. c the third index is the synthesis index which assumes the c values 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 all the arrays are two c dimensional. c c idvw the first dimension of the arrays vt,wt as it appears in c the program that calls vtsec. 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 vt,wt as it appears in c the program that calls vtsec. jdvw must be at least nlon. c c br,bi two or three dimensional arrays (see input parameter nt) c cr,ci that contain the vector spherical harmonic coefficients c of (v,w) as computed by subroutine vhaec. c c mdab the first dimension of the arrays br,bi,cr, and ci as it c appears in the program that calls vtsec. 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 vtsec. ndab must be at c least nlat. c c wvts an array which must be initialized by subroutine vtseci. c once initialized, wvts can be used repeatedly by vtsec c as long as nlon and nlat remain unchanged. wvts must c not be altered between calls of vtsec. c c lwvts the dimension of the array wvts as it appears in the c program that calls vtsec. 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 lwvts must be at least c c 4*nlat*l2+3*max0(l1-2,0)*(nlat+nlat-l1-1)+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 vtsec. 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 c nlat*(2*nt*nlon+max0(6*l2,nlon)) c c if ityp .gt. 2 then lwork must be at least c c l2*(2*nt*nlon+max0(6*nlat,nlon)) c c ************************************************************** c c output parameters c c vt,wt two or three dimensional arrays (see input parameter nt) c in which the derivative of (v,w) with respect to c colatitude theta is stored. vt(i,j),wt(i,j) contain the c derivatives at colatitude theta(i) = (i-1)*pi/(nlat-1) c and longitude phi(j) = (j-1)*2*pi/nlon. the index ranges c are defined above at the input parameter ityp. vt and wt c are computed from the formulas for v and w given in c subroutine vhsec but with vbar and wbar replaced with c their derivatives with respect to colatitude. these c derivatives are denoted by vtbar and wtbar. c c in terms of real variables this expansion takes the form c c for i=1,...,nlat and j=1,...,nlon c c vt(i,j) = the sum from n=1 to n=nlat-1 of c c .5*br(1,n+1)*vtbar(0,n,theta(i)) c c plus the sum from m=1 to m=mmax-1 of the sum from n=m to c n=nlat-1 of the real part of c c (br(m+1,n+1)*vtbar(m,n,theta(i)) c -ci(m+1,n+1)*wtbar(m,n,theta(i)))*cos(m*phi(j)) c -(bi(m+1,n+1)*vtbar(m,n,theta(i)) c +cr(m+1,n+1)*wtbar(m,n,theta(i)))*sin(m*phi(j)) c c and for i=1,...,nlat and j=1,...,nlon c c wt(i,j) = the sum from n=1 to n=nlat-1 of c c -.5*cr(1,n+1)*vtbar(0,n,theta(i)) c c plus the sum from m=1 to m=mmax-1 of the sum from n=m to c n=nlat-1 of the real part of c c -(cr(m+1,n+1)*vtbar(m,n,theta(i)) c +bi(m+1,n+1)*wtbar(m,n,theta(i)))*cos(m*phi(j)) c +(ci(m+1,n+1)*vtbar(m,n,theta(i)) c -br(m+1,n+1)*wtbar(m,n,theta(i)))*sin(m*phi(j)) c c c br(m+1,nlat),bi(m+1,nlat),cr(m+1,nlat), and ci(m+1,nlat) are c assumed zero for m even. 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 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 lwvts c = 10 error in the specification of lwork c c c ******************************************************************* c c subroutine vtseci(nlat,nlon,wvts,lwvts,dwork,ldwork,ierror) c c subroutine vtseci initializes the array wvts which can then be c used repeatedly by subroutine vtsec until nlat or nlon is changed. 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 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 lwvts the dimension of the array wvts as it appears in the c program that calls vtsec. 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 lwvts must be at least c c 4*nlat*l2+3*max0(l1-2,0)*(nlat+nlat-l1-1)+nlon+15 c c c dwork a double precision work array that does not have to be saved. c c ldwork the dimension of the array work as it appears in the c program that calls vtsec. lwork must be at least c 2*(nlat+1) c c ************************************************************** c c output parameters c c wvts an array which is initialized for use by subroutine vtsec. c once initialized, wvts can be used repeatedly by vtsec c as long as nlat or nlon remain unchanged. wvts must not c be altered between calls of vtsec. 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 lwvts c = 4 error in the specification of ldwork c c ********************************************************************** c