关于Fortran计算Cape值错误

Seaside

源文件用的是era5的0.25度的再分析资料，程序里一直到计算抬升凝结高度(LCL)处都是正常的，计算自由对流高度的地方（代码标红的地方），输出的lfcp和lfct都是0，导致下一步计算EL也都出问题，我实在是看不懂这个计算，还请用过的大神可以帮忙指点一下，有什么需要注意的地方吗？
是需要考虑资料缺测吗，如果要剔除缺测值的话，用什么方法比较好，直接跳过计算吗？因为我把前期计算的量都输出了看了一下，没有NAN也没有无穷大，所以即使是缺测的话应该也不会出现所有值都是0的情况吧，看了两天也没看明白到底怎么改
还有另一个子程序计算CAPE的代码（见附件），但是只能放一维的高度场进去算，我如果想算一个三维空间场我要怎么call这个子函数呀，是把z放最内层循环，用一个一维数组去读取他的高度场吗，因为我看到的子程序都是对某一个具体的点去计算，我不太明白这种读取一维的数组的子函数要怎么去写代码，有人用过类似的子函数吗，可以教教我吗，感激不尽！
或者NCL有没有计算CAPE的函数，我是再分析资料，不是wrf资料，我看到wrf可以直接计算，不知道再分析资料有没有类似的函数呢？
下面是程序源代码
···················································································································
program main
implicit none
integer,parameter :: x=201,y=161,z=27,nt=32
real,parameter :: g=9.8,a=7.5,b=237.3,rd=287
integer i,j,k,it
real tk(x,y,z)       !温度
real f(x,y,z)      !相对湿度
real q(x,y,z)      !比湿
real citase(x,y,z)  !假相当位温
real tdk(x,y,z)     !露点
real p(z),t(x,y,z),td(x,y,z),fp(x,y,z),test(x,y)
real ts,tc(x,y),thse(x,y),ty,ec(x,y),qc(x,y),pc(x,y),tt,pp,total(x,y),cape,e,ps,et,qt,ete,qte
real tpp,tppv,tce,tde,tenv,elt(x,y),lfct(x,y),elp(x,y),lfcp(x,y),ey,tdy,py,l,m,n,tota
data p/1000, 975, 950, 925, 900, 875, 850, 825, 800, 775, 750, 700, 650, 600, 550, 500, 450, 400, 350, 300, 250, 225, 200, 175, 150, 125, 100/
!-------------------------read
open(11,file='rhtk.grd',form='binary')
open(18,file='cape.grd',form='binary')
do it=1,nt
read(11) (((tk(i,j,k),i=1,x),j=1,y),k=1,z)
read(11) (((f(i,j,k),i=1,x),j=1,y),k=1,z)
print*,"1"
!----------------------------------cal thse & td
   do k=1,z
        do j=1,y
                do i=1,x
                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))
                enddo
        enddo
enddo
print*,"3"
do k=1,z
do j=1,y
do i=1,x
                t(i,j,k)=tk(i,j,k)-273.16
                td(i,j,k)=tdk(i,j,k)-273.16
                fp(i,j,k)=f(i,j,k)   
enddo
enddo
enddo
print*,"4"
!----------------cal LCL p.t.thse 抬升凝结高度
do j=1,y
do i=1,x
l=p(1)
m=t(i,j,1)
n=td(i,j,1)
tc(i,j)=m-0.976*(m-n)/(0.976-0.000833*((237.3+n)**2)/(273.16+n))  !LCL温度
pc(i,j)=l*((273.16+tc(i,j))/(273.16+m))**3.50
if(tc(i,j)>=0.0)then
ec(i,j)=6.11*exp((17.67*tc(i,j))/(b+tc(i,j))) !水汽压
else
ec(i,j)=6.11*exp((21.87*tc(i,j))/(265.5+tc(i,j)))
endif
qc(i,j)=0.622*ec(i,j)/(pc(i,j)-0.378*ec(i,j)) ! q=0.622*E/(p-0.378E)   
thse(i,j)=(tc(i,j)+273.16)*exp(0.28586*LOG(1000.0/pc(i,j))+2500.0*(qc(i,j)/(338.52-0.24*(tc(i,j)+273.16)+1.24*tc(i,j))))
enddo
enddo
print*,"6"
!---------------------------cal LFC 自由对流高度
do j=1,y
do i=1,x
loop1: do py=pc(i,j),50,-0.5    !pc抬升凝结高度汽压
        do k=1,z-1
        if(py>p(k+1).and.py<p(k))then
        ty=(t(i,j,k)-t(i,j,k+1))*(py-p(k))/(p(k)-p(k+1))+t(i,j,k)
        endif
        enddo
