program rd0840alpha c c This program will read ds084.0 data files retrieved using "t2d0840" program. c c DEC Alpha FORTRAN version: To compile, f77 rd0840alpha.f c C The converted grids are in arrays. c parameter (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 infunit, idum, itwo(2000) character one*8000, inf*72, two*8000 c character nmclabel*32 real*4 fhour, wave1(maxz2) real*4 wave2(maxuu2) integer*4 iyear,imonth,iday,ihour c integer*4 rcsize(mxrc),p(12) c c The (idim) grids of row one of final (idim,jdim,12) grids are at the North Pole, c i.e. (1, 1,?) is at 0E,90N or North Pole c (1,jdim,?) is at 0E,90S or South Pole c real*4 z(idim,jdim,12), u(idim,jdim,12), v(idim,jdim,12), & rh(idim,jdim,6), vt(idim,jdim,12) real*4 grdv(idim*jdim), sinclt(jdim) equivalence (itwo,two) 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/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/ 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 c sinclt(1)=1. ! added 93oct14 sinclt(jdim)=1. ! added 93oct14 c c open the file c infunit=1 write(*,*) 'enter input file? ' read(*,'(a)') inf open(unit=infunit,file=inf,form='unformatted', & access='sequential',status='unknown') c nfile=0 10 continue c c read the first record c read(infunit,end=90,err=9999) (one(i:i),i=1,rcsize(1)) call ebcasc(one(1:32),nmclabel,rcsize(1)) write(*,*) nmclabel c c read the second record c Note: Although DOC says there are 20 bytes, we found its acutal c length is 24 bytes. c read(infunit) (one(i:i),i=1,rcsize(2)) call swap4(one,two,rcsize(2)) call gbyte(itwo(1),idum,0,32) call unp360(idum,fhour,1) call gbyte(itwo(1),iyear,32,32) call gbyte(itwo(1),imonth,64,32) call gbyte(itwo(1),iday,96,32) call gbyte(itwo(1),ihour,128,32) write(*,*) fhour,iyear,imonth,iday,ihour c c read the Z,U,V for all 12 pressure levels c do k=3,38 zmax=-99999. zmin=99999. umax=-99999. umin=99999. vmax=-99999. vmin=99999. kp=k/3 kind=mod(k,3) read(infunit) (one(i:i),i=1,rcsize(k)) if(kind.eq.0) then call swap4(one,two,rcsize(k)) call unp360(itwo,wave1,maxz2) call sphert(idir,grdv,wave1,mlt,fac, & idim,jdim,mx31-1,iromb) do j=1,jdim do i=1,idim z(i,j,kp)=grdv((j-1)*idim+i) if(z(i,j,kp).gt.zmax) zmax=z(i,j,kp) if(z(i,j,kp).lt.zmin) zmin=z(i,j,kp) enddo enddo c write(*,991) ' Z',(z(ii,36,kp),ii=10,15) 991 format(1x,a2,6E12.4) write(*,*) ' z max,min: ',zmax,zmin,' at ',p(kp),'mb' elseif(kind.eq.1) then call swap4(one,two,rcsize(k)) call unp360(itwo,wave2,maxuv2) c write(*,*) 'UU 1-10 ',(wave2(ii),ii=501,510) call sphert(idir,grdv,wave2,mlt,fac,idim, & jdim,mx32-1,iromb) do j=1,jdim do i=1,idim u(i,j,kp)=grdv((j-1)*idim+i)/sinclt(j) ! changed 93oct14 if(u(i,j,kp).gt.umax) umax=u(i,j,kp) if(u(i,j,kp).lt.umin) umin=u(i,j,kp) enddo enddo c write(*,991) ' U',(u(ii,36,kp),ii=10,15) write(*,*) ' u max,min: ',umax,umin,' at ',p(kp),'mb' else call swap4(one,two,rcsize(k)) call unp360(itwo,wave2,maxuv2) call sphert(idir,grdv,wave2,mlt,fac,idim, & jdim,mx32-1,iromb) do j=1,jdim do i=1,idim v(i,j,kp)=grdv((j-1)*idim+i)/sinclt(j) ! changed 93oct14 if(v(i,j,kp).gt.vmax) vmax=v(i,j,kp) if(v(i,j,kp).lt.vmin) vmin=v(i,j,kp) enddo enddo c write(*,991) ' V',(v(ii,36,kp),ii=10,15) write(*,*) ' v max,min: ',vmax,vmin,' at ',p(kp),'mb' endif enddo c c RH records at 6 pressure levels c do k=39,44 rhmax=-99999. rhmin=99999. kp=k-38 read(infunit) (one(i:i),i=1,rcsize(k)) call swap4(one,two,rcsize(k)) call unp360(itwo,wave1,maxz2) call sphert(idir,grdv,wave1,mlt,fac,idim, & jdim,mx31-1,iromb) do j=1,jdim do i=1,idim rh(i,j,kp)=grdv((j-1)*idim+i) if(rh(i,j,kp).gt.rhmax) rhmax=rh(i,j,kp) if(rh(i,j,kp).lt.rhmin) rhmin=rh(i,j,kp) enddo enddo c write(*,991) 'RH',(rh(ii,36,kp),ii=10,15) write(*,*) 'rh max,min: ',rhmax,rhmin,' at ',p(kp),'mb' enddo c c VT records: 12 pressure levels c do k=45,56 vtmax=-99999. vtmin=99999. kp=k-44 read(infunit) (one(i:i),i=1,rcsize(k)) call swap4(one,two,rcsize(k)) call unp360(itwo,wave1,maxz2) call sphert(idir,grdv,wave1,mlt,fac,idim, & jdim,mx31-1,iromb) do j=1,jdim do i=1,idim vt(i,j,kp)=grdv((j-1)*idim+i) if(vt(i,j,kp).gt.vtmax) vtmax=vt(i,j,kp) if(vt(i,j,kp).lt.vtmin) vtmin=vt(i,j,kp) enddo enddo c write(*,991) 'VT',(vt(ii,36,kp),ii=10,15) write(*,*) 'vt max,min: ',vtmax,vtmin,' at ',p(kp),'mb' enddo write(*,*) ' DONE one forecast hour.... ' c go to 10 c c close the file c 90 continue write(*,*) ' finish reading input file....' close(infunit) c c to clear the SUN IEEE computational flags c c idum=ieee_flags("clear","exception","all",outstr) c 9998 stop'Normal Termination' c 9999 write(*,*) '**** ERROR during read, stop' call exit(0) 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*(*) integer*4 num 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 The subroutines below are from NMC 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 DIMENSION GRID(1),WAVE(1) DIMENSION FAC(1) C DIMENSION PNM (LIMWVX) C DIMENSION SYMM(LIMWVX) DIMENSION COSCLT(LIMGAU) DIMENSION IFAX(10),DTRIGS(1000) C DIMENSION WORK(MAXWK) DIMENSION SWRK(LIMWV8) DIMENSION 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) c write (*, 9980) imax, ifax 9980 format ("calling sfax with imax: ", I8, /, "ifax: ",10I7) 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)) 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 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 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 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 c write (*, 9901) ifax 9901 format (" before setf991, ifax:", 10I7) 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)) c write (*, 9902) ifax 9902 format (" after setf991, ifax:", 10I7) 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*4 PNMX(1) REAL*4 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*4 PNMX(1) REAL*4 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 RETURN C END SUBROUTINE LGNDR1(COA,MFP,ALP,DALP,IROMB) C IMPLICIT REAL*8(A-H,O-Z) c DIMENSION ALP(*),DALP(1) DIMENSION 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) DIMENSION 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 SUBROUTINE BSSLZ1(BES,N) C C IMPLICIT REAL*8(A-H,O-Z) DIMENSION BES(N) DIMENSION 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) DIMENSION 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 A(N),WORK(N),TRIGS(N),IFAX(10) DIMENSION 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) DIMENSION A(N),WORK(N),TRIGS(N) DIMENSION 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) DIMENSION WORK(N),A(N),TRIGS(N) DIMENSION 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) DIMENSION A(N),B(N),C(N),D(N),TRIGS(N) DIMENSION 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 DATA(1),NN(1),IFACT(32),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 SUBROUTINE GBYTE (IN,IOUT,ISKIP,NBYTE) CALL GBYTES (IN,IOUT,ISKIP,NBYTE,0,1) RETURN END SUBROUTINE SBYTE (IOUT,IN,ISKIP,NBYTE) CALL SBYTES (IOUT,IN,ISKIP,NBYTE,0,1) RETURN END SUBROUTINE GBYTES (IN,IOUT,ISKIP,NBYTE,NSKIP,N) C Get bytes - unpack bits: Extract arbitrary size values from a C packed bit string, right justifying each value in the unpacked C array. DIMENSION IN(*), IOUT(*) C IN = packed array input C IO = unpacked array output C ISKIP = initial number of bits to skip C NBYTE = number of bits to take C NSKIP = additional number of bits to skip on each iteration C N = number of iterations C************************************** MACHINE SPECIFIC CHANGES START HERE C Machine dependent information required: C LMWD = Number of bits in a word on this machine C MASKS = Set of word masks where the first element has only the C right most bit set to 1, the second has the two, ... C LEFTSH = Shift left bits in word M to the by N bits C RGHTSH = Shift right C OR = Logical OR (add) on this machine. C AND = Logical AND (multiply) on this machine C This is for Sun UNIX Fortran, DEC Alpha, and RS6000 PARAMETER (LMWD=32) DIMENSION MASKS(LMWD) SAVE MASKS DATA MASKS /'1'X,'3'X,'7'X,'F'X, '1F'X,'3F'X,'7F'X,'FF'X, +'1FF'X,'3FF'X,'7FF'X,'FFF'X, '1FFF'X,'3FFF'X,'7FFF'X,'FFFF'X, +'1FFFF'X, '3FFFF'X, '7FFFF'X, 'FFFFF'X, +'1FFFFF'X, '3FFFFF'X, '7FFFFF'X, 'FFFFFF'X, +'1FFFFFF'X, '3FFFFFF'X, '7FFFFFF'X, 'FFFFFFF'X, +'1FFFFFFF'X, '3FFFFFFF'X, '7FFFFFFF'X, 'FFFFFFFF'X/ C +'1FFFFFFFF'X, '3FFFFFFFF'X, '7FFFFFFFF'X, 'FFFFFFFFF'X, C +'1FFFFFFFFF'X, '3FFFFFFFFF'X, '7FFFFFFFFF'X, 'FFFFFFFFFF'X, C +'1FFFFFFFFFF'X, '3FFFFFFFFFF'X, '7FFFFFFFFFF'X, 'FFFFFFFFFFF'X, C +'1FFFFFFFFFFF'X,'3FFFFFFFFFFF'X,'7FFFFFFFFFFF'X,'FFFFFFFFFFFF'X, C +'1FFFFFFFFFFFF'X, '3FFFFFFFFFFFF'X, '7FFFFFFFFFFFF'X, C + 'FFFFFFFFFFFFF'X, C +'1FFFFFFFFFFFFF'X, '3FFFFFFFFFFFFF'X, '7FFFFFFFFFFFFF'X, C 'FFFFFFFFFFFFFF'X, C +'1FFFFFFFFFFFFFF'X, '3FFFFFFFFFFFFFF'X, '7FFFFFFFFFFFFFF'X, C 'FFFFFFFFFFFFFFF'X, C +'1FFFFFFFFFFFFFFF'X,'3FFFFFFFFFFFFFFF'X,'7FFFFFFFFFFFFFFF'X, C 'FFFFFFFFFFFFFFFF'X/ C IBM PC using Microsoft Fortran uses different syntax: C DATA MASKS/16#1,16#3,16#7,16#F,16#1F,16#3F,16#7F,16#FF, C + 16#1FF,16#3FF,16#7FF,16#FFF,16#1FFF,16#3FFF,16#7FFF,16#FFFF, C + 16#1FFFF,16#3FFFF,16#7FFFF,16#FFFFF,16#1FFFFF,16#3FFFFF, C + 16#7FFFFF,16#FFFFFF,16#1FFFFFF,16#3FFFFFF,16#7FFFFFF,16#FFFFFFF, C + 16#1FFFFFFF,16#3FFFFFFF,16#7FFFFFFF,16#FFFFFFFF/ INTEGER RGHTSH, OR, AND LEFTSH(M,N) = ISHFT(M,N) RGHTSH(M,N) = ISHFT(M,-N) C OR(M,N) = M.OR.N C AND(M,N) = M.AND.N C OR(M,N) = IOR(M,N) C AND(M,N) = IAND(M,N) C************************************** MACHINE SPECIFIC CHANGES END HERE C History: written by Robert C. Gammill, jul 1972. C NBYTE must be less than or equal to LMWD ICON = LMWD-NBYTE IF (ICON.LT.0) RETURN MASK = MASKS (NBYTE) C INDEX = number of words into IN before the next "byte" appears C II = number of bits the "byte" is from the left side of the word C ISTEP = number of bits from the start of one "byte" to the next C IWORDS = number of words to skip from one "byte" to the next C IBITS = number of bits to skip after skipping IWORDS C MOVER = number of bits to the right, a byte must be moved to be C right adjusted INDEX = ISKIP/LMWD II = MOD (ISKIP,LMWD) ISTEP = NBYTE+NSKIP IWORDS= ISTEP/LMWD IBITS = MOD (ISTEP,LMWD) DO 6 I=1,N MOVER = ICON-II IF (MOVER) 2,3,4 C The "byte" is split across a word break. 2 MOVEL = -MOVER MOVER = LMWD-MOVEL NP1 = LEFTSH (IN(INDEX+1),MOVEL) NP2 = RGHTSH (IN(INDEX+2),MOVER) IOUT(I) = AND (OR (NP1,NP2) , MASK) GO TO 5 C The "byte" is already right adjusted. 3 IOUT(I) = AND (IN (INDEX+1) , MASK) GO TO 5 C Right adjust the "byte". 4 IOUT(I) = AND (RGHTSH (IN (INDEX+1),MOVER) , MASK) 5 II = II+IBITS INDEX = INDEX+IWORDS IF (II .LT. LMWD) GO TO 6 II = II-LMWD INDEX = INDEX+1 6 CONTINUE RETURN END SUBROUTINE SBYTES (IOUT,IN,ISKIP,NBYTE,NSKIP,N) C Store bytes - pack bits: Put arbitrary size values into a C packed bit string, taking the low order bits from each value C in the unpacked array. DIMENSION IN(*), IOUT(*) C IOUT = packed array output C IN = unpacked array input C ISKIP = initial number of bits to skip C NBYTE = number of bits to pack C NSKIP = additional number of bits to skip on each iteration C N = number of iterations C************************************** MACHINE SPECIFIC CHANGES START HERE C Machine dependent information required: C LMWD = Number of bits in a word on this machine C MASKS = Set of word masks where the first element has only the C right most bit set to 1, the second has the two, ... C LEFTSH = Shift left bits in word M to the by N bits C RGHTSH = Shift right C OR = Logical OR (add) on this machine C AND = Logical AND (multiply) on this machine C NOT = Logical NOT (negation) on this machine C This is for Sun UNIX Fortran PARAMETER (LMWD=32) DIMENSION MASKS(LMWD) SAVE MASKS DATA MASKS /'1'X,'3'X,'7'X,'F'X, '1F'X,'3F'X,'7F'X,'FF'X, +'1FF'X,'3FF'X,'7FF'X,'FFF'X, '1FFF'X,'3FFF'X,'7FFF'X,'FFFF'X, +'1FFFF'X, '3FFFF'X, '7FFFF'X, 'FFFFF'X, +'1FFFFF'X, '3FFFFF'X, '7FFFFF'X, 'FFFFFF'X, +'1FFFFFF'X, '3FFFFFF'X, '7FFFFFF'X, 'FFFFFFF'X, +'1FFFFFFF'X, '3FFFFFFF'X, '7FFFFFFF'X, 'FFFFFFFF'X/ INTEGER RGHTSH, OR, AND LEFTSH(M,N) = ISHFT(M,N) RGHTSH(M,N) = ISHFT(M,-N) C OR(M,N) = M.OR.N C AND(M,N) = M.AND.N C OR(M,N) = IOR(M,N) C AND(M,N) = IAND(M,N) C NOT(M) = .NOT.M C*********************************************************************** C NBYTE must be less than or equal to LMWD ICON = LMWD-NBYTE IF (ICON .LT. 0) RETURN MASK = MASKS(NBYTE) C INDEX = number of words into IOUT the next "byte" is to be stored C II = number of bits in from the left side of the word to store it C ISTEP = number of bits from the start of one "byte" to the next C IWORDS = number of words to skip from one "byte" to the next C IBITS = number of bits to skip after skipping IWORDS C MOVER = number of bits to the right, a byte must be moved to be C right adjusted INDEX = ISKIP/LMWD II = MOD(ISKIP,LMWD) ISTEP = NBYTE+NSKIP IWORDS = ISTEP/LMWD IBITS = MOD(ISTEP,LMWD) DO 6 I=1,N J = AND (MASK,IN(I)) MOVEL = ICON-II IF (MOVEL) 2,3,4 C The "byte" is to be split across a word break 2 MSK = MASKS (NBYTE+MOVEL) IOUT(INDEX+1) = OR (AND(NOT(MSK),IOUT(INDEX+1)),RGHTSH(J,-MOVEL)) ITEMP = AND (MASKS(LMWD+MOVEL),IOUT(INDEX+2)) IOUT(INDEX+2) = OR(ITEMP,LEFTSH(J,LMWD+MOVEL)) GO TO 5 C The "byte" is to be stored right-adjusted 3 IOUT(INDEX+1) = OR ( AND (NOT(MASK),IOUT(INDEX+1)) , J) GO TO 5 C The "byte" is to be stored in middle of word, so shift left. 4 MSK = LEFTSH(MASK,MOVEL) IOUT(INDEX+1) = OR(AND(NOT(MSK),IOUT(INDEX+1)),LEFTSH(J,MOVEL)) 5 II = II+IBITS INDEX = INDEX+IWORDS IF (II .LT. LMWD) GO TO 6 II = II-LMWD INDEX = INDEX+1 6 CONTINUE RETURN END subroutine swap4(in,io,nn) c swaps bytes in groups of 4 to compensate for byte swapping within c words which occurs on DEC (VAX) and PC machines. c c in - input array to be swapped c io - ouput array with bytes swapped c nn - number of bytes to be swapped character*1 in(1),io(1),ih do 10 i=1,nn,4 ih=in(i) io(i)=in(i+3) io(i+3)=ih ih=in(i+1) io(i+1)=in(i+2) io(i+2)=ih 10 continue return end subroutine swap2(in,io,nn) c swaps bytes in groups of 2 to compensate for byte swapping within c words which occurs on DEC (VAX) and PC machines. c c in - input array to be swapped c io - ouput array with bytes swapped c nn - number of bytes to be swapped character*1 in(1),io(1),ih do 10 i=1,nn,2 ih=in(i) io(i)=in(i+1) io(i+1)=ih 10 continue return end subroutine filt(m,n,in) character*1 m(1),n(1),l,u data l/32/,u/127/ do 10 i=1,in n(i)=m(i) if(n(i).lt. l) n(i)=l if(n(i).gt. u) n(i)=u 10 continue return end