C------------------------------------------------------------------------------C C R.J.Purser, National Centers for Environmental Prediction, Washington D.C. C C wd23jp@sun1.wwb.noaa.gov 1996 C C SUBROUTINE XMTOXC C C Transform from map coordinates to cartesian coordinates on unit-sphere C C C C --> XM: coordinates (2 components) in [-1.,1.] in specified map panel C C <-- XC: cartesian coordinates (3 components) of corresponding point C C <-- DXCDXM: Jacobian of the transformation (a (3*2)-matrix) C C --> IPAN: index of the specified map panel (between 1 and 6) C C------------------------------------------------------------------------------C SUBROUTINE XMTOXC(XM,XC,DXCDXM,IPAN) PARAMETER(N=64,NP=N+1) ccc PARAMETER(N=128,NP=N+1) !!!!!!DRAGAN COMMON/QPAN/X(3,-1:NP,-1:NP) DIMENSION XM(2),XC(3),DXCDXM(3,2) *,RPAN(3,3,6) DATA RPAN/0.,1.,0., 0.,0.,1., 1., 0.,0. * ,-1.,0.,0., 0.,0.,1., 0., 1.,0. * ,0.,-1.,0., 0.,0.,1., -1., 0.,0. * ,1., 0.,0., 0.,0.,1., 0.,-1.,0. * ,0., 1.,0., -1.,0.,0., 0., 0.,1. * ,0., 1.,0., 1.,0.,0., 0.,0.,-1./ NV=N+3 NM=N-1 AX=ABS(XM(1))*N AY=ABS(XM(2))*N IX2=AX IF(IX2.GE.N)IX2=NM IX1=IX2-1 RX2=AX-IX2 RX3=RX2-1. IY2=AY IF(IY2.GE.N)IY2=NM IY1=IY2-1 RY2=AY-IY2 RY3=RY2-1. W22=RX3*RY3 W32=-RX2*RY3 W23=-RX3*RY2 W33=RX2*RY2 DO I=1,3 XC(I)=0. DXCDXM(I,1)=0. DXCDXM(I,2)=0. ENDDO C ACCUMULATE FOUR CONTRIBUTIONS FROM THE 3*3-TEMPLATE INTERPOLATIONS C CENTERED AT THIS ELEMENT'S CORNERS... CALL QIN3(X(1,IX1,IY1),RX2,RY2,W22,RY3,RX3 *,XC,DXCDXM(1,1),DXCDXM(1,2),NV) CALL QIN3(X(1,IX2,IY1),RX3,RY2,W32,-RY3,-RX2 *,XC,DXCDXM(1,1),DXCDXM(1,2),NV) CALL QIN3(X(1,IX1,IY2),RX2,RY3,W23,-RY2,-RX3 *,XC,DXCDXM(1,1),DXCDXM(1,2),NV) IF(IX2.NE.NM.OR.IY2.NE.NM)THEN CALL QIN3(X(1,IX2,IY2),RX3,RY3,W33,RY2,RX2 * ,XC,DXCDXM(1,1),DXCDXM(1,2),NV) ELSE C ... BUT INVOKE SPECIAL "FRACTIONAL POWER" CODE WHEN ONE OF THESE ELEMENT- C CORNER POINTS IS THE CORNER OF THE MAP PANEL ITSELF: CALL FIN3(X(1,IX2,IY2),RX3,RY3,W33,RY2,RX2 * ,XC,DXCDXM(1,1),DXCDXM(1,2),NV) ENDIF DO I=1,3 DXCDXM(I,1)=DXCDXM(I,1)*N DXCDXM(I,2)=DXCDXM(I,2)*N ENDDO C NORMALIZATION AND ORTHOGONALIZATION T1=0. T2=0. S=0. DO I=1,3 T1=T1+XC(I)*DXCDXM(I,1) T2=T2+XC(I)*DXCDXM(I,2) S=S+XC(I)**2 ENDDO S=1./SQRT(S) T1=T1*S T2=T2*S DO I=1,3 XC(I)=XC(I)*S DXCDXM(I,1)=DXCDXM(I,1)-T1*XC(I) DXCDXM(I,2)=DXCDXM(I,2)-T2*XC(I) ENDDO IF(XM(1).LT.0.)THEN XC(1)=-XC(1) DXCDXM(2,1)=-DXCDXM(2,1) DXCDXM(3,1)=-DXCDXM(3,1) DXCDXM(1,2)=-DXCDXM(1,2) ENDIF IF(XM(2).LT.0.)THEN XC(2)=-XC(2) DXCDXM(2,1)=-DXCDXM(2,1) DXCDXM(1,2)=-DXCDXM(1,2) DXCDXM(3,2)=-DXCDXM(3,2) ENDIF CALL ROTATE(RPAN(1,1,IPAN),XC) CALL ROTATE(RPAN(1,1,IPAN),DXCDXM(1,1)) CALL ROTATE(RPAN(1,1,IPAN),DXCDXM(1,2)) RETURN END C------------------------------------------------------------------------------C C R.J.Purser, National Centers for Environmental Prediction, Washington D.C. C C wd23jp@sun1.wwb.noaa.gov 1996 C C SUBROUTINE XCTOXM C C Transform from map coordinates to cartesian coordinates on unit-sphere C C C C --> XC: cartesian coordinates (3) of point on unit sphere C C <-- XM: map coordinates (2) of this point C C <-- DXMDXC: generalized Jacobian of the transformation (a (2*3)-matrix) C C <-- IPAN: index of the map panel (between 1 and 6) containing the point C C------------------------------------------------------------------------------C SUBROUTINE XCTOXM(XC,XM,DXMDXC,IPAN) DIMENSION XC(3),XM(2),DXMDXC(2,3),DXCDXM(3,2),G2(2,2),XCT(3) *,RPAN(3,3,6) DATA RPAN/0.,1.,0., 0.,0.,1., 1., 0.,0. * ,-1.,0.,0., 0.,0.,1., 0., 1.,0. * ,0.,-1.,0., 0.,0.,1., -1., 0.,0. * ,1., 0.,0., 0.,0.,1., 0.,-1.,0. * ,0., 1.,0., -1.,0.,0., 0., 0.,1. * ,0., 1.,0., 1.,0.,0., 0.,0.,-1./ DATA NIT/20/,SMAX/1.E-12/ AXC=ABS(XC(1)) AYC=ABS(XC(2)) AZC=ABS(XC(3)) IF(AXC.GT.AYC)THEN IF(AXC.GT.AZC)THEN IF(XC(1).GT.0.)GOTO 1 ! |X| GOTO 3 ! biggest ENDIF ELSEIF(AYC.GT.AZC)THEN IF(XC(2).GT.0.)GOTO 2 ! |Y| GOTO 4 ! biggest ENDIF IF(XC(3).GT.0.)GOTO 5 ! |Z| GOTO 6 ! biggest 1 IPAN=1 GOTO 7 2 IPAN=2 GOTO 7 3 IPAN=3 GOTO 7 4 IPAN=4 GOTO 7 5 IPAN=5 GOTO 7 6 IPAN=6 7 CALL MULTM(RPAN(1,1,IPAN),XC,XCT,3,3,1,3,3,3) D=1./XCT(3) XM(1)=XCT(1)*D XM(2)=XCT(2)*D CALL XMTOXC(XM,XCT,DXCDXM,IPAN) DO IT=1,NIT CALL MULTM(DXCDXM,DXCDXM,G2,2,3,2,3,3,2) DO J=1,3 DO I=1,2 DXMDXC(I,J)=DXCDXM(J,I) ENDDO ENDDO DET=G2(1,1)*G2(2,2)-G2(1,2)*G2(2,1) IF(ABS(DET).LE.1.E-14)THEN C DON'T TRY TO GET JACOBIAN AT THE CORNER SINGULARITIES WHERE DET VANISHES!!! CALL ZERM(DXMDXC,2,3,2) RETURN ENDIF CALL LINVMM(G2,DXMDXC,2,3,2,2) S=0. DO I=1,3 XCT(I)=XCT(I)-XC(I) S=S+XCT(I)**2 ENDDO CALL MSBMM(DXMDXC,XCT,XM,2,3,1,2,3,2) IF(S.LE.SMAX)RETURN CALL XMTOXC(XM,XCT,DXCDXM,IPAN) ENDDO STOP' ITERATIONS IN XCTOXM INSUFFICIENT FOR CONVERGENCE' END SUBROUTINE RELAXA(XO,XB,XC,XD,XF,XG,ORC,AODD) IMPLICIT REAL*8(A-H,O-Z) DIMENSION XO(3),XA(3),XB(3),XC(3),XD(3),XE(3),XF(3),XG(3),XH(3) *,JOFI(3) DATA JOFI/3,2,1/ DO I=1,3 J=JOFI(I) XA(I)=XC(J) XE(I)=XF(J) XH(I)=XG(J) ENDDO CALL RELAX(XO,XA,XB,XC,XD,XE,XF,XG,XH,ORC,AODD) RETURN END SUBROUTINE RELAXB(XO,XA,XC,XD,XG,XH,ORC,AODD) IMPLICIT REAL*8(A-H,O-Z) DIMENSION XO(3),XA(3),XB(3),XC(3),XD(3),XE(3),XF(3),XG(3),XH(3) *,JOFI(3) DATA JOFI/1,3,2/ DO I=1,3 J=JOFI(I) XB(I)=XD(J) XE(I)=XH(J) XF(I)=XG(J) ENDDO CALL RELAX(XO,XA,XB,XC,XD,XE,XF,XG,XH,ORC,AODD) RETURN END SUBROUTINE RELAXC(XO,XA,XB,XD,XE,XH,ORC,AODD) IMPLICIT REAL*8(A-H,O-Z) DIMENSION XO(3),XA(3),XB(3),XC(3),XD(3),XE(3),XF(3),XG(3),XH(3) *,LOFI(3) DATA LOFI/-1,1,1/ DO I=1,3 L=LOFI(I) XC(I)=XA(I)*L XF(I)=XE(I)*L XG(I)=XH(I)*L ENDDO CALL RELAX(XO,XA,XB,XC,XD,XE,XF,XG,XH,ORC,AODD) RETURN END SUBROUTINE RELAXD(XO,XA,XB,XC,XE,XF,ORC,AODD) IMPLICIT REAL*8(A-H,O-Z) DIMENSION XO(3),XA(3),XB(3),XC(3),XD(3),XE(3),XF(3),XG(3),XH(3) *,LOFI(3) DATA LOFI/1,-1,1/ DO I=1,3 L=LOFI(I) XD(I)=XB(I)*L XG(I)=XF(I)*L XH(I)=XE(I)*L ENDDO CALL RELAX(XO,XA,XB,XC,XD,XE,XF,XG,XH,ORC,AODD) RETURN END SUBROUTINE RELAX(X,XA,XB,XC,XD,XE,XF,XG,XH,ORC,AODD) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (DELTA=.01D0,DELTAI=1.D0/DELTA) DIMENSION 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) SX=0. SY=0. DO I=1,3 U(I)=(XA(I)-XC(I))*.5 V(I)=(XB(I)-XD(I))*.5 SX=SX+X(I)*U(I) SY=SY+X(I)*V(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 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*.5 DYY=(DYY-DY)*DELTAI DETI=ORC/(DXX*DYY-DXY*DXY) SX=(DYY*DX-DXY*DY)*DETI SY=(DXX*DY-DXY*DX)*DETI S=0. DO I=1,3 X(I)=X(I)-U(I)*SX-V(I)*SY S=S+X(I)**2 ENDDO S=1./SQRT(S) DO I=1,3 X(I)=X(I)*S ENDDO RETURN END SUBROUTINE REL(X,XA,XB,XC,XD,XE,XF,XG,XH,U,V,AODD,DX,DY) C C XF XB XE C QF QE C XC X XA C QG QH C XG XD XH C IMPLICIT REAL*8(A-H,O-Z) 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 QIN3(V,RX2,RY2,W,DWX,DWY,VC,DVCX,DVCY,NV) C APPLY QUADRATIC INTERPOLATION AT NON-CORNER ELEMENT OF TABLE DIMENSION V(3,NV,*),VC(3),DVCX(3),DVCY(3),PX(3),PY(3) *,DPX(3),DPY(3),YV(3),YDVX(3),YDVY(3),XV(3,3),XDVX(3,3) RX1=RX2+1. RX3=RX2-1. PX(1)= (RX2*RX3)*.5 PX(2)=-(RX1*RX3) PX(3)= (RX1*RX2)*.5 DPX(1)= (RX2+RX3)*.5 DPX(2)=-(RX1+RX3) DPX(3)= (RX1+RX2)*.5 RY1=RY2+1. RY3=RY2-1. PY(1)= (RY2*RY3)*.5 PY(2)=-(RY1*RY3) PY(3)= (RY1*RY2)*.5 DPY(1)= (RY2+RY3)*.5 DPY(2)=-(RY1+RY3) DPY(3)= (RY1+RY2)*.5 DO I=1,3 YV(I)=0. YDVX(I)=0. YDVY(I)=0. ENDDO DO IY=1,3 DO I=1,3 XV(I,IY)=0. XDVX(I,IY)=0. ENDDO DO IX=1,3 DO I=1,3 XV(I,IY)=XV(I,IY)+PX(IX)*V(I,IX,IY) XDVX(I,IY)=XDVX(I,IY)+DPX(IX)*V(I,IX,IY) ENDDO ENDDO DO I=1,3 YV(I)=YV(I)+PY(IY)*XV(I,IY) YDVY(I)=YDVY(I)+DPY(IY)*XV(I,IY) YDVX(I)=YDVX(I)+PY(IY)*XDVX(I,IY) ENDDO ENDDO DO I=1,3 VC(I)=VC(I)+W*YV(I) DVCX(I)=DVCX(I)+W*YDVX(I)+DWX*YV(I) DVCY(I)=DVCY(I)+W*YDVY(I)+DWY*YV(I) ENDDO RETURN END SUBROUTINE INFIN3 C INITIALIZE PARAMETERS NEEDED BY FIN3 IMPLICIT COMPLEX(C) COMMON/CORNEL/C3O8,R4O3,DCO(7,7) DIMENSION IPAR(7) DATA IPAR/1,1,-1,1,-1,1,1/ R1O3=1./3 C3O8=CMPLX(0.,1.) C3O8=-1./CSQRT(C3O8) R4O3=4./3 CALL ZERM(DCO,7,7,7) RHO=2.**R1O3 RHOS=RHO*RHO S=SQRT(.75) SIG=R1O3/(RHO+2.) RHOP=RHO+1. RHOM=RHO-1. DCO(1,1)=1. DCO(2,2)=4.*SIG DCO(2,3)=SIG*RHOS/2 DCO(2,4)=-2.*SIG DCO(2,5)=-SIG*RHOS DCO(3,3)=SIG*S*RHOS DCO(3,4)=4.*SIG*S DCO(4,1)=-15.*SIG DCO(4,2)=4.*SIG*RHOP DCO(4,3)=-SIG*(RHOS-1.) DCO(4,4)=2.*SIG*(2.-RHO) DCO(4,5)=SIG*(2.*RHOS+1.) DCO(5,3)=2.*SIG*RHOS*S DCO(5,4)=-4.*SIG*RHO*S DCO(6,1)=DCO(4,1) DCO(6,2)=-4.*SIG*RHOM DCO(6,3)=RHOS/6 DCO(6,4)=2./3 DCO(6,5)=-SIG*(2.*RHOS-1.) DCO(7,1)=-RHO*(2.-RHO)/4. DCO(7,2)=SIG*RHO DCO(7,3)=-SIG/2 DCO(7,4)=DCO(7,2) DCO(7,5)=DCO(7,3) DO I=2,7 DCO(I,6)=DCO(I,4)*IPAR(I) DCO(I,7)=DCO(I,3)*IPAR(I) ENDDO RETURN END SUBROUTINE FIN3(V,XZ,YZ,W,DWX,DWY,VC,DVCX,DVCY,NV) C APPLY FRACTIONAL-POWER TECHNIQUE IN CORNER ELEMENT INTERPOLATION IMPLICIT COMPLEX(C) COMMON/CORNEL/C3O8,R4O3,DCO(7,7) DIMENSION V(3,NV,*),VC(3),DVCX(3),DVCY(3) DIMENSION P(7),PX(7),PY(7) *,VCT(3),VCTX(3),VCTY(3),VN(7,3),DV(7,3) DATA P/7*0./,PX/7*0./PY/7*0./ IF(XZ.EQ.0..AND.YZ.EQ.0.)THEN DO I=1,3 VC(I)=VC(I)+W*V(I,2,2) DVCX(I)=DVCX(I)+DWX*V(I,2,2) DVCY(I)=DVCY(I)+DWY*V(I,2,2) ENDDO RETURN ENDIF ! PRINT'('' ENTERING NONTRIVIAL FIN3'')' WR4O3=W*R4O3 CZ=CMPLX(XZ,YZ) C=-(CZ*C3O8)**R4O3 CR=C/CZ WR=REAL(CR)*WR4O3 WQ=AIMAG(CR)*WR4O3 X=REAL(C) Y=AIMAG(C) DO I=1,3 VN(1,I)=V(I,2,2) VN(2,I)=V(I,3,2) VN(3,I)=V(I,1,3) VN(4,I)=V(I,1,2) VN(5,I)=V(I,1,1) VN(6,I)=V(I,2,1) VN(7,I)=V(I,3,1) CALL MULMM(DCO,VN(1,I),DV(1,I),7,7,1,7,7,7) ENDDO P(1)=1. P(2)=X P(3)=Y P(4)=.5*X**2 P(5)=X*Y P(6)=.5*Y**2 P(7)=X**3-3.*X*Y**2 PX(2)=1. PX(4)=X PX(5)=Y PX(7)=3.*(X**2-Y**2) PY(3)=1. PY(5)=X PY(6)=Y PY(7)=-6*X*Y DO I=1,3 VCT(I)=0. VCTX(I)=0. VCTY(I)=0. DO J=1,7 VCT(I)=VCT(I)+P(J)*DV(J,I) VCTX(I)=VCTX(I)+PX(J)*DV(J,I) VCTY(I)=VCTY(I)+PY(J)*DV(J,I) ENDDO VC(I)=VC(I)+W*VCT(I) DVCX(I)=DVCX(I)+WR*VCTX(I)+WQ*VCTY(I)+DWX*VCT(I) DVCY(I)=DVCY(I)-WQ*VCTX(I)+WR*VCTY(I)+DWY*VCT(I) ENDDO RETURN END SUBROUTINE ROTATE(R3,XX) DIMENSION R3(3,3),XX(3),VV(3) CALL MULMM(R3,XX,VV,3,3,1,3,3,3) DO I=1,3 XX(I)=VV(I) ENDDO RETURN END