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 gradgc.f c c this file includes documentation and code for c subroutine gradgc i c c ... files which must be loaded with gradgc.f c c sphcom.f, hrfft.f, shagc.f,vhsgc.f c c subroutine gradgc(nlat,nlon,isym,nt,v,w,idvw,jdvw,a,b,mdab,ndab, c + wvhsgc,lvhsgc,work,lwork,ierror) c c given the scalar spherical harmonic coefficients a and b, precomputed c by subroutine shagc for a scalar field sf, subroutine gradgc computes c an irrotational vector field (v,w) such that c c gradient(sf) = (v,w). c c v is the colatitudinal and w is the east longitudinal component c of the gradient. i.e., c c v(i,j) = d(sf(i,j))/dtheta c c and c c w(i,j) = 1/sint*d(sf(i,j))/dlambda c c at the gaussian colatitude point theta(i) (see nlat as input c parameter) and longitude lambda(j) = (j-1)*2*pi/nlon where c where sint = sin(theta(i)). required associated legendre polynomials c are recomputed rather than stored as they are in subroutine gradgs. this c saves storage (compare lsav with lsav in gradgs) but increases c computational requirements. 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 than c 3. 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 shagc to compute the arrays a and b from the c scalar field sf. isym determines whether (v,w) are c computed on the full or half sphere as follows: c c = 0 c c sf 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 sf is antisymmetric about the equator. in this case w is c antisymmetric and v is symmetric 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 = 2 c c sf is symmetric 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 nt is the number of scalar and vector fields. some c computational efficiency is obtained for multiple fields. c the arrays a,b,v, and w can be three dimensional corresponding c to an indexed multiple array sf. in this case, multiple c vector synthesis will be performed to compute each vector c field. the third index for a,b,v, and w 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 c parameters is simplified by assuming that nt=1 or that a,b,v, c and w are two dimensional arrays. c c idvw the first dimension of the arrays v,w as it appears in c the program that calls gradgc. 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 gradgc. 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 scalar field array sf as computed by subroutine shagc. c *** a,b must be computed by shagc prior to calling gradgc. c c mdab the first dimension of the arrays a and b as it appears in c the program that calls gradgc (and shagc). 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 gradgc (and shagc). ndab must be at c least nlat. c c c wvhsgc an array which must be initialized by subroutine vhsgci. c once initialized, c wvhsgc can be used repeatedly by gradgc as long as nlon c and nlat remain unchanged. wvhsgc must not be altered c between calls of gradgc. c c c lvhsgc the dimension of the array wvhsgc as it appears in the c program that calls gradgc. Let 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 lvhsgc 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 gradgc. 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 if isym = 0 then lwork must be at least c c nlat*(2*nt*nlon+max0(6*l2,nlon)) + nlat*(2*l1*nt+1) c c if isym = 1 or 2 then lwork must be at least c c l2*(2*nt*nlon+max0(6*nlat,nlon)) + nlat*(2*l1*nt+1) c c c c ************************************************************** c c output parameters c c c v,w two or three dimensional arrays (see input parameter nt) that c contain an irrotational vector field such that the gradient of c the scalar field sf is (v,w). w(i,j) is the east longitude c component and v(i,j) is the colatitudinal component of velocity c at gaussian colatitude and longitude lambda(j) = (j-1)*2*pi/nlon c the indices for v and w are defined at the input parameter c isym. the vorticity of (v,w) is zero. note that any nonzero c vector field on the sphere will be multiple valued at the poles c [reference swarztrauber]. 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 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 lvhsgc c = 10 error in the specification of lwork c ********************************************************************** c c