- 积分
- 101
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2022-9-9
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
各位同仁，我不是气象专业背景的，但是文章的一个部分需要做moisture budget，因此需要计算vertically integrated moisture flux convergence (VIMFC)。数据的精度是逐小时的，老板的要求做到多年平均之后降水趋势项，也就是moisture storage为零。
参考了https://mp.weixin.qq.com/s/1Pco6yXf1v-91Dse_0gdyQ等代码之后，我目前的代码如下，但是有些问题：（1）模型本身采用atmosphere_hybrid_sigma_pressure_coordinate，能不能直接在这个气压水平上进行积分；还是说需要重新插值，用于插值的大气压力水平应该如何选择，才能规避掉青藏高原那边的问题。（2）这个究竟要做到多久的时间段，才能看到moisture storage近乎零。（3）以及其他任何关于原理、单位的问题，谢谢大家的任何建议！
;*************************************************
; mfc_div_1.ncl
;
; Concepts illustrated:
; - Read daily mean wind components, humidity and sfc. pressure
; from different files
; - Reorder the input (N==>S) grid order to (S==>N) via NCL syntax ::-1
; - Calculate mass weightined layer thickness [units="kg/m2"]
; - Calculate moisture flux [uq, vq]
; - Calculate moisture flux divergence using spherical harmonics
; - Integrate the moisture flux divergence using mass weighting
; - Plot a number of quantities
;*************************************************
;---Calculate the Horizontal Moisture Flux Convergence [MFC]
;*************************************************
;---High frequency source data: hourly/3hr/6hr/12hr/daily .... NOT monthly values
;---References:
;---http://www.cgd.ucar.edu/cas/catalog/newbudgets/
;---http://tornado.sfsu.edu/geosciences/classes/e260/AtmosphericRivers/Moisture%20Flux.pdf
;---https://www.spc.noaa.gov/publications/banacos/mfc-sls.pdf
;===================================================================
; Data Source: ESRL Physical Sciences Division
; https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.html
; NCEP Reanalysis data provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA,
; from their Web site at https://www.esrl.noaa.gov/psd/
;===================================================================
;;;先计算散度再积分
ptop = 0. ; 'shum' upper level
ptop@units = "hPa"
g = 9.80665 ; m/s2
pbot = 110000. ;to ensure integral from psfc
;---ESRL: CDC data
dir = ""
fils = systemfunc("cd "+dir+" ; ls *.cam.h1.*.nc") ; all files beginning with 'ACCESS_SRF'
f1 = addfiles(dir+fils, "r") ; read variables from all files
a = ListCount(f1)
ref = addfile("508-02.cam.h0.2055-09.nc","r")
P0mb = 1000.
lat = ref->lat
;lev_p = (/0,50,100,150,200,250,300,350,400,450,500,550,600,650,700,750,800,850,900,950,1000/)
lev_p = (/0,50,100,150,200,250,300,350,400,450,500,550/)
lev_p@units = "hPa" ; required for vinth2p
lev_p@positive = "down"
print(a)
opt = True
opt@nval_crit = 10 ; require at least 10 values
ndim = 0
do m = 0,a-1
;uu := calculate_daily_values (f1[m]->U, "avg", ndim, opt)
;vv := calculate_daily_values (f1[m]->V, "avg", ndim, opt)
;qq := calculate_daily_values (f1[m]->Q, "avg", ndim, opt)
;ps := calculate_daily_values (f1[m]->PS, "avg", ndim, opt)
;PRECC := calculate_daily_values (f1[m]->PRECC, "avg", ndim, opt)
;PRECL := calculate_daily_values (f1[m]->PRECL, "avg", ndim, opt)
;EVAP := calculate_daily_values (f1[m]->QFLX, "avg", ndim, opt)
uu := f1[m]->U
vv := f1[m]->V
qq := f1[m]->Q
ps := f1[m]->PS
PRECC := f1[m]->PRECC
PRECL := f1[m]->PRECL
EVAP := f1[m]->QFLX
hyam := f1[m]->hyam
hybm := f1[m]->hybm
intyp_m=1 ;1--linear,2--log
kxtrp = False ;Logical. False => no extrapolation when the pressure level is outside of the range of psfc.
u = vinth2p (uu, hyam,hybm, lev_p ,ps, intyp_m, P0mb, 1, kxtrp )
v = vinth2p (vv, hyam,hybm, lev_p ,ps, intyp_m, P0mb, 1, kxtrp )
q = vinth2p (qq, hyam,hybm, lev_p ,ps, intyp_m, P0mb, 1, kxtrp )
;if (m .eq. 0) then
;print(u)
;end if
u!1 = "lev"
u!2 = "lat"
u!3 = "lon"
u&lev = lev_p
u&lat = ref->U&lat
u&lon = ref->U&lon
copy_VarCoords(u,v)
copy_VarCoords(u,q)
ps!1 = "lat"
ps!2 = "lon"
ps&lat = ref->U&lat
ps&lon = ref->U&lon
copy_VarCoords(ps,PRECC)
copy_VarCoords(ps,PRECL)
copy_VarCoords(ps,EVAP)
;---Vertical levels
ptop = ptop*100.
ptop@units = "Pa"
plev = lev_p ; hPa
plev = plev*100 ; [0,...,100000] Pa [kg/(m-s2)]
plev@units = "Pa"
;---Change [kg/kg] to [g/kg]; not necessary: but common units for q
q = q*1000
q@units = "g/kg"
;---Layer thickness: ; Pa=>[kg/(m-s2)], (time,level,lat,lon)
;---Mass weighting: (dp/g) => [Pa/(m/s2)] => (Pa-s2)/m => [kg/(m-s2)][s2/m] => (kg/m2)
;---Reference: http://www.cgd.ucar.edu/cas/catalog/newbudgets/
dp = dpres_plevel_Wrap(plev, ps, ptop, 0) ; Pa; layar thickness
dpg = dp/g
dpg@long_name = "Layer Mass Weighting"
dpg@units = "kg/m2" ; dp/g, Pa/(m s-2), reduce to kg m-2
;---Moisture flux components at each pressure level
uq = u*q ; (:,:,:,:)
uq@long_name = "Zonal Moisture Flux [uq]"
uq@units = "["+u@units+"]["+q@units+"]" ; [m/s][g/kg]
copy_VarCoords(u,uq) ; (time,level,lat,lon)
vq = v*q ; (:,:,:,:)
vq@long_name = "Meridional Moisture Flux [vq]"
vq@units = "["+v@units+"]["+q@units+"]"
copy_VarCoords(v,vq) ; (time,level,lat,lon)
;---Integrated mass weighted moisture flux components
linlog = 1
uq_dpg = uq*dpg ; mass weighted 'uq'; [m/s][g/kg][kg/m2]=>[m/s][g/kg]
iuq = dim_sum_n(uq_dpg, 1)
;iuq = vibeta(plev(::-1),uq(time|:,lat|:,lon|:,lev|::-1),linlog,ps,pbot,ptop)/g
iuq@long_name = "Integrated Zonal UQ [uq*dpg]"
iuq@LONG_NAME = "Sum: Mass Weighted Integrated Zonal Moisture Flux [uq*dpg]"
iuq@units = "[m/s][g/kg]"
;copy_VarCoords(u, iuq); (time,lat,lon)
copy_VarCoords(ps, iuq); (time,lat,lon)
delete(uq_dpg)
vq_dpg = vq*dpg ; mass weighted 'vq'; [m/s][g/kg][kg/m2]=>[m/s][g/kg]
ivq = dim_sum_n(vq_dpg, 1)
;ivq = vibeta(plev(::-1),vq(time|:,lat|:,lon|:,lev|::-1),linlog,ps,pbot,ptop)/g ;;;;units are all Pa.
ivq@long_name = "Integrated Meridional VQ [vq*dpg]"
ivq@LONG_NAME = "Sum: Mass Weighted Integrated Meridional Moisture Flux [vq*dpg]"
ivq@units = "[m/s][g/kg]"
;copy_VarCoords(v, ivq); (time,lat,lon)
copy_VarCoords(ps, ivq)
delete(vq_dpg)
;---Divergence of moisture flux: uv2dvF => global 'fixed' rectilinear grid
duvq = uv2dvF_Wrap(uq, vq) ; (time,level,lat,lon)
duvq@long_name = "Divergence of Moisture Flux"
duvq@units = "g/(kg-s)" ; (1/m)*[(m/s)(g/kg)] => [g/(kg-s)]
;---Mass weighted integration [sum] of the divergence of moisture flux
duvq_dpg = duvq*dpg ; [g/(kg-s)][kg/m2] => [g/(m2-s)]
iduvq = dim_sum_n(duvq_dpg, 1)
;iduvq = dim_sum_n(duvq_dpg, 0)
iduvq@long_name = "Integrated Mass Wgt MFC"
iduvq@LONG_NAME = "Integrated Mass Weighted Moisture Flux Convergence"
iduvq@units = "g/(m2-s)"
;copy_VarCoords(u, iduvq) ; (time,lat,lon)
copy_VarCoords(ps, iduvq)
delete(duvq_dpg)
iduvq = 24.*3600.*iduvq/1000.
iduvq@units = "mm/day"
copy_VarCoords(u(:,0,:,:), iduvq) ; (time,lat,lon)
VIMFC = iduvq ; keep meta data
VIMFC = -VIMFC ; Note the preceding -1 [negative precedes integration]
VIMFC@long_name = "VIMFC"
copy_VarCoords(u(:,0,:,:), VIMFC)
;*************************************************
; Calculate divergence: Use Wrap to include meta data
; Calculate divergent wind components; used for graphics
;*************************************************
fname = "/iodati/mgr04_users/yhzhou/sppfinal/508-03-new-new-new/508-03.moisture.h1." + m + ".nc"
output = addfile(fname, "c")
output->VIMFC = VIMFC
output->PRECC = PRECC
output->PRECL = PRECL
output->EVAP = EVAP
; 删除创建的变量以释放内存
delete(PRECL)
delete(PRECC)
delete(EVAP)
delete(u)
delete(v)
delete(q)
end do
|
|
窥视卡
雷达卡
发表于 2024-4-19 15:05:49
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主
