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 ivrtgs.f c c this file includes documentation and code for c subroutine ivrtgs c c ... files which must be loaded with ivrtgs.f c c sphcom.f, hrfft.f, vhsgs.f,shags.f, gaqd.f c c c subroutine ivrtgs(nlat,nlon,isym,nt,v,w,idvw,jdvw,a,b,mdab,ndab, c + wvhsgs,lvhsgs,work,lwork,pertrb,ierror) c c given the scalar spherical harmonic coefficients a and b, precomputed c by subroutine shags for a scalar array vt, subroutine ivrtgs computes c a divergence free vector field (v,w) whose vorticity is vt - pertrb. c w is the east longitude component and v is the colatitudinal component. c pertrb is a constant which must be subtracted from vt for (v,w) to c exist (see the description of pertrb below). usually pertrb is zero c or small relative to vt. the divergence of (v,w), as computed by c ivrtgs, is the zero scalar field. v(i,j) and w(i,j) are the c colatitudinal and east longitude velocity components at gaussian c colatitude theta(i) (see nlat as input parameter) and longitude c lambda(j) = (j-1)*2*pi/nlon. the c c vorticity(v(i,j),w(i,j)) c c = [-dv/dlambda + d(sint*w)/dtheta]/sint c c = vort(i,j) - pertrb c c and c c divergence(v(i,j),w(i,j)) c c = [d(sint*v)/dtheta + dw/dlambda]/sint c c = 0.0 c c where sint = sin(theta(i)). required associated legendre polynomials c are stored rather than recomputed as they are in subroutine ivrtgc. 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 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 3. 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 c isym this has the same value as the isym that was input to c subroutine shags to compute the arrays a and b. isym c determines whether (v,w) are computed on the full or half c sphere as follows: c c = 0 c vt is not symmetric about the equator. in this case c the vector field (v,w) is computed on the entire sphere. c i.e., in the arrays v(i,j),w(i,j) for i=1,...,nlat and c j=1,...,nlon. c c = 1 c c vt is symmetric about the equator. in this case w is c antiymmetric and v is symmetric about the equator. v c and w are computed on the northern hemisphere only. i.e., c if nlat is odd they are computed for i=1,...,(nlat+1)/2 c and j=1,...,nlon. if nlat is even they are computed for c i=1,...,nlat/2 and j=1,...,nlon. c c = 2 c c vt is antisymmetric about the equator. in this case w is c symmetric and v is antisymmetric about the equator. w c and v are computed on the northern hemisphere only. i.e., c if nlat is odd they are computed for i=1,...,(nlat+1)/2 c and j=1,...,nlon. if nlat is even they are computed for c i=1,...,nlat/2 and j=1,...,nlon. c c c nt in the program that calls ivrtgs, nt is the number of vorticity c and vector fields. some computational efficiency is obtained c for multiple fields. the arrays a,b,v, and w can be three c dimensional and pertrb can be one dimensional corresponding c to an indexed multiple array vt. in this case, multiple vector c synthesis will be performed to compute each vector field. the c third index for a,b,v,w and first for pertrb is the synthesis c index which assumes the values k=1,...,nt. for a single c synthesis set nt=1. the description of the remaining parameters c is simplified by assuming that nt=1 or that a,b,v,w are two c dimensional and pertrb is a constant. c c idvw the first dimension of the arrays v,w as it appears in c the program that calls ivrtgs. if isym = 0 then idvw c must be at least nlat. if isym = 1 or 2 and nlat is c even then idvw must be at least nlat/2. if isym = 1 or 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 ivrtgs. jdvw must be at least nlon. c c a,b two or three dimensional arrays (see input parameter nt) c that contain scalar spherical harmonic coefficients c of the vorticity array vt as computed by subroutine shags. c *** a,b must be computed by shags prior to calling ivrtgs. c c mdab the first dimension of the arrays a and b as it appears in c the program that calls ivrtgs (and shags). mdab must be at c least min0(nlat,(nlon+2)/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 a and b as it appears in c the program that calls ivrtgs (and shags). ndab must be at c least nlat. c c c wvhsgs an array which must be initialized by subroutine vhsgsi. c once initialized c wvhsgs can be used repeatedly by ivrtgs as long as nlon c and nlat remain unchanged. wvhsgs must not be altered c between calls of ivrtgs. c c c lvhsgs the dimension of the array wvhsgs as it appears in the c program that calls ivrtgs. 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 lvhsgs must be at least c c (l1*l2*(nlat+nlat-l1+1))/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 ivrtgs. define c c l2 = nlat/2 if nlat is even or c l2 = (nlat+1)/2 if nlat is odd c l1 = min0(nlat,nlon/2) if nlon is even or c l1 = min0(nlat,(nlon+1)/2) if nlon is odd c c if isym = 0 then lwork must be at least c c (2*nt+1)*nlat*nlon + nlat*(2*nt*l1+1) c c if isym = 1 or 2 then lwork must be at least c c (2*nt+1)*l2*nlon + nlat*(2*nt*l1+1) c c c ************************************************************** c c output parameters c c c v,w two or three dimensional arrays (see input parameter nt) that c contain a divergence free vector field whose vorticity is c vt - pertrb at the gaussian colatitude point theta(i) c and longitude point lambda(j)=(j-1)*2*pi/nlon. w is the east c longitude component and v is the colatitudinal component. the c indices for v and w are defined at the input parameter isym. c the divergence of (v,w) is the zero scalar field. c c pertrb a nt dimensional array (see input parameter nt and assume nt=1 c for the description that follows). vt - pertrb is a scalar c field which can be the vorticity of a vector field (v,w). c pertrb is related to the scalar harmonic coefficients a,b c of vt (computed by shags) by the formula c c pertrb = a(1,1)/(2.*sqrt(2.)) c c an unperturbed vt can be the vorticity of a vector field c only if a(1,1) is zero. if a(1,1) is nonzero (flagged by c pertrb nonzero) then subtracting pertrb from vt yields a c scalar field for which a(1,1) is zero. 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 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 lvhsgs c = 10 error in the specification of lwork c ********************************************************************** c c