C COMPUTES STREAMFN,VELOCITY POTENTIAL & DIVERGENCE C VIA SPHERICAL HARMONICS (K. MO) AT R-30 RESOLUTION WHICH C IS ABOUT 2.5 DEG LAT/LON ON GAUSSIAN GRID c c c IMPORTANT: Use "-old_rl" option when compiling on SGI c---------------------------------------------------------------------- c c---- Note that pgm "cdas_data_2_nic" reads the strmfn & vel. pot. c fields that are computed in this pgm & written to unit 2 c so if changes are made to the data order in file "cdas_r30", c similar modifications must be done to "cdas_data_2_nic" c c-------------------------------------------------------------------- c c Modified for 4-digit year on 6-18-97 (JJ) c c dimension u(144,73), v(144,73), psi(144,73), chi(144,73), 1 div(144,73), uc(144,73), vc(144,73) c ccc character*18 arg ccc character*12 date_sys c ccc data arg /'date >sys_clk.date'/ ccc data date_sys /'sys_clk.date'/ data isunit1 /21/, isunit2 /22/, icunit /88/, ioctl /60/ c c Note: file "uvwind" created by script "get_200uv.scr" c open (1,file= 1 '/export-8/cacsrv1/cpc/anal/cdas/bulletin/pgms/uvwind', 2 recl=144*73*4,form='unformatted',access='direct') c open (2,file='../ctl/cdas_r30.grads',recl=144*73*4, 1 form='unformatted',access='direct') ccc open (3,file= ccc 1 '/export-8/cacsrv1/cpc/anal/cdas/month/prs/mean/extr.uvth_ieee.cl ccc 2im79-95.data',recl=144*73*4,form='unformatted',access='direct') open (3,file= x'/export-8/cacsrv1/cpc/anal/cdas/bulletin/pgms/uvclim', x recl=144*73*4,form='unformatted',access='direct') c open (4,file='/export-8/cacsrv1/cpc/anal/cdas/month/prs/mean/cdas_ xyear_mon') c read(1,rec=1) u read(1,rec=2) v c call STRSPH(isunit1,isunit2,U,V,PSI,CHI,DIV) write(2,rec=1) psi write(2,rec=2) chi write(2,rec=3) div c........Get date from system clock ccc call system(arg) ccc call gtdate_clk(icunit,date_sys,iyr,mon,ida) c c ccc jmon=mon-1 ccc jyr=iyr ccc if(jmon.eq.0) then ccc jmon=12 ccc jyr=iyr-1 ccc endif c c........ Get date from file "cdas_year_mon" c read(4,*) iyrmon jyr=iyrmon/100 jmon=mod(iyrmon,100) c write(*,*)' using month: ',jmon write(*,*)' for anomaly computations in "spherical_r30.exe"' c----- c...... get climatological means of u & v for anomaly computationc c cccc iurec= 7 + (jmon-1)*32 !rec # for 200 U from "CDDB"-structured CDAS climo cccc ivrec= 15 + (jmon-1)*32 !rec # for 200 V from "CDDB"-structured CDAS climo c iurec=1 ivrec=2 c----- read(3,rec=iurec) uc read(3,rec=ivrec) vc c c...... compute anomalous psi, chi & div c do j=1,73 do i=1,144 u(i,j)=u(i,j)-uc(i,j) v(i,j)=v(i,j)-vc(i,j) enddo enddo c c call STRSPH(isunit1,isunit2,U,V,PSI,CHI,DIV) write(2,rec=4) psi write(2,rec=5) chi write(2,rec=6) div c call updctl(ioctl,jmon,jyr) c stop 'done' end c c-------------------------------------------------------------------- c subroutine gtdate_clk(icunit,date_sys,iyr,mon,ida) character*3 xmon character*12 date_sys character*3 cmon(12) data cmon /'Jan','Feb','Mar','Apr','May','Jun', 1 'Jul','Aug','Sep','Oct','Nov','Dec'/ c c Reads UNIX date information from file "sys_clk.date" that is c created in "Main" from call to "system" c c open (icunit,file=date_sys) read(icunit,105) xmon,ida,iyr 105 format(4x,a3,i3,14x,i4) c do imon=1,12 if(xmon.eq.cmon(imon)) go to 250 enddo write(*,*) 'no match for month: ',xmon stop 'pgm stopped in subroutine "gtdate_clk"' c 250 mon=imon return end c c--------------------------------------------------------- c SUBROUTINE STRSPH(isunit1,isunit2,U,V,PSI,CHI,DIV) C C STREAMFUNCTION ROUTINE COMPUTED USING SPHERICAL HARMONICS DONATED C BY KINGTSE MO (11/19/90) C DIMENSION WAVE(1922),WAVE2(1922), WAVE3(1922) DIMENSION FS(961), COR(73), FFF(961), FAC(961), ALAT(73) DIMENSION PSI(144,73), CHI(144,73), DIV(144,73), GG(144,73) DIMENSION U(144,73), V(144,73), GLOB(144,73) COMMON/WAVY/PNMC(43560,3),PVV(43560) C DATA IOPEN /0/ C IF(IOPEN.NE.0) GO TO 150 C C 1ST TIME IN THIS ROUTINE, SO READ IN SPHERICAL DATA C C READ FROM THE DISK FILE 'NWS.WD52.KMO.PNM' THE HARMONICS C TO R30 open(isunit1,file='../pgms/r30_coefs1', 1 recl=43560*4,form='unformatted',access='direct') open(isunit2,file='../pgms/r30_coefs2', 1 recl=43560*3*4,form='unformatted',access='direct') c READ(ISUNIT1,rec=1) PVV READ(ISUNIT2,rec=1) PNMC IOPEN=1 C 150 RA=6375000. FNO=1 CONST=19.4*20.4 PI=4.D0*ATAN(1.D0) HFPI=PI*0.5 MEND=31 DY=2.*RA*PI/72. OMEGA=7.292E-05 C DO 50 I=1,73 AAT=HFPI-DFLOAT(I-1)*PI/73. ALAT(I)=COS(AAT) COR(I)=2.*OMEGA*SIN(AAT) cccc PRINT 8870,I,COR(I) 8870 FORMAT(1X,I8,1E15.6) 50 CONTINUE C DO 45 K=1,961 FS(K)=0. 45 CONTINUE C NM=0 DO 4711 MM=1,MEND MEND1=MM+MEND-1 DO 4711 NN=MM,MEND1 NM=NM+1 S2=NN*(NN-1)/(19.4*20.4) S4=-(S2*S2) FS(NM)=EXP(S4) 4711 CONTINUE C NM=0 DO 40 MM=1,MEND MEND1=MM+MEND-1 DO 40 NN=MM,MEND1 NM=NM+1 FAC(NM)=MM-1 40 CONTINUE C NM=0 DO 43 MM=1,MEND MEND1=MM+MEND-1 DO 43 NN=MM,MEND1 NM=NM+1 FNM=(NN-1)*(NN) FFF(NM)=0. IF(FNM.NE.0.) FFF(NM)=-RA*RA/FNM 43 CONTINUE C C TO DO CALCULATIONS IN THE SPHERICAL COORDINATES FO=0 DO 70 I=1,144 DO 70 J=1,73 70 GLOB(I,J)=U(I,J)*ALAT(J)/RA C DO 735 I=1,144 GLOB(I,1)=0. 735 GLOB(I,73)=0. C CALL SPHERT(3,GLOB,WAVE,0,FAC,144,73,30,1) C DO 75 I=1,144 DO 75 J=2,72 75 GLOB(I,J)=V(I,J)/(ALAT(J)*RA) C DO 755 I=1,144 GLOB(I,1)=0. 755 GLOB(I,73)=0. C CALL SPHERT(2,GLOB,WAVE2,-1,FAC,144,73,30,1) C DO 961 I=1,1922 WAVE3(I)=WAVE2(I)-WAVE(I) 961 CONTINUE C C CALL SPHERT(1,PSI,WAVE3,1,FFF,144,73,30,1) C C TRANSFORM TO SPHERICAL CORR FOR U AND V C DO 63 I=1,144 DO 63 J=1,73 GLOB(I,J)=V(I,J)*ALAT(J)/RA 63 CONTINUE C DO 93 I=1,144 GLOB(I,1)=0. GLOB(I,73)=0. 93 CONTINUE CALL SPHERT(3,GLOB,WAVE,0,FAC,144,73,30,1) DO 85 I=1,144 DO 85 J=2,72 85 GLOB(I,J)=U(I,J)/(ALAT(J)*RA) DO 855 I=1,144 GLOB(I,1)=0. 855 GLOB(I,73)=0. CALL SPHERT(2,GLOB,WAVE2,-1,FAC,144,73,30,1) DO 86 I=1,1922 WAVE3(I)=WAVE2(I)+WAVE(I) WAVE2(I)=WAVE3(I) 86 CONTINUE C C TRANSFORM BACK TO DIVERGENCE CALL SPHERT(1,DIV,WAVE3,1,FS,144,73,30,1) C C TRANSFORM BACK TO CHI CALL SPHERT(1,GG,WAVE2,1,FFF,144,73,30,1) CALL SPHERT(1,CHI,WAVE3,1,FFF,144,73,30,1) C RETURN END C KINGTSE MO ROUTINES TO ASSIST COMPUTATION OF STRMFN & VEL POT C USING SPHERICAL HARMONICS C C USED BY PGM: MICSTRSP IN END-OF-MONTH GRAPHICS JOBS C C JANOWIAK 11/23/90 C SUBROUTINE SPHERT 1(JDR,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 103 ...... AS IDIR= 3 BUT FOR EQUIDISTANT LAT/LON GRID C 104 ...... AS IDIR= 4 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=32,LIMGAU=80,LROMB=1, 1 LIMWVX=(LIMWV+1)*(LIMWV+2)/2*(1-LROMB)+ 2 (LIMWV+1)*(LIMWV+1)*LROMB, 3 LIMGH =(LIMGAU+1)/2,LIMGHP=LIMGH+1, 4 MAXPNM=LIMWVX*LIMGH, 5 MAXWK =LIMWVX*LIMGHP*2) C IF FOR T40,CHANGE LIMWV=40,LIMGAU=80,LROMB=0 C IF FOR R30,CHANGE LIMWV=32, LIMGAU=80, LROMB=1 DIMENSION GRID(1),WAVE(1) DIMENSION FAC(1),IDC(3) C COMMON/SWK/ WORK(MAXWK) COMMON/WAVY/PNMC(MAXPNM,3),PVV(MAXPNM) DIMENSION SYMM(LIMWVX) DIMENSION COSCLT(LIMGAU) ,PNM(MAXPNM) DIMENSION IFAX(10),DTRIGS(1000) C DIMENSION DDW(LIMGAU,6) C DATA MFP1/0/,IRMB1/-999/ DATA MFP2/0/,IRMB2/-999/ DATA IFP/0/,JFP/0/,IDFP/0/ DATA IXTRA/2/ DATA ICALL/0/ IDC(1)=-101 IDC(2)=101 IDC(3)=104 IDIR=IDC(JDR) 8001 FORMAT(5X,'IMAX,JMAX,MAXWV,IROMB=',4I8) C IMAXP=IMAX+IXTRA JMAXHF=(JMAX+1)/2 IJHMAX=IMAX*JMAXHF IPJHMX=IMAXP*JMAXHF IJMAX=IMAX*JMAX IPJMAX=IMAXP*JMAX MEND1=MAXWV+1 MEND1T=MEND1*2 NMMAX=(MEND1+1)*MEND1/2 IF(IROMB.EQ.1) NMMAX=MEND1*MEND1 NMMAXT=NMMAX*2 C 8002 FORMAT(5X,'IMAXP,JMAXHF,IJHMAX,IPJHMX,IJMAX,IPJMAX',6I8) 8003 FORMAT(5X,'MEND1,MEND1T,NMMAX,NMMAXT',4I8) C MAXREQ=NMMAX*JMAXHF IF(MAXREQ.GT.MAXPNM) THEN WRITE(6,8005) 8005 FORMAT(5X,'INCREASE LIMWV OR LIMGAU OF SUBROUTINE SPHERT') STOP 'ABEND' ENDIF C MAXREQ=NMMAX*(JMAXHF+1)*2 IF(MAXREQ.GT.MAXWK) THEN WRITE(6,8006) 8006 FORMAT(5X,'INCREASE LIMWV OR LIMGAU OF SUBROUTINE SPHERT') STOP 'ABEND' ENDIF C IF(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(NMMAX*2.GT.IMAXP*JMAX) THEN IXREQ=FLOAT(NMMAX*2)/FLOAT(JMAX)-FLOAT(IMAX)+1 WRITE(6,8008) IXREQ 8008 FORMAT(5X,'INCREASE IXTRA OF SUBROUTINE SPHERT TO ',I8) STOP 'ABEND' ENDIF C IF(IFP.EQ.IMAX) GO TO 1 IMODE=2 IF(MOD(IMAX,2).NE.0) IMODE=1 CALL SFAX(IFAX,IMAX,IMODE) CALL SFTRIG(DTRIGS,IMAX,IMODE) IFP=IMAX 1 CONTINUE 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 IF(IDFP.EQ.IDIR.AND.JFP.EQ.JMAX.AND.MFP2.EQ.MAXWV 1 .AND.IRMB2.EQ.IROMB) GO TO 4 DO 82 K=1,MAXPNM PNM(K)=PNMC(K,JDR) IF(MAXWV.EQ.31) PNM(K)=PVV(K) 82 CONTINUE IRMB2=IROMB JFP=JMAX MFP2=MAXWV IDFP=IDIR 4 CONTINUE C 9 CONTINUE JKEND=JMAX IPJKED=IPJMAX JMH=JMAX/2 C C END OF CONSTANTS C IF(IDIR.LT.0) GO TO 1000 C C FOREWARD TRANSFORM (GRID TO WAVE) C C GLOBAL CASE DO 100 K=1,1 IX1=IJMAX *(K-1) IY1=IPJMAX*(K-1) DO 100 J=1,JMAXHF IX2 =IMAX *(J-1)+IX1 IY2 =IMAXP*(J-1)+IY1 IX2R=IMAX*(JMAX-J)+IX1 IY2H=IMAXP*(J-1)+IY1+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(MOD(IMAX,2).EQ.0) THEN CALL SFT991(WORK,WORK(IPJKED+1),DTRIGS,IFAX,1,IMAXP,IMAX,JKEND,-1, 1 DDW(1,1),DDW(1,2),DDW(1,3),DDW(1,4),DDW(1,5),DDW(1,6)) ELSE CALL FOURT(WORK,IMAX,1,-1,0,WORK(IPJKED+1)) ENDIF C 115 CONTINUE IDPNM=MOD(IDIR,10) C DO 110 K=1,1 C IX1=IPJMAX*(K-1)-IMAXP IY1=NMMAXT*(K-1) NM=0 DO 160 MM=1,MEND1 IX2=2*(MM-1)+IX1 NEND1=MEND1 IF(IROMB.EQ.1) NEND1=MM+MEND1-1 DO 160 NN=MM,NEND1 NM=NM+1 IY2=2*(NM-1)+IY1+IPJKED IZ1=JMAXHF*(NM-1) MODNM=MOD(NN-MM,2) IF(IDPNM.EQ.2.OR.IDPNM.EQ.4) MODNM=1-MODNM IF(MODNM.EQ.1) GO TO 180 C C SYMMETRIC PART DO 170 II=1,2 IX3=IX2+II WORK(IY2+II)=0. DO 170 J=1,JMAXHF WORK(IY2+II)=WORK(IY2+II)+WORK(IMAXP*J+IX3)*PNM(IZ1+J) 170 CONTINUE GO TO 160 C C ANTISYMMETRIC PART 180 CONTINUE DO 190 II=1,2 IX3=IX2+II+IPJHMX WORK(IY2+II)=0. DO 190 J=1,JMH WORK(IY2+II)=WORK(IY2+II)+WORK(IMAXP*J+IX3)*PNM(IZ1+J) 190 CONTINUE 160 CONTINUE C 110 CONTINUE C NMKEND=NMMAXT IF(MLT.NE.0) GO TO 121 DO 120 NMK=1,NMKEND WAVE(NMK)=WORK(IPJKED+NMK) 120 CONTINUE RETURN C 121 CONTINUE IF(MLT.EQ.-1) GO TO 122 DO 123 K=1,1 IX1=NMMAXT*(K-1)+IPJKED DO 123 NM=1,NMMAX WAVE(NM*2-1)=WORK(IX1+NM*2-1)*FAC(NM) WAVE(NM*2)=WORK(IX1+NM*2)*FAC(NM) 123 CONTINUE RETURN C 122 CONTINUE C MULTIPLY I*FAC DO 124 K=1,1 IX1=NMMAXT*(K-1)+IPJKED C REAL PART DO 126 NM=1,NMMAX WAVE(2*NM-1)=-FAC(NM)*WORK(2*NM+IX1) 126 CONTINUE C IMAGINARY PART DO 127 NM=1,NMMAX WAVE(2*NM )= FAC(NM)*WORK(2*NM+IX1-1) 127 CONTINUE 124 CONTINUE RETURN C C END OF GLOBAL CASE FOREWARD TRANSFORM C 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,IPJKED WORK(IJK)=0. 205 CONTINUE C DO 200 K=1,1 C IX1=IPJMAX*(K-1) IY1=NMFT*(K-1) DO 250 J=1,JMAXHF IZ1=NMMAX*(J-1) IX2 =IX1+IMAXP*(J-1) IX2R=IX1+IMAXP*(JMAX-J) DO 250 II=1,2 IX3 =IX2 +II-2 IX3R=IX2R+II-2 IY2 =IY1 +II-2 C IF(MLT.NE.0) GO TO 261 DO 260 NM=1,NMF WORK(IPJKED+NM)=PNM(IZ1+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(IPJKED+NM)=PNM(IZ1+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(IPJKED+NM)=PNM(IZ1+NM)*(-WAVE(NM*2+IY2+1)*FAC(NM)) 264 CONTINUE GO TO 269 C IMAGINARY PART 266 CONTINUE DO 267 NM=1,NMF WORK(IPJKED+NM)=PNM(IZ1+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(IPJKED+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(IPJKED+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(IPJKED+NM)*SYMM(NM) 490 CONTINUE 495 CONTINUE C 470 CONTINUE C 250 CONTINUE 200 CONTINUE C CALL SFT991(WORK,WORK(IPJKED+1),DTRIGS,IFAX,1,IMAXP,IMAX,JKEND,+1, 1 DDW(1,1),DDW(1,2),DDW(1,3),DDW(1,4),DDW(1,5),DDW(1,6)) C DO 290 K=1,1 IX1=IJMAX*(K-1) IY1=IPJMAX*(K-1) DO 290 J=1,JMAX IX2=IMAX *(J-1)+IX1 IY2=IMAXP*(J-1)+IY1 DO 290 I=1,IMAX GRID(IX2+I)=WORK(IY2+I) 290 CONTINUE C 8888 FORMAT(20A4) RETURN END C DEV. DIV. ROUTINES TO ASSIST COMPUTATION OF STRMFN & VEL POT C USING SPHERICAL HARMONICS C C USED BY PGM: MICSTRSP IN END-OF-MONTH GRAPHICS JOBS C C JANOWIAK 11/23/90 C 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 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 C DEV. DIV. ROUTINES TO ASSIST COMPUTATION OF STRMFN & VEL POT 00000010 C USING SPHERICAL HARMONICS 00000020 C 00000030 C USED BY PGM: MICSTRSP IN END-OF-MONTH GRAPHICS JOBS 00000040 C 00000050 C JANOWIAK 11/23/90 00000060 C 00000070 SUBROUTINE FOURT(DATA,NN,NDIM,ISIGN,IFORM,WORK) 00000080 C 00000090 DIMENSION DATA(1),NN(1),IFACT(32),WORK(1) 00000100 TWOPI=6.283185307 00000110 RTHLF=.7071067812 00000120 NP0=0 00000130 NPREV=0 00000140 IF(NDIM-1)920,1,1 00000150 1 NTOT=2 00000160 DO 2 IDIM=1,NDIM 00000170 IF(NN(IDIM))920,920,2 00000180 2 NTOT=NTOT*NN(IDIM) 00000190 C 00000200 C MAIN LOOP FOR EACH DIMENSION 00000210 C 00000220 NP1=2 00000230 DO 910 IDIM=1,NDIM 00000240 N=NN(IDIM) 00000250 NP2=NP1*N 00000260 IF(N-1)920,900,5 00000270 C 00000280 C IS N A POWER OF TWO AND IF NOT, WHAT ARE ITS FACTORS 00000290 C 00000300 5 M=N 00000310 NTWO=NP1 00000320 IF=1 00000330 IDIV=2 00000340 10 IQUOT=M/IDIV 00000350 IREM=M-IDIV*IQUOT 00000360 IF(IQUOT-IDIV)50,11,11 00000370 11 IF(IREM)20,12,20 00000380 12 NTWO=NTWO+NTWO 00000390 IFACT(IF)=IDIV 00000400 IF=IF+1 00000410 M=IQUOT 00000420 GO TO 10 00000430 20 IDIV=3 00000440 INON2=IF 00000450 30 IQUOT=M/IDIV 00000460 IREM=M-IDIV*IQUOT 00000470 IF(IQUOT-IDIV)60,31,31 00000480 31 IF(IREM)40,32,40 00000490 32 IFACT(IF)=IDIV 00000500 IF=IF+1 00000510 M=IQUOT 00000520 GO TO 30 00000530 40 IDIV=IDIV+2 00000540 GO TO 30 00000550 50 INON2=IF 00000560 IF(IREM)60,51,60 00000570 51 NTWO=NTWO+NTWO 00000580 GO TO 70 00000590 60 IFACT(IF)=M 00000600 70 NON2P=NP2/NTWO 00000610 C 00000620 C SEPARATE FOUR CASES-- 00000630 C 1. COMPLEX TRANSFORM 00000640 C 2. REAL TRANSFORM FOR THE 2ND, 3RD, ETC. DIMENSION. METHOD-- 00000650 C TRANSFORM HALF THE DATA, SUPPLYING THE OTHER HALF BY CON- 00000660 C JUGATE SYMMETRY. 00000670 C 3. REAL TRANSFORM FOR THE 1ST DIMENSION, N ODD. METHOD-- 00000680 C SET THE IMAGINARY PARTS TO ZERO. 00000690 C 4. REAL TRANSFORM FOR THE 1ST DIMENSION, N EVEN. METHOD-- 00000700 C TRANSFORM A COMPLEX ARRAY OF LENGTH N/2 WHOSE REAL PARTS 00000710 C ARE THE EVEN NUMBERED REAL VALUES AND WHOSE IMAGINARY PARTS 00000720 C ARE THE ODD NUMBERED REAL VALUES. SEPARATE AND SUPPLY 00000730 C THE SECOND HALF BY CONJUGATE SYMMETRY. 00000740 C 00000750 IFMIN = 1 00000760 I1RNG = NP1 00000770 IF (IFORM .LE. 0 .AND. IDIM .LT. 4) GO TO 71 00000780 ICASE = 1 00000790 GO TO 100 00000800 C 00000810 71 IF (IDIM .LE. 1) GO TO 72 00000820 ICASE = 2 00000830 I1RNG = NP0 * (1 + NPREV / 2) 00000840 GO TO 100 00000850 C 00000860 72 IF (NTWO .GT. NP1) GO TO 73 00000870 ICASE = 3 00000880 GO TO 100 00000890 C 00000900 73 ICASE = 4 00000910 IFMIN = 2 00000920 NTWO = NTWO / 2 00000930 N = N / 2 00000940 NP2 = NP2 / 2 00000950 NTOT = NTOT / 2 00000960 C 00000970 I = - 1 00000980 DO 80 J = 1, NTOT 00000990 I = I + 2 00001000 DATA(J) = DATA(I) 00001010 80 CONTINUE 00001020 C 00001030 C SHUFFLE DATA BY BIT REVERSAL, SINCE N=2**K. AS THE SHUFFLING 00001040 C CAN BE DONE BY SIMPLE INTERCHANGE, NO WORKING ARRAY IS NEEDED 00001050 C 00001060 100 IF(NON2P-1)101,101,200 00001070 101 NP2HF=NP2/2 00001080 J=1 00001090 DO 150 I2=1,NP2,NP1 00001100 IF(J-I2)121,130,130 00001110 121 I1MAX=I2+NP1-2 00001120 DO 125 I1=I2,I1MAX,2 00001130 DO 125 I3=I1,NTOT,NP2 00001140 J3=J+I3-I2 00001150 TEMPR=DATA(I3) 00001160 TEMPI=DATA(I3+1) 00001170 DATA(I3)=DATA(J3) 00001180 DATA(I3+1)=DATA(J3+1) 00001190 DATA(J3)=TEMPR 00001200 125 DATA(J3+1)=TEMPI 00001210 130 M=NP2HF 00001220 140 IF(J-M)150,150,141 00001230 141 J=J-M 00001240 M=M/2 00001250 IF(M-NP1)150,140,140 00001260 150 J=J+M 00001270 GO TO 300 00001280 C 00001290 C SHUFFLE DATA BY DIGIT REVERSAL FOR GENERAL N 00001300 C 00001310 200 NWORK=2*N 00001320 DO 270 I1=1,NP1,2 00001330 DO 270 I3=I1,NTOT,NP2 00001340 J=I3 00001350 DO 260 I=1,NWORK,2 00001360 IF(ICASE-3)210,220,210 00001370 210 WORK(I)=DATA(J) 00001380 WORK(I+1)=DATA(J+1) 00001390 GO TO 240 00001400 220 WORK(I)=DATA(J) 00001410 WORK(I+1)=0. 00001420 240 IFP2=NP2 00001430 IF=IFMIN 00001440 250 IFP1=IFP2/IFACT(IF) 00001450 J=J+IFP1 00001460 IF(J-I3-IFP2)260,255,255 00001470 255 J=J-IFP2 00001480 IFP2=IFP1 00001490 IF=IF+1 00001500 IF(IFP2-NP1)260,260,250 00001510 260 CONTINUE 00001520 I2MAX=I3+NP2-NP1 00001530 I=1 00001540 DO 270 I2=I3,I2MAX,NP1 00001550 DATA(I2)=WORK(I) 00001560 DATA(I2+1)=WORK(I+1) 00001570 270 I=I+2 00001580 C 00001590 C MAIN LOOP FOR FACTORS OF TWO. 00001600 C W=EXP(ISIGN*2*PI*SQRT(-1)*M/(4*MMAX)). CHECK FOR W=ISIGN*SQRT(-1)00001610 C AND REPEAT FOR W=W*(1+ISIGN*SQRT(-1))/SQRT(2). 00001620 C 00001630 300 IF(NTWO-NP1)600,600,305 00001640 305 NP1TW=NP1+NP1 00001650 IPAR=NTWO/NP1 00001660 310 IF(IPAR-2)350,330,320 00001670 320 IPAR=IPAR/4 00001680 GO TO 310 00001690 330 DO 340 I1=1,I1RNG,2 00001700 DO 340 K1=I1,NTOT,NP1TW 00001710 K2=K1+NP1 00001720 TEMPR=DATA(K2) 00001730 TEMPI=DATA(K2+1) 00001740 DATA(K2)=DATA(K1)-TEMPR 00001750 DATA(K2+1)=DATA(K1+1)-TEMPI 00001760 DATA(K1)=DATA(K1)+TEMPR 00001770 340 DATA(K1+1)=DATA(K1+1)+TEMPI 00001780 350 MMAX=NP1 00001790 360 IF(MMAX-NTWO/2)370,600,600 00001800 370 LMAX=MAX0(NP1TW,MMAX/2) 00001810 DO 570 L=NP1,LMAX,NP1TW 00001820 M=L 00001830 IF(MMAX-NP1)420,420,380 00001840 380 THETA=-TWOPI*FLOAT(L)/FLOAT(4*MMAX) 00001850 IF(ISIGN)400,390,390 00001860 390 THETA=-THETA 00001870 400 WR=COS(THETA) 00001880 WI=SIN(THETA) 00001890 410 W2R=WR*WR-WI*WI 00001900 W2I=2.*WR*WI 00001910 W3R=W2R*WR-W2I*WI 00001920 W3I=W2R*WI+W2I*WR 00001930 420 DO 530 I1=1,I1RNG,2 00001940 KMIN=I1+IPAR*M 00001950 IF(MMAX-NP1)430,430,440 00001960 430 KMIN=I1 00001970 440 KDIF=IPAR*MMAX 00001980 450 KSTEP=4*KDIF 00001990 IF(KSTEP-NTWO)460,460,530 00002000 460 DO 520 K1=KMIN,NTOT,KSTEP 00002010 K2=K1+KDIF 00002020 K3=K2+KDIF 00002030 K4=K3+KDIF 00002040 IF(MMAX-NP1)470,470,480 00002050 470 U1R=DATA(K1)+DATA(K2) 00002060 U1I=DATA(K1+1)+DATA(K2+1) 00002070 U2R=DATA(K3)+DATA(K4) 00002080 U2I=DATA(K3+1)+DATA(K4+1) 00002090 U3R=DATA(K1)-DATA(K2) 00002100 U3I=DATA(K1+1)-DATA(K2+1) 00002110 IF(ISIGN)471,472,472 00002120 471 U4R=DATA(K3+1)-DATA(K4+1) 00002130 U4I=DATA(K4)-DATA(K3) 00002140 GO TO 510 00002150 472 U4R=DATA(K4+1)-DATA(K3+1) 00002160 U4I=DATA(K3)-DATA(K4) 00002170 GO TO 510 00002180 480 T2R=W2R*DATA(K2)-W2I*DATA(K2+1) 00002190 T2I=W2R*DATA(K2+1)+W2I*DATA(K2) 00002200 T3R=WR*DATA(K3)-WI*DATA(K3+1) 00002210 T3I=WR*DATA(K3+1)+WI*DATA(K3) 00002220 T4R=W3R*DATA(K4)-W3I*DATA(K4+1) 00002230 T4I=W3R*DATA(K4+1)+W3I*DATA(K4) 00002240 U1R=DATA(K1)+T2R 00002250 U1I=DATA(K1+1)+T2I 00002260 U2R=T3R+T4R 00002270 U2I=T3I+T4I 00002280 U3R=DATA(K1)-T2R 00002290 U3I=DATA(K1+1)-T2I 00002300 IF(ISIGN)490,500,500 00002310 490 U4R=T3I-T4I 00002320 U4I=T4R-T3R 00002330 GO TO 510 00002340 500 U4R=T4I-T3I 00002350 U4I=T3R-T4R 00002360 510 DATA(K1)=U1R+U2R 00002370 DATA(K1+1)=U1I+U2I 00002380 DATA(K2)=U3R+U4R 00002390 DATA(K2+1)=U3I+U4I 00002400 DATA(K3)=U1R-U2R 00002410 DATA(K3+1)=U1I-U2I 00002420 DATA(K4)=U3R-U4R 00002430 520 DATA(K4+1)=U3I-U4I 00002440 KDIF=KSTEP 00002450 KMIN=4*(KMIN-I1)+I1 00002460 GO TO 450 00002470 530 CONTINUE 00002480 M=M+LMAX 00002490 IF(M-MMAX)540,540,570 00002500 540 IF(ISIGN)550,560,560 00002510 550 TEMPR=WR 00002520 WR=(WR+WI)*RTHLF 00002530 WI=(WI-TEMPR)*RTHLF 00002540 GO TO 410 00002550 560 TEMPR=WR 00002560 WR=(WR-WI)*RTHLF 00002570 WI=(TEMPR+WI)*RTHLF 00002580 GO TO 410 00002590 570 CONTINUE 00002600 IPAR=3-IPAR 00002610 MMAX=MMAX+MMAX 00002620 GO TO 360 00002630 C 00002640 C MAIN LOOP FOR FACTORS NOT EQUAL TO TWO. 00002650 C W=EXP(ISIGN*2*PI*SQRT(-1)*(J1+J2-I3-1)/IFP2) 00002660 C 00002670 600 IF(NON2P-1)700,700,601 00002680 601 IFP1=NTWO 00002690 IF=INON2 00002700 610 IFP2=IFACT(IF)*IFP1 00002710 THETA=-TWOPI/FLOAT(IFACT(IF)) 00002720 IF(ISIGN)612,611,611 00002730 611 THETA=-THETA 00002740 612 WSTPR=COS(THETA) 00002750 WSTPI=SIN(THETA) 00002760 DO 650 J1=1,IFP1,NP1 00002770 THETM=-TWOPI*FLOAT(J1-1)/FLOAT(IFP2) 00002780 IF(ISIGN)614,613,613 00002790 613 THETM=-THETM 00002800 614 WMINR=COS(THETM) 00002810 WMINI=SIN(THETM) 00002820 I1MAX=J1+I1RNG-2 00002830 DO 650 I1=J1,I1MAX,2 00002840 DO 650 I3=I1,NTOT,NP2 00002850 I=1 00002860 WR=WMINR 00002870 WI=WMINI 00002880 J2MAX=I3+IFP2-IFP1 00002890 DO 640 J2=I3,J2MAX,IFP1 00002900 TWOWR=WR+WR 00002910 J3MAX=J2+NP2-IFP2 00002920 DO 630 J3=J2,J3MAX,IFP2 00002930 JMIN=J3-J2+I3 00002940 J=JMIN+IFP2-IFP1 00002950 SR=DATA(J) 00002960 SI=DATA(J+1) 00002970 OLDSR=0. 00002980 OLDSI=0. 00002990 J=J-IFP1 00003000 620 STMPR=SR 00003010 STMPI=SI 00003020 SR=TWOWR*SR-OLDSR+DATA(J) 00003030 SI=TWOWR*SI-OLDSI+DATA(J+1) 00003040 OLDSR=STMPR 00003050 OLDSI=STMPI 00003060 J=J-IFP1 00003070 IF(J-JMIN)621,621,620 00003080 621 WORK(I)=WR*SR-WI*SI-OLDSR+DATA(J) 00003090 WORK(I+1)=WI*SR+WR*SI-OLDSI+DATA(J+1) 00003100 630 I=I+2 00003110 WTEMP=WR*WSTPI 00003120 WR=WR*WSTPR-WI*WSTPI 00003130 640 WI=WI*WSTPR+WTEMP 00003140 I=1 00003150 DO 650 J2=I3,J2MAX,IFP1 00003160 J3MAX=J2+NP2-IFP2 00003170 DO 650 J3=J2,J3MAX,IFP2 00003180 DATA(J3)=WORK(I) 00003190 DATA(J3+1)=WORK(I+1) 00003200 650 I=I+2 00003210 IF=IF+1 00003220 IFP1=IFP2 00003230 IF(IFP1-NP2)610,700,700 00003240 C 00003250 C COMPLETE A REAL TRANSFORM IN THE 1ST DIMENSION, N EVEN, BY CON- 00003260 C JUGATE SYMMETRIES. 00003270 C 00003280 700 GO TO (900,800,900,701),ICASE 00003290 701 NHALF=N 00003300 N=N+N 00003310 THETA=-TWOPI/FLOAT(N) 00003320 IF(ISIGN)703,702,702 00003330 702 THETA=-THETA 00003340 703 WSTPR=COS(THETA) 00003350 WSTPI=SIN(THETA) 00003360 WR=WSTPR 00003370 WI=WSTPI 00003380 IMIN=3 00003390 JMIN=2*NHALF-1 00003400 GO TO 725 00003410 710 J=JMIN 00003420 DO 720 I=IMIN,NTOT,NP2 00003430 SUMR=(DATA(I)+DATA(J))/2. 00003440 SUMI=(DATA(I+1)+DATA(J+1))/2. 00003450 DIFR=(DATA(I)-DATA(J))/2. 00003460 DIFI=(DATA(I+1)-DATA(J+1))/2. 00003470 TEMPR=WR*SUMI+WI*DIFR 00003480 TEMPI=WI*SUMI-WR*DIFR 00003490 DATA(I)=SUMR+TEMPR 00003500 DATA(I+1)=DIFI+TEMPI 00003510 DATA(J)=SUMR-TEMPR 00003520 DATA(J+1)=-DIFI+TEMPI 00003530 720 J=J+NP2 00003540 IMIN=IMIN+2 00003550 JMIN=JMIN-2 00003560 WTEMP=WR*WSTPI 00003570 WR=WR*WSTPR-WI*WSTPI 00003580 WI=WI*WSTPR+WTEMP 00003590 725 IF(IMIN-JMIN)710,730,740 00003600 730 IF(ISIGN)731,740,740 00003610 731 DO 735 I=IMIN,NTOT,NP2 00003620 735 DATA(I+1)=-DATA(I+1) 00003630 740 NP2=NP2+NP2 00003640 NTOT=NTOT+NTOT 00003650 J=NTOT+1 00003660 IMAX=NTOT/2+1 00003670 745 IMIN=IMAX-2*NHALF 00003680 I=IMIN 00003690 GO TO 755 00003700 750 DATA(J)=DATA(I) 00003710 DATA(J+1)=-DATA(I+1) 00003720 755 I=I+2 00003730 J=J-2 00003740 IF(I-IMAX)750,760,760 00003750 760 DATA(J)=DATA(IMIN)-DATA(IMIN+1) 00003760 DATA(J+1)=0. 00003770 IF(I-J)770,780,780 00003780 765 DATA(J)=DATA(I) 00003790 DATA(J+1)=DATA(I+1) 00003800 770 I=I-2 00003810 J=J-2 00003820 IF(I-IMIN)775,775,765 00003830 775 DATA(J)=DATA(IMIN)+DATA(IMIN+1) 00003840 DATA(J+1)=0. 00003850 IMAX=IMIN 00003860 GO TO 745 00003870 780 DATA(1)=DATA(1)+DATA(2) 00003880 DATA(2)=0. 00003890 GO TO 900 00003900 C 00003910 C COMPLETE A REAL TRANSFORM FOR THE 2ND, 3RD, ETC. DIMENSION BY 00003920 C CONJUGATE SYMMETRIES. 00003930 C 00003940 800 IF(I1RNG-NP1)805,900,900 00003950 805 DO 860 I3=1,NTOT,NP2 00003960 I2MAX=I3+NP2-NP1 00003970 DO 860 I2=I3,I2MAX,NP1 00003980 IMAX=I2+NP1-2 00003990 IMIN=I2+I1RNG 00004000 JMAX=2*I3+NP1-IMIN 00004010 IF(I2-I3)820,820,810 00004020 810 JMAX=JMAX+NP2 00004030 820 IF(IDIM-2)850,850,830 00004040 830 J=JMAX+NP0 00004050 DO 840 I=IMIN,IMAX,2 00004060 DATA(I)=DATA(J) 00004070 DATA(I+1)=-DATA(J+1) 00004080 840 J=J-2 00004090 850 J=JMAX 00004100 DO 860 I=IMIN,IMAX,NP0 00004110 DATA(I)=DATA(J) 00004120 DATA(I+1)=-DATA(J+1) 00004130 860 J=J-NP0 00004140 C 00004150 C END OF LOOP ON EACH DIMENSION 00004160 C 00004170 900 NP0=NP1 00004180 NP1=NP2 00004190 910 NPREV=N 00004200 920 IF(ISIGN.EQ.1) RETURN 00004210 NTOT=2 00004220 DO 930 IDIM=1,NDIM 00004230 930 NTOT=NTOT*NN(IDIM) 00004240 NTOTHF=NTOT/2 00004250 DO 940 ITOT=1,NTOT 00004260 940 DATA(ITOT)=DATA(ITOT)/NTOTHF 00004270 RETURN 00004280 END 00004290 c c c subroutine updctl(ioctl,imon,iyr) c c Updates 'TDEF' parameter in '.ctl' file for GrADS c character*80 ctl(22) character*3 cmon(12) character*4 icyr c data cmon /'jan', 'feb', 'mar', 'apr', 'may', 'jun', 1 'jul', 'aug', 'sep', 'oct', 'nov', 'dec'/ c write(icyr,'(i4)') iyr c open (ioctl,file= 1 '/export-8/cacsrv1/cpc/anal/cdas/bulletin/ctl/r30_cdas.ctl') do line=1,22 read(ioctl,'(a80)') ctl(line) enddo rewind ioctl c c c.... Modify "ctl" files c c ctl(14)(15:17)=cmon(imon) ctl(14)(18:21)=icyr c write(ioctl,'(a80)') ctl c write(*,*) '".ctl" files updated as follows:' write(*,105) ctl 105 format(a80) return end