! WD23JP ! ***************** ! * SOCT2.FOR * ! * PURSER 1996 * ! ***************** ! !------------------------------------------------------------------------------C ! R.J.Purser, NCEP, Washington D.C. 1996 C ! SUBROUTINE RELSUBP C ! Perform one sweep of the relaxations of the variational problem within the C ! first octant (0 to 45 degrees from x axis), or subvert the code to compute C ! residuals instead. Everything done in double precision. C ! C ! <-> X estimated cartesian locations for each octogon grid point C ! <-> FORX forcings associated with locations X C ! <-> DELX residuals associated with X C ! --> NA steps between octagon corner and center line C ! --> N steps between octagon center and edge C ! --> ORC overrelaxation coefficient (typically 1.2 or 1.3) C ! --> AODD homogeneity parameter normalized by step size C ! <-- ACC accumulated rms norm C ! --> MODE mode of application: 0 => standard relaxation iteration C ! 1 => get residuals, put into DELX C ! 2 => adjust forcings, put into FORX C !------------------------------------------------------------------------------C SUBROUTINE RELSUBP(X,FORX,DELX,NA,N,ORC,AODD,ACC,MODE) DIMENSION X(3,-1:N+1,-1:N+1),DELX(3,-1:N+1,-1:N+1) *,FORX(3,-1:N+1,-1:N+1),XCOR0(3,3),XCOR1(3,3) NM=N-1 NP=N+1 NA4=NA+N NA2=(NA4+1)/2 NAM=NA-1 NAP=NA+1 NACC=0 ACC=0. ! STANDARD RELAXATIONS APPLIED TO LOWER-BAND OF PRINCIPAL SUB-PANEL.. DO IY=0,NA-1 IYM=IY-1 IYP=IY+1 DO IX=IY,N IXM=IX-1 IXP=IX+1 CALL RELAX(FORX(1,IX,IY),DELX(1,IX,IY) * ,X(1,IX,IY),X(1,IXP,IY),X(1,IX,IYP),X(1,IXM,IY) * ,X(1,IX,IYM),X(1,IXP,IYP),X(1,IXM,IYP),X(1,IXM,IYM) * ,X(1,IXP,IYM),ORC,AODD,ACC,NACC,MODE) ENDDO ENDDO ! ..APPLIED TO MIDLINE.. IY=NA IYM=IY-1 IYP=IY+1 DO IX=IY,NM IXM=IX-1 IXP=IX+1 CALL RELAX(FORX(1,IX,IY),DELX(1,IX,IY) * ,X(1,IX,IY),X(1,IXP,IY),X(1,IX,IYP),X(1,IXM,IY) * ,X(1,IX,IYM),X(1,IXP,IYP),X(1,IXM,IYP),X(1,IXM,IYM) * ,X(1,IXP,IYM),ORC,AODD,ACC,NACC,MODE) ENDDO ! ..APPLIED TO UPPER BAND: DO IY=NA+1,NA2 IYM=IY-1 IYP=IY+1 DO IX=IY,NA4-IY IXM=IX-1 IXP=IX+1 CALL RELAX(FORX(1,IX,IY),DELX(1,IX,IY) * ,X(1,IX,IY),X(1,IXP,IY),X(1,IX,IYP),X(1,IXM,IY) * ,X(1,IX,IYM),X(1,IXP,IYP),X(1,IXM,IYP),X(1,IXM,IYM) * ,X(1,IXP,IYM),ORC,AODD,ACC,NACC,MODE) ENDDO ENDDO DO I=1,3 XCOR0(I,1)=X(I,N,NAM) XCOR1(I,1)=X(I,NM,NAM) XCOR0(I,3)=X(I,NM,NA) XCOR1(I,3)=X(I,NM,NAP) XCOR0(I,2)=X(I,NP,NA) XCOR1(I,2)=X(I,NP,NAM) ENDDO CALL RELV(FORX(1,N,NA),DELX(1,N,NA) *,X(1,N,NA),XCOR0,XCOR1,3,ORC,AODD,ACC,NACC,MODE) IF(MODE.EQ.0)THEN ACC=ACC/NACC ACC=SQRT(ACC) CALL SYMMET(X,NA,N) ELSEIF(MODE.EQ.1)THEN CALL SYMMET(DELX,NA,N) ELSE CALL SYMMET(FORX,NA,N) ENDIF RETURN END !------------------------------------------------------------------------------C ! R.J.Purser, NCEP, Washington D.C. 1996 C ! SUBROUTINE RELAX C ! Relaxation (MODE.eq.0), or residual determination (MODE.eq.1), or C ! adjustment of forcing (MODE.eq.2), at a generic grid point. C ! C ! <-> FORX,DELX,X: vectors of forcing, residual and solution at this point C ! C ! XF XB XE: C ! --> XC XA: surrounding solution vectors in the arrangement shown C ! XG XD XH: C ! C ! --> ORC: overrelaxation cooeficient C ! --> AODD: homogeneity parameter rescaled to grid units C ! --> ACC: accumulated squared residual C ! --> NACC: counter of contributors to ACC C ! --> MODE: mode of application C !------------------------------------------------------------------------------C SUBROUTINE RELAX(FORX,DELX,X *,XA,XB,XC,XD,XE,XF,XG,XH,ORC,AODD,ACC,NACC,MODE) PARAMETER (DELTA=.3E-4,DELTAI=1./DELTA) DIMENSION FORX(3),DELX(3),X(3),XA(3),XB(3),XC(3),XD(3) *,XE(3),XF(3),XG(3),XH(3),XU(3),XV(3),U(3),V(3),XS(3) SX=0. SY=0. DO I=1,3 U(I)=(XA(I)-XC(I))*.5D0 V(I)=(XB(I)-XD(I))*.5D0 SX=SX+X(I)*U(I) SY=SY+X(I)*V(I) XS(I)=X(I) ENDDO !misha ! print *,'================= XC =',XC ! print *,'================= XD =',XD !misha DO I=1,3 U(I)=U(I)-X(I)*SX V(I)=V(I)-X(I)*SY ENDDO SX=0. SY=0. DO I=1,3 XU(I)=X(I)+U(I)*DELTA XV(I)=X(I)+V(I)*DELTA SX=SX+XU(I)**2 SY=SY+XV(I)**2 ENDDO SX=1./SQRT(SX) SY=1./SQRT(SY) DO I=1,3 XU(I)=XU(I)*SX XV(I)=XV(I)*SY ENDDO CALL REL(X ,XA,XB,XC,XD,XE,XF,XG,XH,U,V,AODD,DX,DY) CALL REL(XU,XA,XB,XC,XD,XE,XF,XG,XH,U,V,AODD,DXX,DYX) CALL REL(XV,XA,XB,XC,XD,XE,XF,XG,XH,U,V,AODD,DXY,DYY) DXX=(DXX-DX)*DELTAI DXY=(DXY-DY+DYX-DX)*DELTAI*.5D0 DYY=(DYY-DY)*DELTAI DET=DXX*DYY-DXY*DXY ! IF(DET.EQ.0.D0)CALL WHAT(X,XA,XB,XC,XD,XE,XF,XG,XH) DETI=1./DET SX=(DYY*DX-DXY*DY)*DETI SY=(DXX*DY-DXY*DX)*DETI IF(MODE.EQ.0)THEN ! RELAXATION WITH LAGRANGE MULTIPLIER DELX S=0. DO I=1,3 X(I)=X(I)+(FORX(I)-U(I)*SX-V(I)*SY)*ORC S=S+X(I)**2 ENDDO S=1./SQRT(S) DO I=1,3 X(I)=X(I)*S ENDDO DO I=1,3 ACC=ACC+(XS(I)-X(I))**2 ENDDO ELSEIF(MODE.EQ.1)THEN ! DIAGNOSE RESIDUAL WITHOUT CORRECTING DO I=1,3 DELX(I)=(FORX(I)-U(I)*SX-V(I)*SY) ACC=ACC+DELX(I)**2 ENDDO ELSEIF(MODE.EQ.2)THEN ! RESET LAGRANGE MULTIPLIER DO I=1,3 ACC=ACC+FORX(I)**2 FORX(I)=FORX(I)+U(I)*SX+V(I)*SY ENDDO ELSE STOP'INVALID MODE IN RELAX' ENDIF NACC=NACC+1 RETURN END !------------------------------------------------------------------------------C ! R.J.Purser, NCEP, Washington D.C. 1996 C ! SUBROUTINE RELV C ! Relaxation (MODE.eq.0), or residual determination (MODE.eq.1), or C ! adjustment of forcing (MODE.eq.2), at a vertex of the octagon C ! C ! <-> FORX,DELX,X: vectors of forcing, residual and solution at this point C ! --> XCOR0: X at the NC nearest surrounding grid points C ! --> XCOR1: X at the NC next nearest surrounding points C ! --> NC: number of panels coming together at this corner (=3 for octagon) C ! --> ORC: over relaxation coefficient C ! --> AODD: homogeneity parameter rescaled to grid units C ! --> ACC: accumulated squared residual C ! --> NACC: counter of contributors to ACC C ! --> MODE: mode of application C !------------------------------------------------------------------------------C SUBROUTINE RELV(FORX,DELX,X *,XCOR0,XCOR1,NC,ORC,AODD,ACC,NACC,MODE) PARAMETER (DELTA=.3E-4,DELTAI=1./DELTA) DIMENSION FORX(3),DELX(3),X(3),XCOR0(3,3),XCOR1(3,3) *,XU(3),XV(3),U(3),V(3),XS(3) SX=0. SY=0. DO I=1,3 U(I)=XCOR0(I,1)-X(I) V(I)=XCOR1(I,1)-XCOR0(I,1) SX=SX+X(I)*U(I) SY=SY+X(I)*V(I) XS(I)=X(I) ENDDO DO I=1,3 U(I)=U(I)-X(I)*SX V(I)=V(I)-X(I)*SY ENDDO SX=0. SY=0. DO I=1,3 XU(I)=X(I)+U(I)*DELTA XV(I)=X(I)+V(I)*DELTA SX=SX+XU(I)**2 SY=SY+XV(I)**2 ENDDO SX=1./SQRT(SX) SY=1./SQRT(SY) DO I=1,3 XU(I)=XU(I)*SX XV(I)=XV(I)*SY ENDDO CALL DELV(X,XCOR0,XCOR1,NC *,U,V,AODD,DX,DY) CALL DELV(XU,XCOR0,XCOR1,NC *,U,V,AODD,DXX,DYX) CALL DELV(XV,XCOR0,XCOR1,NC *,U,V,AODD,DXY,DYY) DXX=(DXX-DX)*DELTAI DXY=(DXY-DY+DYX-DX)*DELTAI*.5D0 DYY=(DYY-DY)*DELTAI DETI=1./(DXX*DYY-DXY*DXY) SX=(DYY*DX-DXY*DY)*DETI SY=(DXX*DY-DXY*DX)*DETI IF(MODE.EQ.0)THEN ! RELAXATION WITH LAGRANGE MULTIPLIER DELX S=0. DO I=1,3 X(I)=X(I)+(FORX(I)-U(I)*SX-V(I)*SY)*ORC S=S+X(I)**2 ENDDO S=1./SQRT(S) DO I=1,3 X(I)=X(I)*S ENDDO DO I=1,3 ACC=ACC+(XS(I)-X(I))**2 ENDDO ELSEIF(MODE.EQ.1)THEN ! DIAGNOSE RESIDUAL WITHOUT CORRECTING DO I=1,3 DELX(I)=(FORX(I)-U(I)*SX-V(I)*SY) ACC=ACC+DELX(I)**2 ENDDO ELSEIF(MODE.EQ.2)THEN ! RESET LAGRANGE MULTIPLIER DO I=1,3 ACC=ACC+FORX(I)**2 FORX(I)=FORX(I)+U(I)*SX+V(I)*SY ENDDO ELSE STOP'INVALID MODE IN RELAX' ENDIF NACC=NACC+1 RETURN END SUBROUTINE DELV(X,XCOR0,XCOR1,NC,U,V,AODD,DX,DY) DIMENSION X(3),XCOR0(3,*),XCOR1(3,*),U(3),V(3) DX=0. DY=0. DO IC=1,NC ICP=IC+1 IF(ICP.GT.NC)ICP=1 UE=0. UAB=0. VAB=0. VE=0. SAB=0. SE=0. DO I=1,3 UI=U(I) VI=V(I) XI=X(I) D1=XCOR1(I,IC)-XI D2=XCOR0(I,ICP)-XCOR1(I,IC) D3=XCOR1(I,IC)-XI DX=DX+UI*D3 DY=DY+VI*D3 UE=UE+UI*D1 UAB=UAB+UI*D2 VE=VE+VI*D1 VAB=VAB+VI*D2 SAB=SAB+D2*D2 SE=SE+D1*D2 ENDDO DX=DX+AODD*(UE*SAB-UAB*SE) DY=DY+AODD*(VE*SAB-VAB*SE) ENDDO RETURN END SUBROUTINE REL(X,XA,XB,XC,XD,XE,XF,XG,XH,U,V,AODD,DX,DY) ! ! XF XB XE ! QF QE ! XC X XA ! QG QH ! XG XD XH ! DIMENSION X(3),XA(3),XB(3),XC(3),XD(3),XE(3),XF(3),XG(3),XH(3) *,U(3),V(3) DX=0. DY=0. UAB=0. UBC=0. UCD=0. UDA=0. UE=0. UF=0. UG=0. UH=0. VAB=0. VBC=0. VCD=0. VDA=0. VE=0. VF=0. VG=0. VH=0. SAB=0. SBC=0. SCD=0. SDA=0. SE=0. SF=0. SG=0. SH=0. DO I=1,3 UI=U(I) VI=V(I) XI=X(I) D1=XE(I)-XI D2=XB(I)-XA(I) D3=XA(I)-XI UE=UE+UI*D1 UAB=UAB+UI*D2 DX=DX+UI*D3 VE=VE+VI*D1 VAB=VAB+VI*D2 DY=DY+VI*D3 SAB=SAB+D2*D2 SE=SE+D1*D2 D1=XF(I)-XI D2=XC(I)-XB(I) D3=XB(I)-XI UF=UF+UI*D1 UBC=UBC+UI*D2 DX=DX+UI*D3 VF=VF+VI*D1 VBC=VBC+VI*D2 DY=DY+VI*D3 SBC=SBC+D2*D2 SF=SF+D1*D2 D1=XG(I)-XI D2=XD(I)-XC(I) D3=XC(I)-XI UG=UG+UI*D1 UCD=UCD+UI*D2 DX=DX+UI*D3 VG=VG+VI*D1 VCD=VCD+VI*D2 DY=DY+VI*D3 SCD=SCD+D2*D2 SG=SG+D1*D2 D1=XH(I)-XI D2=XA(I)-XD(I) D3=XD(I)-XI UH=UH+UI*D1 UDA=UDA+UI*D2 DX=DX+UI*D3 VH=VH+VI*D1 VDA=VDA+VI*D2 DY=DY+VI*D3 SDA=SDA+D2*D2 SH=SH+D1*D2 ENDDO DX=DX+AODD*(UE*SAB+UF*SBC+UG*SCD+UH*SDA * -UAB*SE-UBC*SF-UCD*SG-UDA*SH) DY=DY+AODD*(VE*SAB+VF*SBC+VG*SCD+VH*SDA * -VAB*SE-VBC*SF-VCD*SG-VDA*SH) RETURN END SUBROUTINE SYMMET(X,NA,N) DIMENSION X(3,-1:N+1,-1:N+1) DIMENSION LOFIA(3),LOFIC(3),LOFID(3),LOFIE(3) *,JOFIF(3) DATA LOFIA/1,1,-1/,LOFIC/-1,1,1/,LOFID/1,-1,1/,LOFIE/1,1,-1/ DATA JOFIF/2,1,3/ NM=N-1 NP=N+1 NA4=NA+N NA2=(NA4+1)/2 ! SYMMETRIZE ON BOUNDARY OF SUB-PANEL: X(1,0,0)=0. X(2,0,0)=0. DO IX=1,N X(2,IX,0)=0. ENDDO DO IY=0,NA X(3,N,IY)=0. ENDDO DO IY=NA+1,NA2 IX=NA4-IY X(3,IX,IY)=0. ENDDO DO IY=1,NA2 IX=IY S=(X(1,IX,IY)+X(2,IX,IY))/2.D0 X(1,IX,IY)=S X(2,IX,IY)=S ENDDO ! SYMMETRIZE JUST OUTSIDE BOUNDARY OF SUB-PANEL: DO I=1,3 L=LOFIC(I) DO IY=0,1 X(I,-1,IY)=L*X(I,1,IY) ENDDO L=LOFIA(I) DO IY=0,NA X(I,NP,IY)=L*X(I,NM,IY) ENDDO J=JOFIF(I) DO IY=1,NA2 IX=IY-1 X(I,IX,IY)=X(J,IX+1,IY-1) ENDDO DO IY=2,NA2+1 IX=IY-2 X(I,IX,IY)=X(J,IX+2,IY-2) ENDDO L=LOFID(I) DO IX=-1,NP X(I,IX,-1)=L*X(I,IX,1) ENDDO ENDDO DO I=1,3 L=LOFIE(I) DO IY=NA+1,NA2+1 IX=NA4-IY+1 X(I,IX,IY)=L*X(I,IX-1,IY-1) ENDDO DO IY=NA+1,NA2+1 IX=NA4-IY+2 X(I,IX,IY)=L*X(I,IX-2,IY-2) ENDDO ENDDO RETURN END SUBROUTINE ADDXX2(X,X2,NA,N) DIMENSION X(3,-1:N+1,-1:N+1),X2(3,-1:N*2+1,-1:N*2+1) NA4=NA+N NA2=(NA4+1)/2 N2=N*2 N2A=NA*2 IY=0 IY2=IY*2 IY2M=IY2-1 IY2P=IY2+1 DO IX=0,N IX2=IX*2 IX2M=IX2-1 IX2P=IX2+1 DO I=1,3 R=X(I,IX,IY) X2(I,IX2,IY2)=X2(I,IX2,IY2)+R R=.5D0*R X2(I,IX2P,IY2)=X2(I,IX2P,IY2)+R X2(I,IX2,IY2P)=X2(I,IX2,IY2P)+R X2(I,IX2M,IY2)=X2(I,IX2M,IY2)+R X2(I,IX2,IY2M)=X2(I,IX2,IY2M)+R R=R*.5D0 X2(I,IX2P,IY2P)=X2(I,IX2P,IY2P)+R X2(I,IX2M,IY2P)=X2(I,IX2M,IY2P)+R X2(I,IX2M,IY2M)=X2(I,IX2M,IY2M)+R X2(I,IX2P,IY2M)=X2(I,IX2P,IY2M)+R ENDDO ENDDO DO IY=1,NA IY2=IY*2 IY2M=IY2-1 IY2P=IY2+1 DO IX=IY-1,N IX2=IX*2 IX2M=IX2-1 IX2P=IX2+1 DO I=1,3 R=X(I,IX,IY) X2(I,IX2,IY2)=X2(I,IX2,IY2)+R R=.5D0*R X2(I,IX2P,IY2)=X2(I,IX2P,IY2)+R X2(I,IX2,IY2P)=X2(I,IX2,IY2P)+R X2(I,IX2M,IY2)=X2(I,IX2M,IY2)+R X2(I,IX2,IY2M)=X2(I,IX2,IY2M)+R R=R*.5D0 X2(I,IX2P,IY2P)=X2(I,IX2P,IY2P)+R X2(I,IX2M,IY2P)=X2(I,IX2M,IY2P)+R X2(I,IX2M,IY2M)=X2(I,IX2M,IY2M)+R X2(I,IX2P,IY2M)=X2(I,IX2P,IY2M)+R ENDDO ENDDO ENDDO DO IY=NA+1,NA2+1 IY2=IY*2 IY2M=IY2-1 IY2P=IY2+1 DO IX=IY-1,NA4-IY+1 IX2=IX*2 IX2M=IX2-1 IX2P=IX2+1 DO I=1,3 R=X(I,IX,IY) X2(I,IX2,IY2)=X2(I,IX2,IY2)+R R=.5D0*R X2(I,IX2P,IY2)=X2(I,IX2P,IY2)+R X2(I,IX2,IY2P)=X2(I,IX2,IY2P)+R X2(I,IX2M,IY2)=X2(I,IX2M,IY2)+R X2(I,IX2,IY2M)=X2(I,IX2,IY2M)+R R=R*.5D0 X2(I,IX2P,IY2P)=X2(I,IX2P,IY2P)+R X2(I,IX2M,IY2P)=X2(I,IX2M,IY2P)+R X2(I,IX2M,IY2M)=X2(I,IX2M,IY2M)+R X2(I,IX2P,IY2M)=X2(I,IX2P,IY2M)+R ENDDO ENDDO ENDDO CALL FNORM(X2,N2A,N2) CALL SYMMET(X2,N2A,N2) RETURN END SUBROUTINE ADDX2X(X2,X,NA,N) DIMENSION X(3,-1:N+1,-1:N+1),X2(3,-1:N*2+1,-1:N*2+1) NA4=NA+N NA2=(NA4+1)/2 N2=N*2 N2A=NA*2 N2M=N2-1 N2P=N2+1 N2AM=N2A-1 N2AP=N2A+1 DO IY=0,NA-1 IY2=IY*2 IY2M=IY2-1 IY2P=IY2+1 DO IX=IY,N IX2=IX*2 IX2M=IX2-1 IX2P=IX2+1 DO I=1,3 X(I,IX,IY)=X(I,IX,IY)+X2(I,IX2,IY2)+.5D0*( * X2(I,IX2P,IY2)+X2(I,IX2,IY2P)+X2(I,IX2M,IY2)+X2(I,IX2,IY2M)+.5*( * X2(I,IX2P,IY2P)+X2(I,IX2M,IY2P)+X2(I,IX2M,IY2M)+X2(I,IX2P,IY2M) * )) ENDDO ENDDO ENDDO IY=NA IY2=IY*2 IY2M=IY2-1 IY2P=IY2+1 DO IX=IY,N-1 IX2=IX*2 IX2M=IX2-1 IX2P=IX2+1 DO I=1,3 X(I,IX,IY)=X(I,IX,IY)+X2(I,IX2,IY2) +.5D0*( * X2(I,IX2P,IY2)+X2(I,IX2,IY2P)+X2(I,IX2M,IY2)+X2(I,IX2,IY2M)+.5*( *X2(I,IX2P,IY2P)+X2(I,IX2M,IY2P)+X2(I,IX2M,IY2M)+X2(I,IX2P,IY2M))) ENDDO ENDDO DO IY=NA+1,NA2 IY2=IY*2 IY2M=IY2-1 IY2P=IY2+1 DO IX=IY,NA4-IY IX2=IX*2 IX2M=IX2-1 IX2P=IX2+1 DO I=1,3 X(I,IX,IY)=X(I,IX,IY)+X2(I,IX2,IY2) +.5D0*( * X2(I,IX2P,IY2)+X2(I,IX2,IY2P)+X2(I,IX2M,IY2)+X2(I,IX2,IY2M)+.5*( *X2(I,IX2P,IY2P)+X2(I,IX2M,IY2P)+X2(I,IX2M,IY2M)+X2(I,IX2P,IY2M))) ENDDO ENDDO ENDDO DO I=1,3 X(I,N,NA)=X(I,N,NA)+X2(I,N2,N2A) * +.5D0*(X2(I,N2,N2AM)+X2(I,N2P,N2A)+X2(I,N2M,N2A) * +.5D0*(X2(I,N2M,N2AM)+X2(I,N2P,N2AM)+X2(I,N2M,N2AP))) ENDDO RETURN END SUBROUTINE SETX2X(X2,X,NA,N) DIMENSION X(3,-1:N+1,-1:N+1),X2(3,-1:N*2+1,-1:N*2+1) NA4=NA+N NA2=(NA4+1)/2 DO IY=0,NA IY2=IY*2 DO IX=IY,N IX2=IX*2 DO I=1,3 X(I,IX,IY)=X2(I,IX2,IY2) ENDDO ENDDO ENDDO DO IY=NA+1,NA2 IY2=IY*2 DO IX=IY,NA4-IY IX2=IX*2 DO I=1,3 X(I,IX,IY)=X2(I,IX2,IY2) ENDDO ENDDO ENDDO CALL SYMMET(X,NA,N) RETURN END SUBROUTINE DIFXX2(X,X2,NA,N) DIMENSION X(3,-1:N+1,-1:N+1),X2(3,-1:N*2+1,-1:N*2+1) NA4=NA+N NA2=(NA4+1)/2 DO IY=0,NA IYM=IY-1 IYP=IY+1 IY2=IY*2 DO IX=0,N IXM=IX-1 IXP=IX+1 IX2=IX*2 DO I=1,3 X(I,IX,IY)=X(I,IX,IY)-X2(I,IX2,IY2) ENDDO ENDDO ENDDO DO IY=NA+1,NA2 IYM=IY-1 IYP=IY+1 IY2=IY*2 DO IX=IY,NA4-IY IXM=IX-1 IXP=IX+1 IX2=IX*2 DO I=1,3 X(I,IX,IY)=X(I,IX,IY)-X2(I,IX2,IY2) ENDDO ENDDO ENDDO CALL SYMMET(X,NA,N) RETURN END SUBROUTINE ZERV(A,N) DIMENSION A(*) DO I=1,N A(I)=0. ENDDO RETURN END SUBROUTINE FNORM(X,NA,N) DIMENSION X(3,-1:N+1,-1:N+1) NA4=NA+N NA2=(NA4+1)/2 DO IY=0,NA DO IX=IY,N CALL NORM3(X(1,IX,IY)) ENDDO ENDDO DO IY=NA+1,NA2 DO IX=IY,NA4-IY CALL NORM3(X(1,IX,IY)) ENDDO ENDDO RETURN END SUBROUTINE NORM3(X) ! NORMALIZE THE 3-VECTOR: DIMENSION X(3) S=0. DO I=1,3 S=S+X(I)*X(I) ENDDO S=1./SQRT(S) DO I=1,3 X(I)=X(I)*S ENDDO RETURN END