call ttt(py,thse(i,j),-50.,tc(i,j),ts)
if(abs(ts-ty)<0.01)then
exit loop1
endif
enddo loop1
        if(py<p(z))then
        cape=0.
         else
lfcp(i,j)=py
lfct(i,j)=ty
         endif
write(*,*)"lfcp is:",lfcp(i,j),"lfct is:",lfct(i,j)
enddo
enddo
!---------------cal EL 平衡高度
do j=1,y
do i=1,x
loop2: do py=lfcp(i,j)-5,50,-0.5
        do k=1,z-1
        if((py-p(k))*(py-p(k+1))<=0)then
        ty=(t(i,j,k)-t(i,j,k+1))*(py-p(k+1))/(p(k)-p(k+1))+t(i,j,k+1)
        endif
   
        enddo
call ttt(py,thse(i,j),-50.,lfct(i,j),ts)
if(abs(ts-ty)<0.01)then
exit loop2
endif
enddo loop2

        if(py<p(z))then
     elp(i,j)=p(z)
         elt(i,j)=t(i,j,z)
        print *,"EL higher than the highest level"
else
elp(i,j)=py
elt(i,j)=ty
         endif
write(*,*)"elp is:",elp(i,j),"elt is:",elt(i,j)
enddo
enddo
!----------cal CAPE
do j=1,y
do i=1,x
tota=0.0
do pp=lfcp(i,j),elp(i,j),-0.5
        do k=1,z-1
        if((pp-p(k))*(pp-p(k+1))<=0)then
        tt=(t(i,j,k)-t(i,j,k+1))*(pp-p(k+1))/(p(k)-p(k+1))+t(i,j,k+1)
        endif
        enddo
        call ttp(pp,thse,tt,tpp)
  tppv=(tpp+273.16)*(1+0.378*(6.11*10**(a*tpp/(b+tpp)))/pp)
  call tenvironment(pp,z,t,td,p,tce,tde)
  tenv=(tce+273.16)*(1+0.378*(6.11*10**(a*tde/(b+tde)))/pp)
  cape=(rd*(tppv-tenv)/pp)*0.5
  tota=tota+cape
end do
total(i,j)=tota
!100 write(*,*)"total_cape=",total(21,2),"J/KG"
   
  print*,total(i,j)
enddo
enddo
write(18)((total(i,j),i=1,x),j=1,y)
enddo
close(18)
close(11)
end

!---------------cal thse&td
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 30
      T0=273.16
        b=4.9283
        a=6764.9
      ET=6.11*((T0/td)**b)*exp(a/T0-a/td)
      qs=0.622*ET/(p-0.378*ET)
20     continue
30      return
end
!-------------------根据气压求湿绝热上升气块温度
subroutine ttt(pp,thse,elt,lfct,tt)
real pp,e,q,tt,lfct,a,b,elt,ths
a=7.5
b=237.3
do tt=elt,lfct,0.5
if(tt>=0.0)then
e=6.11*exp(17.67*tt/(b+tt)) !&Euml;&reg;&AElig;&ucirc;&Ntilde;&sup1;
else
e=6.11*exp(21.87*tt/(265.5+tt))
endif
q=0.622*e/(pp-0.378*e)
  ths=(tt+273.16)*exp(0.28586*LOG(1000.0/pp)+2500*q/(338.52-0.24*(273.16+tt)+1.24*tt))
   ths=ths-273.16
  if((abs(ths-thse)<0.1))         exit
