之前一直使用的從師兄師姐那里繼承來的Fortran腳本計(jì)算考廉，并且需要自行將原始再分析資料轉(zhuǎn)化為二進(jìn)制文件惯殊，過于繁瑣
之前在NCL官網(wǎng)有看見過寫好的基于NCEP資料計(jì)算的腳本
NCL: Divergent and Rotational Winds; Mass Flux; Yanai Heat (ucar.edu)
也在某公眾號(hào)上看過python版本的，但是博主對(duì)比之后會(huì)有細(xì)微差異
https://mp.weixin.qq.com/s/6H6TOT_OxQOYku9cbHolXA
自己試了一下NCL的版本纳本，原腳本是寫給日分辨率資料使用的窍蓝，我用來計(jì)算月均值貌似也可以，目前看著運(yùn)行沒啥異常繁成。
計(jì)算主要分為兩步吓笙，先計(jì)算q1q2，為局地變化項(xiàng)巾腕，平流項(xiàng)面睛，垂直對(duì)流項(xiàng)之和，然后分別乘以Cp和L
但是由于本人天原和動(dòng)力氣象都沒學(xué)好尊搬，且后來一直做環(huán)境方向叁鉴，所以對(duì)這個(gè)腳本還存在很多疑惑的點(diǎn)，先記錄下來以后有機(jī)會(huì)一個(gè)一個(gè)解決毁嗦。
1.一般文獻(xiàn)里Q1Q2的單位都是w/m2亲茅，但是腳本計(jì)算的逐層的Q1Q2單位不一致，那要是畫剖面怎么辦呢狗准？
并且，似乎茵肃，原作者也不清楚該怎么對(duì)逐層數(shù)據(jù)的單位進(jìn)行轉(zhuǎn)換
2. 對(duì)q2積分完之后所有格點(diǎn)的值都變成了空值腔长，不知道是不是因?yàn)橛玫脑轮担坷碚撋喜粦?yīng)該啊验残，畢竟q1算的也沒問題（我也不知道是不是真的沒問題）
3. 為什么垂直積分的時(shí)候沒有用到地表氣壓捞附？
4.對(duì)Q1Q2整層積分之后，單位經(jīng)過自己的換算都變成了w/m2
先注釋掉開頭自定義函數(shù)中對(duì)q1q2單位向day的轉(zhuǎn)化，直接以原始的s進(jìn)行計(jì)算
并且在主程序中添加對(duì)整層積分完q1q2向Q1Q2的計(jì)算
最后一點(diǎn)鸟召，一般我們使用的都是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????\
??, rankp, ranku, rankv, rankH, rankT, ranko??\
??, Cp, Lc, dTdt, ss, opt_adv, long_name, units?\
??, gridType, advT, q1, Q1, dHdt, advH, loq, q2, Q2
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]
;*******************************************
?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("-----")
;****************************************
;---Compute static stability [K/Pa] <==>?or "K-m-s2/kg"
;****************************************
?ss?= static_stability (p , T, 1, 0)
?printVarSummary(ss)
?printMinMax(ss,0)
?print("-----")
;****************************************
;---Compute advection term: spherical harmonics
;---U*(dT/dx) + V*(dT/dy):?[m/s][K/m] -> [K/s]
;****************************************
?opt_adv= 0
?long_name = "temp advection"
?units?= "K/s"
?gridType?= 1
?advT??= advect_variable(u,v,T,gridType,long_name,units,opt_adv)
;****************************************
;---Net Heat
;****************************************
?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("-----")
;****************************************
;---Apparent Heat Source
;****************************************
?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("-----")
;*******************************************
;---Compute local dH/dt?
;*******************************************
?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("-----")
;****************************************
;---Compute advection term: global fixed grid: spherical harmonics
;---U*(dH/dlon) + V*(dH/dlat)
;****************************************
?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)
;****************************************
;---Compute vertical movement of H
;****************************************
?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("-----")
;****************************************
;---Apparent Moisture Sink
;****************************************
?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("-----")
;****************************************
;---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"
?return( [/q1, q2, Q1, Q2/] )
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"
?netCDF?= True?????????????; Write netCDF
;---Variable and file handling
?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,:)???????; %?
;---Need specific humidity
?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)
;---Convert "hours since ..." to "seconds since ..."
?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)
;---Convert hPa -> Pa for function
?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]
;********************************************
;---Vertical integration
; J[oule]??kg-m2/s2 = N-m = Pa/m3 = W-s??; energy????
;********************************************
?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("-----")
;***********************************************
;---Save to a netcdf file
;***********************************************
?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
; ;***********************************************
; ;---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)
