Ran Xue

···

undef("def_dAdt")
;输入位涡的气候态和异常场
;数据：time lat lon 
function def_dAdt(pv_cli[*][*][*]:numeric,pv_ano[*][*][*]:numeric,youwantlevel)

local pi,dName,lon,lat,level,gradpv,pv_grad_lat,pv_grad_lon,A,opt,dt,t,gradpv2,dQdy,time

begin

pi      = atan(1.0)*4.
dName   = getvardims(pv_cli)

if (any(ismissing(dName(1:)))) then
print("fatal:  def_DADt: requires that all the rightmost dimensions be named")
        exit
end if

lat     = pv_cli&$dName(1)$
lon     = pv_cli&$dName(2)$
; cosine
coslat  = cos(lat(:)*pi/180.)
coslat@_FillValue = default_fillvalue(typeof(coslat))
coslat  = (/where(coslat.le.0.,coslat@_FillValue,coslat)/)
; 1-D -> 3-D
coslattmp   = conform_dims(dimsizes(pv_cli),coslat,1)
;------------------------------------------------------

; dQdx    = center_finite_diff_n(pv_cli,lon*pi/180., False,0,2)
dQdy    = center_finite_diff_n(pv_cli,lat*pi/180., False,0,1)
copy_VarMeta(pv_cli, dQdy)

if (any(dQdy.eq.0)) then
    print("your data: dQdy have zero!!!!")
end if

;------------------------------------------------------
gradpv       = grad_latlon_cfd(pv_cli, pv_cli&latitude, pv_cli&longitude, True, False)
pv_grad_lat  = gradpv[0]    
pv_grad_lon  = gradpv[1]

copy_VarMeta(pv_cli, pv_grad_lat)
copy_VarMeta(pv_cli, pv_grad_lon)

printVarSummary(pv_grad_lat)
printMinMax(pv_grad_lat, 0)

if (any(pv_grad_lat.eq.0)) then
    print("your data: pv_grad_lat have zero!!!!")
end if

if (any(pv_grad_lon.eq.0)) then
    print("your data: pv_grad_lat have zero!!!!")
end if


pv_grad_lat =where(pv_grad_lat.eq.0, pv_grad_lat@_FillValue, pv_grad_lat)
pv_grad_lon =where(pv_grad_lon.eq.0, pv_grad_lon@_FillValue, pv_grad_lon)
dQdy        =where(dQdy.eq.0,               dQdy@_FillValue, dQdy)

;------填补缺失值
guess     = 1                ; use zonal means
is_cyclic = True             ; cyclic [global]
nscan     = 1500             ; usually much less than this
eps       = 1.e-2            ; variable dependent
relc      = 0.6              ; relaxation coefficient
opt       = 0                ; not used

; poisson_grid_fill( pv_grad_lat, is_cyclic, guess, nscan, eps, relc, opt)
; poisson_grid_fill( pv_grad_lon, is_cyclic, guess, nscan, eps, relc, opt)
; poisson_grid_fill( dQdy, is_cyclic, guess, nscan, eps, relc, opt)

;------------------------------------------------------

re      = 6.37*10^6

gradpv2=sqrt(pv_grad_lat^2+pv_grad_lon^2)
copy_VarMeta(pv_cli, gradpv2)


if (any(gradpv2.eq.0)) then
    print("your data: gradpv2 have zero!!!!")
end if
gradpv2=where(gradpv2.eq.0, gradpv2@_FillValue, gradpv2)
;poisson_grid_fill( gradpv2, is_cyclic, guess, nscan, eps, relc, opt)
print("-----gradpv2-----")
printMinMax(gradpv2, 0)

;-----3种方法计算A
;A=(youwantlevel*pv_ano^2)/(2.*1000.*gradpv2)

A=(youwantlevel*pv_ano^2*re*coslattmp^2)/(2.*1000.*dQdy)

;A=(1.225*pv_ano^2)/(2.*dQdy)

copy_VarMeta(pv_cli, A)
printVarSummary(A)


print("-----pv_ano-----")
printMinMax(pv_ano, 0)
print("-----pv_ano^2-----")
printMinMax(pv_ano^2, 0)
print("-----dQdy-----")
printMinMax(dQdy^2, 0)
print("-----A-----")
printMinMax(A, 0)

;utc     =cd_calendar(pv_cli&$dName(0)$, 2)
;day     = tointeger(utc(:,2))
;day=ispan(0, 60, 1)
;print(utc)

;---Convert "hours since ..." to "seconds since ..."

time    = pv_cli&time   
print(time@units);hours since 1900-01-01 00:00:00

dt = 24.*3600.
t = time*3600.


if(any(ismissing(A))) then
    print("Your data: A is missing.")
end if
poisson_grid_fill( A, is_cyclic, guess, nscan, eps, relc, opt)

dAdt    = tofloat(center_finite_diff_n(A,t,False,0,0));最后一个是时间维 单位s
copy_VarMeta(pv_cli, dAdt)
printVarSummary(dAdt)
printMinMax(dAdt, 0)


return(dAdt)

end

···