! WD23JP ! ***************** ! * SOCT1.FOR * ! * PURSER 1996 * ! ***************** ! !------------------------------------------------------------------------------C ! R.J.Purser, NCEP, Washington D.C. 1996 C ! SUBROUTINE INSOCT1 C ! Get conformal solution, then refine it to get "smoothed" grid solution C ! If using multigrid method, set parameter NGEN to positive number of C ! "generations" of grids. Otherwise, revert to SOR and ensure that enough C ! iterations (NITSOR), independent of the m.g. algorithm, are given. C ! C ! --> HOM homogeneity parameter C ! --- RA,QA work arrays of size NFT C ! --> NFT suitable period (eg 128) for FFT C ! --- JUMBLE integer work arrays of size NFT C ! --- WORK real work array of size NFT C ! --> TFILE name of file into which table is put (CHARACTER*39) C ! --> NCY no. of multigrid (m.g.) cycles (eg 15 for m.g., 0 for S.O.R.) C ! --> NIT number of SOR iterations at each stage of m.g. (eg 3) C ! --> NITSOR number of final SOR iterations (10 for m.g., 1000 for SOR) C !------------------------------------------------------------------------------C SUBROUTINE INSOCT1(HOM,RA,QA,NFT,JUMBLE,WORK,TFILE *,NCY,NIT,NITSOR) PARAMETER(NGEN=6,NRAT1=1,NRAT2=3 *,NXD=9*(NRAT2**2*(1+4*(4**NGEN-1)/3)+6*NRAT2*(1+2*(2**NGEN-1)) * +9*(NGEN+1))) PARAMETER(N=NRAT2*2**NGEN,NP=N+1) COMPLEX CI,CIR,CIQ CHARACTER*39 TFILE COMMON/QPAN/X(3,-1:NP,-1:NP) COMMON/CSTOCT1/A1(40),A2(40),B1(40),B2(40),A,B,CX,CY *,SSR,N1,ST,STI,CI,CIR,CIQ DIMENSION RA(0:*),QA(0:*),JUMBLE(NFT),WORK(NFT) !jim DIMENSION RA(0:*),QA(0:*),JUMBLE(N),WORK(N) DIMENSION XD(NXD) A=FLOAT(NRAT1)/NRAT2 N1=20 CALL INDOCT(RA,QA,NFT,JUMBLE,WORK) CALL IROCTA(XD,NIT,NITSOR,NCY,HOM,TFILE) RETURN END !------------------------------------------------------------------------------C ! R.J.Purser, NCEP, Washington D.C. 1996 C ! SUBROUTINE INDOCT C ! Initialize the Taylor series for the conformal octagon solution C ! C ! --- RA,QA arrays of size N used for real & imag components of FFT C ! --> N suitable size (eg N=128) for one cycle of data in the FFT C ! --- JUMBLE integer work array of size N for FFT C ! --- WORK real work array of size N for FFT C !------------------------------------------------------------------------------C SUBROUTINE INDOCT(RA,QA,N,JUMBLE,WORK) COMPLEX CI,CIR,CIQ,CIA,CIAP,CIRS2,Z,W,CANG,CDANG1,CDANG2 COMMON/CSTOCT1/A1(40),A2(40),B1(40),B2(40),A,B,CX,CY *,SSR,N1,ST,STI,CI,CIR,CIQ DIMENSION RA(0:*),QA(0:*),JUMBLE(N),WORK(N) ST=4./3. STI=3./4. CI=(0.,1.) CIR=(-CI)**.75 CIQ=CI**.25 NH=N/2 R2=SQRT(2.) ORC=.5 SW1=SQRT(1.+A*A) SW2=AMIN1(A*2,R2*(1.-A)) SSR=(SW1/SW2)**2 RS1=.95*SW1 RS2=.95*SW2 RS3=RS1**4 RS4=RS2**ST DO I=1,N1 A1(I)=0. ENDDO A1(1)= 1.05060 A1(2)= -.172629 A1(3)= .125716 A1(4)= .00015657 A1(5)= -.0017866 A1(6)= -.0031168 DO I=1,N1 A2(I)=0. ENDDO A2(1)= .688411 A2(2)=-.180924 A2(3)=-.186593 A2(4)= .024811 A2(5)=-.044888 A2(6)=-.049190 B=-.147179 JUMBLE(1)=0 PIO4=ATAN(1.) CDANG1=CI*PIO4/NH CDANG2=-3.*CDANG1 CIA=CI*A CIAP=CIA+1. CIRS2=CI*RS2 DO KIT=1,50 DO I=0,NH CANG=CDANG1*I Z=RS1*CEXP(CANG) !mish Z=RS1*EXP(CANG) CALL TOCT(Z,W) W=W**4 U=REAL(W) V=AIMAG(W) RA(I)=U QA(I)=-V RA(N-I)=U QA(N-I)=V ENDDO QA(0)=0. QA(NH)=0. CALL CFFT(RA,QA,N,1.,WORK,JUMBLE) RS3P=1. DO I=1,N1 RS3P=RS3P*RS3 A1(I)=A1(I) +ORC*(RA(I)/RS3P-A1(I)) ENDDO DO I=0,NH CANG=CDANG2*I Z=CIAP-CIRS2*CEXP(CANG) CALL TOCT(Z,W) W=CI*(W-1.)/(W+1.) U=REAL(W) V=AIMAG(W) RA(I)=U QA(I)=V RA(N-I)=U QA(N-I)=-V ENDDO QA(0)=0. QA(NH)=0. CALL CFFT(RA,QA,N,1.,WORK,JUMBLE) B=B+ORC*(RA(0)-B) RS4P=1. DO I=1,N1 RS4P=RS4P*RS4 A2(I)=A2(I) +ORC*(RA(I)/RS4P-A2(I)) ENDDO ENDDO W=B W=(W+CI)/(CI-W) CX=REAL(W) CY=AIMAG(W) CALL CINVRT(A1,B1,N1,RA) CALL CINVRT(A2,B2,N1,RA) PRINT'('' IMAGE OF OCTAGON-VERTEX ON UNIT-CIRCLE:'')' WRITE(6,62)CX,CY PRINT'('' TAYLOR COEFFICIENTS, SERIES A1,A2,B1,B2:'')' WRITE(6,63)B DO I=1,N1 WRITE(6,64)I,A1(I),A2(I),B1(I),B2(I) ENDDO 62 FORMAT(2(1X,E12.6)) 63 FORMAT(4X,'0',13X,4(1X,E12.6)) 64 FORMAT(1X,I4,4(1X,E12.6)) RETURN END !------------------------------------------------------------------------------C ! R.J.Purser, NCEP, Washington D.C. 1996 C ! SUBROUTINE IROCTA C ! Solve the variational problem for positions of reference grid of principal C ! quadrant of standard smoothed octagon for the given homogeneity parameter. C ! Solution method is either nonlinear multigrid (NGEN.GT.0, NCY.GT.0) or else C ! reduces to SOR when NGEN=0. C ! C ! --- XD double precision work array C ! --> NIT number of SOR iterations used at each stage of multigrid scheme C ! --> NITSOR number of SOR iterations used at the end to polish final solution C ! --> NCY number of multigrid "V-cycles" (ignored when NGEN=0) C ! --> A single precision homogeneity parameter in the variational problem C ! --> TFILE character*39 name of file containing final single precision table C !------------------------------------------------------------------------------C SUBROUTINE IROCTA(XD,NIT,NITSOR,NCY,A,TFILE) CHARACTER*39 TFILE REAL A,X PARAMETER(NGEN=6,NRAT1=1,NRAT2=3,N=NRAT2*2**NGEN,NP=N+1) COMMON/QPAN/X(3,-1:NP,-1:NP) DIMENSION XD(0:*),NG(0:10),NXG(0:10),AODDG(0:10) *,DG(0:10),LX(0:10+1),LF(0:10),LD(0:10) ORC=1.2 NG(0)=NRAT2 DO IGEN=1,NGEN IGENM=IGEN-1 NG(IGEN)=2*NG(IGENM) ENDDO LX(0)=0 DO IGEN=0,NGEN DG(IGEN)=1./NG(IGEN) AODDG(IGEN)=A/(8.*DG(IGEN)**2) NXG(IGEN)=3*(NG(IGEN)+3)**2 LF(IGEN)=LX(IGEN)+NXG(IGEN) LD(IGEN)=LF(IGEN)+NXG(IGEN) LX(IGEN+1)=LD(IGEN)+NXG(IGEN) ENDDO LXI=LX(NGEN) NGI=NG(NGEN) NGA=(NRAT1*NGI)/NRAT2 LFI=LF(NGEN) LDI=LD(NGEN) AODDGI=AODDG(NGEN) ! INITIALIZE FINE GRID ESTIMATE USING CONFORMAL SOLUTION: CALL INSUBP(XD(LXI),NGA,NGI) !misha ! goto 123 !misha ! PERFORM NIT PRELIMINARY ITERATIONS OF RELAXATION (SMOOTHING) ON FINE GRID: DO IT=1,NIT CALL RELSUBP(XD(LXI),XD(LFI),XD(LDI),NGA,NGI,ORC,AODDGI,ACC,0) ENDDO PRINT'('' LEVEL '',I1,'' ACC= '',E12.6)',NGEN,ACC ! PERFORM NCY CYCLES OF NONLINEAR MULTIGRID WITH NGEN REFINEMENTS: DO ICY=1,NCY ! GO FROM FINER TO COARSER GRID REPRESENTATIONS: DO IGEN=NGEN-1,0,-1 IGENP=IGEN+1 LXIP=LX(IGENP) LFIP=LF(IGENP) LDIP=LD(IGENP) NGIP=NG(IGENP) NGAP=(NRAT1*NGIP)/NRAT2 AODDGIP=AODDG(IGENP) LXI=LX(IGEN) LFI=LF(IGEN) LDI=LD(IGEN) NGI=NG(IGEN) NGA=(NRAT1*NGI)/NRAT2 NXGI=NXG(IGEN) AODDGI=AODDG(IGEN) ! EVALUATE FIELD OF RESIDUALS (D) AT FINER GRID: CALL RELSUBP(XD(LXIP),XD(LFIP),XD(LDIP),NGAP,NGIP,1.0 * ,AODDGIP,ACC,1) CALL ZERV(XD(LFI),NXGI) ! COMBINE FINER GRID RESIDUALS (D) INTO A COARSER GRID REPRESENTATION (F): CALL ZERV(XD(LFI),NXGI) CALL ADDX2X(XD(LDIP),XD(LFI),NGA,NGI) ! SET COARSER GRID ESTIMATE (X) TO BE CONSISTENT WITH FINER GRID ESTIMATE (X): CALL SETX2X(XD(LXIP),XD(LXI),NGA,NGI) ! SET COARSER GRID FORCING (F) TO MISMATCH BETWEEN FINER AND COARSER RESIDUALS: CALL RELSUBP(XD(LXI),XD(LFI),XD(LDI),NGA,NGI,1.0,AODDGI,ACC,2) ! PERFORM NIT RELAXATION ITERATIONS WITH THIS FORCING (F) AT COARSER GRID: DO IT=1,NIT CALL RELSUBP(XD(LXI),XD(LFI),XD(LDI),NGA,NGI,ORC,AODDGI,ACC,0) ENDDO PRINT'('' LEVEL '',I1,'' ACC= '',E12.6)',IGEN,ACC ENDDO ! GO FROM COARSER TO FINER GRID REPRESENTATIONS: DO IGEN=0,NGEN-1 IGENP=IGEN+1 LXIP=LX(IGENP) LFIP=LF(IGENP) LDIP=LD(IGENP) NGIP=NG(IGENP) NGAP=(NRAT1*NGIP)/NRAT2 AODDGIP=AODDG(IGENP) LXI=LX(IGEN) LFI=LF(IGEN) LDI=LD(IGEN) NGI=NG(IGEN) NGA=(NRAT1*NGI)/NRAT2 NXGI=NXG(IGEN) AODDGI=AODDG(IGEN) ! FIND DIFFERENCE (ON COARSER GRID) BETWEEN COARSER AND FINER GRID ESTIMATES X: CALL DIFXX2(XD(LXI),XD(LXIP),NGA,NGI) ! INTERPOLATE AND ADD THIS DIFFERENCE TO FINER GRID ESTIMATE (X): CALL ADDXX2(XD(LXI),XD(LXIP),NGA,NGI) ! PERFORM NIT RELAXATION ITERATIONS AT FINER GRID TO SMOOTH INTERPOLANTS: DO IT=1,NIT CALL RELSUBP(XD(LXIP),XD(LFIP),XD(LDIP),NGAP,NGIP,ORC * ,AODDGIP,ACC,0) ENDDO PRINT'('' LEVEL '',I1,'' ACC= '',E12.6)',IGENP,ACC ENDDO PRINT'('' CYCLE '',I4,'' COMPLETE'')',ICY ENDDO ! MULTIGRID CYCLES COMPLETE. NOW POLISH RESULT WITH EXTRA S.O.R. ITERATIONS: LXI=LX(NGEN) LFI=LF(NGEN) LDI=LD(NGEN) NGI=NG(NGEN) NGA=(NRAT1*NGI)/NRAT2 NXGI=NXG(NGEN) AODDGI=AODDG(NGEN) DO IT=1,NITSOR CALL RELSUBP(XD(LXI),XD(LFI),XD(LDI),NGA,NGI,ORC * ,AODDGI,ACC,0) PRINT'('' LEVEL '',I1,'' ACC= '',E12.6)',NGEN,ACC ENDDO ! COPY FINE GRID SOLUTION TO SINGLE PRECISION TABLE AND STORE FOR POSTERITY: !misha 123 CALL QTABLE(XD(LXI),X,NGA,NGI) ! CALL QTABLE(XD(LXI),X,NGA,NGI) !misha !jim OPEN(UNIT=9,FILE=TFILE,STATUS='UNKNOWN',FORM='BINARY') !misha OPEN(UNIT=9,FILE=TFILE,STATUS='UNKNOWN',FORM='unformatted') !misha WRITE(9)X CLOSE(UNIT=9) !test ! do i=-1,np ! do j=-1,np ! print *,x(1,i,j), x(2,i,j), x(3,i,j) ! end do ! end do !test ! TASK COMPLETE; NOW GRAB ANOTHER CUP OF COFFEE. RETURN END !------------------------------------------------------------------------------C ! R.J.Purser, NCEP, Washington D.C. 1996 C ! SUBROUTINE QTABLE C ! Symmetrize and copy results at finest grid of multigrid solution to a C ! single precision table of reference grid positions. The reference grid C ! covers one quadrant of the square containing the standard octagon, C ! plus an extra line of data points all round. Location (0,0) is the center C ! of the square containing the octogon, (N,N) is the corner of the quadrant C ! of the square, (N,NA) and (NA,N) are the two corners of the octagon in the C ! quadrant. C ! C ! <-> XD double precision array, symmetrized on output, of table data C ! <-- X single precision array of symmetrized table of cartesianlocations C ! --> NA number of grid steps between octagon corner and center line C ! --> N number of grid steps from center to side of square C !------------------------------------------------------------------------------C SUBROUTINE QTABLE(XD,X,NA,N) DIMENSION XD(3,-1:N+1,-1:N+1),X(3,-1:N+1,-1:N+1),LOFIE(3) *,JOFIF(3) DATA LOFIE/1,1,-1/,JOFIF/2,1,3/ ! SYMMETRIZE ACROSS E-EDGE NP=N+1 NA4=NA+N DO I=1,3 L=LOFIE(I) DO IY=NA+1,NP JX=NA4-IY DO IX=JX+2,NP JY=NA4-IX XD(I,IX,IY)=XD(I,JX,JY)*L ENDDO ENDDO ENDDO ! COPY TO X WITH SYMMETRIZATION ACROSS THE DIAGONAL, X=Y: DO I=1,3 J=JOFIF(I) DO IY=-1,NP DO IX=-1,IY X(I,IX,IY)=XD(J,IY,IX) X(J,IY,IX)=XD(J,IY,IX) ENDDO ENDDO ENDDO RETURN END !------------------------------------------------------------------------------C ! R.J.Purser, NCEP, Washington D.C. 1996 C ! SUBROUTINE INSUBP C ! Initialize octagon-quadrant cartesian locations to those given by the C ! conformal mapping solution C ! C ! <-- X double precision array of locations of ocatgon quadrant & margin C ! --> MA steps between octagon corner and center line C ! --> M1 steps between edge and center line C !------------------------------------------------------------------------------C SUBROUTINE INSUBP(X,MA,M1) ! INITIALIZATION OF GRID VALUES USING CONFORMAL SOLUTION: COMPLEX Z,W DIMENSION X(3,-1:M1+1,-1:M1+1) D=1./M1 DO IY=0,M1 YA=IY*D DO IX=0,MIN(M1,M1+MA-IY) XA=IX*D Z=XA+YA*(0,1) !mish ! Print *, "BEFORE ix=",ix ! Print *, "BEFORE z=",z ! Print *, "BEFORE xa=",xa ! Print *, "BEFORE ya=",ya !mish CALL TOCT(Z,W) !mish ! Print *, "AFTER ix =",ix !mish XW=REAL(W) YW=AIMAG(W) XX=XW**2+YW**2 S=1./(1.+XX) X(1,IX,IY)=2.*XW*S X(2,IX,IY)=2.*YW*S X(3,IX,IY)=(1.-XX)*S ENDDO ENDDO CALL SYMMET(X(1,-1,-1),MA,M1) !mish ! Print *, "AFTER" ! stop !mish RETURN END !------------------------------------------------------------------------------C ! R.J.Purser, National Meteorological Center, Washington D.C. 1995 C ! SUBROUTINE TOCT C ! Transform from complez-Z in the standard unit-octagon to complex-W in the C ! unit-circle C !------------------------------------------------------------------------------C SUBROUTINE TOCT(Z,W) COMPLEX CI,CIR,CIQ,Z,W,ZT LOGICAL KX,KY,KXY,KX1,KXY1,KS COMMON/CSTOCT1/A1(40),A2(40),B1(40),B2(40),A,B,CX,CY *,SSR,N1,ST,STI,CI,CIR,CIQ X=REAL(Z) Y=AIMAG(Z) if(z.eq.0.) then W=X+Y*(0,1) return end if if(x.eq.1.) x=0.9999 if(y.eq.1.) y=0.9999 KX=X.LT.0. IF(KX)X=-X KY=Y.LT.0. IF(KY)Y=-Y KXY=Y.GT.X IF(KXY)THEN T=X X=Y Y=T ENDIF KX1=X.GT.1. KXY1=X+Y.GT.1.+A IF(KX1)THEN X=2.-X ELSEIF(KXY1)THEN T=X X=1.+A-Y Y=1.+A-T ENDIF DD1=X**2+Y**2 DD2=(1.-X)**2+(A-Y)**2 ZT=X+Y*(0,1) KS=DD1.LT.SSR*DD2 !misha ! print *, "come to break1, ks=", ks ! print *, " , x= ",x ! print *, " , y= ",y !misha IF(KS)THEN ZT=ZT**4 CALL TAY(ZT,A1,N1,W) W=CIQ*(W/CI)**.25 ELSE ZT=A+CI*ZT-CI ZT=ZT**ST CALL TAY(ZT,A2,N1,W) W=W+B W=(W+CI)/(CI-W) ENDIF IF(KX1.OR.KXY1)W=W/CABS(W)**2 X=REAL(W) Y=AIMAG(W) IF(KXY)THEN T=X X=Y Y=T ENDIF IF(KY)Y=-Y IF(KX)X=-X W=X+Y*(0,1) !misha ! print *, " , x= ",x ! print *, " , y= ",y !misha RETURN END !------------------------------------------------------------------------------C ! R.J.Purser, National Meteorological Center, Washington D.C. 1994 C ! SUBROUTINE TAY C ! Evaluate the complex function W of Z whose real C ! Taylor series coefficients are RA. C ! C ! --> Z function argument (complex) C ! --> RA Taylor coefficients (real) C ! --> N number of coeffients (starting with the linear term) C ! <-- W Taylor-series approximation of the function (complex) C !------------------------------------------------------------------------------C SUBROUTINE TAY(Z,RA,N,W) COMPLEX Z,W DIMENSION RA(*) W=0. DO I=N,1,-1 W=(W+RA(I))*Z ENDDO RETURN END !------------------------------------------------------------------------------C ! R.J.Purser, National Meteorological Center, Washington D.C. 1994 C ! SUBROUTINE CINVRT C ! Compute the Taylor series coefficients Z for the functional-inverse of C ! the function whose Taylor series coefficients are W, for the case where C ! the constant coefficients of both Z and W are 0 C ! C ! --> W Taylor coefficients of original function (starting with linear) C ! <-- Z Taylor coefficients of inverse function C ! --> M number of Taylor series coefficients computed C ! --- WORK workspace array consisting of at least 3*M elements C !------------------------------------------------------------------------------C SUBROUTINE CINVRT(W,Z,M,WORK) DIMENSION W(*),Z(*),WORK(M,3) W1=W(1) DO I=1,M WORK(I,1)=W(I) WORK(I,2)=0. ENDDO WORK(1,2)=1. DO J=1,M ZJ=WORK(J,2)/WORK(J,1) Z(J)=ZJ DO I=J,M WORK(I,2)=WORK(I,2)-ZJ*WORK(I,1) ENDDO CALL CONV(WORK(1,1),W,WORK(1,3),M) DO I=1,M WORK(I,1)=WORK(I,3) ENDDO ENDDO RETURN END !------------------------------------------------------------------------------C ! R.J.Purser, National Meteorological Center, Washington D.C. 1994 C ! SUBROUTINE CONV C ! Convolve double-precision series A with B to form C, up to M terms C ! starting with element 1 C ! C ! A,B --> inputs (convolution factors) C ! C <-- output (convolution product) C ! M --> number of elements of A, B, C C !------------------------------------------------------------------------------C SUBROUTINE CONV(A,B,C,M) DIMENSION A(*),B(*),C(*) DO K=1,M C(K)=0. ENDDO DO I=1,M DO J=1,M-I K=I+J C(K)=C(K)+A(I)*B(J) ENDDO ENDDO RETURN END