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 ivlapec.f c c this file includes documentation and code for c subroutine ivlapec c c ... files which must be loaded with ivlapec.f c c sphcom.f, hrfft.f, vhaec.f, vhsec.f c c c c subroutine ivlapec(nlat,nlon,ityp,nt,v,w,idvw,jdvw,br,bi,cr,ci, c +mdbc,ndbc,wvhsec,lvhsec,work,lwork,ierror) c c c subroutine ivlapec computes a the vector field (v,w) whose vector c laplacian is (vlap,wlap). w and wlap are east longitudinal c components of the vectors. v and vlap are colatitudinal components c of the vectors. br,bi,cr, and ci are the vector harmonic coefficients c of (vlap,wlap). these must be precomputed by vhaec and are input c parameters to ivlapec. (v,w) have the same symmetry or lack of c symmetry about the about the equator as (vlap,wlap). the input c parameters ityp,nt,mdbc,ndbc must have the same values used by c vhaec to compute br,bi,cr, and ci for (vlap,wlap). c 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 longitude 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 this parameter should have the same value input to subroutine c vhaec to compute the coefficients br,bi,cr, and ci for the c vector field (vlap,wlap). ityp is set as follows: c c = 0 no symmetries exist in (vlap,wlap) about the equator. (v,w) c is computed and stored on the entire sphere in the arrays c arrays v(i,j) and w(i,j) for i=1,...,nlat and j=1,...,nlon. c c = 1 no symmetries exist in (vlap,wlap) about the equator. (v,w) c is computed and stored on the entire sphere in the arrays c v(i,j) and w(i,j) for i=1,...,nlat and j=1,...,nlon. the c vorticity of (vlap,wlap) is zero so the coefficients cr and c ci are zero and are not used. the vorticity of (v,w) is c also zero. c c c = 2 no symmetries exist in (vlap,wlap) about the equator. (v,w) c is computed and stored on the entire sphere in the arrays c w(i,j) and v(i,j) for i=1,...,nlat and j=1,...,nlon. the c divergence of (vlap,wlap) is zero so the coefficients br and c bi are zero and are not used. the divergence of (v,w) is c also zero. c c = 3 wlap is antisymmetric and vlap is symmetric about the c equator. consequently w is antisymmetric and v is symmetric. c (v,w) is computed and stored on the northern hemisphere c only. if nlat is odd, storage is in the arrays v(i,j), c w(i,j) for i=1,...,(nlat+1)/2 and j=1,...,nlon. if nlat c is even, storage is in the arrays v(i,j),w(i,j) for c i=1,...,nlat/2 and j=1,...,nlon. c c = 4 wlap is antisymmetric and vlap is symmetric about the c equator. consequently w is antisymmetric and v is symmetric. c (v,w) is computed and stored on the northern hemisphere c only. if nlat is odd, storage is in the arrays v(i,j), c w(i,j) for i=1,...,(nlat+1)/2 and j=1,...,nlon. if nlat c is even, storage is in the arrays v(i,j),w(i,j) for c i=1,...,nlat/2 and j=1,...,nlon. the vorticity of (vlap, c wlap) is zero so the coefficients cr,ci are zero and c are not used. the vorticity of (v,w) is also zero. c c = 5 wlap is antisymmetric and vlap is symmetric about the c equator. consequently w is antisymmetric and v is symmetric. c (v,w) is computed and stored on the northern hemisphere c only. if nlat is odd, storage is in the arrays w(i,j), c v(i,j) for i=1,...,(nlat+1)/2 and j=1,...,nlon. if nlat c is even, storage is in the arrays w(i,j),v(i,j) for c i=1,...,nlat/2 and j=1,...,nlon. the divergence of (vlap, c wlap) is zero so the coefficients br,bi are zero and c are not used. the divergence of (v,w) is also zero. c c c = 6 wlap is symmetric and vlap is antisymmetric about the c equator. consequently w is symmetric and v is antisymmetric. c (v,w) is computed and stored on the northern hemisphere c only. if nlat is odd, storage is in the arrays w(i,j), c v(i,j) for i=1,...,(nlat+1)/2 and j=1,...,nlon. if nlat c is even, storage is in the arrays w(i,j),v(i,j) for c i=1,...,nlat/2 and j=1,...,nlon. c c = 7 wlap is symmetric and vlap is antisymmetric about the c equator. consequently w is symmetric and v is antisymmetric. c (v,w) is computed and stored on the northern hemisphere c only. if nlat is odd, storage is in the arrays w(i,j), c v(i,j) for i=1,...,(nlat+1)/2 and j=1,...,nlon. if nlat c is even, storage is in the arrays w(i,j),v(i,j) for c i=1,...,nlat/2 and j=1,...,nlon. the vorticity of (vlap, c wlap) is zero so the coefficients cr,ci are zero and c are not used. the vorticity of (v,w) is also zero. c c = 8 wlap is symmetric and vlap is antisymmetric about the c equator. consequently w is symmetric and v is antisymmetric. c (v,w) is computed and stored on the northern hemisphere c only. if nlat is odd, storage is in the arrays w(i,j), c v(i,j) for i=1,...,(nlat+1)/2 and j=1,...,nlon. if nlat c is even, storage is in the arrays w(i,j),v(i,j) for c i=1,...,nlat/2 and j=1,...,nlon. the divergence of (vlap, c wlap) is zero so the coefficients br,bi are zero and c are not used. the divergence of (v,w) is also zero. c c c nt nt is the number of vector fields (vlap,wlap). some computational c efficiency is obtained for multiple fields. in the program c that calls ivlapec, the arrays v,w,br,bi,cr and ci can be c three dimensional corresponding to an indexed multiple vector c field. in this case multiple vector synthesis will be performed c to compute the (v,w) for each field (vlap,wlap). the third c index is the synthesis index which assumes the values k=1,...,nt. c for a single synthesis set nt=1. the description of the c remaining parameters is simplified by assuming that nt=1 or c that all arrays are two dimensional. c c idvw the first dimension of the arrays w and v as it appears in c the program that calls ivlapec. if ityp=0,1, or 2 then idvw c must be at least nlat. if ityp > 2 and nlat is even then idvw c must be at least nlat/2. if ityp > 2 and nlat is odd then idvw c must be at least (nlat+1)/2. c c jdvw the second dimension of the arrays w and v as it appears in c the program that calls ivlapec. jdvw must be at least nlon. c c c br,bi two or three dimensional arrays (see input parameter nt) c cr,ci that contain vector spherical harmonic coefficients of the c vector field (vlap,wlap) as computed by subroutine vhaec. c br,bi,cr and ci must be computed by vhaec prior to calling c ivlapec. if ityp=1,4, or 7 then cr,ci are not used and can c be dummy arguments. if ityp=2,5, or 8 then br,bi are not c used and can be dummy arguments. c c mdbc the first dimension of the arrays br,bi,cr and ci as it c appears in the program that calls ivlapec. mdbc must be c at least min0(nlat,nlon/2) if nlon is even or at least c min0(nlat,(nlon+1)/2) if nlon is odd. c c ndbc the second dimension of the arrays br,bi,cr and ci as it c appears in the program that calls ivlapec. ndbc must be at c least nlat. c c wvhsec an array which must be initialized by subroutine vhseci. c once initialized, wvhsec c can be used repeatedly by ivlapec as long as nlat and nlon c remain unchanged. wvhsec must not be altered between calls c of ivlapec. c c lvhsec the dimension of the array wvhsec as it appears in the c program that calls ivlapec. 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 then lvhsec must be at least c c 4*nlat*l2+3*max0(l1-2,0)*(nlat+nlat-l1-1)+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 ivlapec. 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 ityp .le. 2 then c c nlat*(2*nt*nlon+max0(6*l2,nlon)) + nlat*(4*nt*l1+1) c c or if ityp .gt. 2 let c c l2*(2*nt*nlon+max0(6*nlat,nlon)) + nlat*(4*nt*l1+1) c c will suffice as a minimum length for lwork c (see ierror=10 below) c c ************************************************************** c c output parameters c c c v,w two or three dimensional arrays (see input parameter nt) that c contain a vector field whose vector laplacian is (vlap,wlap). c w(i,j) is the east longitude and v(i,j) is the colatitudinal c component of the vector. v(i,j) and w(i,j) are given on the c sphere at the colatitude c c theta(i) = (i-1)*pi/(nlat-1) c c for i=1,...,nlat and east longitude c c lambda(j) = (j-1)*2*pi/nlon c c for j=1,...,nlon. c c let cost and sint be the cosine and sine at colatitude theta. c let d( )/dlambda and d( )/dtheta be the first order partial c derivatives in longitude and colatitude. let sf be either v c or w. define: c c del2s(sf) = [d(sint*d(sf)/dtheta)/dtheta + c 2 2 c d (sf)/dlambda /sint]/sint c c then the vector laplacian of (v,w) in (vlap,wlap) satisfies c c vlap = del2s(v) + (2*cost*dw/dlambda - v)/sint**2 c c and c c wlap = del2s(w) - (2*cost*dv/dlambda + w)/sint**2 c c c ierror a parameter which flags errors in input parameters as follows: c c = 0 no errors detected c c = 1 error in the specification of nlat c c = 2 error in the specification of nlon c c = 3 error in the specification of ityp c c = 4 error in the specification of nt c c = 5 error in the specification of idvw c c = 6 error in the specification of jdvw c c = 7 error in the specification of mdbc c c = 8 error in the specification of ndbc c c = 9 error in the specification of lvhsec c c = 10 error in the specification of lwork c c c ********************************************************************** c c end of documentation for ivlapec c c ********************************************************************** c