二维贴体坐标和常用坐标的转换问题

[渐渐成长]

课题组有一个较老的计算程序，使用的是贴体坐标系。控制方程也需要进行转换在贴体坐标系下相应的形式。如下
里面的alpha和beta等参数定义如下

XI和ETA是贴体坐标系的两个方向，我在求以上jacobi,beta,gama,alpha时遇到一些问题，请给位老看掌掌眼。老程序应用的网格中beta等参数没有负的，我的新网格会出现负的。但是从公式来看，确实存在负数的可能性。不得其解。我的网格是

参数求解的代码如下：

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!     计算JCB,ALPHA,BETA,GAMA,XXI,XETA,YXI,YETA等等
!     
!    在计算的过程中，边界采用前差，内点采用中心差分
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
PROGRAM MAIN
IMPLICIT NONE
INTEGER I,J,INDCOS
INTEGER NIC,NJC
INTEGER,PARAMETER :: IT=206,JT=92
INTEGER LIMIT(IT,JT),KONES(IT,JT)
REAL X(IT,JT),Y(IT,JT),JCB(IT,JT),ALPHA(IT,JT),BETA(IT,JT)
REAL GAMA(IT,JT),XXI(IT,JT),XETA(IT,JT),YXI(IT,JT),YETA(IT,JT)

!!!!!!!!!!!!网格文件的打开!!!!!!!!!!!!!
OPEN(10,FILE='ZhongXinJieMianErWei.DAT',STATUS='OLD')
READ(10,*)NIC,NJC
DO J=1,NJC
	DO I=1,NIC
		READ(10,*)X(I,J),Y(I,J),LIMIT(I,J),KONES(I,J)
	END DO
END DO
CLOSE(10)

!!!!!!!!!!!!!!JCB等的计算!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!拐角点(1,1)
XXI(1,1)=X(2,1)-X(1,1)
XETA(1,1)=X(1,2)-X(1,1)
YXI(1,1)=Y(2,1)-Y(1,1)
YETA(1,1)=Y(1,2)-Y(1,1)
!!!!!拐角点(1,NJC)
XXI(1,NJC)=X(2,NJC)-X(1,NJC)
XETA(1,NJC)=X(1,NJC)-X(1,NJC-1)
YXI(1,NJC)=Y(2,NJC)-Y(1,NJC)
YETA(1,NJC)=Y(1,NJC)-Y(1,NJC-1)
!!!!!拐角点(NIC,NJC)
XXI(NIC,NJC)=X(NIC,NJC)-X(NIC-1,NJC)
XETA(NIC,NJC)=X(NIC,NJC)-X(NIC,NJC-1)
YXI(NIC,NJC)=Y(NIC,NJC)-Y(NIC-1,NJC)
YETA(NIC,NJC)=Y(NIC,NJC)-Y(NIC,NJC-1)
!!!!!拐角点(NIC,1)
XXI(NIC,1)=X(NIC,1)-X(NIC-1,1)
XETA(NIC,1)=X(NIC,2)-X(NIC,1)
YXI(NIC,1)=Y(NIC,1)-Y(NIC-1,1)
YETA(NIC,1)=Y(NIC,2)-Y(NIC,1)

!!!!!左边界点
DO J=2,NJC-1
	XXI(1,J)=X(2,J)-X(1,J)
	YXI(1,J)=Y(2,J)-Y(1,J)
	XETA(1,J)=(X(1,J+1)-X(1,J-1))/2.0
	YETA(1,J)=(Y(1,J+1)-Y(1,J-1))/2.0
END DO
!!!!!上边界点
DO I=2,NIC-1
	XXI(I,NJC)=(X(I+1,NJC)-X(I-1,NJC))/2.0
	YXI(I,NJC)=(Y(I+1,NJC)-Y(I-1,NJC))/2.0
	XETA(I,NJC)=X(I,NJC)-X(I,NJC-1)
	YETA(I,NJC)=Y(I,NJC)-Y(I,NJC-1)
END DO
!!!!!下边界点
DO I=2,NIC-1
	XXI(I,1)=(X(I+1,1)-X(I-1,1))/2.0
	YXI(I,1)=(Y(I+1,1)-Y(I-1,1))/2.0
	XETA(I,1)=X(I,2)-X(I,1)
	YETA(I,1)=Y(I,2)-Y(I,1)
END DO
!!!!!右边界点
DO J=2,NJC-1
	XXI(NIC,J)=X(NIC,J)-X(NIC-1,J)
	YXI(NIC,J)=Y(NIC,J)-Y(NIC-1,J)
	XETA(NIC,J)=(X(NIC,J+1)-X(NIC,J-1))/2.0
	YETA(NIC,J)=(Y(NIC,J+1)-Y(NIC,J-1))/2.0
END DO

!!!!!内部点
DO J=2,NJC-1
	DO I=2,NIC-1
		XXI(I,J)=(X(I+1,J)-X(I-1,J))/2.0
		XETA(I,J)=(X(I,J+1)-X(I,J-1))/2.0
		YXI(I,J)=(Y(I+1,J)-Y(I-1,J))/2.0
		YETA(I,J)=(Y(I,J+1)-Y(I,J-1))/2.0
	END DO
END DO
!!!!!!几何参数的计算
DO J=1,NJC
	DO I=1,NIC
		ALPHA(I,J)=XETA(I,J)*XETA(I,J)+YETA(I,J)*YETA(I,J)
		GAMA(I,J)=XXI(I,J)*XXI(I,J)+YXI(I,J)*YXI(I,J)
		BETA(I,J)=XXI(I,J)*XETA(I,J)+YXI(I,J)*YETA(I,J)
		JCB(I,J)=XXI(I,J)*YETA(I,J)-XETA(I,J)*YXI(I,J)
	END DO
END DO


!!!!!!!!!!!!!!!!!!!!!!!!!!结果的输出!!!!!!!!!!!!!!!!!!!!!!
OPEN(11,FILE='GRIDJCB.DAT',STATUS='REPLACE')
DO J=1,NJC
	DO I=1,NIC
		WRITE(11,*)JCB(I,J),ALPHA(I,J),BETA(I,J),GAMA(I,J),XXI(I,J),&
					XETA(I,J),YXI(I,J),YETA(I,J)
	END DO
END DO
CLOSE(11)


!!!!!!!!!!!!!!!!!!!!!!!!结果检查!!!!!!!!!!!!!!!!!!!!!!!!!
OPEN(12,FILE='CheckJCB.DAT',STATUS='REPLACE')
DO J=1,NJC
	DO I=1,NIC
		IF(JCB(I,J).LT.0)WRITE(12,*)"JCB",I,J,JCB(I,J)
		IF(ALPHA(I,J).LT.0)WRITE(12,*)"ALPHA",I,J,ALPHA(I,J)
		IF(BETA(I,J).LT.0)WRITE(12,*)"BETA",I,J,BETA(I,J)
		IF(GAMA(I,J).LT.0)WRITE(12,*)"GAMA",I,J,GAMA(I,J)
		IF(XXI(I,J).LT.0)WRITE(12,*)"XXI",I,J,XXI(I,J)
		IF(XETA(I,J).LT.0)WRITE(12,*)"XETA",I,J,XETA(I,J)
		IF(YXI(I,J).LT.0)WRITE(12,*)"YXI",I,J,YXI(I,J)
		IF(YETA(I,J).LT.0)WRITE(12,*)"YETA",I,J,YETA(I,J)
	END DO
END DO
CLOSE(12)
END PROGRAM

[ibelief]

Fortran 代码，令人感动

[李东岳]

Fortran 代码，令人感动