加载中…大气视热源Q1和视水汽汇Q2的计算
(2022-04-11 22:52:55) | 分类： NCL的小tips |
之前一直使用的从师兄师姐那里继承来的Fortran脚本计算，并且需要自行将原始再分析资料转化为二进制文件，过于繁琐
![]( "大气视热源Q1和视水汽汇Q2的计算")
并且在主程序中添加对整层积分完q1q2向Q1Q2的计算
![]( "大气视热源Q1和视水汽汇Q2的计算")
, rankp, ranku, rankv,
rankH, rankT, ranko
\
, Cp, Lc, dTdt, ss,
opt_adv, long_name, units
\
, gridType, advT, q1,
Q1, dHdt, advH, loq, q2, Q2
dTdt = center_finite_diff_n
(T,time,False,0,0) ; 'time' is
'seconds since'
copy_VarCoords(T, dTdt)
dTdt@longname = "Temperature: Local Time
derivative"
dTdt@units = "K/s"
printVarSummary(dTdt)
printMinMax(dTdt,0)
print("-----")
ss = static_stability (p
, T, 1, 0)
printVarSummary(ss)
printMinMax(ss,0)
print("-----")
opt_adv
= 0
long_name = "temp advection"
units
= "K/s"
gridType =
1
advT
=
advect_variable(u,v,T,gridType,long_name,units,opt_adv)
q1 = dTdt -
(omega*ss-advT)
; term_1 - term_2
q1@long_name = "q1: Net Heat Flux"
q1@units
= "K/s"
copy_VarCoords(T,q1)
printVarSummary(q1)
printMinMax(q1,0)
print("-----")
Cp
= 1004.64
Cp@long_name = "specific heat of dry air at
constant pressure"
Cp@units
= "J/(K-kg)"
; [kg-m2/s2]/(K-kg) => [kg-m2]/[K-kg-s2] =>
m2/[K-s2]
Q1 =
Cp*q1 ; [J/(K-kg)][K/s] -> [J/(kg-s)] ->
[(kg-m2/s2)/(kg-s)) -> [ m2/s3 ]
; [J/(K-kg)][K/s] ->
[(J/s)(1/kg)] -> W/kg
???
Q1@long_name = "Q1: Total Diabatic Heating
as the Apparent Heat Source"
Q1@units
= "m2/s3"
copy_VarCoords(T,Q1)
printVarSummary(Q1)
printMinMax(Q1,0)
print("-----")
dHdt = center_finite_diff_n
(H,time,False,0,0)
copy_VarCoords(H, dHdt)
dHdt@longname = "Specific Humidity: Local
Time derivative"
dHdt@units = "g/(kg-s)"
;
(g/kg)/s
printVarSummary(dHdt)
printMinMax(dHdt,0)
print("-----")
long_name = "moisture advection"
units
= "g/(kg-s)"
; (m/s)*(g/kg)*(1/m)
advH
=
advect_variable(u,v,H,gridType,long_name,units,opt_adv)
dHdp = center_finite_diff_n
(H,p,False,1,npr)
copy_VarCoords(H, dHdp)
dHdp@longname = "Specific Humidity: Local
Vertical Derivative"
dHdp@units = "g/(kg-Pa)"
; (g/kg)/Pa
printVarSummary(dHdp)
printMinMax(dHdp,0)
print("-----")
dHdp =
omega*dHdp
; overwrite ...
convenience
dHdp@longname = "Specific Humidity:
Vertical Moisture Transport"
dHdp@units
= "g/(kg-s)"
;
[(Pa/s)(g/kg)/Pa)]
printVarSummary(dHdp)
printMinMax(dHdp,0)
print("-----")
q2 =
-(dHdt +advH +dHdp)
q2@long_name = "q2: moisture
sink" ;
"?apparent? drying rate"
q2@units
= "g/(kg-s)"
copy_VarCoords(T,q2)
printVarSummary(q2)
Lc
= 2.26e6 ;
[J/kg]=[m2/s2] Latent Heat of
Condensation/Vaporization
Lc@long_name = "Latent Heat of
Condensation/Vaporization"
Lc@units
= "J/kg" ;
==> "m2/s2"
Q2 =
Lc*q2
;
(J/kg)(g/(kg-s))->(m2/s2)(g/(kg-s)) ->[(g/kg)][m2/s3]
; J[oule]-> kg-m2/s2 = N-m =
Pa/m3 = W-s
; energy
Q2@long_name = "Q2: Apparent Moisture
Sink"
Q2@units
=
"(g/kg)[m2/s3]"
copy_VarCoords(T,Q2)
printVarSummary(Q2)
printMinMax(Q2,0)
print("-----")
return( [/q1, q2, Q1, Q2/] )
netCDF =
True
; Write
netCDF
diri =
"/public/home/sunxiaoyun/datadir/ncep_monthly/pressure_new/"
f1
=
addfile(diri+"air.mon.mean.nc","r")
; daily mean data
f2
=
addfile(diri+"omega.mon.mean.nc","r")
f3
= addfile(diri+"uwnd.mon.mean.nc","r")
f4
= addfile(diri+"vwnd.mon.mean.nc","r")
f5
= addfile(diri+"rhum.mon.mean.nc","r")
; convenience: make all:
S->N
temp =
f1->air(:,0:7,::-1,:)
; degK
omega =
f2->omega(:,0:7,::-1,:)
;
Pascal/s
uwnd =
f3->uwnd(:,0:7,::-1,:)
; m/s
vwnd =
f4->vwnd(:,0:7,::-1,:)
; m/s
rhum =
f5->rhum(:,:,::-1,:)
; %
p
=
f1->level(0:7)
; hPa
[*]
p4d
= conform(temp, p, 1)
printVarSummary(p)
printVarSummary(p4d)
printVarSummary(rhum)
q
= mixhum_ptrh (p4d, temp,
rhum,-2)
;
g/kg
delete( [/p4d, rhum/] )
q@long_name = "specific humidity"
q@units = "g/kg"
copy_VarCoords(temp, q)
printVarSummary(q)
printMinMax(q, 0)
time =
f5->time
; "hours
since 1800-1-1"
time =
time*3600
time@units = "seconds
since 1800-1-1 00:00:0.0"
printVarSummary(time)
print("---")
t
= time
; Lyndz'
name
ymdh = cd_calendar(time,-3)
p
= p*100
;
Pa [100000,...,10000]
p@units = "Pa"
p!0
= "p"
p&p
=
p
; not necessary
printVarSummary(p)
print("---")
Cp
= 1004.64
; specific heat of dry air at
constant pressure
Cp@units =
"J/(K*kg)"
Lc
= 2.26e6 ;
[J/kg]=[m2/s2] Latent Heat of
Condensation/Vaporization
Lc@long_name = "Latent Heat of
Condensation/Vaporization"
Lc@units
= "J/kg" ;
==> "m2/s2"
npr = 1
opt = True
q1q2 =
Q1Q2_yanai(time,p,uwnd,vwnd,q,temp,omega,npr,opt)
print(q1q2)
q1 =
q1q2[0]
q2 =
q1q2[1]
Q1 =
q1q2[2]
Q2 =
q1q2[3]
ptop =
30000.0
; Pa
pbot = 100000.0
vopt =
1
g
= 9.80665
; [m/s2] gravity at 45 deg lat used by the
WMO
dp =
dpres_plevel_Wrap(p,pbot,ptop,0)
dpg =
dp/g
; Pa/(m/s2)=> (Pa-s2)/m
dpg@long_name = "Layer Mass
Weighting"
dpg@units
= "kg/m2"
; dp/g
=> Pa/(m/s2) => [kg/(m-s2)][m/s2] reduce to
(kg/m2)
;
Pa (s2/m)
=> [kg/(m-s2)][s2/m]=>[kg/m2]
dpg :=
conform(q1,dpg,1)
; reassign [convenience]
q1int = wgt_vertical_n(q1,dpg,vopt,1) ;
SUM[q1*dpg] => (K/s)(kg/m2)
q1int@long_name = "Heat Source: vertically
integrated"
q1int@units
=
"(K/s)(kg/m2)"
; (K/s)(kg/m2) => (kg-K)/(m2-s)
printVarSummary(q1int)
copy_VarCoords(temp(:,0,:,:),q1int)
printMinMax(q1int,0)
print("-----")
q2int = wgt_vertical_n(q2,dpg,vopt,1)
q2int@long_name = "Moisture Sink:
vertically integrated"
q2int@units
= "g/(s-m2)"
copy_VarCoords(temp(:,0,:,:),q2int)
printMinMax(q2int,0)
print("-----")
Q1int
= Cp*q1int
; (K/s)(kg/m2)*(J/(K*kg)) => J/(s*m2) =>
w/m2
Q1int@long_name = "Q1: Total Diabatic
Heating as the Apparent Heat Source"
Q1int@units
=
"w/m2"
copy_VarCoords(temp(:,0,:,:),Q1int)
printMinMax(Q1int,0)
print("-----")
Q2int
= Lc*q2int*1000
;
g/(s-m2)*(J/kg) => J/(1000*s*m2) => w/m2/1000
Q2int@long_name = "Q2: Apparent Moisture
Sink"
Q2int@units
= "w/m2"
copy_VarCoords(temp(:,0,:,:),Q2int)
printMinMax(Q2int,0)
print("-----")
if (netCDF)
diro
= "./"
filo
= "YANAI.apparent_heat_ncep_mon.nc"
ptho
= diro+filo
system("/bin/rm -f "+ptho)
ncdf
= addfile(ptho,"c")
fAtt
= True
fAtt@title
= "Apparent Heat Source based
on Yanai et al. 1973"
fAtt@source_name = "NCEP-NCAR
Reanalysis 2"
fAtt@source_URL =
"https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis2.html"
fAtt@source
= "NOAA/OAR/ESRL PSD, Boulder, Colorado,
USA"
fAtt@Conventions = "None"
fAtt@creation_date = systemfunc("date")
fileattdef(ncdf,fAtt)
; copy file attributes
filedimdef(ncdf,"time",-1,True) ; make time an UNLIMITED
dimension
ncdf->Q1 = Q1
ncdf->Q2 = Q2
ncdf->Q1int = Q1int
ncdf->Q2int = Q2int
end if
之前在NCL官网有看见过写好的基于NCEP资料计算的脚本
也在某公众号上看过python版本的，但是博主对比之后会有细微差异
https://mp.weixin.qq.com/s/6H6TOT_OxQOYku9cbHolXA
自己试了一下NCL的版本，原脚本是写给日分辨率资料使用的，我用来计算月均值貌似也可以，目前看着运行没啥异常。
计算主要分为两步，先计算q1q2，为距地变化项，平流项，垂直对流项之和，然后分别乘以Cp和L
但是由于本人天原和动力气象都没学好，且后来一直做环境方向，所以对这个脚本还存在很多疑惑的点，先记录下来以后有机会一个一个解决。
1.一般文献里Q1Q2的单位都是w/m2，但是脚本计算的逐层的Q1Q2单位不一致，那要是画剖面怎么办呢？
并且，似乎，原作者也不清楚该怎么对逐层数据的单位进行转换
2.
对q2积分完之后所有格点的值都变成了空值，不知道是不是因为用的月值？理论上不应该啊，毕竟q1算的也没问题（我也不知道是不是真的没问题）
3. 为什么垂直积分的时候没有用到地表气压？
4.对Q1Q2整层积分之后，单位经过自己的换算都变成了w/m2
先注释掉开头自定义函数中对q1q2单位向day的转化，直接以原始的s进行计算
并且在主程序中添加对整层积分完q1q2向Q1Q2的计算
最后一点，一般我们使用的都是Q1Q2，不知道原作者为什么输出了q1q2，顺带把输出变量也给改了。
Code:
undef("Q1Q2_yanai.ncl")
function
Q1Q2_yanai(time[*]:numeric,p,u,v,H,T,omega,npr[1]:integer,opt[1]:logical)
;+++++++++++++++++++++++++++++++++++++++++++++++++++++
; NOT SUPPORTED: NOT FULLY TESTED : WORK in
PROGRESS
;
**CHECK
UNITS**
;+++++++++++++++++++++++++++++++++++++++++++++++++++++
; Nomenclature
; time - "seconds since
..."
; p
- Pressure [Pa]
; u,v -
zonal, meridional wind components[m/s]
; H
- specific humidity [g/kg]
; T
- temperature [K] or [C]
; omega - vertical
velocity [Pa/s]
; npr -
demension number corresponding to pressure dimension
; opt -
set to False
;+++++++++++++++++++++++++++++++++++++++++++++++++++++
;; Q1 = Cp*dTdt -
Cp*(omega*ss-advT)
;; Q1 = Cp*[ dTdt-
omega*ss-advT ]
;; q1 = dTdt-
omega*ss-advT
;; Q1 = Cp*[ q1
]
;;
;; Q1 = Total diabatic
heating [including radiation, latent heating, and
;;
sfc heat flux) and sub-grid scale heat flux
convergence
;;---
;; q2 = -(dHdt +advH
+dHdp)
;; Q2 = Lc*[ q2 ]
;;
;; Q2 = the latent heating
due to condensation or evaporation processes
;;
and subgrid-scale moisture flux
convergences,
;++++++++++++++++++++++++++++++++++++++++++++
; References:
;
;
https://renqlsysu.github.io/2019/02/01/apparent_heat_source/
;
; Yanai, M. (1961):
; A detailed analysis of typhoon formation.
; J. Meteor. Soc. Japan , 39 , 187–214
;
; Yanai et al (1973):
; Determination of bulk properties of
tropical cloud clusters
; from large-scale heat
and moisture budgets.
; J.
Atmos. Sci. , 30 ,
611–627,
;
https://doi.org/10.1175/1520-0469(1973)030<0611:DOBPOT>2.0.CO;2
;
; Yanai, M and T.Tomita (1998):
;
https://pdfs.semanticscholar.org/fb57/a6a59cc4a684194b5e622ea6f875d0b4439a.pdf
;********************************************
; Diabatic heating in the atmosphere is a combined consequence
of
; radiative fluxes, phase changes of water substance, and
turbulent
; flux of sensible heat from the earth's
surface.
;
; In the tropics, it is the major driving force of the
atmospheric circulation.
; It responds to the vertical gradient of diabatic
heating.
;********************************************
local dimp, dimu, dimv, dimH, dimT, dimo
\
begin
;;
Use for error testing
;;dimp
= dimsizes(p)
;;dimu
= dimsizes(u)
;;dimv
= dimsizes(v)
;;dimH
= dimsizes(H)
;;dimT
= dimsizes(T)
;;dimo
= dimsizes(omega)
;;rankp
= dimsizes(dimp)
;;ranku
= dimsizes(dimu)
;;rankv
= dimsizes(dimv)
;;rankH
= dimsizes(dimH)
;;rankT
= dimsizes(dimT)
;;ranko
= dimsizes(dimo)
;*******************************************
;---Compute local dT/dt [K/s]
;*******************************************
;****************************************
;---Compute static stability [K/Pa]
<==> or "K-m-s2/kg"
;****************************************
;****************************************
;---Compute advection term: spherical harmonics
;---U*(dT/dx) + V*(dT/dy): [m/s][K/m] ->
[K/s]
;****************************************
;****************************************
;---Net Heat
;****************************************
;****************************************
;---Apparent Heat Source
;****************************************
;*******************************************
;---Compute local dH/dt
;*******************************************
;****************************************
;---Compute advection term: global fixed grid: spherical
harmonics
;---U*(dH/dlon) + V*(dH/dlat)
;****************************************
;****************************************
;---Compute vertical movement of H
;****************************************
;****************************************
;---Apparent Moisture Sink
;****************************************
;****************************************
;---Make q1, q2 to per day
;****************************************
; q1 = q1*86400
; q2 = q2*86400
; q1@units
= "K/day"
; q2@units
= "g/(kg-day)"
;****************************************
;---Make Q1, Q2 to per W/m2
;****************************************
;;Q1 = Q1*?????
;;Q2 = Q2*?????
;;Q1@units = "W/m2"
;;Q2@units = "W/m2"
end
;==============================================================================
;
MAIN
;==============================================================================
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/csm/contributed.ncl"
;---Variable and file handling
;---Need specific humidity
;---Convert "hours since ..." to "seconds since ..."
;---Convert hPa -> Pa for function
;++++++++++++++++++++++++++++++++++++++++++++++++++++++
;********************************************
;---Vertical integration
; J[oule]
kg-m2/s2 = N-m = Pa/m3 = W-s
;
energy
;********************************************
;***********************************************
;---Save to a netcdf file
;***********************************************
; ;***********************************************
; ;---Plot
; ;***********************************************
; nt
= 4
; YMDH =
ymdh(nt)
; LEVP =
600
; opt = True
; opt@PrintStat =
True
; stat_q1 =
stat_dispersion(q1(nt,{LEVP},:,:), opt )
; stat_q2 =
stat_dispersion(q2(nt,{LEVP},:,:), opt )
; stat_q1i=
stat_dispersion(q1int(nt,:,:), opt )
; stat_q2i=
stat_dispersion(q2int(nt,:,:), opt )
; plot =
new(2,graphic)
; wks
=
gsn_open_wks("png","q2q1_yanai")
; send
graphics to PNG file
; ;--- mfc_adv and mfc_con at a specified pressure level
; res
=
True
; plot mods desired
;
res@gsnDraw
= False
; don't draw yet
;
res@gsnFrame
=
False
; don't
advance frame yet
;
res@cnFillOn
=
True
; turn on color
;
res@cnLinesOn
= False
; turn off contour lines
;
res@cnLineLabelsOn
= False
; turn off contour lines
; ;res@cnFillPalette
= "ViBlGrWhYeOrRe" ; set
White-in-Middle color map
;
res@cnFillPalette
= "amwg256"
; set
White-in-Middle color map
; ;res@cnFillMode
=
"RasterFill"
;
res@mpFillOn
=
False
; turn off
map fill
; ;;res@lbLabelBarOn
= False
; turn off individual cb's
;
res@lbOrientation
= "Vertical"
;
; Use a common scale
; res@cnLevelSelectionMode
= "ManualLevels"; manual set levels so lb consistent
;
res@cnMaxLevelValF
=
5.0
; max level
;
res@cnMinLevelValF
=
-res@cnMaxLevelValF
; min level
;
res@cnLevelSpacingF
=
0.20
; contour interval
;
res@gsnCenterString
= LEVP+"hPa"
; plot(0) =
gsn_csm_contour_map(wks,q1(nt,{LEVP},:,:),res)
;
res@cnMaxLevelValF
=
5.0
; max level
;
res@cnMinLevelValF
=
-res@cnMaxLevelValF
; min level
;
res@cnLevelSpacingF
=
0.20
; contour interval
; plot(1) =
gsn_csm_contour_map(wks,q2(nt,{LEVP},:,:),res)
; resP
= True
; modify
the panel plot
;
resP@gsnMaximize
= True
; resP@gsnPanelMainString
:= YMDH
; ;;resP@gsnPanelLabelBar
= True
; add common colorbar
;
gsn_panel(wks,plot,(/2,1/),resP)
; now draw as one plot
;
delete(res@gsnCenterString) ; not used for this
plot
;
res@cnMaxLevelValF
=
14000.0
; max level
;
res@cnMinLevelValF
=
-res@cnMaxLevelValF
; min level
;
res@cnLevelSpacingF
=
500.0
; contour interval
; plot(0) =
gsn_csm_contour_map(wks,q1int(nt,:,:),res)
;
res@cnMaxLevelValF
=
7000.0 ; max
level
;
res@cnMinLevelValF
=
-res@cnMaxLevelValF
; min level
;
res@cnLevelSpacingF
=
500.0
; contour interval
; plot(1) =
gsn_csm_contour_map(wks,q2int(nt,:,:),res)
;
gsn_panel(wks,plot,(/2,1/),resP)
; now draw as one plot
; ;---Cross section
; rescx
= True
; plot mods desired
;
rescx@gsnMaximize
= True
; LAT =
10
; if (LAT.ge.0) then
;
rescx@tiMainString
= YMDH+": "+LAT+"N"
; else
;
rescx@tiMainString
= YMDH+": "+ABS(LAT)+"S"
; end if
;
rescx@cnLevelSelectionMode =
"ManualLevels"
; manual contour levels
;
rescx@cnLinesOn
= False
; turn off contour
lines
;
rescx@cnLineLabelsOn
=
False
;
rescx@cnFillOn
=
True
; turn on
color fill
;
rescx@cnFillPalette
=
"ViBlGrWhYeOrRe" ; set White-in-Middle color map
;
rescx@cnFillPalette
=
"amwg256"
; set White-in-Middle color map
;
rescx@cnMaxLevelValF
= 5.0
; max level
;
rescx@cnMinLevelValF
=
-rescx@cnMaxLevelValF ; min level
;
rescx@cnLevelSpacingF
=
0.25
; contour
interval
; plt_q1 =
gsn_csm_pres_hgt(wks,q1(nt,{1000:300},{LAT},:),rescx)
;
rescx@cnMaxLevelValF
= 5.0
; max level
;
rescx@cnMinLevelValF
=
-rescx@cnMaxLevelValF ; min level
;
rescx@cnLevelSpacingF
=
0.25
; contour
interval
; plt_q2 =
gsn_csm_pres_hgt(wks,q2(nt,{1000:300},{LAT},:),rescx)
喜欢
0
赠金笔