subroutine readspec(datanm,datag,what,itime_index,nyear,lfrst, * iggflag) c c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c ! ! c ! The names of the data files are hard-wired into the code. ! c ! ! c ! You'll need, the netCDF and UDUNITS libraries from UniData. ! c ! ! c ! netCDF can be found at: ! c ! http://www.unidata.ucar.edu/packages/netcdf/index.html ! c ! ! c ! UDUNITS can be found at: ! c ! http://www.unidata.ucar.edu/packages/udunits/index.html ! c ! ! c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c c==> reads spectral coefficient files from nmc reanalysis. c spectral coefficients returned in complex array datanm c (dimension 2016x28). c if iggflag = 1, data also returned on 192x94 gaussian grid c in array datag (dimension 192x94x28). c c input variables: c what is a character*4 field identifier. c allowable values are: 'vort','div ','shum','vair', c 'lnps' and 'orog' (vorticity, divergence, specific c humidity, virtual temperature, natural log of surface c pressure (in centibars) and orographic height (meters)). c c NOTE: if what = 'lnps' or 'orog', values are returned c in (datanm(nm,1),nm=1,2016) and datag(i,j,1). Other 27 levels of c these arrays are unchanged. c c itime_index is an integer time index denoting the c record to read from the netCDF file. c c nyear is an integer denoting the year of data desired c (1984, 1989 etc) c c lfrst is a logical value that must be set to .TRUE. on first c call so spherical harmonic transform arrays are initialized. c lfrst is set to .FALSE. internally on first call. c c iggflag = 1 means gaussian grid data is returned. c NOTE: the gaussian grid point j=1 is the point closest to c the north pole. i=1 is the greenwich meridian. c parameter (nlon = 192) parameter (nlat = 94) parameter (mwaves = 63) parameter (nmdim = (mwaves+1)*mwaves/2) parameter (levels = 28) c real datag(nlon,nlat,levels) complex datanm(nmdim,levels),scrm(mwaves,nlat) logical lfrst character*4 what c real trig((3*nlon/2)+1),aa(nlat*(nlon+2)),cc(nlat*(nlon+1)) double precision workdp(nlat,2) real worksp(4*nlat*(nlat+1)+2) integer ifax(13),nmodes(mwaves),indxm(nmdim),indxn(nmdim) real pnm(nmdim,nlat),hnm(nmdim,nlat), * gwrc(nlat),gaulats(nlat),weights(nlat), * del(nmdim),deli(nmdim) save trig,aa,cc,ifax,nmodes,indxm,indxn,pnm,hnm,gwrc, * gaulats,weights,del,deli c if (lfrst) then call sphinit(6.3712e6,nlon,nlat,mwaves,nmdim, * workdp,worksp, * pnm,hnm,gwrc,gaulats,weights, * del,deli,nmodes,indxm,indxn) call set99(trig,ifax,nlon) lfrst = .FALSE. end if c if (what .eq. 'vort' .or. what .eq. 'div ' .or. * what .eq. 'shum' .or. what .eq. 'vair') then call specread(datanm,what,itime_index,nyear) if (iggflag .eq. 1) then do k=1,levels call spharm(datag(1,1,k),datanm(1,k),scrm, * pnm,weights,nmodes, * aa,cc,ifax,trig, * nlon,nlat,mwaves,nmdim,-1) enddo end if else if (what .eq. 'lnps') then call specreadps(datanm(1,1),itime_index,nyear) if (iggflag .eq. 1) then call spharm(datag(1,1,1),datanm(1,1),scrm, * pnm,weights,nmodes, * aa,cc,ifax,trig, * nlon,nlat,mwaves,nmdim,-1) end if else if (what .eq. 'orog') then call specreadzs(datanm(1,1)) if (iggflag .eq. 1) then call spharm(datag(1,1,1),datanm(1,1),scrm, * pnm,weights,nmodes, * aa,cc,ifax,trig, * nlon,nlat,mwaves,nmdim,-1) end if else write(6,*) 'invalid field identifier!' write(6,*) 'allowed are:' write(6,*) 'vort,div ,shum,vair,lnps,orog' stop end if c return end c subroutine getuvspec(vrtnm,divnm,ug,vg,itime_index,nyear,lfrst) c c==> reads spectral coefficients of vorticity and divergence c from nmc reanalysis, converts data to u*cos(lat) and c v*cos(lat) on 192x94 gaussian grid. c c on output, vrtnm(2016,28), divnm(2016,28) contain c spectral coefficients of vorticity and divergence. c ug(192,94,28) and vg(192,94,28) contain u*cos(lat) c and v*cos(lat) on gaussian grid. c NOTE: the gaussian grid point j=1 is the point closest to c the north pole. i=1 is the greenwich meridian. c c itime_index is an integer time index denoting the c record to read from the netCDF file. c c nyear is an integer denoting the year of data desired c (1984, 1989 etc) c c lfrst is a logical value that must be set to .TRUE. on first c call so spherical harmonic transform arrays are initialized. c lfrst is set to .FALSE. internally on first call. c c parameter (nlon = 192) parameter (nlat = 94) parameter (mwaves = 63) parameter (nmdim = (mwaves+1)*mwaves/2) parameter (levels = 28) c real ug(nlon,nlat,levels),vg(nlon,nlat,levels) complex vrtnm(nmdim,levels),divnm(nmdim,levels), * um(mwaves,nlat),vm(mwaves,nlat) logical lfrst c real trig((3*nlon/2)+1),aa(nlat*(nlon+2)),cc(nlat*(nlon+1)) double precision workdp(nlat,2) real worksp(4*nlat*(nlat+1)+2) integer ifax(13),nmodes(mwaves),indxm(nmdim),indxn(nmdim) real pnm(nmdim,nlat),hnm(nmdim,nlat), * gwrc(nlat),gaulats(nlat),weights(nlat), * del(nmdim),deli(nmdim) save trig,aa,cc,ifax,nmodes,indxm,indxn,pnm,hnm,gwrc, * gaulats,weights,del,deli c if (lfrst) then call sphinit(6.3712e6,nlon,nlat,mwaves,nmdim, * workdp,worksp, * pnm,hnm,gwrc,gaulats,weights, * del,deli,nmodes,indxm,indxn) call set99(trig,ifax,nlon) LFRST = .FALSE. end if c call specread(vrtnm,'vort',itime_index,nyear) call specread(divnm,'div ',itime_index,nyear) do k=1,levels call getuv(vrtnm(1,k),divnm(1,k),um,vm, * ug(1,1,k),vg(1,1,k),6.3712e6, * nlon,nlat,mwaves,nmdim, * aa,cc,ifax,trig, * pnm,hnm,deli,nmodes) enddo c return end c subroutine specread(datanm,what,time_index,nyear) c c==> reads spectral sigma files for what = 'vort', 'div', c 'shum' and 'vair'. c Input: what (an character*4 denoting the name of the variable), c time_index (an integer denoting the output time c to read from the file) and nyear (an integer denoting c the year). c datanm is the output complex 2016x28 array of spectral coefficients. c INCLUDE '/usr/local/netcdf/include/netcdf2.3.inc' c parameter (mwaves = 63) parameter (nmdim = (mwaves+1)*mwaves/2) parameter (nmdim2 = 2*nmdim) parameter (levels = 28) parameter (levels2 = 56) c complex datanm(nmdim, levels) c character*4 what character*80 ncfilename c real mean(levels2), add_offset(levels), scale_factor(levels) integer*2 cscaled(2*nmdim) c integer time_index,rcode C integer ncid, specid, meanid, add_offsetid, scale_factorid c integer cstart(3), ccount(3) data cstart /1,1,1/ data ccount /nmdim2, 1, 1/ integer mstart(2), mcount(2) data mstart /1,1/ data mcount /levels2,1/ c integer astart(2), acount(2) data astart /1,1/ data acount /levels,1/ mstart(2) = time_index astart(2) = time_index cstart(3) = time_index len_what = index(what,' ') - 1 if (len_what .eq. -1) len_what = 4 if ( len_what .eq. 3 ) then write(ncfilename, 10) what(1:len_what), nyear - 1900 10 format('/dmdata/Datasets.public/nmc.reanalysis/spectral/', * a3,'.spec.',i2,'.nc') else if ( len_what .eq. 4 ) then write(ncfilename, 20) what(1:len_what), nyear - 1900 20 format('/dmdata/Datasets.public/nmc.reanalysis/spectral/', * a4,'.spec.',i2,'.nc') end if ncid = ncopn(ncfilename, NCREAD, RCODE) specid = ncvid(NCID,what,RCODE) meanid = ncvid(ncid, 'mean', RCODE) add_offsetid = ncvid(ncid, 'add_offset', RCODE) scale_factorid = ncvid(ncid, 'scale_factor', RCODE) call ncvgt(ncid, meanid, mstart, $ mcount, mean, RCODE) call ncvgt(ncid, add_offsetid, astart, $ acount, add_offset, RCODE) call ncvgt(ncid, scale_factorid, astart, $ acount, scale_factor, RCODE) l = 1 do ii = 1, levels cstart(2) = ii call ncvgt(ncid, specid, cstart, $ ccount, cscaled, RCODE) k = 3 do i = 1, nmdim if( i .eq. 1 ) then datanm(i,ii) = cmplx(mean(l), mean(l+1)) l = l + 2 else unpackr = scale_factor(ii)* $ float(cscaled(k)) + add_offset(ii) unpackc = scale_factor(ii)* $ float(cscaled(k+1)) + add_offset(ii) datanm(i,ii) = cmplx(unpackr, $ unpackc ) k = k + 2 endif enddo enddo call ncclos(ncid, RCODE) return end c subroutine specreadps(datanm,time_index,nyear) c c==> reads spectral sigma files for natural log of sfce pressure. c Input: time_index (an integer denoting the output time c to read from the file) and nyear (an integer denoting c the year). c datanm is the output complex 2016 element array of spectral coefficients. c INCLUDE '/usr/local/netcdf/include/netcdf2.3.inc' c parameter (mwaves = 63) parameter (nmdim = (mwaves+1)*mwaves/2) parameter (nmdim2 = 2*nmdim) c complex datanm(nmdim) c character*4 what character*80 ncfilename c real mean(2), add_offset, scale_factor integer*2 cscaled(2*nmdim) c integer time_index,rcode C integer ncid, specid, meanid, add_offsetid, scale_factorid c integer cstart(2), ccount(2) data cstart /1,1/ data ccount /nmdim2, 1/ integer mstart(2), mcount(2) data mstart /1,1/ data mcount /2,1/ integer astart data astart /1/ mstart(2) = time_index astart = time_index cstart(2) = time_index what = 'pres' write(ncfilename, 10) nyear - 1900 10 format('/dmdata/Datasets.public/nmc.reanalysis/spectral/', * 'pres.nlog.sfc.spec.', * i2,'.nc') ncid = ncopn(ncfilename, NCREAD, RCODE) specid = ncvid(NCID,what,RCODE) meanid = ncvid(ncid, 'mean', RCODE) add_offsetid = ncvid(ncid, 'add_offset', RCODE) scale_factorid = ncvid(ncid, 'scale_factor', RCODE) call ncvgt(ncid, meanid, mstart, $ mcount, mean, RCODE) call ncvgt1(ncid, add_offsetid, astart, $ add_offset, RCODE) call ncvgt1(ncid, scale_factorid, astart, $ scale_factor, RCODE) call ncvgt(ncid, specid, cstart, $ ccount, cscaled, RCODE) k = 3 do i = 1, nmdim if( i .eq. 1 ) then datanm(i) = cmplx(mean(1), mean(2)) else unpackr = scale_factor* $ float(cscaled(k)) + add_offset unpackc = scale_factor* $ float(cscaled(k+1)) + add_offset datanm(i) = cmplx(unpackr, $ unpackc ) k = k + 2 endif enddo call ncclos(ncid, RCODE) return end c subroutine specreadzs(datanm) c c==> reads the spectral orography file. c datanm is the output complex 2016 element array of spectral coefficients. c parameter (mwaves = 63) parameter (nmdim = (mwaves+1)*mwaves/2) c character*80 input_file complex datanm(nmdim) character*4 what integer istart(1),icount(1) c INCLUDE '/usr/local/netcdf/include/netcdf2.3.inc' c input_file='/dmdata/Datasets.public/nmc.reanalysis/'// * 'spectral/orog.spec.nc' ncid = ncopn(input_file, NCREAD, IERCODE) what = 'orog' ispecid = ncvid(NCID,what,IERCODE) istart(1) = 1 icount(1) = 2*nmdim call ncvgt(ncid, ispecid, istart, $ icount, datanm, IERCODE) call ncclos(ncid,IERCODE) c return end c subroutine sphinit(re,nlons,nlats,mwaves,nmdim, * workdp,worksp, * pnm,hnm,gwrc,gaulats,weights, * del,deli,nmodes,indxm,indxn) c double precision workdp(nlats,2) real worksp(4*nlats*(nlats+1)+2) c real pnm(nmdim,nlats), * hnm(nmdim,nlats), * gwrc(nlats),gaulats(nlats), * weights(nlats), * del(nmdim),deli(nmdim) integer nmodes(mwaves),indxm(nmdim), * indxn(nmdim) c pi = 4.*atan(1.0) C C==> COMPUTE USEFUL INDICES: C INDXM(NM) CONTAINS THE VALUES OF THE FOURIER WAVENUMBERS. C INDXN(NM) CONTAINS THE VALUES OF THE N INDEX C NMODES(I) CONTAINS THE NUMBER OF N MODES FOR EACH M C DEL(NM) CONTAINS N*(N+1). c deli(nm) is the inverse of del(nm). C NMSTRT = 0 DO 200 I=1,MWAVES M = I - 1 NMODES(I) = MWAVES - M DO 100 J=1,NMODES(I) NM = NMSTRT + J N = M + J - 1 INDXM(NM) = M INDXN(NM) = N DEL(NM) = N*(N+1) if (n .ne. 0) deli(nm) = 1./del(nm) 100 CONTINUE NMSTRT = NMSTRT + NMODES(I) 200 CONTINUE deli(1) = 0.0 C C==> CALCULATE GAUSSIAN WEIGHTS AND SIN(LATITUDES). C call gaqd(nlats,workdp(1,1),workdp(1,2), * worksp,4*nlats*(nlats+1)+2,ierror) if (ierror .ne. 0) then write(6,*) ' failure in GAQD! IERROR = ',ierror stop end if c c==> convert gaulats from colatitude to sin(latitude), going from c north pole to south pole. c do 250 j=1,nlats gaulats(j) = dcos(workdp(j,1)) weights(j) = workdp(j,2) 250 continue c DO 300 JLAT=1,NLATS C C==> FILL GWRC(J) WITH THE GAUSSIAN WEIGHTS DIVIDED BY THE a times the C cosine SQUARED OF THE LATITUDE FROM THE NORTH POLE TO THE C south pole, INCLUSIVE. C SINLAT = GAULATS(JLAT) COSSQ = 1. - SINLAT**2 GWRC(JLAT) = WEIGHTS(JLAT)/(RE*COSSQ) C C==> COMPUTE ASSOCIATED LEGENDRE POLYNOMIALS AND THEIR DERIVATIVES C FOR EACH LATITUDE (going from north pole to south pole). c THE ASSOCIATED LEGENDRE FUNCTIONS ARE DEFINED: C C PBAR(M,N,THETA) = SQRT((2*N+1)*FACTORIAL(N-M)/(2*FACTORIAL(N+M))) C *SIN(THETA)**M/(2**N*FACTORIAL(N)) TIMES THE C (N+M)TH DERIVATIVE OF (X**2-1)**N WITH RESPECT C TO X=COS(THETA) c c note: theta = 0.5*pi - phi, where phi is latitude and c theta is colatitude. c Therefore, cos(theta) = sin(phi) and sin(theta) = cos(phi). C C where n is degree (subscript), m is order C (superscript), and theta is colatitude in C radians. The functions are orthogonal C polynomials on the surface of the sphere and C are normalized so that the integral of C (PBAR(n,m,theta)**2)*sin(theta) on the interval C theta=0 to theta=pi equals 1. C C NOTE THAT PBAR(0,0,THETA) = SQRT(2)/2, AND C PBAR(1,0,THETA) = 0.5*SQRT(6)*SIN(LAT). C C THE DERIVATIVES OF PBAR (HBAR) ARE DEFINED SUCH THAT C HBAR(M,N,THETA) = -(1-X**2)*D(PBAR)/DX, WHERE X IS C COS OF COLATITUDE (OR SINE OF LATITUDE). C CALL LEGEND(SINLAT,NMODES,MWAVES,NMDIM, * PNM(1,JLAT),HNM(1,JLAT)) 300 CONTINUE c return end c SUBROUTINE DINTQL(NM,N,D,E,Z,IERR) C DINTQL IS A DOUBLE PRECISION MODIFICATION OF INTQL2 OFF C EISPACK TO BE USED BY GAQD IN SPHEREPACK FOR COMPUTING C GAUSSIAN WEIGHTS AND POINTS C INTEGER I,J,K,L,M,N,II,NM,MML,IERR DOUBLE PRECISION D(N),E(N),Z(NM,N) DOUBLE PRECISION B,C,F,G,P,R,S,TST1,TST2,DPYTHA IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C E(N) = 0.0D0 C DO 240 L = 1, N J = 0 C .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... 105 DO 110 M = L, N-1 TST1 = DABS(D(M)) + DABS(D(M+1)) TST2 = TST1 + DABS(E(M)) IF (TST2 .EQ. TST1) GO TO 120 110 CONTINUE M = N C 120 P = D(L) IF (M .EQ. L) GO TO 240 IF (J .EQ. 30) GO TO 1000 J = J + 1 C .......... FORM SHIFT .......... G = (D(L+1) - P) / (2.0D0 * E(L)) R = DPYTHA(G,1.0D0) G = D(M) - P + E(L) / (G + SIGN(R,G)) S = 1.0D0 C = 1.0D0 P = 0.0D0 MML = M - L C .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... DO 200 II = 1, MML I = M - II F = S * E(I) B = C * E(I) R = DPYTHA(F,G) E(I+1) = R IF (R .EQ. 0.0D0) GO TO 210 S = F / R C = G / R G = D(I+1) - P R = (D(I) - G) * S + 2.0D0 * C * B P = S * R D(I+1) = G + P G = C * R - B C .......... FORM VECTOR .......... DO 180 K = 1, N F = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * F Z(K,I) = C * Z(K,I) - S * F 180 CONTINUE C 200 CONTINUE C D(L) = D(L) - P E(L) = G E(M) = 0.0D0 GO TO 105 C .......... RECOVER FROM UNDERFLOW .......... 210 D(I+1) = D(I+1) - P E(M) = 0.0D0 GO TO 105 240 CONTINUE C .......... ORDER EIGENVALUES AND EIGENVECTORS .......... DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .GE. P) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C .......... SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS .......... 1000 IERR = L 1001 RETURN END DOUBLE PRECISION FUNCTION DPYTHA(A,B) DOUBLE PRECISION A,B C DPYTHA IS A DOUBLE PRECISION MODIFICATION OF PYTHAG OFF EISPACK C FOR USE BY DIMTQL C C FINDS SQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW C DOUBLE PRECISION P,R,S,T,U P = DABS(A) IF (DABS(B).GE.DABS(A)) P = DABS(B) C P = AMAX1(DABS(A),DABS(B)) IF (P .EQ. 0.0D0) GO TO 20 R = (DABS(A)/P)**2 IF (DABS(B).LT.DABS(A)) R = (DABS(B)/P)**2 C R = (AMIN1(DABS(A),DABS(B))/P)**2 10 CONTINUE T = 4.0D0 + R IF (T .EQ. 4.0D0) GO TO 20 S = R/T U = 1.0D0 + 2.0D0*S P = U*P R = (S/U)**2 * R GO TO 10 20 DPYTHA = P RETURN END SUBROUTINE DRST(NM,N,W,E,MATZ,Z,IERR) C DRST IS A DOUBLE PRECISION MODIFICATION OF RST OFF EISPACK C TO BE USED TO COMPUTE GAUSSIAN POINTS AND WEIGHTS C INTEGER I,J,N,NM,IERR,MATZ DOUBLE PRECISION W(N),E(N),Z(NM,N) C C .......... FIND BOTH EIGENVALUES AND EIGENVECTORS .......... 20 DO 40 I = 1, N C DO 30 J = 1, N Z(J,I) = 0.0D0 30 CONTINUE C Z(I,I) = 1.0D0 40 CONTINUE C CALL DINTQL(NM,N,W,E,Z,IERR) RETURN END C SUBROUTINE GAQD(NLAT,THETA,WTS,WORK,LWORK,IERROR) C C SUBROUTINE GAQD COMPUTES THE NLAT GAUSSIAN COLATITUDES AND WEIGHTS C IN DOUBLE PRECISION. THE COLATITUDES ARE IN RADIANS AND LIE IN THE C IN THE INTERVAL (0,PI). C C INPUT PARAMETERS C C NLAT THE NUMBER OF GAUSSIAN COLATITUDES IN THE INTERVAL (0,PI) C (BETWEEN THE TWO POLES). NLAT MUST BE GREATER THAN ZERO. C C WORK A TEMPORARY WORK SPACE C C LWORK THE LENGTH OF THE WORK SPACE IN THE ROUTINE CALLING GAQD C LWORK MUST BE AT LEAST 4*NLAT*(NLAT+1)+2. C C OUTPUT PARAMETERS C C THETA A DOUBLE PRECISION VECTOR OF LENGTH NLAT CONTAINING THE C NLAT GAUSSIAN COLATITUDES ON THE SPHERE IN INCREASING RADI C IN THE INTERVAL (O,PI). C C WTS A DOUBLE PRECISION VECTOR OF LENGTH NLAT CONTAINING THE C NLAT GAUSSIAN WEIGHTS. C C IERROR = 0 NO ERRORS C = 1 IF LWORK.LT.4*NLAT*(NLAT+1)+2 C = 2 IF NLAT.LE.0 C = 3 IF UNABLE TO COMPUTE GAUSSIAN POINTS C (FAILURE IN IN EIGENVALUE ROUTINE) C C ***************************************************************** SUBROUTINE GAQD(NLAT,THETA,WTS,WORK,LWORK,IERROR) DIMENSION WORK(LWORK),THETA(NLAT),WTS(NLAT) DOUBLE PRECISION THETA,WTS,X N = NLAT IERROR = 1 C CHECK WORK SPACE LENGTH IF (LWORK.LT.4*N*(N+1)+2) RETURN IERROR = 2 IF (N.LE.0) RETURN IERROR = 0 IF (N.GT.2) THEN C PARTITION WORK SPACE FOR DOUBLE PRECISION EIGENVALUE(VECTOR COMPUT I1 = 1 I2 = I1+2*N I3 = I2+2*N CALL GAQD1(N,THETA,WTS,WORK(I1),WORK(I2),WORK(I3),IERROR) IF (IERROR.NE.0) THEN IERROR = 3 RETURN END IF RETURN ELSE IF (N.EQ.1) THEN WTS(1) = 2.0D0 THETA(1) = DACOS(0.0D0) ELSE IF (N.EQ.2) THEN C COMPUTE WEIGHTS AND POINTS ANALYTICALLY WHEN N=2 WTS(1) = 1.0D0 WTS(2) = 1.0D0 X = DSQRT(1.0D0/3.0D0) THETA(1) = DACOS(X) THETA(2) = DACOS(-X) RETURN END IF END SUBROUTINE GAQD1(N,THETA,WTS,W,E,WRK,IER) DIMENSION THETA(N),WTS(N),W(N),E(N),WRK(1) DOUBLE PRECISION THETA,WTS,TEMP,W,E,WRK C SET SYMMETRIC TRIDIAGNONAL MATRIX SUBDIAGONAL AND DIAGONAL C COEFFICIENTS FOR MATRIX COMING FROM COEFFICIENTS IN THE C RECURSION FORMULA FOR LEGENDRE POLYNOMIALS C A(N)*P(N-1)+B(N)*P(N)+C(N)*P(N+1) = 0. WRK(1)=0.D0 WRK(N+1) = 0.D0 W(1)=0.D0 E(1) = 0.D0 DO 100 J=2,N WRK(J)= (J-1.D0)/DSQRT((2.D0*J-1.D0)*(2.D0*J-3.D0)) WRK(J+N)=0.D0 E(J) = WRK(J) W(J) = 0.D0 100 CONTINUE C COMPUTE EIGENVALUES OF MATRIX MATZ = 1 INDX = 2*N+1 NP1=N+1 CALL DRST(N,N,W,E,MATZ,WRK(INDX),IER) IF (IER.NE.0) RETURN C COMPUTE GAUSSIAN WEIGHTS AND POINTS DO 101 J=1,N THETA(J) = DACOS(W(J)) C SET GAUSSIAN WEIGHTS AS 1ST COMPONENTS OF EIGENVECTORS SQUARED INDX = (J+1)*N+1 WTS(J) = 2.0D0*WRK(INDX)**2 101 CONTINUE C REVERSE ORDER OF GAUSSIAN POINTS TO BE C MONOTONIC INCREASING IN RADIANS N2 = N/2 DO 102 I=1,N2 TEMP = THETA(I) THETA(I) = THETA(N-I+1) THETA(N-I+1) = TEMP 102 CONTINUE RETURN END c SUBROUTINE LEGEND(X,NMODES,MWAVES,NMDIM, * PMN, HMN) C C THIS SUBROUTINE COMPUTES ASSOCIATED LEGENDRE C FUNCTIONS, PMN, AND THE DERIVATIVE QUANTITY C HMN = -(1 - X*X) * D(PMN)/DX AT X = COS( COLATITUDE ) C GIVEN SOME COLATITUDE IN THE DOMAIN. C C ACCURACY: C C THE RECURSION RELATION EMPLOYED IS LESS ACCURATE C NEAR THE POLES FOR M + N .GE. ROUGHLY 75 BECAUSE THE RECURSION C INVOLVES THE DIFFERENCE OF NEARLY EQUAL NUMBERS, SO C THE ASSOCIATED LEGENDRE FUNCTIONS ARE COMPUTED USING DOUBLE C PRECISION ACCUMULATORS. C SEE BELOUSOV, TABLES OF NORMALIZED ASSOCIATED LEGENDRE C POLYNOMIALS, C. 1962, FOR INFORMATION ON THE ACCURATE C COMPUTATION OF ASSOCIATED LEGENDRE FUNCTIONS. C REAL PMN(NMDIM), X HMN(NMDIM) C INTEGER NMODES(MWAVES) C double precision A, B, PROD, SINSQ, X eps,epsnmp1,xd,EPSPMN, PMNJ, PMNJM1, PMNJM2 C C C**** SET PARAMETERS FOR ENTRY INTO THE RECURSIVE FORMULAE. C C xd = dble(x) SINSQ = 1.D0 - XD * XD C A = 1.D0 B = 0.D0 PROD = 1.D0 C C C**** LOOP FOR THE 'M' INDEX. C NMSTRT = 0 DO 70 I = 1, MWAVES C M = (I - 1) NMAX = NMODES(I) C C WHEN M=0 (I=1), STANDARD LEGENDRE POLYNOMIALS ARE C GENERATED. C IF (M .NE. 0) THEN A = A + 2.D0 B = B + 2.D0 PROD = PROD * SINSQ * A / B END IF C C**** GENERATE PMN AND HMN FOR J = 1 AND 2. C PMNJM2 = SQRT(0.5D0 * PROD) NM = NMSTRT + 1 PMN(NM) = PMNJM2 C PMNJM1 = SQRT( DBLE(2 * M + 3) ) * XD * PMNJM2 IF (NM .NE. NMDIM) * PMN(NM+1) = PMNJM1 C NP1 = M + 1 EPSNMP1 = SQRT( DBLE(NP1*NP1 - M*M) / DBLE(4*NP1*NP1 - 1) ) EPSPMN = XD * PMNJM1 - EPSNMP1 * PMNJM2 C HMN(NM) = DBLE(M) * EPSNMP1 * PMNJM1 IF (NM .NE. NMDIM) * HMN(NM+1) = DBLE(M+1) * EPSPMN - * DBLE(M+2) * EPSNMP1 * PMNJM2 c C**** LOOP FOR THE 'N' INDEX. C NOW APPLY THE RECURSION FORMULAE FOR J .GE. 3. C DO 50 J = 3, NMAX N = M + J - 1 NM = NMSTRT + J EPS = SQRT( DBLE(N*N - M*M) / DBLE(4*N*N - 1) ) PMNJ = EPSPMN / EPS PMN(NM) = PMNJ C C COMPUTE EPS * PMN FOR J+1. C EPSPMN = XD * PMNJ - EPS * PMNJM1 C C COMPUTE THE DERIVATIVE. C HMN(NM) = DBLE(N) * EPSPMN - * DBLE(N+1) * EPS * PMNJM1 C PMNJM2 = PMNJM1 PMNJM1 = PMNJ 50 CONTINUE NMSTRT = NMSTRT + NMAX 70 CONTINUE RETURN END c SUBROUTINE FFT99(A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN) C C PURPOSE PERFORMS MULTIPLE FAST FOURIER TRANSFORMS. THIS PACKAGE C WILL PERFORM A NUMBER OF SIMULTANEOUS REAL/HALF-COMPLEX C PERIODIC FOURIER TRANSFORMS OR CORRESPONDING INVERSE C TRANSFORMS, I.E. GIVEN A SET OF REAL DATA VECTORS, THE C PACKAGE RETURNS A SET OF 'HALF-COMPLEX' FOURIER C COEFFICIENT VECTORS, OR VICE VERSA. THE LENGTH OF THE C TRANSFORMS MUST BE AN EVEN NUMBER GREATER THAN 4 THAT HAS C NO OTHER FACTORS EXCEPT POSSIBLY POWERS OF 2, 3, AND 5. C THIS IS AN ALL-FORTRAN VERSION OF A OPTIMIZED ROUTINE C FFT99 WRITTEN FOR XMP/YMPs BY DR. CLIVE TEMPERTON OF C ECMWF. C C THE PACKAGE FFT99F CONTAINS SEVERAL USER-LEVEL ROUTINES: C C SUBROUTINE SET99 C AN INITIALIZATION ROUTINE THAT MUST BE CALLED ONCE C BEFORE A SEQUENCE OF CALLS TO THE FFT ROUTINES C (PROVIDED THAT N IS NOT CHANGED). C C SUBROUTINES FFT99 AND FFT991 C TWO FFT ROUTINES THAT RETURN SLIGHTLY DIFFERENT C ARRANGEMENTS OF THE DATA IN GRIDPOINT SPACE. C C USAGE LET N BE OF THE FORM 2**P * 3**Q * 5**R, WHERE P .GE. 1, C Q .GE. 0, AND R .GE. 0. THEN A TYPICAL SEQUENCE OF C CALLS TO TRANSFORM A GIVEN SET OF REAL VECTORS OF LENGTH C N TO A SET OF 'HALF-COMPLEX' FOURIER COEFFICIENT VECTORS C OF LENGTH N IS C C DIMENSION IFAX(13),TRIGS(3*N/2+1),A(M*(N+2)), C + WORK(M*(N+1)) C C CALL SET99 (TRIGS, IFAX, N) C CALL FFT99 (A,WORK,TRIGS,IFAX,INC,JUMP,N,M,ISIGN) C C SEE THE INDIVIDUAL WRITE-UPS FOR SET99, FFT99, AND C FFT991 BELOW, FOR A DETAILED DESCRIPTION OF THE C ARGUMENTS. C C HISTORY THE PACKAGE WAS WRITTEN BY CLIVE TEMPERTON AT ECMWF IN C NOVEMBER, 1978. IT WAS MODIFIED, DOCUMENTED, AND TESTED C FOR NCAR BY RUSS REW IN SEPTEMBER, 1980. C C----------------------------------------------------------------------- C C SUBROUTINE SET99 (TRIGS, IFAX, N) C C PURPOSE A SET-UP ROUTINE FOR FFT99 AND FFT991. IT NEED ONLY BE C CALLED ONCE BEFORE A SEQUENCE OF CALLS TO THE FFT C ROUTINES (PROVIDED THAT N IS NOT CHANGED). C C ARGUMENT IFAX(13),TRIGS(3*N/2+1) C DIMENSIONS C C ARGUMENTS C C ON INPUT TRIGS C A FLOATING POINT ARRAY OF DIMENSION 3*N/2 IF N/2 IS C EVEN, OR 3*N/2+1 IF N/2 IS ODD. C C IFAX C AN INTEGER ARRAY. THE NUMBER OF ELEMENTS ACTUALLY USED C WILL DEPEND ON THE FACTORIZATION OF N. DIMENSIONING C IFAX FOR 13 SUFFICES FOR ALL N LESS THAN A MILLION. C C N C AN EVEN NUMBER GREATER THAN 4 THAT HAS NO PRIME FACTOR C GREATER THAN 5. N IS THE LENGTH OF THE TRANSFORMS (SEE C THE DOCUMENTATION FOR FFT99 AND FFT991 FOR THE C DEFINITIONS OF THE TRANSFORMS). C C ON OUTPUT IFAX C CONTAINS THE FACTORIZATION OF N/2. IFAX(1) IS THE C NUMBER OF FACTORS, AND THE FACTORS THEMSELVES ARE STORED C IN IFAX(2),IFAX(3),... IF SET99 IS CALLED WITH N ODD, C OR IF N HAS ANY PRIME FACTORS GREATER THAN 5, IFAX(1) C IS SET TO -99. C C TRIGS C AN ARRAY OF TRIGONOMETRIC FUNCTION VALUES SUBSEQUENTLY C USED BY THE FFT ROUTINES. C C----------------------------------------------------------------------- C C SUBROUTINE FFT991 (A,WORK,TRIGS,IFAX,INC,JUMP,N,M,ISIGN) C AND C SUBROUTINE FFT99 (A,WORK,TRIGS,IFAX,INC,JUMP,N,M,ISIGN) C C PURPOSE PERFORM A NUMBER OF SIMULTANEOUS REAL/HALF-COMPLEX C PERIODIC FOURIER TRANSFORMS OR CORRESPONDING INVERSE C TRANSFORMS, USING ORDINARY SPATIAL ORDER OF GRIDPOINT C VALUES (FFT991) OR EXPLICIT CYCLIC CONTINUITY IN THE C GRIDPOINT VALUES (FFT99). GIVEN A SET C OF REAL DATA VECTORS, THE PACKAGE RETURNS A SET OF C 'HALF-COMPLEX' FOURIER COEFFICIENT VECTORS, OR VICE C VERSA. THE LENGTH OF THE TRANSFORMS MUST BE AN EVEN C NUMBER THAT HAS NO OTHER FACTORS EXCEPT POSSIBLY POWERS C OF 2, 3, AND 5. THIS IS AN ALL-FORTRAN VERSION OF C OPTIMIZED ROUTINE FFT991 WRITTEN FOR XMP/YMPs BY C DR. CLIVE TEMPERTON OF ECMWF. C C ARGUMENT A(M*(N+2)), WORK(M*(N+1)), TRIGS(3*N/2+1), IFAX(13) C DIMENSIONS C C ARGUMENTS C C ON INPUT A C AN ARRAY OF LENGTH M*(N+2) CONTAINING THE INPUT DATA C OR COEFFICIENT VECTORS. THIS ARRAY IS OVERWRITTEN BY C THE RESULTS. C C WORK C A WORK ARRAY OF DIMENSION M*(N+1) C C TRIGS C AN ARRAY SET UP BY SET99, WHICH MUST BE CALLED FIRST. C C IFAX C AN ARRAY SET UP BY SET99, WHICH MUST BE CALLED FIRST. C C INC C THE INCREMENT (IN WORDS) BETWEEN SUCCESSIVE ELEMENTS OF C EACH DATA OR COEFFICIENT VECTOR (E.G. INC=1 FOR C CONSECUTIVELY STORED DATA). C C JUMP C THE INCREMENT (IN WORDS) BETWEEN THE FIRST ELEMENTS OF C SUCCESSIVE DATA OR COEFFICIENT VECTORS. ON CRAYS, C TRY TO ARRANGE DATA SO THAT JUMP IS NOT A MULTIPLE OF 8 C (TO AVOID MEMORY BANK CONFLICTS). FOR CLARIFICATION OF C INC AND JUMP, SEE THE EXAMPLES BELOW. C C N C THE LENGTH OF EACH TRANSFORM (SEE DEFINITION OF C TRANSFORMS, BELOW). C C M C THE NUMBER OF TRANSFORMS TO BE DONE SIMULTANEOUSLY. C C ISIGN C = +1 FOR A TRANSFORM FROM FOURIER COEFFICIENTS TO C GRIDPOINT VALUES. C = -1 FOR A TRANSFORM FROM GRIDPOINT VALUES TO FOURIER C COEFFICIENTS. C C ON OUTPUT A C IF ISIGN = +1, AND M COEFFICIENT VECTORS ARE SUPPLIED C EACH CONTAINING THE SEQUENCE: C C A(0),B(0),A(1),B(1),...,A(N/2),B(N/2) (N+2 VALUES) C C THEN THE RESULT CONSISTS OF M DATA VECTORS EACH C CONTAINING THE CORRESPONDING N+2 GRIDPOINT VALUES: C C FOR FFT991, X(0), X(1), X(2),...,X(N-1),0,0. C FOR FFT99, X(N-1),X(0),X(1),X(2),...,X(N-1),X(0). C (EXPLICIT CYCLIC CONTINUITY) C C WHEN ISIGN = +1, THE TRANSFORM IS DEFINED BY: C X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N)) C WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K) C AND I=SQRT (-1) C C IF ISIGN = -1, AND M DATA VECTORS ARE SUPPLIED EACH C CONTAINING A SEQUENCE OF GRIDPOINT VALUES X(J) AS C DEFINED ABOVE, THEN THE RESULT CONSISTS OF M VECTORS C EACH CONTAINING THE CORRESPONDING FOURIER COFFICIENTS C A(K), B(K), 0 .LE. K .LE N/2. C C WHEN ISIGN = -1, THE INVERSE TRANSFORM IS DEFINED BY: C C(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*EXP(-2*I*J*K*PI/N)) C WHERE C(K)=A(K)+I*B(K) AND I=SQRT(-1) C C A CALL WITH ISIGN=+1 FOLLOWED BY A CALL WITH ISIGN=-1 C (OR VICE VERSA) RETURNS THE ORIGINAL DATA. C C NOTE: THE FACT THAT THE GRIDPOINT VALUES X(J) ARE REAL C IMPLIES THAT B(0)=B(N/2)=0. FOR A CALL WITH ISIGN=+1, C IT IS NOT ACTUALLY NECESSARY TO SUPPLY THESE ZEROS. C C EXAMPLES GIVEN 19 DATA VECTORS EACH OF LENGTH 64 (+2 FOR EXPLICIT C CYCLIC CONTINUITY), COMPUTE THE CORRESPONDING VECTORS OF C FOURIER COEFFICIENTS. THE DATA MAY, FOR EXAMPLE, BE C ARRANGED LIKE THIS: C C FIRST DATA A(1)= . . . A(66)= A(70) C VECTOR X(63) X(0) X(1) X(2) ... X(63) X(0) (4 EMPTY LOCATIONS) C C SECOND DATA A(71)= . . . A(140) C VECTOR X(63) X(0) X(1) X(2) ... X(63) X(0) (4 EMPTY LOCATIONS) C C AND SO ON. HERE INC=1, JUMP=70, N=64, M=19, ISIGN=-1, C AND FFT99 SHOULD BE USED (BECAUSE OF THE EXPLICIT CYCLIC C CONTINUITY). C C ALTERNATIVELY THE DATA MAY BE ARRANGED LIKE THIS: C C FIRST SECOND LAST C DATA DATA DATA C VECTOR VECTOR VECTOR C C A(1)= A(2)= A(19)= C C X(63) X(63) . . . X(63) C A(20)= X(0) X(0) . . . X(0) C A(39)= X(1) X(1) . . . X(1) C . . . C . . . C . . . C C IN WHICH CASE WE HAVE INC=19, JUMP=1, AND THE REMAINING C PARAMETERS ARE THE SAME AS BEFORE. IN EITHER CASE, EACH C COEFFICIENT VECTOR OVERWRITES THE CORRESPONDING INPUT C DATA VECTOR. C C----------------------------------------------------------------------- DIMENSION A(N),WORK(N),TRIGS(N),IFAX(1) C C SUBROUTINE "FFT99" - MULTIPLE FAST REAL PERIODIC TRANSFORM C CORRESPONDING TO OLD SCALAR ROUTINE FFT9 C PROCEDURE USED TO CONVERT TO HALF-LENGTH COMPLEX TRANSFORM C IS GIVEN BY COOLEY, LEWIS AND WELCH (J. SOUND VIB., VOL. 12 C (1970), 315-337) C C A IS THE ARRAY CONTAINING INPUT AND OUTPUT DATA C WORK IS AN AREA OF SIZE (N+1)*LOT C TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES C IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N/2 C INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR' C (E.G. INC=1 FOR CONSECUTIVELY STORED DATA) C JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR C N IS THE LENGTH OF THE DATA VECTORS C LOT IS THE NUMBER OF DATA VECTORS C ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT C = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL C C ORDERING OF COEFFICIENTS: C A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2) C WHERE B(0)=B(N/2)=0; (N+2) LOCATIONS REQUIRED C C ORDERING OF DATA: C X(N-1),X(0),X(1),X(2),...,X(N),X(0) C I.E. EXPLICIT CYCLIC CONTINUITY; (N+2) LOCATIONS REQUIRED C C VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS IN C PARALLEL C C *** N.B. N IS ASSUMED TO BE AN EVEN NUMBER C C DEFINITION OF TRANSFORMS: C ------------------------- C C ISIGN=+1: X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N)) C WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K) C C ISIGN=-1: A(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N)) C B(K)=-(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N)) C C C C NFAX=IFAX(1) NX=N+1 NH=N/2 INK=INC+INC IF (ISIGN.EQ.+1) GO TO 30 C C IF NECESSARY, TRANSFER DATA TO WORK AREA IGO=50 IF (MOD(NFAX,2).EQ.1) GOTO 40 IBASE=INC+1 JBASE=1 DO 20 L=1,LOT I=IBASE J=JBASE CDIR$ IVDEP DO 10 M=1,N WORK(J)=A(I) I=I+INC J=J+1 10 CONTINUE IBASE=IBASE+JUMP JBASE=JBASE+NX 20 CONTINUE C IGO=60 GO TO 40 C C PREPROCESSING (ISIGN=+1) C ------------------------ C 30 CONTINUE CALL FFT99A(A,WORK,TRIGS,INC,JUMP,N,LOT) IGO=60 C C COMPLEX TRANSFORM C ----------------- C 40 CONTINUE IA=INC+1 LA=1 DO 80 K=1,NFAX IF (IGO.EQ.60) GO TO 60 50 CONTINUE CALL VPASSM(A(IA),A(IA+INC),WORK(1),WORK(2),TRIGS, * INK,2,JUMP,NX,LOT,NH,IFAX(K+1),LA) IGO=60 GO TO 70 60 CONTINUE CALL VPASSM(WORK(1),WORK(2),A(IA),A(IA+INC),TRIGS, * 2,INK,NX,JUMP,LOT,NH,IFAX(K+1),LA) IGO=50 70 CONTINUE LA=LA*IFAX(K+1) 80 CONTINUE C IF (ISIGN.EQ.-1) GO TO 130 C C IF NECESSARY, TRANSFER DATA FROM WORK AREA IF (MOD(NFAX,2).EQ.1) GO TO 110 IBASE=1 JBASE=IA DO 100 L=1,LOT I=IBASE J=JBASE CDIR$ IVDEP DO 90 M=1,N A(J)=WORK(I) I=I+1 J=J+INC 90 CONTINUE IBASE=IBASE+NX JBASE=JBASE+JUMP 100 CONTINUE C C FILL IN CYCLIC BOUNDARY POINTS 110 CONTINUE IA=1 IB=N*INC+1 CDIR$ IVDEP DO 120 L=1,LOT A(IA)=A(IB) A(IB+INC)=A(IA+INC) IA=IA+JUMP IB=IB+JUMP 120 CONTINUE GO TO 140 C C POSTPROCESSING (ISIGN=-1): C -------------------------- C 130 CONTINUE CALL FFT99B(WORK,A,TRIGS,INC,JUMP,N,LOT) C 140 CONTINUE RETURN END SUBROUTINE FFT99A(A,WORK,TRIGS,INC,JUMP,N,LOT) DIMENSION A(N),WORK(N),TRIGS(N) C C SUBROUTINE FFT99A - PREPROCESSING STEP FOR FFT99, ISIGN=+1 C (SPECTRAL TO GRIDPOINT TRANSFORM) C NH=N/2 NX=N+1 INK=INC+INC C C A(0) AND A(N/2) IA=1 IB=N*INC+1 JA=1 JB=2 CDIR$ IVDEP DO 10 L=1,LOT WORK(JA)=A(IA)+A(IB) WORK(JB)=A(IA)-A(IB) IA=IA+JUMP IB=IB+JUMP JA=JA+NX JB=JB+NX 10 CONTINUE C C REMAINING WAVENUMBERS IABASE=2*INC+1 IBBASE=(N-2)*INC+1 JABASE=3 JBBASE=N-1 C DO 30 K=3,NH,2 IA=IABASE IB=IBBASE JA=JABASE JB=JBBASE C=TRIGS(N+K) S=TRIGS(N+K+1) CDIR$ IVDEP DO 20 L=1,LOT WORK(JA)=(A(IA)+A(IB))- * (S*(A(IA)-A(IB))+C*(A(IA+INC)+A(IB+INC))) WORK(JB)=(A(IA)+A(IB))+ * (S*(A(IA)-A(IB))+C*(A(IA+INC)+A(IB+INC))) WORK(JA+1)=(C*(A(IA)-A(IB))-S*(A(IA+INC)+A(IB+INC)))+ * (A(IA+INC)-A(IB+INC)) WORK(JB+1)=(C*(A(IA)-A(IB))-S*(A(IA+INC)+A(IB+INC)))- * (A(IA+INC)-A(IB+INC)) IA=IA+JUMP IB=IB+JUMP JA=JA+NX JB=JB+NX 20 CONTINUE IABASE=IABASE+INK IBBASE=IBBASE-INK JABASE=JABASE+2 JBBASE=JBBASE-2 30 CONTINUE C IF (IABASE.NE.IBBASE) GO TO 50 C WAVENUMBER N/4 (IF IT EXISTS) IA=IABASE JA=JABASE CDIR$ IVDEP DO 40 L=1,LOT WORK(JA)=2.0*A(IA) WORK(JA+1)=-2.0*A(IA+INC) IA=IA+JUMP JA=JA+NX 40 CONTINUE C 50 CONTINUE RETURN END SUBROUTINE FFT99B(WORK,A,TRIGS,INC,JUMP,N,LOT) DIMENSION WORK(N),A(N),TRIGS(N) C C SUBROUTINE FFT99B - POSTPROCESSING STEP FOR FFT99, ISIGN=-1 C (GRIDPOINT TO SPECTRAL TRANSFORM) C NH=N/2 NX=N+1 INK=INC+INC C C A(0) AND A(N/2) SCALE=1.0/FLOAT(N) IA=1 IB=2 JA=1 JB=N*INC+1 CDIR$ IVDEP DO 10 L=1,LOT A(JA)=SCALE*(WORK(IA)+WORK(IB)) A(JB)=SCALE*(WORK(IA)-WORK(IB)) A(JA+INC)=0.0 A(JB+INC)=0.0 IA=IA+NX IB=IB+NX JA=JA+JUMP JB=JB+JUMP 10 CONTINUE C C REMAINING WAVENUMBERS SCALE=0.5*SCALE IABASE=3 IBBASE=N-1 JABASE=2*INC+1 JBBASE=(N-2)*INC+1 C DO 30 K=3,NH,2 IA=IABASE IB=IBBASE JA=JABASE JB=JBBASE C=TRIGS(N+K) S=TRIGS(N+K+1) CDIR$ IVDEP DO 20 L=1,LOT A(JA)=SCALE*((WORK(IA)+WORK(IB)) * +(C*(WORK(IA+1)+WORK(IB+1))+S*(WORK(IA)-WORK(IB)))) A(JB)=SCALE*((WORK(IA)+WORK(IB)) * -(C*(WORK(IA+1)+WORK(IB+1))+S*(WORK(IA)-WORK(IB)))) A(JA+INC)=SCALE*((C*(WORK(IA)-WORK(IB))-S*(WORK(IA+1)+WORK(IB+1))) * +(WORK(IB+1)-WORK(IA+1))) A(JB+INC)=SCALE*((C*(WORK(IA)-WORK(IB))-S*(WORK(IA+1)+WORK(IB+1))) * -(WORK(IB+1)-WORK(IA+1))) IA=IA+NX IB=IB+NX JA=JA+JUMP JB=JB+JUMP 20 CONTINUE IABASE=IABASE+2 IBBASE=IBBASE-2 JABASE=JABASE+INK JBBASE=JBBASE-INK 30 CONTINUE C IF (IABASE.NE.IBBASE) GO TO 50 C WAVENUMBER N/4 (IF IT EXISTS) IA=IABASE JA=JABASE SCALE=2.0*SCALE CDIR$ IVDEP DO 40 L=1,LOT A(JA)=SCALE*WORK(IA) A(JA+INC)=-SCALE*WORK(IA+1) IA=IA+NX JA=JA+JUMP 40 CONTINUE C 50 CONTINUE RETURN END SUBROUTINE FFT991(A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN) DIMENSION A(N),WORK(N),TRIGS(N),IFAX(1) C C SUBROUTINE "FFT991" - MULTIPLE REAL/HALF-COMPLEX PERIODIC C FAST FOURIER TRANSFORM C C SAME AS FFT99 EXCEPT THAT ORDERING OF DATA CORRESPONDS TO C THAT IN MRFFT2 C C PROCEDURE USED TO CONVERT TO HALF-LENGTH COMPLEX TRANSFORM C IS GIVEN BY COOLEY, LEWIS AND WELCH (J. SOUND VIB., VOL. 12 C (1970), 315-337) C C A IS THE ARRAY CONTAINING INPUT AND OUTPUT DATA C WORK IS AN AREA OF SIZE (N+1)*LOT C TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES C IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N/2 C INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR' C (E.G. INC=1 FOR CONSECUTIVELY STORED DATA) C JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR C N IS THE LENGTH OF THE DATA VECTORS C LOT IS THE NUMBER OF DATA VECTORS C ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT C = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL C C ORDERING OF COEFFICIENTS: C A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2) C WHERE B(0)=B(N/2)=0; (N+2) LOCATIONS REQUIRED C C ORDERING OF DATA: C X(0),X(1),X(2),...,X(N-1) C C VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS IN C PARALLEL C C *** N.B. N IS ASSUMED TO BE AN EVEN NUMBER C C DEFINITION OF TRANSFORMS: C ------------------------- C C ISIGN=+1: X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N)) C WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K) C C ISIGN=-1: A(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N)) C B(K)=-(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N)) C C C NFAX=IFAX(1) NX=N+1 NH=N/2 INK=INC+INC IF (ISIGN.EQ.+1) GO TO 30 C C IF NECESSARY, TRANSFER DATA TO WORK AREA IGO=50 IF (MOD(NFAX,2).EQ.1) GOTO 40 IBASE=1 JBASE=1 DO 20 L=1,LOT I=IBASE J=JBASE CDIR$ IVDEP DO 10 M=1,N WORK(J)=A(I) I=I+INC J=J+1 10 CONTINUE IBASE=IBASE+JUMP JBASE=JBASE+NX 20 CONTINUE C IGO=60 GO TO 40 C C PREPROCESSING (ISIGN=+1) C ------------------------ C 30 CONTINUE CALL FFT99A(A,WORK,TRIGS,INC,JUMP,N,LOT) IGO=60 C C COMPLEX TRANSFORM C ----------------- C 40 CONTINUE IA=1 LA=1 DO 80 K=1,NFAX IF (IGO.EQ.60) GO TO 60 50 CONTINUE CALL VPASSM(A(IA),A(IA+INC),WORK(1),WORK(2),TRIGS, * INK,2,JUMP,NX,LOT,NH,IFAX(K+1),LA) IGO=60 GO TO 70 60 CONTINUE CALL VPASSM(WORK(1),WORK(2),A(IA),A(IA+INC),TRIGS, * 2,INK,NX,JUMP,LOT,NH,IFAX(K+1),LA) IGO=50 70 CONTINUE LA=LA*IFAX(K+1) 80 CONTINUE C IF (ISIGN.EQ.-1) GO TO 130 C C IF NECESSARY, TRANSFER DATA FROM WORK AREA IF (MOD(NFAX,2).EQ.1) GO TO 110 IBASE=1 JBASE=1 DO 100 L=1,LOT I=IBASE J=JBASE CDIR$ IVDEP DO 90 M=1,N A(J)=WORK(I) I=I+1 J=J+INC 90 CONTINUE IBASE=IBASE+NX JBASE=JBASE+JUMP 100 CONTINUE C C FILL IN ZEROS AT END 110 CONTINUE IB=N*INC+1 CDIR$ IVDEP DO 120 L=1,LOT A(IB)=0.0 A(IB+INC)=0.0 IB=IB+JUMP 120 CONTINUE GO TO 140 C C POSTPROCESSING (ISIGN=-1): C -------------------------- C 130 CONTINUE CALL FFT99B(WORK,A,TRIGS,INC,JUMP,N,LOT) C 140 CONTINUE RETURN END SUBROUTINE SET99 (TRIGS, IFAX, N) DIMENSION IFAX(13),TRIGS(1) C C MODE 3 IS USED FOR REAL/HALF-COMPLEX TRANSFORMS. IT IS POSSIBLE C TO DO COMPLEX/COMPLEX TRANSFORMS WITH OTHER VALUES OF MODE, BUT C DOCUMENTATION OF THE DETAILS WERE NOT AVAILABLE WHEN THIS ROUTINE C WAS WRITTEN. C DATA MODE /3/ CALL FAX (IFAX, N, MODE) I = IFAX(1) IF (IFAX(I+1) .GT. 5 .OR. N .LE. 4) IFAX(1) = -99 IF (IFAX(1) .LE. 0 ) THEN WRITE(6,*) ' SET99 -- INVALID N' STOP'SET99' ENDIF CALL FFTRIG (TRIGS, N, MODE) RETURN END SUBROUTINE FAX(IFAX,N,MODE) DIMENSION IFAX(10) NN=N IF (IABS(MODE).EQ.1) GO TO 10 IF (IABS(MODE).EQ.8) GO TO 10 NN=N/2 IF ((NN+NN).EQ.N) GO TO 10 IFAX(1)=-99 RETURN 10 K=1 C TEST FOR FACTORS OF 4 20 IF (MOD(NN,4).NE.0) GO TO 30 K=K+1 IFAX(K)=4 NN=NN/4 IF (NN.EQ.1) GO TO 80 GO TO 20 C TEST FOR EXTRA FACTOR OF 2 30 IF (MOD(NN,2).NE.0) GO TO 40 K=K+1 IFAX(K)=2 NN=NN/2 IF (NN.EQ.1) GO TO 80 C TEST FOR FACTORS OF 3 40 IF (MOD(NN,3).NE.0) GO TO 50 K=K+1 IFAX(K)=3 NN=NN/3 IF (NN.EQ.1) GO TO 80 GO TO 40 C NOW FIND REMAINING FACTORS 50 L=5 INC=2 C INC ALTERNATELY TAKES ON VALUES 2 AND 4 60 IF (MOD(NN,L).NE.0) GO TO 70 K=K+1 IFAX(K)=L NN=NN/L IF (NN.EQ.1) GO TO 80 GO TO 60 70 L=L+INC INC=6-INC GO TO 60 80 IFAX(1)=K-1 C IFAX(1) CONTAINS NUMBER OF FACTORS NFAX=IFAX(1) C SORT FACTORS INTO ASCENDING ORDER IF (NFAX.EQ.1) GO TO 110 DO 100 II=2,NFAX ISTOP=NFAX+2-II DO 90 I=2,ISTOP IF (IFAX(I+1).GE.IFAX(I)) GO TO 90 ITEM=IFAX(I) IFAX(I)=IFAX(I+1) IFAX(I+1)=ITEM 90 CONTINUE 100 CONTINUE 110 CONTINUE RETURN END SUBROUTINE FFTRIG(TRIGS,N,MODE) DIMENSION TRIGS(1) PI=2.0*ASIN(1.0) IMODE=IABS(MODE) NN=N IF (IMODE.GT.1.AND.IMODE.LT.6) NN=N/2 DEL=(PI+PI)/FLOAT(NN) L=NN+NN DO 10 I=1,L,2 ANGLE=0.5*FLOAT(I-1)*DEL TRIGS(I)=COS(ANGLE) TRIGS(I+1)=SIN(ANGLE) 10 CONTINUE IF (IMODE.EQ.1) RETURN IF (IMODE.EQ.8) RETURN DEL=0.5*DEL NH=(NN+1)/2 L=NH+NH LA=NN+NN DO 20 I=1,L,2 ANGLE=0.5*FLOAT(I-1)*DEL TRIGS(LA+I)=COS(ANGLE) TRIGS(LA+I+1)=SIN(ANGLE) 20 CONTINUE IF (IMODE.LE.3) RETURN DEL=0.5*DEL LA=LA+NN IF (MODE.EQ.5) GO TO 40 DO 30 I=2,NN ANGLE=FLOAT(I-1)*DEL TRIGS(LA+I)=2.0*SIN(ANGLE) 30 CONTINUE RETURN 40 CONTINUE DEL=0.5*DEL DO 50 I=2,N ANGLE=FLOAT(I-1)*DEL TRIGS(LA+I)=SIN(ANGLE) 50 CONTINUE RETURN END SUBROUTINE VPASSM(A,B,C,D,TRIGS,INC1,INC2,INC3,INC4,LOT,N,IFAC,LA) DIMENSION A(N),B(N),C(N),D(N),TRIGS(N) C C SUBROUTINE "VPASSM" - MULTIPLE VERSION OF "VPASSA" C PERFORMS ONE PASS THROUGH DATA C AS PART OF MULTIPLE COMPLEX FFT ROUTINE C A IS FIRST REAL INPUT VECTOR C B IS FIRST IMAGINARY INPUT VECTOR C C IS FIRST REAL OUTPUT VECTOR C D IS FIRST IMAGINARY OUTPUT VECTOR C TRIGS IS PRECALCULATED TABLE OF SINES " COSINES C INC1 IS ADDRESSING INCREMENT FOR A AND B C INC2 IS ADDRESSING INCREMENT FOR C AND D C INC3 IS ADDRESSING INCREMENT BETWEEN A"S & B"S C INC4 IS ADDRESSING INCREMENT BETWEEN C"S & D"S C LOT IS THE NUMBER OF VECTORS C N IS LENGTH OF VECTORS C IFAC IS CURRENT FACTOR OF N C LA IS PRODUCT OF PREVIOUS FACTORS C DATA SIN36/0.587785252292473/,COS36/0.809016994374947/, * SIN72/0.951056516295154/,COS72/0.309016994374947/, * SIN60/0.866025403784437/ C M=N/IFAC IINK=M*INC1 JINK=LA*INC2 JUMP=(IFAC-1)*JINK IBASE=0 JBASE=0 IGO=IFAC-1 IF (IGO.GT.4) RETURN GO TO (10,50,90,130),IGO C C CODING FOR FACTOR 2 C 10 IA=1 JA=1 IB=IA+IINK JB=JA+JINK DO 20 L=1,LA I=IBASE J=JBASE CDIR$ IVDEP DO 15 IJK=1,LOT C(JA+J)=A(IA+I)+A(IB+I) D(JA+J)=B(IA+I)+B(IB+I) C(JB+J)=A(IA+I)-A(IB+I) D(JB+J)=B(IA+I)-B(IB+I) I=I+INC3 J=J+INC4 15 CONTINUE IBASE=IBASE+INC1 JBASE=JBASE+INC2 20 CONTINUE IF (LA.EQ.M) RETURN LA1=LA+1 JBASE=JBASE+JUMP DO 40 K=LA1,M,LA KB=K+K-2 C1=TRIGS(KB+1) S1=TRIGS(KB+2) DO 30 L=1,LA I=IBASE J=JBASE CDIR$ IVDEP DO 25 IJK=1,LOT C(JA+J)=A(IA+I)+A(IB+I) D(JA+J)=B(IA+I)+B(IB+I) C(JB+J)=C1*(A(IA+I)-A(IB+I))-S1*(B(IA+I)-B(IB+I)) D(JB+J)=S1*(A(IA+I)-A(IB+I))+C1*(B(IA+I)-B(IB+I)) I=I+INC3 J=J+INC4 25 CONTINUE IBASE=IBASE+INC1 JBASE=JBASE+INC2 30 CONTINUE JBASE=JBASE+JUMP 40 CONTINUE RETURN C C CODING FOR FACTOR 3 C 50 IA=1 JA=1 IB=IA+IINK JB=JA+JINK IC=IB+IINK JC=JB+JINK DO 60 L=1,LA I=IBASE J=JBASE CDIR$ IVDEP DO 55 IJK=1,LOT C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I)) D(JA+J)=B(IA+I)+(B(IB+I)+B(IC+I)) C(JB+J)=(A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I))) C(JC+J)=(A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I))) D(JB+J)=(B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I))) D(JC+J)=(B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I))) I=I+INC3 J=J+INC4 55 CONTINUE IBASE=IBASE+INC1 JBASE=JBASE+INC2 60 CONTINUE IF (LA.EQ.M) RETURN LA1=LA+1 JBASE=JBASE+JUMP DO 80 K=LA1,M,LA KB=K+K-2 KC=KB+KB C1=TRIGS(KB+1) S1=TRIGS(KB+2) C2=TRIGS(KC+1) S2=TRIGS(KC+2) DO 70 L=1,LA I=IBASE J=JBASE CDIR$ IVDEP DO 65 IJK=1,LOT C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I)) D(JA+J)=B(IA+I)+(B(IB+I)+B(IC+I)) C(JB+J)= * C1*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I)))) * -S1*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I)))) D(JB+J)= * S1*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I)))) * +C1*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I)))) C(JC+J)= * C2*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I)))) * -S2*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I)))) D(JC+J)= * S2*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I)))) * +C2*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I)))) I=I+INC3 J=J+INC4 65 CONTINUE IBASE=IBASE+INC1 JBASE=JBASE+INC2 70 CONTINUE JBASE=JBASE+JUMP 80 CONTINUE RETURN C C CODING FOR FACTOR 4 C 90 IA=1 JA=1 IB=IA+IINK JB=JA+JINK IC=IB+IINK JC=JB+JINK ID=IC+IINK JD=JC+JINK DO 100 L=1,LA I=IBASE J=JBASE CDIR$ IVDEP DO 95 IJK=1,LOT C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I)) C(JC+J)=(A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)) D(JA+J)=(B(IA+I)+B(IC+I))+(B(IB+I)+B(ID+I)) D(JC+J)=(B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I)) C(JB+J)=(A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I)) C(JD+J)=(A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I)) D(JB+J)=(B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I)) D(JD+J)=(B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I)) I=I+INC3 J=J+INC4 95 CONTINUE IBASE=IBASE+INC1 JBASE=JBASE+INC2 100 CONTINUE IF (LA.EQ.M) RETURN LA1=LA+1 JBASE=JBASE+JUMP DO 120 K=LA1,M,LA KB=K+K-2 KC=KB+KB KD=KC+KB C1=TRIGS(KB+1) S1=TRIGS(KB+2) C2=TRIGS(KC+1) S2=TRIGS(KC+2) C3=TRIGS(KD+1) S3=TRIGS(KD+2) DO 110 L=1,LA I=IBASE J=JBASE CDIR$ IVDEP DO 105 IJK=1,LOT C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I)) D(JA+J)=(B(IA+I)+B(IC+I))+(B(IB+I)+B(ID+I)) C(JC+J)= * C2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))) * -S2*((B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I))) D(JC+J)= * S2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))) * +C2*((B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I))) C(JB+J)= * C1*((A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I))) * -S1*((B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I))) D(JB+J)= * S1*((A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I))) * +C1*((B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I))) C(JD+J)= * C3*((A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I))) * -S3*((B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I))) D(JD+J)= * S3*((A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I))) * +C3*((B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I))) I=I+INC3 J=J+INC4 105 CONTINUE IBASE=IBASE+INC1 JBASE=JBASE+INC2 110 CONTINUE JBASE=JBASE+JUMP 120 CONTINUE RETURN C C CODING FOR FACTOR 5 C 130 IA=1 JA=1 IB=IA+IINK JB=JA+JINK IC=IB+IINK JC=JB+JINK ID=IC+IINK JD=JC+JINK IE=ID+IINK JE=JD+JINK DO 140 L=1,LA I=IBASE J=JBASE CDIR$ IVDEP DO 135 IJK=1,LOT C(JA+J)=A(IA+I)+(A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I)) D(JA+J)=B(IA+I)+(B(IB+I)+B(IE+I))+(B(IC+I)+B(ID+I)) C(JB+J)=(A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) * -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))) C(JE+J)=(A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) * +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))) D(JB+J)=(B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) * +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))) D(JE+J)=(B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) * -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))) C(JC+J)=(A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) * -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))) C(JD+J)=(A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) * +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))) D(JC+J)=(B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) * +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))) D(JD+J)=(B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) * -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))) I=I+INC3 J=J+INC4 135 CONTINUE IBASE=IBASE+INC1 JBASE=JBASE+INC2 140 CONTINUE IF (LA.EQ.M) RETURN LA1=LA+1 JBASE=JBASE+JUMP DO 160 K=LA1,M,LA KB=K+K-2 KC=KB+KB KD=KC+KB KE=KD+KB C1=TRIGS(KB+1) S1=TRIGS(KB+2) C2=TRIGS(KC+1) S2=TRIGS(KC+2) C3=TRIGS(KD+1) S3=TRIGS(KD+2) C4=TRIGS(KE+1) S4=TRIGS(KE+2) DO 150 L=1,LA I=IBASE J=JBASE CDIR$ IVDEP DO 145 IJK=1,LOT C(JA+J)=A(IA+I)+(A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I)) D(JA+J)=B(IA+I)+(B(IB+I)+B(IE+I))+(B(IC+I)+B(ID+I)) C(JB+J)= * C1*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) * -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))) * -S1*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) * +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))) D(JB+J)= * S1*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) * -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))) * +C1*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) * +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))) C(JE+J)= * C4*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) * +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))) * -S4*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) * -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))) D(JE+J)= * S4*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) * +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))) * +C4*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) * -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))) C(JC+J)= * C2*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) * -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))) * -S2*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) * +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))) D(JC+J)= * S2*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) * -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))) * +C2*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) * +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))) C(JD+J)= * C3*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) * +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))) * -S3*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) * -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))) D(JD+J)= * S3*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) * +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))) * +C3*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) * -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))) I=I+INC3 J=J+INC4 145 CONTINUE IBASE=IBASE+INC1 JBASE=JBASE+INC2 150 CONTINUE JBASE=JBASE+JUMP 160 CONTINUE RETURN END c subroutine spharm(ugrid,anm,am, * pnm,weights,nmodes, * a,c,ifax,trig, * nlons,nlats,mwaves,nmdim,idir) c REAL UGRID(NLONS,nlats),pnm(nmdim,nlats),weights(nlats) COMPLEX ANM(nmdim), AM(MWAVES,nlats) integer ifax(13),nmodes(mwaves) real trig((3*nlons/2)+1), * a(nlats*(nlons+2)),c(nlats*(nlons+1)) c IF (IDIR .eq. +1) then c C==> GRID SPACE TO SPECTRAL SPACE TRANSFORMATION C FIRST, INITIALIZE ARRAY. C DO 100 NM = 1, NMDIM ANM(NM) = CMPLX(0.0,0.0) 100 CONTINUE c c==> perform ffts on each latitude. C call rfft(ugrid,am, * a,c,ifax,trig, * nlons,nlats,mwaves,1) C C==> SUM OVER ALL GAUSSIAN LATITUDES FOR EACH MODE AND EACH WAVE TO C OBTAIN THE TRANSFORMED VARIABLE IN SPECTRAL SPACE. C do 300 j=1,nlats NMSTRT = 0 DO 300 m = 1, MWAVES DO 200 n = 1, NModes(m) NM = NMSTRT + n anm(NM)=anm(NM)+pnm(nm,j)*weights(j)*am(m,j) 200 CONTINUE nmstrt = nmstrt + nmodes(m) 300 continue c C==> SPECTRAL SPACE TO GRID SPACE TRANSFORMATION. C else if (idir .eq. -1) then c DO 700 J = 1, NLATS C C==> INVERSE LEGENDRE TRANSFORM TO GET VALUES OF THE ZONAL FOURIER C TRANSFORM AT LATITUDE j. C C==> SUM THE VARIOUS MERIDIONAL MODES TO BUILD THE FOURIER SERIES C COEFFICIENT FOR ZONAL WAVENUMBER M=I-1 AT the GIVEN LATITUDE. C NMSTRT = 0 dO 700 m = 1, MWAVES am(m,j) = cmplx(0., 0.) DO 500 n = 1, NModes(m) NM = NMSTRT + n am(m,j) = am(m,j) + anm(NM) * Pnm(NM,j) 500 CONTINUE NMSTRT = NMSTRT + NModes(m) 700 CONTINUE C C==> FOURIER TRANSFORM TO COMPUTE THE VALUES OF THE VARIABLE IN GRID C SPACE at THE J-TH LATITUDE. C call rfft(UGRID,AM, * a,c,ifax,trig, * nlons,nlats,mwaves,-1) c else write(6,*) 'error in spharm: idir must be -1 or +1!' write(6,*) 'execution terminated in subroutine spharm' end if c RETURN END c subroutine rfft(data,coeff, * a,c,ifax,trig, * nlons,nlats,mwaves,idir) c real data(nlons,nlats) complex coeff(mwaves,nlats) c real trig((3*nlons/2)+1), * a(nlats*(nlons+2)),c(nlats*(nlons+1)) integer ifax(13) c c---------------------- c==> forward transform. c---------------------- if (idir .eq. +1) then c c==> copy the data into the work array. c transforms are computed pairwise using a complex fft. c n = 0 do 10 j=1,nlats do 10 i=1,nlons+2 n = n + 1 a(n) = 0.0 if (i .le. nlons) then a(n) = data(i,j) end if 10 continue c call fft991(a,c,trig,ifax,1,nlons+2,nlons,nlats,-1) c n = -1 do 15 j=1,nlats do 15 m=1,(nlons/2)+1 n = n + 2 if (m .le. mwaves) then coeff(m,j) = cmplx(a(n),a(n+1)) end if 15 continue c---------------------- c==> inverse transform. c---------------------- else if (idir .eq. -1) then c do n = 1,nlats*(nlons+2) a(n) = 0.0 enddo n = -1 do 30 j=1,nlats do 30 m=1,(nlons/2)+1 n = n + 2 if (m .le. mwaves) then a(n) = real(coeff(m,j)) a(n+1) = aimag(coeff(m,j)) end if 30 continue c call fft991(a,c,trig,ifax,1,nlons+2,nlons,nlats,1) c n = 0 do 40 j=1,nlats do 40 i=1,nlons+2 n = n + 1 if (i .le. nlons) then data(i,j) = a(n) end if 40 continue c else write(6,*) ' idir must be +1 or -1 in RFFT!' write(6,*) ' execution terminated.' end if c return end c subroutine getuv(vrtnm,divnm,um,vm,ug,vg,re, * nlons,nlats,mwaves,nmdim, * aa,cc,ifax,trig, * pnm,hnm,deli,nmodes) c real ug(nlons,nlats),vg(nlons,nlats) complex vrtnm(nmdim),divnm(nmdim), * um(mwaves,nlats),vm(mwaves,nlats) c real trig((3*nlons/2)+1), * aa(nlats*(nlons+2)),cc(nlats*(nlons+1)) integer ifax(13) c real pnm(nmdim,nlats), * hnm(nmdim,nlats), * deli(nmdim) integer nmodes(mwaves) c do 200 j=1,nlats nmstrt = 0 c do 200 m=1,mwaves rm = m-1 um(m,j) = cmplx(0.,0.) vm(m,j) = cmplx(0.,0.) c do 100 n=1,nmodes(m) nm = nmstrt + n um(m,j) = um(m,j) - re*deli(nm)*( * cmplx(0.,rm)*divnm(nm)*pnm(nm,j) + * vrtnm(nm)*hnm(nm,j) ) vm(m,j) = vm(m,j) - re*deli(nm)*( * cmplx(0.,rm)*vrtnm(nm)*pnm(nm,j) - * divnm(nm)*hnm(nm,j) ) 100 continue c nmstrt = nmstrt + nmodes(m) c 200 continue c call rfft(ug,um, * aa,cc,ifax,trig, * nlons,nlats,mwaves,-1) call rfft(vg,vm, * aa,cc,ifax,trig, * nlons,nlats,mwaves,-1) c return end