📄 cape.f90
字号:
!---------------主程序开始
program main
implicit none
!integer,parameter :: x=61,y=34,z=21
integer,parameter :: x=81,y=71,z=21
real,parameter :: T0=273.16
integer i,j,k
real tk(x,y,z),f(x,y,z),h(x,y,z),u(x,y,z),v(x,y,z)
real citase(x,y,z),q(x,y,z),tdk(x,y,z)
real cape1(x,y),cin1(x,y),li1(x,y),bcape1(x,y),bic1(x,y),dcape1(x,y),ssi1(x,y)
real p(z),t(z),td(z),fp(z)
real te(120),ta(120)
data p/1000, 975, 950, 925, 900, 850, 800, 750, 700, 650, 600, 550, 500, 450, 400, 350, 300, 250, 200, 150, 100/
real CAPE,PLFC,PE,CIN,LI,PC,TC,BCAPE,BIC,DCAPEX,TCON,CCL,H0,shr,SSI
!-----------------------------------------------
!-------------------------从源文件读入数据
open(1,file='D:\learn\work\diagnosis\myNeed\0126_00\0126_00.dat',status='unknown',form='unformatted',access='direct',recl=x*y*z)
!open(1,file='D:\temp\need\ht1.dat',status='unknown',form='unformatted',access='direct',recl=x*y*z)
read(1,rec=1) (((tk(i,j,k),i=1,x),j=1,y),k=1,z)
read(1,rec=2) (((h(i,j,k),i=1,x),j=1,y),k=1,z)
read(1,rec=3) (((u(i,j,k),i=1,x),j=1,y),k=1,z)
read(1,rec=4) (((v(i,j,k),i=1,x),j=1,y),k=1,z)
read(1,rec=5) (((f(i,j,k),i=1,x),j=1,y),k=1,z)
close(1)
!-----------------------------------计算 thse & td
do i=1,x
do j=1,y
do k=1,z
call tem_se(f(i,j,k),tk(i,j,k),p(k),tdk(i,j,k),q(i,j,k),citase(i,j,k))
print*, tk(i,j,k),h(i,j,k),u(i,j,k),f(i,j,k),tdk(i,j,k)
pause
enddo
enddo
enddo
!--------------------calculate CAPE,CIN,LI,SSI
do i=1,x
do j=1,y
do k=1,z
t(k)=tk(i,j,k)-T0
td(k)=tdk(i,j,k)-T0
enddo
CALL CAL_CAPE(P,T,TD,z,CAPE,PLFC,PE,CIN,LI,PC,TC)
call CAL_BCAPE(P,T,TD,z,BCAPE)
call CAL_BIC(P,T,TD,z,BIC)
call CAL_DCAPE(P,T,TD,z,DCAPEX)
H0=3600
call CAL_SSI(p,H0,H(i,j,:),U(i,j,:),V(i,j,:),z,CAPE,shr,SSI)
cape1(i,j)=cape
cin1(i,j)=cin
li1(i,j)=li
bcape1(i,j)=bcape
bic1(i,j)=bic
dcape1(i,j)=dcapex
ssi1(i,j)=ssi
print*, cape1(i,j),cin1(i,j),li1(i,j),bcape1(i,j),bic1(i,j),dcape1(i,j),ssi1(i,j)
! pause
enddo
enddo
!---------------把诊断量写入目标文件
open(3,file='D:\learn\work\diagnosis\myNeed\0126_00\cape0126.dat',status='replace',form='unformatted',access='direct',recl=x*y)
!open(3,file='D:\temp\need\ht1_1.dat',status='replace',form='unformatted',access='direct',recl=x*y)
write(3,rec=1) ((cape1(i,j),i=1,x),j=1,y)
write(3,rec=2) ((cin1(i,j),i=1,x),j=1,y)
write(3,rec=3) ((li1(i,j),i=1,x),j=1,y)
write(3,rec=4) ((bcape1(i,j),i=1,x),j=1,y)
write(3,rec=5) ((bic1(i,j),i=1,x),j=1,y)
write(3,rec=6) ((dcape1(i,j),i=1,x),j=1,y)
write(3,rec=7) ((ssi1(i,j),i=1,x),j=1,y)
do k=1,z
write(3,rec=7+k) ((citase(i,j,k),i=1,x),j=1,y)
write(3,rec=7+z+k) ((q(i,j,k),i=1,x),j=1,y)
enddo
close(3)
stop
end program
!---------------计算假相当位温和露点的子程序
SUBROUTINE tem_se(f,t,p,td,q,se)
T0=273.16
d=4.9283
c=6764.9
ET=6.11*((T0/t)**d)*exp(c/T0-c/t)
qs=0.622*ET/(p-0.378*ET)
q=qs*f/100
call tem_td(q,t,p,td)
c=0.28586*alog((1000.0/p))
d=q/(338.52-0.24*t+1.24*td)
se=t*exp((c+2500*d))
return
end
SUBROUTINE tem_td(q,t,p,td)
real q,t,p,td,qs,ET
td=t
qs=q+1
do 20, while(q.le.qs)
td=td-0.02
if (td<200) goto 10
T0=273.16
d=4.9283
c=6764.9
ET=6.11*((T0/td)**d)*exp(c/T0-c/td)
qs=0.622*ET/(p-0.378*ET)
20 continue
10 return
end
!-----------------------主程序结束
SUBROUTINE CAL_CAPE(P,T,TD,LEV,CAPE,PLFC,PE,CIN,LI,PC,TC)
! THIS SUBROUTINE CALCULATE CONVECTIVE AVAILABLE POTENTIAL ENERGY
! AND IT IS USED WITH CAL_TC_PC,CAL_TA,CAL_TE,CAL_ENGY,
! IN WHICH SUBROUTINE CAL_TA AND CAL_TC_PC ARE USED WITH CAL_THSE
! INPUT T ======> TEMPERATURE (CELSIUS DEGREE)
! TD =====> DEW POINT (CELSIUS DEGREE)
! P ======> PRESSURE (HPA)
! LEV ====> LEVEL
! OUTPUT CAPE ===> EQUIVALENT POTENTIAL TEMPERATURE (J/KG)
! PLFC ===> FREE CONVECTION LEVEL
! PE =====> EMBALANCE LEVEL
! CIN ====> CONVECTIVE INHIBITION
! LI =====> LIFTING INDEX
! TEMP ARRAY
! TE =====> ENVIRONMENT TEMPERATURE (KELVIN)
! TA =====> PARCEL TEMPERATURE (KELVIN)
! TC
DIMENSION P(LEV),T(LEV),TD(LEV)
DIMENSION TE(120),TA(120)
REAL CAPE,TC,PC,THSE,PC0,TC0,TD0,THSE0
REAL PLFC,PE,LI,CIN
INTEGER LEV
! ASSUMING THE BOTTOM POINT AS THE LIFTING LAYER
THSE=220.
PC=P(1)
TC=T(1)
P0=P(1)
T0=T(1)
TD0=TD(1)
CALL CAL_TC_PC(P0,T0,TD0,PC0,TC0,THSE0)
PC=PC0
TC=TC0
CALL CAL_TE(P,T,TE,LEV)
CALL CAL_TA(P0,T0,PC,TC,THSE0,TA)
CALL CAL_PLFC_AND_PE(P,TE,TA,LEV,PLFC,PE)
CALL CAL_VIR_ENGY(P,T,TD,LEV,TE,TA,CAPE)
! CALL CAL_CIN1(plfc,P0,TE,TA,CIN)
CALL CAL_CIN(P0,TE,TA,CIN)
CALL CAL_LI(P0,TE,TA,LI)
! RETURN
END
SUBROUTINE CAL_TC_PC(P,T,TD,PC,TC,THSE)
! CALCULATE LIFTING CONDENSATION LAYER AND LIFTING CONDENSATION TEMPERATURE
! AND ALSO PSEUDO-POTENTIAL TEMPERATURE AT CONDESAATING POINT
! INPUT P ========> LIFTING PRESSURE (HPA)
! T ========> LIFTING TEMPERATURE (CELSIUS DEGREE)
! TD =======> DEW POINT (CELSIUS DEGREE)
! OUTPUT
! PC =======> CONDESATION PRESSURE
! TC =======> CONDENSATION TEMPERATURE
! THSE =====> PSEUDO EQUIVALENT POTENTION TEMPERATURE (KELVIN)
REAL P,T,TD,PC,TC,THSE
TC=T-(T-TD)*0.976/(0.976-0.000833*(237.3+TD)**2/(273.15+TD))
PC=P*((273.15+TC)/(273.15+T))**3.5
CALL CAL_THSE(PC,TC,THSE)
RETURN
END
SUBROUTINE CAL_THSE(P,T,THSE)
! 计算假相当位温，不考虑 0 度以下假相当位温计算公式的差别
REAL E,A,B,P,T,THSE
A=7.5*T/(237.3+T)
E=6.11*10.**A
B=1000.0/(P-E)
THSE=(273.16+T)*B**0.286*EXP((2.5*10.**6-2368.*T)*0.622*E/(1004.*(273.16+T)*(P-E)))
RETURN
END
SUBROUTINE CAL_THSE1(P,T,THSE)
! 本子程序考虑 0 度以下，假相当位温计算公式不同
REAL E,A,B,P,T,THSE
IF(T.GE.0)A=7.5*T/(237.3+T)
IF(T.LT.0)A=9.5*T/(265.5+T)
E=6.11*10.**A
B=1000.0/(P-E)
IF(T.GE.0) THSE=(273.16+T)*B**0.286*EXP((597.3+0.5*T)*0.622*E/(0.24*(273.16+T)*(P-E)))
IF(T.LT.0) THSE=(273.16+T)*B**0.286*EXP((667.0+0.5*T)*0.622*E/(0.24*(273.16+T)*(P-E)))
RETURN
END
SUBROUTINE CAL_TE(PP,TP,TE,KMAX)
! p ,T ,TE,LEV
! ********** CALCULATE TEMPERTURE AT EACH LEVEL *********
! 计算层结每隔 10HPA 的温度（对数气压线性内插）
! PARAMETER (KMAX=30)
DIMENSION PP(KMAX),TP(KMAX),TE(110)
REAL P_L
INTEGER I_NL,IN
DO K=1,110
TE(K)=-60.
ENDDO
! WRITE(*,*)(PP(K),K=1,KMAX)
! WRITE(*,*)(TP(K),K=1,KMAX)
! PAUSE 31
I_NL=KMAX
I=1
IN=1
101 P_L=INT(PP(I)/10.)*10.0
P_T=INT(PP(I+1)/10.)*10.+10.
10 TE(IN)=TP(I)+(TP(I)-TP(I+1))/ALOG(PP(I)/PP(I+1))*ALOG(P_L/PP(I))
IN=IN+1
P_L=P_L-10.
IF(P_T.GT.PP(I))IN=IN-1
IF(P_L.GE.P_T)GOTO 10
I=I+1
! IF(PP(I+1).GT.50.OR.(I+1).LT.(I_NL))THEN
IF((I+1).LE.(I_NL))THEN
GOTO 101
ELSE
GOTO 20
ENDIF
20 CONTINUE
IF((PP(I_NL).GE.100.).AND.(INT(PP(I_NL)/10.)*10.).EQ.(PP(I_NL))) TE(IN)=TP(I_NL)
RETURN
END
SUBROUTINE CAL_TA(PP0,TP0,PC,TC,THSE,TA)
! 计算状态曲线各整层高度（间隔 10HPA）上的温度
! INPUT PP0 ======> PRESSURE AT THE LIFTING LEVEL (HPA)
! TP0 ======> TEMPERATURE AT THE LIFTING LEVEL (CELSIUS DEGREE)
! PC =======> LIFTING CONDENSATION LEVEL (HPA)
! TC =======> LIFTING CONDESATION TEMPERATURE (CELSIUS DEGREE)
! THESE ====> PSEUDO EQUIVALENT POTENTIAL TEMPERATURE (KELVIN)
! OUTPUT ARRAY
! TA =======> PARCEL TEMPERATURE AT EACH LEVEL (CELSIUS DEGREE)
REAL PP0,TP0,PC,TC,THSE,THSE1,THSE2,P_L,T1,T2,RD,CPD
INTEGER KPC
DIMENSION TA(110)
RD=287.04
CPD=1005.
DO K=1,110
TA(K)=-130.
ENDDO
!C ************** DRY ADIABATIC PROCESS BEGIN
KPC=INT(PP0/10)-INT(PC/10)
P_L=INT(PP0/10.)*10.
THD=(TP0+273.15)*(1000./PP0)**(RD/CPD)
DO I=1,KPC
TA(I)=THD*(P_L/1000.)**(RD/CPD)-273.15
P_L=P_L-10.
ENDDO
! DO I=1,KPC
! TA(I)=TP0+(TC-TP0)/ALOG(PC/PP0)*ALOG(P_L/PP0)
! P_L=P_L-10.
! ENDDO
! ************** DRY ADIABATIC PROCESS END
IN=KPC
T1=TC
P1=PC
! WRITE(*,*)' IN = ',IN,' TC = ',TC,' PC = ',PC
KZERO=1
! KZERO IS A JUDGEMENT FOR TEMPERATURE
!*************** MOIST ADIABATIC PROCESS BEGIN
10 CALL CAL_THSE(P_L,T1,THSE1)
CALL CAL_THSE(P_L,T1-0.1,THSE2)
T2=T1-(THSE1-THSE)/(THSE1-THSE2)*0.1
! NOTICE DIFFERENCE OF THSE WHEN TEMPERERTURE OBOVE OR BELEOW ZERO
IF(T2.LT.0.AND.KZERO.EQ.1)THEN
CALL CAL_THSE(P_L,T2,THSE)
KZERO=0
ENDIF
IN=IN+1
TA(IN)=T2
T1=T2
P_L=P_L-10.
! WRITE(*,*)' T1 = ',T1
! PAUSE
IF(P_L.GE.90.)GOTO 10
! *************** MOIST ADIABATIC PROCESS END
RETURN
END
SUBROUTINE CAL_PLFC_AND_PE(P,TE,TA,LEV,PLFC,PE)
! THIS SUBROUTINE CALCULATE LAYER OF FREE CONVECTION AND EMBALANCE LAYER
! NOTICE THIS SUBROUTINE IS NOT INDEPENDENT!!!
! INPUT TE =====> ENVIRONMENT TEMPERATURE (KELVIN)
! TA =====> PARCEL TEMPERATURE (KELVIN)
! P =====> PRESSURE (HPA)
! LEV ====> NUMBER OF LAYER
! OUTPUT PLFC ===> FREE CONVECTION LAYER (HPA)
! PE =====> EMBALANCE PRESSURE (HPA)
REAL PLFC,PE
DIMENSION TE(110),TA(110),P(LEV)
INTEGER LEV,KLFC
REAL DELT1,DELT2
! WRITE(*,'(10F8.2)')(TA(K),K=1,10),(TE(K),K=1,10)
KLFC=0
! IF(TA(1).GE.TE(1).AND.TA(2).GT.TE(2))THEN
IF(TA(1).GT.TE(1).OR.(TA(1).EQ.TE(1).AND.TA(2).GT.TE(2)))THEN
PLFC=P(1)
ELSEIF((TA(1).EQ.TE(1).AND.TA(2).LE.TE(2)).OR.TA(2).LE.TE(2))THEN
P_L=INT(P(1)/10.)*10
P_T=P_L-10
DO K=2,90
IF(TA(K).GT.TE(K))THEN
KLFC=K
P_B=P_T+10
DELT1=TE(K-1)-TA(K-1)
DELT2=TA(K)-TE(K)
PLFC=(P_B*DELT2+P_T*DELT1)/(DELT1+DELT2)
GOTO 290
ENDIF
P_T=P_T-10
ENDDO
IF(KLFC.EQ.0)PLFC=99999.9
ENDIF
290 CONTINUE
NB=INT(P(1)/10.)
NT=INT((P(LEV)-0.001)/10.)+1
NE=NB-NT+1
KPE=NE
IF(INT(P(LEV))/10..EQ.P(LEV)/10.)THEN
P_T=P(LEV)
ELSE
P_T=INT(P(LEV)/10.)*10+10.
ENDIF
IF(TE(NE).LT.TA(NE))THEN
PE=99999.9
ELSE
DO K=NE-1,1,-1
IF(TE(K).LT.TA(K))THEN
KPE=K
P_B=P_T+10
DELT1=TE(K+1)-TA(K+1)
DELT2=TA(K)-TE(K)
PE=(P_B*DELT2+P_T*DELT1)/(DELT1+DELT2)
GOTO 390
ENDIF
P_T=P_T+10
ENDDO
IF(KPE.EQ.NE)PE=99999.9
ENDIF
390 CONTINUE
RETURN
END
SUBROUTINE CAL_VIR_ENGY(P,T,TD,LEV,TE,TA,CAPEV)
!
! NOTICE IT IS SUPPORTED BY CAL_TDE,CAL_Q AND CAL_QS
! INPUT T ======> TEMPERATURE (CELSIUS DEGREE)
! TD =====> DEW POINT (CELSIUS DEGREE)
! P ======> PRESSURE (HPA)
! LEV ====> LEVEL
! TE =====> ENVIRONMENT TEMPERATURE (CELSIUS DEGREE)
! TA =====> PARCEL TEMPERATURE (CELSIUS DEGREE)
! OUTPUT
! CAPEV ==> CONVECTIVE AVAILABLE POTENTIAL ENERGY (J/KG)
! TEMP ARRAY
! PDS ====> PRESSURE (HPA)
! QE =====> SPECIFIC HUMICITY
! QSE ====> SATURATE SPECIFIC HUMIDITY
! TDE ====> DEW POINT (CELSIUS DEGREE)
REAL CAPEV,R,P_L
DIMENSION TE(110),TA(110),TDE(110),P(LEV),T(LEV),TD(LEV)
DIMENSION PDS(110),QE(110),QSE(110)
DO K=1,LEV
IF(TD(K).GT.90)TD(K)=T(K)-30.
ENDDO
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -