用NCL编写计算CTT的程序

denglin19901030

由于本机没有安装6.3.0，所以想仿照气象家园中的有关6.3.0里的CTT计算贴，用NCL 编写一个脚本出来，但是编完之后运行总是报错，请各位大神帮忙瞅瞅哪里有问题呢？以下是根据NCL6.3.0里计算CTT的帖子编的：load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
;*************************************************

begin

eps = 0.622
ussalr = 0.0065
rgas = 287.04
grav  = 9.81
abscoefi = 0.272
abscoef  = 0.145
celkel = 273.15
pf=new((/34,270,270/),float)
nc_files=systemfunc("ls wrfout_d02*")
n=dimsizes(nc_files)

do it=0,0
f=addfile (nc_files(it) , "r")
time1=chartostring(f->Times)
print(time1)
mdims = getfilevardimsizes(f,"P")
nd = dimsizes(mdims)
qcw     = f->QCLOUD(0,:,:,:)*1000                             ; get qc January
qvp     = f->QVAPOR(0,:,:,:)*1000
qci     = f->QICE(0,:,:,:)*1000
P     = f->P(0,:,:,:)
PB    = f->PB(0,:,:,:)
T     = f->T(0,:,:,:)
PHB   = f->PHB(0,:,:,:)
PH    = f->PH(0,:,:,:)
ter   = f->HGT(0,:,:)
pres  = (P+PB)
T     = T+300
tk    =wrf_tk(pres,T)
pres  = pres*0.01
stag_ght  = PH
stag_ght  = (PHB+PH)/9.81
ght       = wrf_user_unstagger(stag_ght,PHB@stagger)

;calculate the surface pressure
do i=0,269
do j=0,269
ratmix = 0.001*qvp(0,j,i)
arg1   = eps+ratmix
arg2   = eps*(1.+ratmix)
vt     = tk(0,j,i)*arg1/arg2
agl_hgt = ght(33,j,i)-ter(j,i)
arg1  = -grav/(rgas*ussalr)
pf(33,j,i)=pres(0,j,i)*(vt/(vt+ussalr*(agl_hgt)))^(arg1)
end do
end do

do i=0,269
do j=0,269
do k=0,32
ripk = 33-k+1
pf(k,j,i)=0.5*(pres(ripk-1,j,i)+pres(ripk-2,j,i))
end do
end do
end do

do i=0,269
do j=0,269
opdepthd = 0
do k=0,33
opdepthu=opdepthd
ripk=33-k
if(k.eq.0)then
dp=200.*(pf(0,j,i)-pres(33,j,i))
else
dp=100.*(pf(k,j,i)-pf(k-1,j,i))
end if
if(qci(k,j,i).eq.0)then
if(tk(k,j,i).lt.celkel)then
opdepthd = opdepthu+abscoefi*qcw(k,j,i)*dp/grav
else
opdepthd = opdepthu+abscoef*qcw(k,j,i)*dp/grav
end if
else
opdepthd=opdepthd+(abscoef*qcw(ripk,j,i)+abscoefi*qci(ripk,j,i)*dp/grav)
end if

if(opdepthd.lt.1.and.k.lt.33)then
continue
else if(opdepthd.lt.1.and.k.eq.33)then
prsctt=pres(0,j,i)
else
fac=(1.-opdepthu)/(opdepthd-opdepthu)
prsctt=pf(k-1,j,i)+fac*(pf(k,j,i)-pf(k-1,j,i))
compara=(/pres(0,j,i),max(pres(33,j,i)),prsctt/)
prsctt=min(compara)
end if
break
end do

do k=1,33
ripk=33-k
p1 = prs(ripk+1,j,i)
p2 = prs(ripk,j,i)
if(prsctt.ge.p1.and.prsctt.le.p2)then
fac=(prsctt-p1)/(p2-p1)
arg1=fac*(tk(ripk,j,i)-tk(ripk+1))-celkel
ctt(j,i)=tk(ripk+1,j,i)+arg1
break
end if
end do
end do
end do
end do
end


附件为气象家园@diandianncl6.3.0画ctt总是错误终于找到原因了？！！！所提供的脚本

denglin19901030

忘记贴错误图了

denglin19901030

denglin19901030 发表于 2015-10-23 20:18
忘记贴错误图了

问题已经解决，是没有加入END IF 的问题，NCL和Fortran的差异，但是算出来的值依旧有问题，查看了WRFUserARW_ctt.ncl,发现里面调用的语句和wrf_fctt有差异，是当中还缺少了哪个转换脚本吗？？