百帕/秒。
2.3.5子程序
2.3.5.1运动学法求OMEGA的子程序
SUBROUTINE YUNDF(U,V,P,DL,DM,F0,L,M,N,OMEGA) REAL(8),DIMENSION(L,M,N)::U,V,DIV,OMEGA,OMEGA1 REAL(8),DIMENSION(L,M,N)::OMT REAL(8),DIMENSION(N)::P
!OMT为由热力学法算出的顶层OMEGA值
REAL(8)::F0,DL,DM
OMT=0 !这里取OMT=0，如有资料可用热力学法计算 OMEGA1=0
CALL DIVER(U,V,P,DL,DM,F0,L,M,N,DIV) DO I=1,L DO J=1,M DO K=2,N
OMEGA1(I,J,K)=OMEGA1(I,J,K-1)-(DIV(I,J,K-1)+DIV(I,J,K))/2*(P(K)-P(K-1))
END DO DO K=1,N
OMEGA(I,J,K)=OMEGA1(I,J,K)-K*(K-1)/(N*(N-1))*(OMEGA1(I,J,N)-OMT(I,J,N)) END DO
END DO END DO END
2.3.5.2求解准地转OMEGA方程的子程序
SUBROUTINE OME(U,V,T,P,L,M,N,F0,DL,DM,OMEGA) !求解准地转OMEGA方程 REAL(8),DIMENSION(L,M,N)::U,V,FAI,KAP,SIGMA REAL(8),DIMENSION(L,M,5)::T,TP
REAL(8),DIMENSION(L,M,N)::OMEGA1,OMEGA2,OMEGA,OMEGA0 REAL(8),DIMENSION(L,M,N)::QA1,QA2,QA3,QB1,QB2,QB3 REAL(8),DIMENSION(N)::P,SIGMAV REAL(8),DIMENSION(M)::FCO
REAL(8)::EPS,CP,RD,F0,DL,DM,ALF EPS=1.E-5 !迭代要求的精度 ALF=1.44 RD=287 CP=1005 OMEGA1=0
OMEGA2=0
CALL LHS(U,V,P,L,M,N,F0,DL,DM,FAI,KAP) DO J=1,M !科氏参数
FCO(J)=2*(7.29E-5)*SIN(F0+(J-1)*DM) END DO
!中间变量
! 计算静力稳定度
CALL AP1(T,L,M,N,P,TP) DO K=1,N DO J=1,M DO I=1,L
SIGMA(I,J,K)=-(RD/P(K)*(TP(I,J,K)-RD/CP*T(I,J,K)/P(K))) END DO END DO END DO DO K=1,N DO I=1,L
DO J=1,M SIGMAV(K)=SIGMAV(K)+SIGMA(I,J,K)
SIGMAV(K)=0
END DO END DO END DO DO K=1,N DO J=1,M
DO I=1,L
IF(SIGMA(I,J,K)<=0)THEN
SIGMA(I,J,K)=SIGMAV(K) END IF END DO END DO !
END DO
计算右边第一项
CALL LAPLACE(FAI,L,M,N,F0,DL,DM,QA1) DO K=1,N DO J=1,M
DO I=1,L
QA1(I,J,K)=QA1(I,J,K)+FCO(J) END DO END DO
END DO
CALL JACOB(FAI,QA1,DL,DM,F0,L,M,N,QA2) CALL AP1(QA2,L,M,N,P,QA3) QA3=QA3*FCO(M/2) N0=0 DO
N0=N0+1
OMEGA0=OMEGA1
CALL ERROR(OMEGA1,QA3,P,SIGMA,L,M,N,F0,FCO,DL,DM,ALF) IF(MAXVAL(ABS(OMEGA1-OMEGA0)) ! END DO open(2,file='omega1.dat') WRITE(2,'( 计算右边第二项 CALL AP1(FAI,L,M,N,P,QB1) CALL JACOB(FAI,QB1,DL,DM,F0,L,M,N,QB2) CALL LAPLACE(QB2,L,M,N,F0,DL,DM,QB3) QB3=-FCO(M/2)*QB3 N1=0 DO N1=N1+1 OMEGA0=OMEGA2 CALL ERROR(OMEGA2,QB3,P,SIGMA,L,M,N,F0,FCO,DL,DM,ALF) IF(MAXVAL(ABS(OMEGA2-OMEGA0)) open(3,file='omega2.dat') WRITE(3,'( OMEGA=OMEGA1+OMEGA2 END SUBROUTINE AP1(A,L,M,N,P,FF) REAL(8),DIMENSION(L,M,N)::A REAL(8),DIMENSION(L,M,N)::FF REAL(8),DIMENSION(N)::P DO J=1,M DO I=1,L DO K=2,N-1 FF(I,J,K)=((P(K)-P(K-1))*A(I,J,K+1)+(P(K+1)+P(K-1)-2*P(K))*A(I,J,K) & !计算垂直方向一次不等距差分 -(P(K+1)-P(K))*A(I,J,K-1))/(2*(P(K)-P(K-1))*(P(K+1)-P(K))) END DO END DO END DO END SUBROUTINE LAPLACE(A,L,M,N,F0,DL,DM,QA) REAL(8),DIMENSION(L,M,N)::A,QA REAL(8)::F0,DL,DM,AA AA=6.371E6 L1=L-1 M1=M-1 DO K=1,N DO I=2,L1 DO J=2,M1 QA(I,J,K)=(A(I+1,J,K)+A(I-1,J,K)-2*A(I,J,K))/(DL*COS(F0+(J-1) & *DM))**2+(A(I,J+1,K)+A(I,J-1,K)-2*A(I,J,K))/(DM*DM) & -TAN(F0+(J-1)*DM)*(A(I,J+1,K)-A(I,J-1,K))/(2*DM) END DO END DO END DO QA=QA/(AA*AA) CALL BOUND(QA,L,M,N) END SUBROUTINE JACOB(A,B,DL,DM,F0,L,M,N,JC) !采用两种JACOB差分方法的平均 REAL(8),DIMENSION(L,M,N)::A,B,JC REAL(8)::F0,DL,DM,J1,J2 REAL(8)::AA=6.371E6 L1=L-1 M1=M-1 DO K=1,N DO J=2,M1 DO I=2,L1 J1=-((B(I+1,J,K)-B(I-1,J,K))*(A(I,J+1,K)-A(I,J-1,K)) & -(B(I,J+1,K)-B(I,J-1,K))*(A(I+1,J,K)-A(I-1,J,K))) & /(AA*AA*COS(F0+(J-1)*DM)*4*DL*DM) J2=-(B(I+1,J+1,K)*(A(I,J+1,K)-A(I+1,J,K)) & -B(I-1,J-1,K)*(A(I-1,J,K)-A(I,J-1,K)) & -B(I-1,J,K)*(A(I,J+1,K)-A(I-1,J,K)) & +B(I+1,J-1,K)*(A(I+1,J,K)-A(I,J-1,K))) & /(AA*AA*COS(F0+(J-1)*DM)*4*DL*DM) JC(I,J,K)=0.6*J1+0.4*J2 END DO END DO END DO CALL BOUND(JC,L,M,N) END SUBROUTINE ERROR(A,Q,P,SI,L,M,N,F0,FCO,DL,DM,ALF) REAL(8),DIMENSION(L,M,N)::A,Q,SI REAL(8),DIMENSION(N)::P REAL(8),DIMENSION(M)::FCO REAL(8)::F0,DL,DM,AA,EPS,ALF,R EPS=1.E3 AA=6.371E6 L1=L-1 M1=M-1 N1=N-1 DO K=2,N1 DO I=2,L1 DO J=2,M1 R=((A(I+1,J,K)+A(I-1,J,K))/(DL*COS(F0+(J-1)*DM))**2 & +(A(I,J+1,K)+A(I,J-1,K))/DM**2-TAN(F0+(J-1)*DM)* & (A(I,J+1,K)-A(I,J-1,K))/(2*DM)+(AA*FCO(M/2))**2* & ((P(K)-P(K-1))*A(I,J,K+1)+(P(K+1)-P(K))*A(I,J,K-1)) & /(SI(I,J,K)*(P(K)-P(K-1))**2*(P(K+1)-P(K))) & -AA*AA*Q(I,J,K)/SI(I,J,K))/(2/(DL*COS(F0+(J-1)*DM))**2 & +2/DM**2+(AA*FCO(M/2))**2*(P(K+1)-P(K-1)) & /(SI(I,J,K)*(P(K)-P(K-1))**2*(P(K+1)-P(K)))) A(I,J,K)=ALF*R+(1-ALF)*A(I,J,K) END DO END DO END DO END