end do
end
!-----------由气压求环境温度---------------------------------
subroutine tenvironment(pp,iz,tc,td,p,tce,tde)
integer iz,i
real pp,tc(iz),td(iz),p(iz),tce,tde
   do i=1,iz-1
       if((pp-p(i))*(pp-p(i+1))<=0) then
        tce=((tc(i)-tc(i+1))*(pp-p(i+1)))/(p(i)-p(i+1))+tc(i+1)
        tde=((td(i)-td(i+1))*(pp-p(i+1)))/(p(i)-p(i+1))+td(i+1)
       endif
   end do
  return         
end
!-------------------求某处气压湿绝热上升气块温度
subroutine ttp(pp,thse,tt,tpp)
real pp,e,q,tt,lfct,a,b,elt,ths
a=7.5
b=237.3
do tpp=tt-50.,tt+50.,0.5
if(tt>=0.0)then
e=6.11*exp(17.67*tt/(b+tt)) !水汽压
else
e=6.11*exp(21.87*tt/(265.5+tt))
endif
q=0.622*e/(pp-0.378*e)
  ths=(tt+273.16)*exp(0.28586*LOG(1000.0/pp)+2500*q/(338.52-0.24*(273.16+tt)+1.24*tt))
   ths=ths-273.16
  if((abs(ths-thse)<0.1))         exit
end do
end

Seaside

有大神可以帮帮忙看一下吗

TSY11235

最近也在计算CAPE值，可否加QQ交流？

Seaside

TSY11235 发表于 2020-7-22 15:07
最近也在计算CAPE值，可否加QQ交流？

我后来就没用这个程序了，实在没看懂，直接用fnl现成的cape了

TSY11235

Seaside 发表于 2020-7-27 15:24
我后来就没用这个程序了，实在没看懂，直接用fnl现成的cape了

我也舍弃了fortran，直接用python计算的。可以调用wrf-python或者metpy库进行计算。
https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.html#module-metpy.calc
https://wrf-python.readthedocs.io/en/latest/diagnostics.html

yahooman

TSY11235 发表于 2020-7-27 22:27
我也舍弃了fortran，直接用python计算的。可以调用wrf-python或者metpy库进行计算。
https://unidata.gi ...

python计算的速度可以吗？

Seaside

TSY11235 发表于 2020-7-27 22:27
我也舍弃了fortran，直接用python计算的。可以调用wrf-python或者metpy库进行计算。
https://unidata.gi ...

谢谢！ 我后面要用的话我再研究一下

candyzh

TSY11235 发表于 2020-7-27 22:27
我也舍弃了fortran，直接用python计算的。可以调用wrf-python或者metpy库进行计算。
https://unidata.gi ...

我也用metpy计算过cape值， 就是不知道准不准确耶。  还有对流抑制能量CIN我发现metpy和micaps计算的结果正负号相反， 不知道为什么会出现这种情况。

TSY11235

yahooman 发表于 2020-7-28 19:46
python计算的速度可以吗？

我的数据量较小，没有做过具体测试，您可以对比一下两种方法。但是python的好处就是调用起来比较方便。

TSY11235

candyzh 发表于 2020-8-5 02:09
我也用metpy计算过cape值， 就是不知道准不准确耶。  还有对流抑制能量CIN我发现metpy和micaps计算的结果 ...

目前我使用的是wrf-python，直接调用的。还没有验证过值的准确性。这是个需要注意的问题。猜想结果正负号应该是根据物理意义来的，去掉负号表示能量值的绝对值，加上负号表示“抑制”？不确定我表述清楚没，也许需要看看这两个库的源码。