program read_ds0840_sun c c c SUN FORTRAN version c parameter (mxsize=512,mxrc=56,mx31=31,mx32=32) parameter (idim=144,jdim=73) parameter (maxz=mx31*mx31,maxuv=mx32*mx31) parameter (maxz2=2*maxz, maxuv2=2*maxuv) parameter (maxuu2=2*mx32*mx32) integer*4 strmopen, strmclose, strmpos integer*4 infunit, iclose, idum external strmopen !$pragma C ( strmopen) external strmclose !$pragma C ( strmclose) external strmread !$pragma C ( strmread) external strmpos !$pragma C ( strmpos) character one*32767, outstr*80 c character nmclabel*32 real*4 fhour,wave1(maxz2), wave2(maxuu2),gd1(idim,jdim) real*8 sinclt(jdim) real*8 grdv(idim*jdim),agrd integer*4 iyear,imonth,iday,ihour c integer*4 rcsize(mxrc),p(mxrc) c data rcsize/ 32, 24,7688,7936,7936,7688,7936,7936,7688,7936, & 7936,7688,7936,7936,7688,7936,7936,7688,7936,7936, & 7688,7936,7936,7688,7936,7936,7688,7936,7936,7688, & 7936,7936,7688,7936,7936,7688,7936,7936,7688,7688, & 7688,7688,7688,7688,7688,7688,7688,7688,7688,7688, & 7688,7688,7688,7688,7688,7688/ data p/-999,-999,1000,1000,1000,850,850,850,700,700,700, & 500, 500, 500, 400, 400,400,300,300,300,250,250,250, & 200, 200, 200, 150, 150,150,100,100,100, 70, 70, 70, & 50, 50, 50,1000, 850,700,500,400,300, & 1000, 850, 700, 500, 400,300,250,200,150,100, 70, 50/ data er/6371.0E3/,omega/7.292E-5/,idir/-101/,mlt/0/,fac/1.0/ save c iromb=1 eriv=1./er ersq=er*er ersqiv=1./ersq pi=4.*atan(1.) rad=pi/180. twomg=2.*omega do j=1,jdim sinclt(j)=sin(float(j-1)*pi/float(jdim-1)) enddo sinclt(1)=1. sinclt(jdim)=1. c c open the file c infunit=strmopen() c c get one character at a time c nfile=0 c c process one file c lastlen=0 10 continue c c find MRF first c c call strmread(infunit, one, mxsize, lenrd) c if(one(1:3).eq.'EOF') go to 90 c lastlen=lastlen+mxsize c do i=1,mxsize c if(one(i:i+2).eq.'MRF') go to 15 c enddo c go to 10 c 15 continue c idum = strmpos(infunit,lastlen, 0) do nrc=1,mxrc zmax=-99999. zmin=99999. umax=-99999. umin=99999. vmax=-99999. vmin=99999. rhmax=-99999. rhmin=99999. vtmax=-99999. vtmin=99999. c do kk=1,maxz2 wave1(kk)=0.0 enddo do kk=1,maxuu2 wave2(kk)=0.0 enddo if(rcsize(nrc).gt.mxsize) then mread=rcsize(nrc)/mxsize mleft=rcsize(nrc)-mread*mxsize do ijk=1,mread kbeg=(ijk-1)*mxsize+1 kend=(ijk-1)*mxsize+mxsize call strmread(infunit,one(kbeg:kend),mxsize,lenrd) enddo kbeg=kend+1 kend=kbeg+mleft call strmread(infunit,one(kbeg:kend),mleft,lenrd) else call strmread(infunit, one, rcsize(nrc), lenrd) endif c write(*,*) ' nrc = ',nrc if(one(1:3).eq.'EOF') go to 90 lastlen=lastlen+rcsize(nrc) if(nrc.eq.1) then call ebcasc(one,nmclabel,rcsize(1)) write(*,*) nmclabel elseif(nrc.eq.2) then call gbyte(one,idum,0,32) call ibm2real(idum,fhour) call gbyte(one,iyear,32,32) call gbyte(one,imonth,64,32) call gbyte(one,iday,96,32) call gbyte(one,ihour,128,32) write(*,*) fhour,iyear,imonth,iday,ihour elseif(nrc.gt.2 .and. nrc.le.38) then !z,u,v kp=(nrc-3)/3+1 write(*,*) ' kp = ',kp idum=mod(nrc-2,3) if(idum.eq.1) then !z c call nibm2real(one,wave1,maxz2) c write(*,*) ' Z 1-10 ',(wave1(ii),ii=1,maxz2) call unp360(one,wave1,maxz2) c if(nrc.eq.3) then c open(unit=30,file='Z1.wave') c write(30,1233) (wave1(kb),kb=1,maxz2) 1233 format(10E12.4) c close(unit=30) c endif c write(*,*) ' Z 1-10 ',(wave1(ii),ii=1,10) call sphert(idir,grdv,wave1,mlt,fac, & idim,jdim,30,iromb) c write(*,991) ' Z',(grdv(ii),ii=4000,4005) 991 format(1x,a2,6E12.4) do j=1,jdim do i=1,idim agrd=grdv((j-1)*idim+i) gd1(i,j)=agrd if(agrd.gt.zmax) then zmax=agrd kzmx=(j-1)*idim+i endif if(agrd.lt.zmin) then zmin=agrd kzmn=(j-1)*idim+i endif enddo enddo write(*,*) ' Z max at ',zmax,kzmx write(*,*) ' Z min at ',zmin,kzmn c if(p(nrc).eq.500) then c open(unit=23,file='Z.dat',access='direct',recl=144*4) c do j=jdim,1,-1 c write(23,rec=jdim-j+1) (gd1(ka,j),ka=1,idim) c enddo c close(unit=23) c endif elseif(idum.eq.2) then !u call unp360(one,wave2,maxuv2) c if(nrc.eq.4) then c open(unit=30,file='U1.wave') c write(30,1233) (wave2(kb),kb=1,maxuv2) c close(unit=30) c endif write(*,*) 'U wave2 at end ',(wave2(ii),ii=maxuv2-9,maxuv2) call sphert(idir,grdv,wave2,mlt,fac, & idim,jdim,31,iromb) write(*,991) ' U',(grdv(ii),ii=4000,4005) do j=1,jdim do i=1,idim agrd=grdv((j-1)*idim+i)/sinclt(j) ! changed 93oct14 gd1(i,j)=agrd if(agrd.gt.umax) then umax=agrd kumxi=i kumxj=j endif if(agrd.lt.umin) then umin=agrd kumni=i kumnj=j endif enddo enddo write(*,*) ' U max at ',umax,kumxi,kumxj write(*,*) ' U min at ',umin,kumni,kumnj c if(p(nrc).eq.500) then c open(unit=23,file='U.dat',access='direct',recl=144*4) c do j=jdim,1,-1 c write(23,rec=jdim-j+1) (gd1(ka,j),ka=1,idim) c enddo c close(unit=23) c endif else !v call unp360(one,wave2,maxuv2) c if(nrc.eq.5) then c open(unit=30,file='V1.wave') c write(30,1233) (wave2(kb),kb=1,maxuv2) c close(unit=30) c endif c write(*,*) 'V wave2 at end ',(wave2(ii),ii=maxuv2-9,maxuv2) call sphert(idir,grdv,wave2,mlt,fac, & idim,jdim,31,iromb) write(*,991) ' V',(grdv(ii),ii=4000,4005) do j=1,jdim do i=1,idim agrd=grdv((j-1)*idim+i)/sinclt(j) ! changed 93oct14 gd1(i,j)=agrd if(agrd.gt.vmax) then vmax=agrd kvmxi=i kvmxj=j endif if(agrd.lt.vmin) then vmin=agrd kvmni=i kvmnj=j endif enddo enddo write(*,*) ' V max at ',vmax,kvmxi,kvmxj write(*,*) ' V min at ',vmin,kvmni,kvmnj c if(p(nrc).eq.500) then c open(unit=23,file='V.dat',access='direct',recl=144*4) c do j=jdim,1,-1 c write(23,rec=jdim-j+1) (gd1(ka,j),ka=1,idim) c enddo c close(unit=23) c endif endif elseif(nrc.gt.38 .and. nrc.le.44) then ! rh call unp360(one,wave1,maxz2) c if(nrc.eq.39) then c open(unit=30,file='RH1.wave') c write(30,1233) (wave1(kb),kb=1,maxz2) c close(unit=30) c endif write(*,*) 'RH wave1 at end ',(wave1(ii),ii=maxz2-9,maxz2) call sphert(idir,grdv,wave1,mlt,fac, & idim,jdim,30,iromb) write(*,991) 'RH',(grdv(ii),ii=4000,4005) do j=1,jdim do i=1,idim agrd=grdv((j-1)*idim+i) gd1(i,j)=agrd if(agrd.gt.rhmax) then rhmax=agrd krhmxi=i krhmxj=j endif if(agrd.lt.rhmin) then rhmin=agrd krhmni=i krhmnj=j endif enddo enddo write(*,*) ' RH max at ',rhmax,krhmxi,krhmxj write(*,*) ' RH min at ',rhmin,krhmni,krhmnj c if(p(nrc).eq.700) then c open(unit=23,file='RH.dat',access='direct',recl=144*4) c do j=jdim,1,-1 c write(23,rec=jdim-j+1) (gd1(ka,j),ka=1,idim) c enddo c close(unit=23) c endif else ! vt call unp360(one,wave1,maxz2) c if(nrc.eq.45) then c open(unit=30,file='VT1.wave') c write(30,1233) (wave1(kb),kb=1,maxz2) c close(unit=30) c endif write(*,*) 'VT wave1 at end ',(wave1(ii),ii=maxz2-9,maxz2) call sphert(idir,grdv,wave1,mlt,fac, & idim,jdim,30,iromb) write(*,991) 'VT',(grdv(ii),ii=4000,4005) do j=1,jdim do i=1,idim agrd=grdv((j-1)*idim+i) gd1(i,j)=agrd if(agrd.gt.vtmax) then vtmax=agrd kvtmxi=i kvtmxj=j endif if(agrd.lt.vtmin) then vtmin=agrd kvtmni=i kvtmnj=j endif enddo enddo write(*,*) ' VT max at ',vtmax,kvtmxi,kvtmxj write(*,*) ' VT min at ',vtmin,kvtmni,kvtmnj c if(p(nrc).eq.500) then c open(unit=23,file='vt.dat',access='direct',recl=144*4) c do j=jdim,1,-1 c write(23,rec=jdim-j+1) (gd1(ka,j),ka=1,idim) c enddo c close(unit=23) c endif endif c c reposition c idum = strmpos(infunit,lastlen, 0) c write(*,*) ' reposition at ',lastlen enddo c go to 10 c c close the file c 90 continue write(*,*) ' finish reading input file....' iclose=strmclose(infunit) idum=ieee_flags("clear","exception","all",outstr) 9998 stop'Normal Termination' c 9999 write(*,*) '**** ERROR during read, stop' call exit(0) end c c subroutine ibm2real(bit32, rval) c c bit32 is the input 32-bit word c c rval is the real number returned c integer isign, ira, irb real rval, rb character bit32*4 c c get the sign bit c call gbyte(bit32,isign,0,1) c c get the 7-bit characteristic c call gbyte(bit32,ira,1,7) c c get the 24-bit mantissa c call gbyte(bit32,irb,8,24) c rb=irb rval=(-1)**isign * rb * 16**(ira-64) / 2**24 c write(*,*) ' ibm2real: ',isign,ira,irb,rval return end c c subroutine nibm2real(in,ro,n) character*(*) in dimension ro(1) iof=0 do i=1,n call gbyte(in,isign,iof,1) call gbyte(in,iexpn,iof+1,7) call gbyte(in,ifrac,iof+8,24) frac=ifrac rval=(-1)**isign * frac * 16**(iexpn-64) / 2**24 ro(i)=rval iof=iof+32 enddo return end c c subroutine ebcasc(nc,nd,num) c c convert num characters from ebcdic to ascii, nc is ebcdic input string c and nd is ascii output string, nc and nd may be the same string. c c the conversion table corresponds to the NCAR import/export conversion c and should give the same results. c character nc*(*),nd*(*) dimension ntab(256) data ntab/ a 000,001,002,003,156,009,134,127,151,141,142,011,012,013,014,015, b 016,017,018,019,157,133,008,135,024,025,146,143,028,029,030,031, c 128,129,130,131,132,010,023,027,136,137,138,139,140,005,006,007, d 144,145,022,147,148,149,150,004,152,153,154,155,020,021,158,026, e 032,160,161,162,163,164,165,166,167,168,213,046,060,040,043,124, f 038,169,170,171,172,173,174,175,176,177,033,036,042,041,059,094, g 045,047,178,179,180,181,182,183,184,185,229,044,037,095,062,063, h 186,187,188,189,190,191,192,193,194,096,058,035,064,039,061,034, i 195,097,098,099,100,101,102,103,104,105,196,197,198,199,200,201, j 202,106,107,108,109,110,111,112,113,114,203,204,205,206,207,208, k 209,126,115,116,117,118,119,120,121,122,210,211,212,091,214,215, l 216,217,218,219,220,221,222,223,224,225,226,227,228,093,230,231, m 123,065,066,067,068,069,070,071,072,073,232,233,234,235,236,237, n 125,074,075,076,077,078,079,080,081,082,238,239,240,241,242,243, o 092,159,083,084,085,086,087,088,089,090,244,245,246,247,248,249, p 048,049,050,051,052,053,054,055,056,057,250,251,252,253,254,255/ do 20 i=1,num nd(i:i)=char(ntab(ichar(nc(i:i))+1)) 20 continue return end c c C==================================================================== SUBROUTINE SPHERT 1(IDIR,GRID,WAVE,MLT,FAC,IMAX,JMAX,MAXWV,IROMB) C C THIS ROUTINE PERFORMS SPHERICAL TRANSFORM C C IDIR.GT.0 GRID TO WAVE. C IDIR= 1 ...... NORMAL TRANSFORM (GAUSSIAN GRID) C 3 ...... TRANSFORM WITH PNM*COS (GAUSSIAN GRID) C 4 ...... TRANSFORM WITH DPNM*COS (GAUSSIAN GRID) C 101 ...... AS IDIR=1 BUT FOR EQUIDISTANT LAT/LON GRID C 102 ...... AS IDIR=2 BUT FOR EQUIDISTANT LAT/LON GRID C IDIR.LT.0 WAVE TO GRID C IDIR= -1 ...... NORMAL TRANSFORM (GAUSSIAN GRID) C -2 ...... TRANSFORM WITH DPNM (GAUSSIAN GRID) C -101 ...... AS IDIR=-1 BUT FOR EQUIDISTANT LAT/LON GRID C -102 ...... AS IDIR=-2 BUT FOR EQUIDISTANT LAT/LON GRID C C 1 MULTIPLY 'FAC(NM)' TO WAVE VALUES C MLT = 0 NO MULTIPLY C -1 MULTIPLY 'I*FAC(NM)' (IMAGINARY) TO WAVE VALUES C C IROMB = 1 FOR RHOMBOIDAL TRUNCATION C = 0 FOR TRIANGULAR TRUNCATION C C NOTE C GRID IS REAL*4 GRID VALUES OF (IMAX,JMAX) C WAVE IS REAL*4 WAVE VALUES OF ( 2,NMMAX) C WORK IS REAL*4 WORKING ARRAY OF (IMAXP,JMAX,2) C C THE FOLLOWING PARAMETER STATEMENT SPECIFIES MAXIMUM POSSIBLE C WAVE TRUNCATION HANDLED IN THIS PROGRAM C C LIMWV ... MAXIMUM POSSIBLE WAVE NUMBER C LIMGAU .. MAXIMUM POSSIBLE GAUSSIAN LATITUDE C LROMB ... =0 TRIANGULAR =1 RHOMBOIDAL TRUNCATION C PARAMETER (LIMWV=150,LIMGAU=500,LROMB=1, 1 LIMWVX=(LIMWV+1)*(LIMWV+2)/2*(1-LROMB)+ 2 (LIMWV+1)*(LIMWV+1)*LROMB, 3 LIMGH =(LIMGAU+1)/2, 4 LIMWVT=LIMWV*2,LIMWVP=LIMWV+1, 5 MAXWK =LIMWVP*LIMGAU*2,LIMWV8=LIMWV*8) C real*8 GRID(1) real*4 WAVE(1) real*8 FAC(1) real*8 PNM (LIMWVX) real*8 SYMM(LIMWVX),DTRIGS(1000) real*8 COSCLT(LIMGAU) DIMENSION IFAX(10) C real*8 WORK(MAXWK) real*8 SWRK(LIMWV8) real*8 DDW(LIMGAU,6) C DIMENSION NFAC(11) C LOGICAL LSFFT save ifax C DATA NFAC/32,64,72,96,128,144,192,216,256,288,384/ DATA MFP1/0/,IRMB1/-999/ DATA MFP2/0/,IRMB2/-999/ DATA IFP/0/,JFP/0/,IDFP/0/ DATA IXTRA/2/ C LSFFT=.FALSE. DO 3 N=1,11 IF(IMAX.EQ.NFAC(N)) LSFFT=.TRUE. 3 CONTINUE C MEND1=MAXWV+1 NMMAX=(MEND1+1)*MEND1/2 IF(IROMB.EQ.1) NMMAX=MEND1*MEND1 NMMAXT=NMMAX*2 IMAXP=IMAX+IXTRA C IF(NMMAXT.GT.IMAXP*JMAX) THEN IXREQ=FLOAT(NMMAXT)/FLOAT(JMAX)-FLOAT(IMAX)+1 WRITE(6,8008) IXREQ 8008 FORMAT(5X,'IXTRA OF SUBROUTINE SPHERT INCREASED TO ',I8) IXTRA=IXREQ ENDIF C IMAXP=IMAX+IXTRA JMAXHF=(JMAX+1)/2 IPJHMX=IMAXP*JMAXHF IPJMAX=IMAXP*JMAX JMH=JMAX/2 C MEND1=MAXWV+1 MEND1T=MEND1*2 C IF(IDIR.GT.0.AND.NMMAX*2.GT.IMAX*JMAX) THEN WRITE(6,8007) 8007 FORMAT(5X,'*** POSSIBLE ARGUMENT LIST ERROR',/,5X, 1 '*** WAVE DATA HAS MORE DEGREE OF FREEDOM THAN GRID DATA') ENDIF C IF(LSFFT) THEN IF(IFP.NE.IMAX) THEN IMODE=2 CALL SFAX(IFAX,IMAX,IMODE) CALL SFTRIG(DTRIGS,IMAX,IMODE) IFP=IMAX ENDIF ENDIF C IF(MFP1.EQ.MAXWV.AND.IRMB1.EQ.IROMB) GO TO 2 NM=0 DO 10 MM=1,MEND1 NEND1=MEND1 IF(IROMB.EQ.1) NEND1=MM+MEND1-1 DO 20 NN=MM,NEND1 NM=NM+1 MODNM=MOD(NN-MM,2) IF(MODNM.EQ.1) GO TO 11 SYMM(NM)=1. GO TO 20 11 CONTINUE SYMM(NM)=-1. 20 CONTINUE 10 CONTINUE MFP1=MAXWV IRMB1=IROMB 2 CONTINUE C C END OF CONSTANTS C IF(IDIR.LT.0) GO TO 1000 C C FOREWARD TRANSFORM (GRID TO WAVE) C DO 100 J=1,JMAXHF IX2 =IMAX *(J-1) IY2 =IMAXP*(J-1) IX2R=IMAX*(JMAX-J) IY2H=IMAXP*(J-1)+IPJHMX DO 105 I=1,IMAX WORK(IY2 +I)=GRID(IX2+I)+GRID(IX2R+I) 105 CONTINUE IF(J.GT.JMH) GO TO 100 DO 106 I=1,IMAX WORK(IY2H+I)=GRID(IX2+I)-GRID(IX2R+I) 106 CONTINUE 100 CONTINUE C IF(LSFFT) THEN CALL SFT991(WORK,WORK(IPJMAX+1),DTRIGS,IFAX,1,IMAXP,IMAX,JMAX,-1, 1 DDW(1,1),DDW(1,2),DDW(1,3),DDW(1,4),DDW(1,5),DDW(1,6)) ELSE IMAXT=IMAX*2 DO 1212 J=1,JMAX IX1=IMAX *(J-1) IX2=IMAXP*(J-1) DO 1214 I=1,IMAX SWRK(2*I-1) = WORK(IX2+I) SWRK(2*I)= 0.0 1214 CONTINUE CALL FOURT(SWRK,IMAX,1,-1,0,WORK(IPJMAX+1)) write(*,*) ' FOURT ',(swrk(iii),iii=1,5) DO 1213 I=1,IMAXP WORK(IX2+I)=SWRK(I) 1213 CONTINUE 1212 CONTINUE ENDIF C IDPNM=MOD(IDIR,10) C NM=0 DO 155 NM=1,NMMAXT WORK(IPJMAX+NM)=0. 155 CONTINUE C DO 170 J=1,JMAXHF C C COMPUTE PNM/ZNM ETC. ON GAUSSIAN OR LAT/LON GRID C IF(IABS(IDIR).LT.100) THEN CALL PNMSGL(COSCLT,PNM,JMAX,MEND1,IDIR,IROMB,J) ELSE CALL ZNMSGL(COSCLT,PNM,JMAX,MEND1,IDIR,IROMB,J) ENDIF C DO 170 II=1,2 C NM=0 DO 160 MM=1,MEND1 IX2=2*(MM-1)-IMAXP NEND1=MEND1 IF(IROMB.EQ.1) NEND1=MM+MEND1-1 DO 160 NN=MM,NEND1 NM=NM+1 IY2=2*(NM-1)+IPJMAX MODNM=MOD(NN-MM,2) IF(IDPNM.EQ.2.OR.IDPNM.EQ.4) MODNM=1-MODNM C C SYMMETRIC PART C IF(MODNM.EQ.0) THEN IX3=IX2+II WORK(IY2+II)=WORK(IY2+II)+WORK(IMAXP*J+IX3)*PNM(NM) ELSE C C ANTISYMMETRIC PART C IF(J.LE.JMH) THEN IX3=IX2+II+IPJHMX WORK(IY2+II)=WORK(IY2+II)+WORK(IMAXP*J+IX3)*PNM(NM) ENDIF ENDIF C 160 CONTINUE 170 CONTINUE C C MULTIPLY REAL OR IMAGINARY CONSTANTS (FUNCTION OF N) C NMKEND=NMMAXT IF(MLT.NE.0) GO TO 121 DO 120 NMK=1,NMKEND WAVE(NMK)=WORK(IPJMAX+NMK) 120 CONTINUE RETURN C 121 CONTINUE IF(MLT.EQ.-1) GO TO 122 DO 123 NM=1,NMMAX WAVE(NM*2-1)=WORK(IPJMAX+NM*2-1)*FAC(NM) WAVE(NM*2)=WORK(IPJMAX+NM*2)*FAC(NM) 123 CONTINUE RETURN C 122 CONTINUE C MULTIPLY I*FAC C REAL PART DO 126 NM=1,NMMAX WAVE(2*NM-1)=-FAC(NM)*WORK(2*NM+IPJMAX) 126 CONTINUE C IMAGINARY PART DO 127 NM=1,NMMAX WAVE(2*NM )= FAC(NM)*WORK(2*NM+IPJMAX-1) 127 CONTINUE RETURN C C END OF GLOBAL CASE FOREWARD TRANSFORM C C----------------------------------------------------------------------- C-------- BACKWARD TRANSFORM ------------------------------------------ C----------------------------------------------------------------------- 1000 CONTINUE C C BACKWARD TRANSFORM (WAVE TO GRID) C NMF=NMMAX NMFT=NMF*2 C C ZERO CLEAR GRIDWK DO 205 IJK=1,IPJMAX WORK(IJK)=0. 205 CONTINUE C DO 250 J=1,JMAXHF C C COMPUTE PNM/ZNM ETC. ON GAUSSIAN GRID C IF(IABS(IDIR).LT.100) THEN CALL PNMSGL(COSCLT,PNM,JMAX,MEND1,IDIR,IROMB,J) ELSE CALL ZNMSGL(COSCLT,PNM,JMAX,MEND1,IDIR,IROMB,J) ENDIF c if(j.eq.2) then write(*,*) ' in ZNMSGL: j=2 ',(pnm(iix),iix=11,15) endif C IX2 =IMAXP*(J-1) IX2R=IMAXP*(JMAX-J) C DO 250 II=1,2 C IX3 =IX2 +II-2 IX3R=IX2R+II-2 IY2 = II-2 C IF(MLT.NE.0) GO TO 261 DO 260 NM=1,NMF WORK(IPJMAX+NM)=PNM(NM)*WAVE(NM*2+IY2) 260 CONTINUE c if(j.eq.2) c & write(*,*) ' in 260 work: j=2 ',(work(iix),iix=11,15) GO TO 269 261 CONTINUE IF(MLT.EQ.-1) GO TO 263 C MULTIPLY 'FAC' DO 262 NM=1,NMF WORK(IPJMAX+NM)=PNM(NM)*WAVE(NM*2+IY2)*FAC(NM) 262 CONTINUE GO TO 269 C MULTIPLY 'I*FAC' 263 CONTINUE IF(II.EQ.2) GO TO 266 C REAL PART DO 264 NM=1,NMF WORK(IPJMAX+NM)=PNM(NM)*(-WAVE(NM*2+IY2+1)*FAC(NM)) 264 CONTINUE GO TO 269 C IMAGINARY PART 266 CONTINUE DO 267 NM=1,NMF WORK(IPJMAX+NM)=PNM(NM)*WAVE(NM*2+IY2-1)*FAC(NM) 267 CONTINUE C 269 CONTINUE NM=0 DO 270 MM=1,MEND1 IX4=IX3+2*MM NEND1=MEND1 IF(IROMB.EQ.1) NEND1=MM+MEND1-1 DO 280 NN=MM,NEND1 NM=NM+1 WORK(IX4)=WORK(IX4)+WORK(IPJMAX+NM) 280 CONTINUE 270 CONTINUE C c if(j.eq.2) c & write(*,*) ' in 260 work: j=2 ',(work(iix),iix=11,15) 275 CONTINUE C GLOBAL CASE IF(J.GT.JMH) GO TO 250 NM=0 DO 470 MM=1,MEND1 IX4R=IX3R+2*MM NEND1=MEND1 IF(IROMB.EQ.1) NEND1=MM+MEND1-1 IF(IDIR.EQ.-2.OR.IDIR.EQ.-102) GO TO 485 DO 480 NN=MM,NEND1 NM=NM+1 WORK(IX4R)=WORK(IX4R)+WORK(IPJMAX+NM)*SYMM(NM) 480 CONTINUE GO TO 495 485 CONTINUE C C CASE OF DPNM DO 490 NN=MM,NEND1 NM=NM+1 WORK(IX4R)=WORK(IX4R)-WORK(IPJMAX+NM)*SYMM(NM) 490 CONTINUE 495 CONTINUE C 470 CONTINUE C 250 CONTINUE C IF(LSFFT) THEN CALL SFT991(WORK,WORK(IPJMAX+1),DTRIGS,IFAX,1,IMAXP,IMAX,JMAX,+1, 1 DDW(1,1),DDW(1,2),DDW(1,3),DDW(1,4),DDW(1,5),DDW(1,6)) DO 290 J=1,JMAX IX2=IMAX *(J-1) IY2=IMAXP*(J-1) DO 290 I=1,IMAX GRID(IX2+I)=WORK(IY2+I) 290 CONTINUE ELSE IMAXT=IMAX*2 DO 2212 J=1,JMAX IX1=IMAX *(J-1) IX2=IMAXP*(J-1) DO 2214 I=1,IMAXP SWRK(I)= WORK(IX2+I) 2214 CONTINUE DO 2215 I=4,IMAXP,2 SWRK(IMAXT-I+3)= WORK(IX2+I-1) SWRK(IMAXT-I+4)=-WORK(IX2+I) 2215 CONTINUE CALL FOURT(SWRK,IMAX,1,1,1,WORK(IPJMAX+1)) II=0 DO 2213 I=1,IMAXT,2 II=II+1 GRID(IX1+II)=SWRK(I) 2213 CONTINUE 2212 CONTINUE ENDIF RETURN END SUBROUTINE PNMSGL(COSCLT,PNMX,JMAX,MEND1,IDIR,IROMB,J) SAVE C C THIS SUBROUTINE COMPUTES LEGENDRE FUNCTION VALUES ON HEMISPHERIC C GAUSSIAN GRIDS AND ASSOCIATED GAUSSIAN QUADRATURE C C IDIR=-1 RETURNS PNM C =-2 DPNM -COS(LAT)*D(PNM)/D(LAT) C = 1 PNMG WEIGHTS*PNM C = 3 PNMGC WEIGHTS*PNM/COS(LAT)**2 C = 4 DPNMGC -WEIGHTS*D(PNM)/D(LAT)/COS(LAT) C C IROMB= 1 RHOMBOIDAL TRUNCATON C = OTHER TRIANGULAR TRUNCATION C PARAMETER (LIMWV=30+1,LROMB=1,LIMGAU=73, 1 LIMWVX=(LIMWV+1)*(LIMWV+2)/2*(1-LROMB)+ 2 (LIMWV+1)*(LIMWV+1)*LROMB) cxxx cxxx 2 (LIMWV+1)*(LIMWV+1)*LROMB+1) cxxx C REAL*8 PNMX(1) REAL*8 COSCLT(1) C REAL*8 PPNM(LIMWVX),HHNM(LIMWVX) REAL*8 GW(LIMGAU),GL(LIMGAU),GWCS(LIMGAU) COMMON/COMPZ/ PPNM,HHNM,GW,GL,GWCS C DATA IFP/0/ C IF(IFP.NE.JMAX) THEN IFP=JMAX NMMAX=MEND1*(MEND1+1)/2 IF(IROMB.EQ.1) NMMAX=MEND1*MEND1 IF(NMMAX.GT.LIMWVX) THEN WRITE(6,1001) 1001 FORMAT(5X,'*** INCREASE LIMWV OF ROUTINE PNMSGL') STOP 'ABEND' ENDIF JMAXHF=(JMAX+1)/2 CALL GAUAW1(GL,GW,JMAX) IF(MOD(JMAX,2).EQ.1) GW(JMAXHF)=GW(JMAXHF)*0.5D0 ENDIF C C START COMPUTATION OF LEGENDRE POLYNOMIAL C COSCLT(J)=(GL(J)) C COSCLT(J)=SNGL(GL(J)) GWCS(J)=GW(J)/(1.D0-GL(J)**2) CALL LGNDR1(GL(J),MEND1,PPNM,HHNM,IROMB) C IF(IDIR.EQ.-1) THEN DO 120 NM=1,NMMAX PNMX(NM)=(PPNM(NM)) C PNMX(NM)=SNGL(PPNM(NM)) 120 CONTINUE C ENDIF IF(IDIR.EQ.-2) THEN DO 130 NM=1,NMMAX PNMX(NM)=(HHNM(NM)) c PNMX(NM)=SNGL(HHNM(NM)) 130 CONTINUE ENDIF C IF(IDIR.EQ.1) THEN DO 140 NM=1,NMMAX PNMX(NM)=(PPNM(NM)*GW(J)) c PNMX(NM)=SNGL(PPNM(NM)*GW(J)) 140 CONTINUE ENDIF C IF(IDIR.EQ.2) THEN C DO 150 NM=1,NMMAX C PNMX(NM)=SNGL(HHNM(NM)*GW(J)) C 150 CONTINUE C ENDIF IF(IDIR.EQ.3) THEN DO 160 NM=1,NMMAX PNMX(NM)=(PPNM(NM)*GWCS(J)) c PNMX(NM)=SNGL(PPNM(NM)*GWCS(J)) 160 CONTINUE ENDIF IF(IDIR.EQ.4) THEN DO 170 NM=1,NMMAX PNMX(NM)=(HHNM(NM)*GWCS(J)) c PNMX(NM)=SNGL(HHNM(NM)*GWCS(J)) 170 CONTINUE ENDIF 110 CONTINUE C RETURN END SUBROUTINE ZNMSGL(COSCLT,PNMX,JMAX,MEND1,IDIR,IROMB,J) C C THIS SUBROUTINE COMPUTES LEGENDRE FUNCTION VALUES ON HEMISPHERIC C EQUIDISTANT GRIDS AND ASSOCIATE EQUIDISTANT QUADRATURE C C RETURNS PNM FOR IDIR=-101 C DPNM FOR IDIR=-102 C cxxxx c PARAMETER (LIMWV=30+1,LROMB=1,LIMGAU=73, c 1 LIMWVX=(LIMWV+2)*(LIMWV+3)/2*(1-LROMB)+ c 2 (LIMWV+2)*(LIMWV+2)*LROMB) cxxxx PARAMETER (LIMWV=30+1,LROMB=1,LIMGAU=73, 1 LIMWVX=(LIMWV+1)*(LIMWV+2)/2*(1-LROMB)+ 2 (LIMWV+1)*(LIMWV+1)*LROMB) C IMPLICIT REAL*8(A-H,O-Z) SAVE REAL*8 PNMX(1) REAL*8 COSCLT(1) cxxx REAL*4 COSCLT(LIMGAU) C COMMON/COMPZ/ PPNM,HHNM,GAUW,GAUL,GWCS REAL*8 PPNM(LIMWVX),HHNM(LIMWVX) REAL*8 GAUW(LIMGAU),GAUL(LIMGAU),GWCS(LIMGAU) c REAL*8 SPC1(LIMGAU),SPC2(LIMGAU),SPC3(LIMGAU) C IF(IDIR.GT.0) THEN PRINT *,'THIS ROUTINE CANNOT BE USED FOR EQUIDISTANT QUADRATURE' STOP 'ABEND' ENDIF C NMMAX=MEND1*(MEND1+1)/2 IF(IROMB.EQ.1) NMMAX=MEND1*MEND1 C IF(NMMAX.GT.LIMWVX) THEN WRITE(6,1001) 1001 FORMAT(5X,'*** INCREASE LIMWV OF ROUTINE ZNMSGL') STOP 'ABEND' ENDIF C JMAXHF=(JMAX+1)/2 PI=4.D0*ATAN(1.D0) HFPI=PI*0.5D0 DLAT=PI/DFLOAT(JMAX-1) C IY1=NMMAX*(J-1) RLAT=HFPI-DFLOAT(J-1)*DLAT SL=SIN(RLAT) COSCLT(J)=SL CALL LGNDR1(SL,MEND1,PPNM,HHNM,IROMB) C IF(IDIR.EQ.-101) THEN DO 20 NM=1,NMMAX PNMX(NM)=(PPNM(NM)) c PNMX(NM)=SNGL(PPNM(NM)) 20 CONTINUE ENDIF IF(IDIR.EQ.-102) THEN DO 25 NM=1,NMMAX PNMX(NM)=(HHNM(NM)) c PNMX(NM)=SNGL(HHNM(NM)) 25 CONTINUE ENDIF C if(j.eq.2) write(*,*) ' in NGMSGL: ',(pnmx(ka),ka=1,10) RETURN C END SUBROUTINE LGNDR1(COA,MFP,ALP,DALP,IROMB) C IMPLICIT REAL*8(A-H,O-Z) c DIMENSION ALP(*),DALP(1) real*8 ALP(1),DALP(1) C C SIA=SQRT(1.D0-COA*COA) LM=2 LMD=1 ALP(1)=SQRT(0.5D0) F1M=SQRT(1.5D0) ALP(2)=F1M*COA DALP(1)=0.D0 DO 1 M1=1,MFP M=M1-1 AM=M A2M=M+M RE1=SQRT(A2M+3.D0) E1=1.D0/RE1 IF(M.EQ.0) GO TO 3 F2M=F1M*SIA/SQRT(A2M) F1M=F2M*RE1 LM=LM+1 ALP(LM)=F2M IF(IROMB.EQ.1) GO TO 2 IF(M1.NE.MFP) GO TO 2 LMD=LMD+1 DALP(LMD)=AM*E1*F1M*COA GO TO 1 2 LM=LM+1 ALP(LM)=F1M*COA LMD=LMD+1 DALP(LMD)=AM*E1*ALP(LM) 3 M2=M+2 MFPX=MFP IF(IROMB.EQ.1) MFPX=MFP+M DO 4 N=M2,MFPX AN=N AN2=N*N E2=SQRT((4.D0*AN2-1.D0)/(AN2-AM*AM)) LM=LM+1 ALP(LM)=E2*(COA*ALP(LM-1)-E1*ALP(LM-2)) E2=1.D0/E2 LMD=LMD+1 DALP(LMD)=(AN-1.D0)*E2*ALP(LM)-AN*E1*ALP(LM-2) E1=E2 4 CONTINUE LM=LM-1 1 CONTINUE C RETURN END SUBROUTINE GAUAW1(A,W,K) C C IMPLICIT REAL*8(A-H,O-Z) real*8 A(K),W(K) C ESP=1.D-14 C=(1.D0-(2.D0/3.14159265358979D0)**2)*0.25D0 FK=K KK=K/2 CALL BSSLZ1(A,KK) DO 30 IS=1,KK XZ=COS(A(IS)/SQRT((FK+0.5D0)**2+C)) ITER=0 10 PKM2=1.D0 PKM1=XZ ITER=ITER+1 IF(ITER.GT.10) GO TO 70 DO 20 N=2,K FN=N PK=((2.D0*FN-1.D0)*XZ*PKM1-(FN-1.D0)*PKM2)/FN PKM2=PKM1 20 PKM1=PK PKM1=PKM2 PKMRK=(FK*(PKM1-XZ*PK))/(1.D0-XZ**2) SP=PK/PKMRK XZ=XZ-SP AVSP=ABS(SP) IF(AVSP.GT.ESP) GO TO 10 A(IS)=XZ W(IS)=(2.D0*(1.D0-XZ**2))/(FK*PKM1)**2 30 CONTINUE IF(K.EQ.KK*2) GO TO 50 A(KK+1)=0.D0 PK=2.D0/FK**2 DO 40 N=2,K,2 FN=N 40 PK=PK*FN**2/(FN-1.D0)**2 W(KK+1)=PK 50 CONTINUE DO 60 N=1,KK L=K+1-N A(L)=-A(N) 60 W(L)=W(N) C PRINT *,'PRINT OUT FROM GAU JMAX=',K C PRINT *,A C PRINT *,W RETURN 70 WRITE(6,6000) 6000 FORMAT(//5X,14HERROR IN GAUAW//) STOP END c SUBROUTINE BSSLZ1(BES,N) IMPLICIT REAL*8(A-H,O-Z) real*8 BES(N) real*8 BZ(50) C DATA PI/3.14159265358979D0/ DATA BZ / 2.4048255577D0, 5.5200781103D0, $ 8.6537279129D0,11.7915344391D0,14.9309177086D0,18.0710639679D0, $ 21.2116366299D0,24.3524715308D0,27.4934791320D0,30.6346064684D0, $ 33.7758202136D0,36.9170983537D0,40.0584257646D0,43.1997917132D0, $ 46.3411883717D0,49.4826098974D0,52.6240518411D0,55.7655107550D0, $ 58.9069839261D0,62.0484691902D0,65.1899648002D0,68.3314693299D0, $ 71.4729816036D0,74.6145006437D0,77.7560256304D0,80.8975558711D0, $ 84.0390907769D0,87.1806298436D0,90.3221726372D0,93.4637187819D0, $ 96.6052679510D0,99.7468198587D0,102.888374254D0,106.029930916D0, $ 109.171489649D0,112.313050280D0,115.454612653D0,118.596176630D0, $ 121.737742088D0,124.879308913D0,128.020877005D0,131.162446275D0, $ 134.304016638D0,137.445588020D0,140.587160352D0,143.728733573D0, $ 146.870307625D0,150.011882457D0,153.153458019D0,156.295034268D0/ NN=N IF(N.LE.50) GO TO 12 BES(50)=BZ(50) DO 5 J=51,N 5 BES(J)=BES(J-1)+PI NN=49 12 DO 15 J=1,NN 15 BES(J)=BZ(J) RETURN END SUBROUTINE SFAX(IFAX,N,MODE) C C DIMENSION IFAX(10) C 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 C 10 CONTINUE K=1 20 CONTINUE 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 30 CONTINUE 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 40 CONTINUE 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 50 CONTINUE L=5 INC=2 60 CONTINUE 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 CONTINUE L=L+INC INC=6-INC GO TO 60 80 CONTINUE IFAX(1)=K-1 NFAX=IFAX(1) 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 C END C SUBROUTINE SFTRIG(TRIGS,N,MODE) C C C* IMPLICIT REAL*8(A-H,O-Z) real*8 TRIGS(1) C C* PI=2.D0*ASIN(1.D0) PI=2. *ASIN(1. ) IMODE=IABS(MODE) NN=N IF(IMODE.GT.1.AND.IMODE.LT.6) NN=N/2 DEL=(PI+PI)/DFLOAT(NN) L=NN+NN DO 10 I=1,L,2 C* ANGLE=0.5D0*DFLOAT(I-1)*DEL ANGLE=0.5 *DFLOAT(I-1)*DEL TRIGS(I)=COS(ANGLE) TRIGS(I+1)=SIN(ANGLE) 10 CONTINUE IF(IMODE.EQ.1) RETURN IF(IMODE.EQ.8) RETURN C* DEL=0.5D0*DEL DEL=0.5 *DEL NH=(NN+1)/2 L=NH+NH LA=NN+NN DO 20 I=1,L,2 C* ANGLE=0.5D0*DFLOAT(I-1)*DEL ANGLE=0.5 *DFLOAT(I-1)*DEL TRIGS(LA+I)=COS(ANGLE) TRIGS(LA+I+1)=SIN(ANGLE) 20 CONTINUE IF(IMODE.LE.3) RETURN C* DEL=0.5D0*DEL DEL=0.5 *DEL LA=LA+NN IF(MODE.EQ.5) GO TO 40 DO 30 I=2,NN ANGLE=DFLOAT(I-1)*DEL C* TRIGS(LA+I)=2.D0*SIN(ANGLE) TRIGS(LA+I)=2. *SIN(ANGLE) 30 CONTINUE RETURN 40 CONTINUE C* DEL=0.5D0*DEL DEL=0.5 *DEL DO 50 I=2,N ANGLE=DFLOAT(I-1)*DEL TRIGS(LA+I)=SIN(ANGLE) 50 CONTINUE RETURN C END C SUBROUTINE SFT991(A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN, $ W1,W2,W3,W4,W5,W6) C C C* IMPLICIT REAL*8(A-H,O-Z) DIMENSION IFAX(10) real*8 A(N),WORK(N),TRIGS(N) real*8 W1(LOT),W2(LOT),W3(LOT),W4(LOT),W5(LOT),W6(LOT) C NFAX=IFAX(1) NX=N+1 NH=N/2 INK=INC+INC IF(ISIGN.EQ.1) GO TO 30 C IGO=50 IF(MOD(NFAX,2).EQ.1) GO TO 40 IBASE=1 JBASE=1 DO 20 L=1,LOT I=IBASE J=JBASE 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 30 CONTINUE CALL SFT99A(A,WORK,TRIGS,INC,JUMP,N,LOT,W1,W2,W3,W4,W5,W6) IGO=60 C 40 CONTINUE IA=1 LA=1 DO 80 K=1,NFAX c write (*, 9950) k, nfax 9950 format ("k is ", I2, " runs to ", I6) c write (*, 9951) ifax 9951 format ("ifax: ", 10I7) IF(IGO.EQ.60) GO TO 60 50 CONTINUE CALL SPASSM(A(IA),A(IA+INC),WORK(1),WORK(2),TRIGS,INK,2,JUMP,NX, $ LOT,NH,IFAX(K+1),LA,W1,W2,W3,W4,W5,W6) IGO=60 GO TO 70 60 CONTINUE CALL SPASSM(WORK(1),WORK(2),A(IA),A(IA+INC),TRIGS,2,INK,NX,JUMP, $ LOT,NH,IFAX(K+1),LA,W1,W2,W3,W4,W5,W6) IGO=50 70 CONTINUE LA=LA*IFAX(K+1) 80 CONTINUE C IF(ISIGN.EQ.-1) GO TO 130 C IF(MOD(NFAX,2).EQ.1) GO TO 110 IBASE=1 JBASE=1 DO 100 L=1,LOT I=IBASE J=JBASE 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 110 CONTINUE RETURN C 130 CONTINUE CALL SFT99B(WORK,A,TRIGS,INC,JUMP,N,LOT,W1,W2,W3,W4,W5,W6) C 140 CONTINUE RETURN C END C SUBROUTINE SFT99A(A,WORK,TRIGS,INC,JUMP,N,LOT,W1,W2,W3,W4,W5,W6) C C C* IMPLICIT REAL*8(A-H,O-Z) real*8 A(N),WORK(N),TRIGS(N) real*8 W1(LOT),W2(LOT),W3(LOT),W4(LOT),W5(LOT),W6(LOT) C* DATA DP2,DM2/2.D0,-2.D0/ DATA DP2,DM2/2. ,-2. / C NH=N/2 NX=N+1 INK=INC+INC C IA=1 IB=N*INC+1 JA=1 JB=2 DO 11 L=1,LOT WORK(JA)=A(IA)+A(IB) IA=IA+JUMP IB=IB+JUMP JA=JA+NX 11 CONTINUE IA=1 IB=N*INC+1 DO 12 L=1,LOT WORK(JB)=A(IA)-A(IB) IA=IA+JUMP IB=IB+JUMP JB=JB+NX 12 CONTINUE C 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) DO 21 L=1,LOT W1(L)=A(IA)+A(IB) W2(L)=A(IA)-A(IB) W3(L)=A(IA+INC)+A(IB+INC) W4(L)=A(IA+INC)-A(IB+INC) IA=IA+JUMP IB=IB+JUMP 21 CONTINUE DO 22 L=1,LOT W5(L)=S*W2(L) W6(L)=C*W2(L) 22 CONTINUE DO 23 L=1,LOT W5(L)=W5(L)+C*W3(L) W6(L)=W6(L)-S*W3(L) 23 CONTINUE DO 24 L=1,LOT WORK(JA)=W1(L)-W5(L) JA=JA+NX 24 CONTINUE DO 25 L=1,LOT WORK(JB)=W1(L)+W5(L) JB=JB+NX 25 CONTINUE JA=JABASE DO 26 L=1,LOT WORK(JA+1)=W6(L)+W4(L) JA=JA+NX 26 CONTINUE JB=JBBASE DO 27 L=1,LOT WORK(JB+1)=W6(L)-W4(L) JB=JB+NX 27 CONTINUE IABASE=IABASE+INK IBBASE=IBBASE-INK JABASE=JABASE+2 JBBASE=JBBASE-2 30 CONTINUE C IF(IABASE.NE.IBBASE) GO TO 50 IA=IABASE JA=JABASE DO 41 L=1,LOT WORK(JA)=DP2*A(IA) IA=IA+JUMP JA=JA+NX 41 CONTINUE IA=IABASE JA=JABASE DO 42 L=1,LOT WORK(JA+1)=DM2*A(IA+INC) IA=IA+JUMP JA=JA+NX 42 CONTINUE C 50 CONTINUE RETURN C END C SUBROUTINE SFT99B(WORK,A,TRIGS,INC,JUMP,N,LOT,W1,W2,W3,W4,W5,W6) C C C* IMPLICIT REAL*8(A-H,O-Z) real*8 WORK(N),A(N),TRIGS(N) real*8 W1(LOT),W2(LOT),W3(LOT),W4(LOT),W5(LOT),W6(LOT) C NH=N/2 NX=N+1 INK=INC+INC C C* SCALE=1.D0/DFLOAT(N) SCALE=1. /DFLOAT(N) IA=1 IB=2 JA=1 JB=N*INC+1 DO 14 L=1,LOT W1(L)=WORK(IA)+WORK(IB) W2(L)=WORK(IA)-WORK(IB) IA=IA+NX IB=IB+NX 14 CONTINUE DO 11 L=1,LOT A(JA)=SCALE*W1(L) JA=JA+JUMP 11 CONTINUE DO 12 L=1,LOT A(JB)=SCALE*W2(L) JB=JB+JUMP 12 CONTINUE JA=1 DO 13 L=1,LOT C* A(JA+INC)=0.D0 A(JA+INC)=0. JA=JA+JUMP 13 CONTINUE C C* SCALE=0.5D0*SCALE 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 )*SCALE S=TRIGS(N+K+1)*SCALE DO 21 L=1,LOT W1(L)=WORK(IA)+WORK(IB) W2(L)=WORK(IA)-WORK(IB) W3(L)=WORK(IA+1)+WORK(IB+1) W4(L)=WORK(IB+1)-WORK(IA+1) IA=IA+NX IB=IB+NX 21 CONTINUE DO 29 L=1,LOT W1(L)=SCALE*W1(L) W4(L)=SCALE*W4(L) 29 CONTINUE DO 22 L=1,LOT W5(L)=C*W3(L) W6(L)=C*W2(L) 22 CONTINUE DO 23 L=1,LOT W5(L)=W5(L)+S*W2(L) W6(L)=W6(L)-S*W3(L) 23 CONTINUE DO 24 L=1,LOT A(JA)=W1(L)+W5(L) JA=JA+JUMP 24 CONTINUE DO 25 L=1,LOT A(JB)=W1(L)-W5(L) JB=JB+JUMP 25 CONTINUE JA=JABASE DO 26 L=1,LOT A(JA+INC)=W6(L)+W4(L) JA=JA+JUMP 26 CONTINUE JB=JBBASE DO 27 L=1,LOT A(JB+INC)=W6(L)-W4(L) JB=JB+JUMP 27 CONTINUE IABASE=IABASE+2 IBBASE=IBBASE-2 JABASE=JABASE+INK JBBASE=JBBASE-INK 30 CONTINUE C IF(IABASE.NE.IBBASE) GO TO 50 IA=IABASE JA=JABASE C* SCALE=2.D0*SCALE SCALE=2. *SCALE SMINS=-SCALE DO 41 L=1,LOT A(JA)=SCALE*WORK(IA) IA=IA+NX JA=JA+JUMP 41 CONTINUE IA=IABASE JA=JABASE DO 42 L=1,LOT A(JA+INC)=SMINS*WORK(IA+1) IA=IA+NX JA=JA+JUMP 42 CONTINUE C 50 CONTINUE RETURN C END C SUBROUTINE SPASSM(A,B,C,D,TRIGS,INC1,INC2,INC3,INC4,LOT,N,IFAC,LA, $ W1,W2,W3,W4,W5,W6) C C C* IMPLICIT REAL*8(A-H,O-Z) real*8 A(N),B(N),C(N),D(N),TRIGS(N) real*8 W1(LOT),W2(LOT),W3(LOT),W4(LOT),W5(LOT),W6(LOT) C* DATA SIN60/0.866025403784437D0/ DATA SIN60/0.8660254 / 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),IGO C 10 CONTINUE IA=1 JA=1 IB=IA+IINK JB=JA+JINK DO 20 L=1,LA I=IBASE J=JBASE DO 15 IJK=1,LOT C(JA+J)=A(IA+I)+A(IB+I) D(JA+J)=B(IA+I)+B(IB+I) I=I+INC3 J=J+INC4 15 CONTINUE I=IBASE J=JBASE DO 16 IJK=1,LOT C(JB+J)=A(IA+I)-A(IB+I) D(JB+J)=B(IA+I)-B(IB+I) I=I+INC3 J=J+INC4 16 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 DO 25 IJK=1,LOT C(JA+J)=A(IA+I)+A(IB+I) D(JA+J)=B(IA+I)+B(IB+I) I=I+INC3 J=J+INC4 25 CONTINUE I=IBASE DO 26 IJK=1,LOT W1(IJK)=A(IA+I)-A(IB+I) W2(IJK)=B(IA+I)-B(IB+I) I=I+INC3 26 CONTINUE J=JBASE DO 27 IJK=1,LOT C(JB+J)=C1*W1(IJK) D(JB+J)=S1*W1(IJK) J=J+INC4 27 CONTINUE J=JBASE DO 28 IJK=1,LOT C(JB+J)=C(JB+J)-S1*W2(IJK) D(JB+J)=D(JB+J)+C1*W2(IJK) J=J+INC4 28 CONTINUE IBASE=IBASE+INC1 JBASE=JBASE+INC2 30 CONTINUE JBASE=JBASE+JUMP 40 CONTINUE RETURN C 50 CONTINUE IA=1 JA=1 IB=IA+IINK JB=JA+JINK IC=IB+IINK JC=JB+JINK DO 60 L=1,LA I=IBASE DO 55 IJK=1,LOT W1(IJK)=A(IB+I)+A(IC+I) W2(IJK)=B(IB+I)+B(IC+I) W3(IJK)=A(IB+I)-A(IC+I) W4(IJK)=B(IB+I)-B(IC+I) I=I+INC3 55 CONTINUE DO 51 IJK=1,LOT W3(IJK)=SIN60*W3(IJK) W4(IJK)=SIN60*W4(IJK) 51 CONTINUE I=IBASE J=JBASE DO 56 IJK=1,LOT C(JA+J)=A(IA+I)+W1(IJK) D(JA+J)=B(IA+I)+W2(IJK) I=I+INC3 J=J+INC4 56 CONTINUE DO 57 IJK=1,LOT C* W1(IJK)=0.5D0*W1(IJK) C* W2(IJK)=0.5D0*W2(IJK) W1(IJK)=0.5 *W1(IJK) W2(IJK)=0.5 *W2(IJK) 57 CONTINUE I=IBASE DO 52 IJK=1,LOT W1(IJK)=A(IA+I)-W1(IJK) W2(IJK)=B(IA+I)-W2(IJK) I=I+INC3 52 CONTINUE J=JBASE DO 58 IJK=1,LOT C(JB+J)=W1(IJK)-W4(IJK) D(JB+J)=W2(IJK)+W3(IJK) J=J+INC4 58 CONTINUE J=JBASE DO 59 IJK=1,LOT C(JC+J)=W1(IJK)+W4(IJK) D(JC+J)=W2(IJK)-W3(IJK) J=J+INC4 59 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 DO 61 IJK=1,LOT W1(IJK)=A(IB+I)+A(IC+I) W2(IJK)=B(IB+I)+B(IC+I) W3(IJK)=A(IB+I)-A(IC+I) W4(IJK)=B(IB+I)-B(IC+I) I=I+INC3 61 CONTINUE DO 71 IJK=1,LOT W3(IJK)=SIN60*W3(IJK) W4(IJK)=SIN60*W4(IJK) 71 CONTINUE I=IBASE J=JBASE DO 62 IJK=1,LOT C(JA+J)=A(IA+I)+W1(IJK) D(JA+J)=B(IA+I)+W2(IJK) I=I+INC3 J=J+INC4 62 CONTINUE DO 63 IJK=1,LOT C* W1(IJK)=0.5D0*W1(IJK) C* W2(IJK)=0.5D0*W2(IJK) W1(IJK)=0.5 *W1(IJK) W2(IJK)=0.5 *W2(IJK) 63 CONTINUE I=IBASE DO 72 IJK=1,LOT W1(IJK)=A(IA+I)-W1(IJK) W2(IJK)=B(IA+I)-W2(IJK) I=I+INC3 72 CONTINUE DO 64 IJK=1,LOT W5(IJK)=W1(IJK)-W4(IJK) W6(IJK)=W2(IJK)+W3(IJK) 64 CONTINUE J=JBASE DO 65 IJK=1,LOT C(JB+J)=C1*W5(IJK) D(JB+J)=S1*W5(IJK) J=J+INC4 65 CONTINUE J=JBASE DO 68 IJK=1,LOT C(JB+J)=C(JB+J)-S1*W6(IJK) D(JB+J)=D(JB+J)+C1*W6(IJK) J=J+INC4 68 CONTINUE DO 66 IJK=1,LOT W5(IJK)=W1(IJK)+W4(IJK) W6(IJK)=W2(IJK)-W3(IJK) 66 CONTINUE J=JBASE DO 67 IJK=1,LOT C(JC+J)=C2*W5(IJK) D(JC+J)=S2*W5(IJK) J=J+INC4 67 CONTINUE J=JBASE DO 69 IJK=1,LOT C(JC+J)=C(JC+J)-S2*W6(IJK) D(JC+J)=D(JC+J)+C2*W6(IJK) J=J+INC4 69 CONTINUE IBASE=IBASE+INC1 JBASE=JBASE+INC2 70 CONTINUE JBASE=JBASE+JUMP 80 CONTINUE RETURN C 90 CONTINUE 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 DO 91 IJK=1,LOT W1(IJK)=A(IA+I)+A(IC+I) W2(IJK)=A(IB+I)+A(ID+I) W3(IJK)=B(IA+I)+B(IC+I) W4(IJK)=B(IB+I)+B(ID+I) I=I+INC3 91 CONTINUE J=JBASE DO 92 IJK=1,LOT C(JA+J)=W1(IJK)+W2(IJK) D(JA+J)=W3(IJK)+W4(IJK) J=J+INC4 92 CONTINUE J=JBASE DO 93 IJK=1,LOT C(JC+J)=W1(IJK)-W2(IJK) D(JC+J)=W3(IJK)-W4(IJK) J=J+INC4 93 CONTINUE I=IBASE DO 94 IJK=1,LOT W1(IJK)=A(IA+I)-A(IC+I) W2(IJK)=A(IB+I)-A(ID+I) W3(IJK)=B(IA+I)-B(IC+I) W4(IJK)=B(IB+I)-B(ID+I) I=I+INC3 94 CONTINUE J=JBASE DO 95 IJK=1,LOT C(JB+J)=W1(IJK)-W4(IJK) D(JB+J)=W3(IJK)+W2(IJK) J=J+INC4 95 CONTINUE J=JBASE DO 96 IJK=1,LOT C(JD+J)=W1(IJK)+W4(IJK) D(JD+J)=W3(IJK)-W2(IJK) J=J+INC4 96 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 DO 101 IJK=1,LOT W1(IJK)=A(IA+I)+A(IC+I) W2(IJK)=A(IB+I)+A(ID+I) W3(IJK)=B(IA+I)+B(IC+I) W4(IJK)=B(IB+I)+B(ID+I) I=I+INC3 101 CONTINUE J=JBASE DO 102 IJK=1,LOT C(JA+J)=W1(IJK)+W2(IJK) D(JA+J)=W3(IJK)+W4(IJK) J=J+INC4 102 CONTINUE DO 103 IJK=1,LOT W1(IJK)=W1(IJK)-W2(IJK) W3(IJK)=W3(IJK)-W4(IJK) 103 CONTINUE J=JBASE DO 104 IJK=1,LOT C(JC+J)=C2*W1(IJK) D(JC+J)=S2*W1(IJK) J=J+INC4 104 CONTINUE J=JBASE DO 105 IJK=1,LOT C(JC+J)=C(JC+J)-S2*W3(IJK) D(JC+J)=D(JC+J)+C2*W3(IJK) J=J+INC4 105 CONTINUE I=IBASE DO 111 IJK=1,LOT W1(IJK)=A(IA+I)-A(IC+I) W2(IJK)=A(IB+I)-A(ID+I) W3(IJK)=B(IA+I)-B(IC+I) W4(IJK)=B(IB+I)-B(ID+I) I=I+INC3 111 CONTINUE DO 112 IJK=1,LOT W5(IJK)=W1(IJK)-W4(IJK) W6(IJK)=W3(IJK)+W2(IJK) 112 CONTINUE J=JBASE DO 113 IJK=1,LOT C(JB+J)=C1*W5(IJK) D(JB+J)=S1*W5(IJK) J=J+INC4 113 CONTINUE J=JBASE DO 117 IJK=1,LOT C(JB+J)=C(JB+J)-S1*W6(IJK) D(JB+J)=D(JB+J)+C1*W6(IJK) J=J+INC4 117 CONTINUE DO 114 IJK=1,LOT W1(IJK)=W1(IJK)+W4(IJK) W3(IJK)=W3(IJK)-W2(IJK) 114 CONTINUE J=JBASE DO 115 IJK=1,LOT C(JD+J)=C3*W1(IJK) D(JD+J)=S3*W1(IJK) J=J+INC4 115 CONTINUE J=JBASE DO 116 IJK=1,LOT C(JD+J)=C(JD+J)-S3*W3(IJK) D(JD+J)=D(JD+J)+C3*W3(IJK) J=J+INC4 116 CONTINUE IBASE=IBASE+INC1 JBASE=JBASE+INC2 110 CONTINUE JBASE=JBASE+JUMP 120 CONTINUE RETURN C C END SUBROUTINE FOURT(DATA,NN,NDIM,ISIGN,IFORM,WORK) C cxxx DIMENSION DATA(*),NN(1),IFACT(32),WORK(1) dimension NN(1),IFACT(32) real*8 DATA(1),WORK(1) TWOPI=6.283185307 RTHLF=.7071067812 NP0=0 NPREV=0 IF(NDIM-1)920,1,1 1 NTOT=2 DO 2 IDIM=1,NDIM IF(NN(IDIM))920,920,2 2 NTOT=NTOT*NN(IDIM) C C MAIN LOOP FOR EACH DIMENSION C NP1=2 DO 910 IDIM=1,NDIM N=NN(IDIM) NP2=NP1*N IF(N-1)920,900,5 C C IS N A POWER OF TWO AND IF NOT, WHAT ARE ITS FACTORS C 5 M=N NTWO=NP1 IF=1 IDIV=2 10 IQUOT=M/IDIV IREM=M-IDIV*IQUOT IF(IQUOT-IDIV)50,11,11 11 IF(IREM)20,12,20 12 NTWO=NTWO+NTWO IFACT(IF)=IDIV IF=IF+1 M=IQUOT GO TO 10 20 IDIV=3 INON2=IF 30 IQUOT=M/IDIV IREM=M-IDIV*IQUOT IF(IQUOT-IDIV)60,31,31 31 IF(IREM)40,32,40 32 IFACT(IF)=IDIV IF=IF+1 M=IQUOT GO TO 30 40 IDIV=IDIV+2 GO TO 30 50 INON2=IF IF(IREM)60,51,60 51 NTWO=NTWO+NTWO GO TO 70 60 IFACT(IF)=M 70 NON2P=NP2/NTWO C C SEPARATE FOUR CASES-- C 1. COMPLEX TRANSFORM C 2. REAL TRANSFORM FOR THE 2ND, 3RD, ETC. DIMENSION. METHOD-- C TRANSFORM HALF THE DATA, SUPPLYING THE OTHER HALF BY CON- C JUGATE SYMMETRY. C 3. REAL TRANSFORM FOR THE 1ST DIMENSION, N ODD. METHOD-- C SET THE IMAGINARY PARTS TO ZERO. C 4. REAL TRANSFORM FOR THE 1ST DIMENSION, N EVEN. METHOD-- C TRANSFORM A COMPLEX ARRAY OF LENGTH N/2 WHOSE REAL PARTS C ARE THE EVEN NUMBERED REAL VALUES AND WHOSE IMAGINARY PARTS C ARE THE ODD NUMBERED REAL VALUES. SEPARATE AND SUPPLY C THE SECOND HALF BY CONJUGATE SYMMETRY. C IFMIN = 1 I1RNG = NP1 IF (IFORM .LE. 0 .AND. IDIM .LT. 4) GO TO 71 ICASE = 1 GO TO 100 C 71 IF (IDIM .LE. 1) GO TO 72 ICASE = 2 I1RNG = NP0 * (1 + NPREV / 2) GO TO 100 C 72 IF (NTWO .GT. NP1) GO TO 73 ICASE = 3 GO TO 100 C 73 ICASE = 4 IFMIN = 2 NTWO = NTWO / 2 N = N / 2 NP2 = NP2 / 2 NTOT = NTOT / 2 C I = - 1 DO 80 J = 1, NTOT I = I + 2 DATA(J) = DATA(I) 80 CONTINUE C C SHUFFLE DATA BY BIT REVERSAL, SINCE N=2**K. AS THE SHUFFLING C CAN BE DONE BY SIMPLE INTERCHANGE, NO WORKING ARRAY IS NEEDED C 100 IF(NON2P-1)101,101,200 101 NP2HF=NP2/2 J=1 DO 150 I2=1,NP2,NP1 IF(J-I2)121,130,130 121 I1MAX=I2+NP1-2 DO 125 I1=I2,I1MAX,2 DO 125 I3=I1,NTOT,NP2 J3=J+I3-I2 TEMPR=DATA(I3) TEMPI=DATA(I3+1) DATA(I3)=DATA(J3) DATA(I3+1)=DATA(J3+1) DATA(J3)=TEMPR 125 DATA(J3+1)=TEMPI 130 M=NP2HF 140 IF(J-M)150,150,141 141 J=J-M M=M/2 IF(M-NP1)150,140,140 150 J=J+M GO TO 300 C C SHUFFLE DATA BY DIGIT REVERSAL FOR GENERAL N C 200 NWORK=2*N DO 270 I1=1,NP1,2 DO 270 I3=I1,NTOT,NP2 J=I3 DO 260 I=1,NWORK,2 IF(ICASE-3)210,220,210 210 WORK(I)=DATA(J) WORK(I+1)=DATA(J+1) GO TO 240 220 WORK(I)=DATA(J) WORK(I+1)=0. 240 IFP2=NP2 IF=IFMIN 250 IFP1=IFP2/IFACT(IF) J=J+IFP1 IF(J-I3-IFP2)260,255,255 255 J=J-IFP2 IFP2=IFP1 IF=IF+1 IF(IFP2-NP1)260,260,250 260 CONTINUE I2MAX=I3+NP2-NP1 I=1 DO 270 I2=I3,I2MAX,NP1 DATA(I2)=WORK(I) DATA(I2+1)=WORK(I+1) 270 I=I+2 C C MAIN LOOP FOR FACTORS OF TWO. C W=EXP(ISIGN*2*PI*SQRT(-1)*M/(4*MMAX)). CHECK FOR W=ISIGN*SQRT(-1) C AND REPEAT FOR W=W*(1+ISIGN*SQRT(-1))/SQRT(2). C 300 IF(NTWO-NP1)600,600,305 305 NP1TW=NP1+NP1 IPAR=NTWO/NP1 310 IF(IPAR-2)350,330,320 320 IPAR=IPAR/4 GO TO 310 330 DO 340 I1=1,I1RNG,2 DO 340 K1=I1,NTOT,NP1TW K2=K1+NP1 TEMPR=DATA(K2) TEMPI=DATA(K2+1) DATA(K2)=DATA(K1)-TEMPR DATA(K2+1)=DATA(K1+1)-TEMPI DATA(K1)=DATA(K1)+TEMPR 340 DATA(K1+1)=DATA(K1+1)+TEMPI 350 MMAX=NP1 360 IF(MMAX-NTWO/2)370,600,600 370 LMAX=MAX0(NP1TW,MMAX/2) DO 570 L=NP1,LMAX,NP1TW M=L IF(MMAX-NP1)420,420,380 380 THETA=-TWOPI*FLOAT(L)/FLOAT(4*MMAX) IF(ISIGN)400,390,390 390 THETA=-THETA 400 WR=COS(THETA) WI=SIN(THETA) 410 W2R=WR*WR-WI*WI W2I=2.*WR*WI W3R=W2R*WR-W2I*WI W3I=W2R*WI+W2I*WR 420 DO 530 I1=1,I1RNG,2 KMIN=I1+IPAR*M IF(MMAX-NP1)430,430,440 430 KMIN=I1 440 KDIF=IPAR*MMAX 450 KSTEP=4*KDIF IF(KSTEP-NTWO)460,460,530 460 DO 520 K1=KMIN,NTOT,KSTEP K2=K1+KDIF K3=K2+KDIF K4=K3+KDIF IF(MMAX-NP1)470,470,480 470 U1R=DATA(K1)+DATA(K2) U1I=DATA(K1+1)+DATA(K2+1) U2R=DATA(K3)+DATA(K4) U2I=DATA(K3+1)+DATA(K4+1) U3R=DATA(K1)-DATA(K2) U3I=DATA(K1+1)-DATA(K2+1) IF(ISIGN)471,472,472 471 U4R=DATA(K3+1)-DATA(K4+1) U4I=DATA(K4)-DATA(K3) GO TO 510 472 U4R=DATA(K4+1)-DATA(K3+1) U4I=DATA(K3)-DATA(K4) GO TO 510 480 T2R=W2R*DATA(K2)-W2I*DATA(K2+1) T2I=W2R*DATA(K2+1)+W2I*DATA(K2) T3R=WR*DATA(K3)-WI*DATA(K3+1) T3I=WR*DATA(K3+1)+WI*DATA(K3) T4R=W3R*DATA(K4)-W3I*DATA(K4+1) T4I=W3R*DATA(K4+1)+W3I*DATA(K4) U1R=DATA(K1)+T2R U1I=DATA(K1+1)+T2I U2R=T3R+T4R U2I=T3I+T4I U3R=DATA(K1)-T2R U3I=DATA(K1+1)-T2I IF(ISIGN)490,500,500 490 U4R=T3I-T4I U4I=T4R-T3R GO TO 510 500 U4R=T4I-T3I U4I=T3R-T4R 510 DATA(K1)=U1R+U2R DATA(K1+1)=U1I+U2I DATA(K2)=U3R+U4R DATA(K2+1)=U3I+U4I DATA(K3)=U1R-U2R DATA(K3+1)=U1I-U2I DATA(K4)=U3R-U4R 520 DATA(K4+1)=U3I-U4I KDIF=KSTEP KMIN=4*(KMIN-I1)+I1 GO TO 450 530 CONTINUE M=M+LMAX IF(M-MMAX)540,540,570 540 IF(ISIGN)550,560,560 550 TEMPR=WR WR=(WR+WI)*RTHLF WI=(WI-TEMPR)*RTHLF GO TO 410 560 TEMPR=WR WR=(WR-WI)*RTHLF WI=(TEMPR+WI)*RTHLF GO TO 410 570 CONTINUE IPAR=3-IPAR MMAX=MMAX+MMAX GO TO 360 C C MAIN LOOP FOR FACTORS NOT EQUAL TO TWO. C W=EXP(ISIGN*2*PI*SQRT(-1)*(J1+J2-I3-1)/IFP2) C 600 IF(NON2P-1)700,700,601 601 IFP1=NTWO IF=INON2 610 IFP2=IFACT(IF)*IFP1 THETA=-TWOPI/FLOAT(IFACT(IF)) IF(ISIGN)612,611,611 611 THETA=-THETA 612 WSTPR=COS(THETA) WSTPI=SIN(THETA) DO 650 J1=1,IFP1,NP1 THETM=-TWOPI*FLOAT(J1-1)/FLOAT(IFP2) IF(ISIGN)614,613,613 613 THETM=-THETM 614 WMINR=COS(THETM) WMINI=SIN(THETM) I1MAX=J1+I1RNG-2 DO 650 I1=J1,I1MAX,2 DO 650 I3=I1,NTOT,NP2 I=1 WR=WMINR WI=WMINI J2MAX=I3+IFP2-IFP1 DO 640 J2=I3,J2MAX,IFP1 TWOWR=WR+WR J3MAX=J2+NP2-IFP2 DO 630 J3=J2,J3MAX,IFP2 JMIN=J3-J2+I3 J=JMIN+IFP2-IFP1 SR=DATA(J) SI=DATA(J+1) OLDSR=0. OLDSI=0. J=J-IFP1 620 STMPR=SR STMPI=SI SR=TWOWR*SR-OLDSR+DATA(J) SI=TWOWR*SI-OLDSI+DATA(J+1) OLDSR=STMPR OLDSI=STMPI J=J-IFP1 IF(J-JMIN)621,621,620 621 WORK(I)=WR*SR-WI*SI-OLDSR+DATA(J) WORK(I+1)=WI*SR+WR*SI-OLDSI+DATA(J+1) 630 I=I+2 WTEMP=WR*WSTPI WR=WR*WSTPR-WI*WSTPI 640 WI=WI*WSTPR+WTEMP I=1 DO 650 J2=I3,J2MAX,IFP1 J3MAX=J2+NP2-IFP2 DO 650 J3=J2,J3MAX,IFP2 DATA(J3)=WORK(I) DATA(J3+1)=WORK(I+1) 650 I=I+2 IF=IF+1 IFP1=IFP2 IF(IFP1-NP2)610,700,700 C C COMPLETE A REAL TRANSFORM IN THE 1ST DIMENSION, N EVEN, BY CON- C JUGATE SYMMETRIES. C 700 GO TO (900,800,900,701),ICASE 701 NHALF=N N=N+N THETA=-TWOPI/FLOAT(N) IF(ISIGN)703,702,702 702 THETA=-THETA 703 WSTPR=COS(THETA) WSTPI=SIN(THETA) WR=WSTPR WI=WSTPI IMIN=3 JMIN=2*NHALF-1 GO TO 725 710 J=JMIN DO 720 I=IMIN,NTOT,NP2 SUMR=(DATA(I)+DATA(J))/2. SUMI=(DATA(I+1)+DATA(J+1))/2. DIFR=(DATA(I)-DATA(J))/2. DIFI=(DATA(I+1)-DATA(J+1))/2. TEMPR=WR*SUMI+WI*DIFR TEMPI=WI*SUMI-WR*DIFR DATA(I)=SUMR+TEMPR DATA(I+1)=DIFI+TEMPI DATA(J)=SUMR-TEMPR DATA(J+1)=-DIFI+TEMPI 720 J=J+NP2 IMIN=IMIN+2 JMIN=JMIN-2 WTEMP=WR*WSTPI WR=WR*WSTPR-WI*WSTPI WI=WI*WSTPR+WTEMP 725 IF(IMIN-JMIN)710,730,740 730 IF(ISIGN)731,740,740 731 DO 735 I=IMIN,NTOT,NP2 735 DATA(I+1)=-DATA(I+1) 740 NP2=NP2+NP2 NTOT=NTOT+NTOT J=NTOT+1 IMAX=NTOT/2+1 745 IMIN=IMAX-2*NHALF I=IMIN GO TO 755 750 DATA(J)=DATA(I) DATA(J+1)=-DATA(I+1) 755 I=I+2 J=J-2 IF(I-IMAX)750,760,760 760 DATA(J)=DATA(IMIN)-DATA(IMIN+1) DATA(J+1)=0. IF(I-J)770,780,780 765 DATA(J)=DATA(I) DATA(J+1)=DATA(I+1) 770 I=I-2 J=J-2 IF(I-IMIN)775,775,765 775 DATA(J)=DATA(IMIN)+DATA(IMIN+1) DATA(J+1)=0. IMAX=IMIN GO TO 745 780 DATA(1)=DATA(1)+DATA(2) DATA(2)=0. GO TO 900 C C COMPLETE A REAL TRANSFORM FOR THE 2ND, 3RD, ETC. DIMENSION BY C CONJUGATE SYMMETRIES. C 800 IF(I1RNG-NP1)805,900,900 805 DO 860 I3=1,NTOT,NP2 I2MAX=I3+NP2-NP1 DO 860 I2=I3,I2MAX,NP1 IMAX=I2+NP1-2 IMIN=I2+I1RNG JMAX=2*I3+NP1-IMIN IF(I2-I3)820,820,810 810 JMAX=JMAX+NP2 820 IF(IDIM-2)850,850,830 830 J=JMAX+NP0 DO 840 I=IMIN,IMAX,2 DATA(I)=DATA(J) DATA(I+1)=-DATA(J+1) 840 J=J-2 850 J=JMAX DO 860 I=IMIN,IMAX,NP0 DATA(I)=DATA(J) DATA(I+1)=-DATA(J+1) 860 J=J-NP0 C C END OF LOOP ON EACH DIMENSION C 900 NP0=NP1 NP1=NP2 910 NPREV=N 920 IF(ISIGN.EQ.1) RETURN NTOT=2 DO 930 IDIM=1,NDIM 930 NTOT=NTOT*NN(IDIM) NTOTHF=NTOT/2 DO 940 ITOT=1,NTOT 940 DATA(ITOT)=DATA(ITOT)/NTOTHF RETURN END c c SUBROUTINE UNP360(IN,IO,N) PARAMETER (IWSZ=32) C FOR EXECUTION ON WORD SIZE IWSZ (>=32) C CONVERTS ARRAYS OF IBM 360/370 WORDS WHICH ARE IN A CONTINUOUS BIT C STREAM INTO CRAY WORDS. C WILL ACCEPT MIXED ARRAYS OF REAL (REAL*4) AND INTEGER (INTEGER*4) C AND RETURNS THE APPROPRIATE CRAY REAL OR FULL WORD INTEGER C VALUES. C IN = ARRAY CONTAINING 360 WORDS AS A CONTINUOUS BIT STREAM. M C IO = ARRAY FOR CRAY OUTPUT WORDS. C N = NUMBER OF VALUES TO CONVERT. c DIMENSION IN(1),IO(1) character*(*) in dimension io(1) IOF=0 DO 20 I=1,N CALL GBYTES(IN,IS,IOF,1,0,1) IOF=IOF+1 CALL GBYTES(IN,IEX,IOF,7,0,1) IOF=IOF+7 CALL GBYTES(IN,IFR,IOF,24,0,1) IF(IEX .EQ. 0 .OR. IEX .EQ. 127) GO TO 10 XV=FLOAT(IFR)*16.**(IEX-70) IF(IS .NE. 0) XV=-XV CALL SBYTES(IO(I),XV,0,IWSZ,0,1) GOTO 20 10 CONTINUE IF(IS .NE. 0) IFR=IFR-16777216 IO(I)=IFR 20 IOF=IOF+24 RETURN END