气象编程 | 大气视热源Q1和视水汽汇Q2程序
大气视热源Q1和视水汽汇Q2的程序提供NCL、FORTRAN、IDL三种选择，推荐使用NCL。
Q1Q2_yanai.ncl: Using high-frequency data, calculate apparent-heat-source (Q1) and apparent-moisture-sink (Q2) quantities as described in:
This example uses daily mean data from NOAA/OAR/ESRL PSD, Boulder, Colorado, USA. Specifically, NCEP_Reanalysis 2 daily mean data spanning 3-9 April 1993.
This is a test script. It is not fully tested. The returned quantities are:
q1 = (dT/dt) - ; K/day
q2 = -[(dH/dt) +V.grad(H) +omega*(dH/dp)] ; g/(kg-Pa)
and
Q1 = Cp*q1 ; apparent Heat Source
Q2 = Lc*q2 ; apparent Moisture Sink
具体NCL代码如下：
undef("Q1Q2_yanai.ncl")
function Q1Q2_yanai(time<li>:numeric,p,u,v,H,T,omega,npr:integer,opt:logical)
;+++++++++++++++++++++++++++++++++++++++++++++++++++++
; NOT SUPPORTED:NOT FULLY TESTED : WORK in PROGRESS
; **CHECK UNITS**
;+++++++++++++++++++++++++++++++++++++++++++++++++++++
; Nomenclature
; time - "seconds since ..."
; p - Pressure
; u,v - zonal, meridional wind components
; H - specific humidity
; T - temperature or
; omega - vertical velocity
; 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:
;
; <a href="https://renqlsysu.github.io/2019/02/01/apparent_heat_source/" target="_blank">https://renqlsysu.github.io/2019/02/01/apparent_heat_source/</a>
;
; 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,
;<a href="https://doi.org/10.1175/1520-0469" target="_blank">https://doi.org/10.1175/1520-0469</a>(1973)030<0611:DOBPOT>2.0.CO;2
;
; Yanai, M and T.Tomita (1998):
; <a href="https://pdfs.semanticscholar.org/fb57/a6a59cc4a684194b5e622ea6f875d0b4439a.pdf" target="_blank">https://pdfs.semanticscholar.org/fb57/a6a59cc4a684194b5e622ea6f875d0b4439a.pdf</a>
;********************************************
; 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
;*******************************************
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 <==>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): ->
;****************************************
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)"; /(K-kg) => / => m2/
Q1= Cp*q1; -> -> [(kg-m2/s2)/(kg-s)) -> [ m2/s3 ]
; -> [(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 ; =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)]
; J->kg-m2/s2 = N-m = Pa/m3 = W-s ; energy
Q2@long_name = "Q2: Apparent Moisture Sink"
Q2@units = "(g/kg)"
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( )
end
;==============================================================================
; MAIN
;==============================================================================
netCDF= True ; Write netCDF
;---Variable and file handling
diri = "/Users/shea/Data/netCDF/"
f1 = addfile(diri+"temp_lag4tolead2.nc","r") ; daily mean data
f2 = addfile(diri+"omega_lag4tolead2.nc","r")
f3 = addfile(diri+"uwnd_lag4tolead2.nc","r")
f4 = addfile(diri+"vwnd_lag4tolead2.nc","r")
f5 = addfile(diri+"rhum_lag4tolead2.nc","r")
; convenience: make all:S->N
temp = f1->air(:,:,::-1,:) ; degK
omega = f2->omega(:,:,::-1,:) ; Pascal/s
uwnd = f3->uwnd(:,:,::-1,:) ; m/s
vwnd = f4->vwnd(:,:,::-1,:) ; m/s
rhum = short2flt(f5->rhum(:,:,::-1,:)) ; %
;---Need specific humidity
p = f1->level ; hPa</li><li>
p4d = conform(temp, p, 1)
q = mixhum_ptrh (p4d, temp, rhum,-2) ; g/kg
delete( )
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)
print(ymdh)
;---Convert hPa -> Pa for function
p = p*100 ; Pa
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)"
npr= 1
opt= True
q1q2 = Q1Q2_yanai(time,p,uwnd,vwnd,q,temp,omega,npr,opt)
print(q1q2)
q1 = q1q2
q2 = q1q2
printVarSummary(q1)
printMinMax(q1,0)
print("================")
printVarSummary(q2)
printMinMax(q2,0)
print("================")
;Q1 = q1q2
;Q2 = q1q2
;********************************************
;---Vertical integration
; J kg-m2/s2 = N-m = Pa/m3 = W-s ; energy
;********************************************
ptop= 30000.0 ; Pa
pbot= 100000.0
vopt= 1
g = 9.80665 ; 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) => reduce to (kg/m2)
; Pa (s2/m) => =>
dpg:= conform(q1,dpg,1) ; reassign
q1int = wgt_vertical_n(q1,dpg,vopt,1) ; SUM=> (K/s)(kg/m2)
q1int@long_name = "Heat Source: vertically integrated"
q1int@units = "(K/day)(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/(day-m2)"
copy_VarCoords(temp(:,0,:,:),q2int)
printMinMax(q2int,0)
print("-----")
;***********************************************
;---Save to a netcdf file
;***********************************************
if (netCDF)
diro = "./"
filo = "YANAI.apparent_heat.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
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:200},{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:200},{LAT},:),rescx)</li>
具体Fortran代码
C***************************************************************
C Compute voticity divegence YANG 1992.10.25
C input UUVVOUTPUT VORT DIV (APPP) *
C***************************************************************
program calsiq
parameter(I0=71,J0=71,K0=21,kk=19,k4=12,df1=1.,df2=1.,f0=-10.)
DIMENSION ww(i0,j0,kk,k4),w(i0,j0,kk),u(i0,j0,kk),v(i0,j0,kk)
DIMENSION uu(i0,j0,kk,k4),TT(i0,j0,kk,k4),HH(i0,j0,kk,k4),P(K0)
& ,t(i0,j0,kk),q(i0,j0,kk),QQ(i0,j0,kk,k4),vv(i0,j0,kk,k4)
DIMENSION PP(KK),siq1(i0,j0,kk,k4),siq2(i0,j0,kk,k4),h(i0,j0,kk)
DIMENSION K88(kk),work(i0,j0,k0,k4),OROG(i0,j0),hs(i0,j0,k4)
DIMENSION rq(i0,j0),tg(i0,j0,kk),th(i0,j0,kk),ft1(i0,j0,kk)
DIMENSION ft(i0,j0,kk),Q1(i0,j0),Q2(i0,j0),ft3(i0,j0,kk)
c COMMON /CONTAA/RG,RD,CP,G,OMG,DLON,DLAT,FO,DF
c & ,DLATD,DLOND,CONV,COV,DPS12
CHARACTER*40,FILEN,FILEN1
DATA P/1000.,975.,950.,925.,900.,850.,800.,750.,700.,650.,600.,
&550.,500.,450.,400.,350.,300.,250.,200.,150.,100/
DATA PP/1000.,950.,900.,850.,800.,750.,700.,650.,600.,
&550.,500.,450.,400.,350.,300.,250.,200.,150.,100/
c DATA K66/1,2,3,4,5,9/
DATA K88/1,3,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21/
C******************************************************************
C Compute VORV DIAV liuhz 1998.10.25 *
C******************************************************************
open(11,file='e:\data\ncep0629\data\orog1.grb'
& ,form='unformatted',access='direct',recl=i0*j0)
read(11,rec=1)
&((orog(I,J),I=1,i0),J=1,j0)
close(11)
open(11,file='e:\data\ncep0629\data\uu0629.grb'
& ,form='unformatted',access='direct',recl=i0*j0*k0*k4)
read(11,rec=1)
&((((work(I,J,k,K2),I=1,i0),J=1,j0),k=1,k0),k2=1,k4)
close(11)
do kt=1,k4
do k=1,kk
k2=k88(k)
do j=1,j0
do i=1,i0
uu(i,j,k,kt)=work(i,j,k2,kt)
enddo
enddo
enddo
enddo
open(11,file='e:\data\ncep0629\data\vv0629.grb'
& ,form='unformatted',access='direct',recl=i0*j0*k0*k4)
read(11,rec=1)
&((((work(I,J,k,K2),I=1,i0),J=1,j0),k=1,k0),k2=1,k4)
close(11)
do kt=1,k4
do k=1,kk
k2=k88(k)
do j=1,j0
do i=1,i0
vv(i,j,k,kt)=work(i,j,k2,kt)
enddo
enddo
enddo
enddo
open(11,file='e:\data\ncep0629\data\tt0629.grb'
& ,form='unformatted',access='direct',recl=i0*j0*k0*k4)
read(11,rec=1)
&((((work(I,J,k,K2),I=1,i0),J=1,j0),k=1,k0),k2=1,k4)
close(11)
do kt=1,k4
do k=1,kk
k2=k88(k)
do j=1,j0
do i=1,i0
tt(i,j,k,kt)=work(i,j,k2,kt)
enddo
enddo
enddo
enddo
open(11,file='e:\data\ncep0629\data\hh0629.grb'
& ,form='unformatted',access='direct',recl=i0*j0*k0*k4)
read(11,rec=1)
&((((work(I,J,k,K2),I=1,i0),J=1,j0),k=1,k0),k2=1,k4)
close(11)
do kt=1,k4
do k=1,kk
k2=k88(k)
do j=1,j0
do i=1,i0
hh(i,j,k,kt)=work(i,j,k2,kt)
enddo
enddo
enddo
enddo
open(11,file='e:\data\ncep0629\data\qq0629.grb'
& ,form='unformatted',access='direct',recl=i0*j0*k0*k4)
read(11,rec=1)
&((((work(I,J,k,K2),I=1,i0),J=1,j0),k=1,k0),k2=1,k4)
close(11)
do kt=1,k4
do k=1,kk
k2=k88(k)
do j=1,j0
do i=1,i0
qq(i,j,k,kt)=work(i,j,k2,kt)
enddo
enddo
enddo
enddo
irec=0
jrec=0
do 100 k=1,k4
do j=1,j0
do i=1,i0
do l=1,kk
h(i,j,l)=hh(i,j,l,k)
u(i,j,l)=uu(i,j,l,k)
v(i,j,l)=vv(i,j,l,k)
t(i,j,l)=tt(i,j,l,k)
q(i,j,l)=qq(i,j,l,k)
enddo
enddo
enddo
OPEN(20,FILE='e:\data\ncep0629\data\q1t.GRB',FORM='UNFORMATTED',
&ACCESS='DIRECT',RECL=i0*j0*(k0-1))
OPEN(21,FILE='e:\data\ncep0629\data\q1v.GRB',FORM='UNFORMATTED',
&ACCESS='DIRECT',RECL=i0*j0*(k0-1))
OPEN(22,FILE='e:\data\ncep0629\data\q1w.GRB',FORM='UNFORMATTED',
&ACCESS='DIRECT',RECL=i0*j0*(k0-1))
OPEN(23,FILE='e:\data\ncep0629\data\w20.GRB',FORM='UNFORMATTED',
&ACCESS='DIRECT',RECL=i0*j0*(k0-1))
OPEN(24,FILE='e:\data\ncep0629\data\q2t.GRB',FORM='UNFORMATTED',
&ACCESS='DIRECT',RECL=i0*j0*(k0-1))
OPEN(25,FILE='e:\data\ncep0629\data\q2v.GRB',FORM='UNFORMATTED',
&ACCESS='DIRECT',RECL=i0*j0*(k0-1))
OPEN(26,FILE='e:\data\ncep0629\data\q2w.GRB',FORM='UNFORMATTED',
&ACCESS='DIRECT',RECL=i0*j0*(k0-1))
OPEN(27,FILE='e:\data\ncep0629\data\q12.GRB',FORM='UNFORMATTED',
&ACCESS='DIRECT',RECL=i0*j0)
DO L=1,kk
do i=1,i0
do j=1,j0
rq(i,j)=u(i,j,l)
enddo
enddo
call smth9(i0,j0,rq,rq)
do i=1,i0
do j=1,j0
u(i,j,l)=rq(i,j)
enddo
enddo
do i=1,i0
do j=1,j0
rq(i,j)=v(i,j,l)
enddo
enddo
call smth9(i0,j0,rq,rq)
do i=1,i0
do j=1,j0
v(i,j,l)=rq(i,j)
enddo
enddo
do i=1,i0
do j=1,j0
rq(i,j)=t(i,j,l)
enddo
enddo
call smth9(i0,j0,rq,rq)
do i=1,i0
do j=1,j0
t(i,j,l)=rq(i,j)
enddo
enddo
do i=1,i0
do j=1,j0
rq(i,j)=q(i,j,l)
enddo
enddo
call smth9(i0,j0,rq,rq)
do i=1,i0
do j=1,j0
q(i,j,l)=rq(i,j)
enddo
enddo
ENDDO
IF(K.GT.1)THEN
DO L=1,kk
do i=1,i0
do j=1,j0
rq(i,j)=tt(i,j,l,k-1)
enddo
enddo
call smth9(i0,j0,rq,rq)
do i=1,i0
do j=1,j0
tg(i,j,l)=rq(i,j)
enddo
enddo
do i=1,i0
do j=1,j0
rq(i,j)=qq(i,j,l,k-1)
enddo
enddo
call smth9(i0,j0,rq,rq)
do i=1,i0
do j=1,j0
ft1(i,j,l)=rq(i,j)
enddo
enddo
ENDDO
ENDIF
IF(K.LT.K4)THEN
DO L=1,kk
do i=1,i0
do j=1,j0
rq(i,j)=tt(i,j,l,K+1)
enddo
enddo
call smth9(i0,j0,rq,rq)
do i=1,i0
do j=1,j0
th(i,j,l)=rq(i,j)
enddo
enddo
do i=1,i0
do j=1,j0
rq(i,j)=qq(i,j,l,K+1)
enddo
enddo
call smth9(i0,j0,rq,rq)
do i=1,i0
do j=1,j0
ft3(i,j,l)=rq(i,j)
enddo
enddo
ENDDO
ENDIF
DT=24*3600.0
C �����¶ȵľֵر仯
IF(K.EQ.1)THEN
DO I=1,i0
DO J=1,j0
DO L=1,kk
TG(I,J,L)=(TH(I,J,L)-T(I,J,L))/DT
fT(I,J,L)=(ft3(I,J,L)-q(I,J,L))/DT
ENDDO
ENDDO
ENDDO
ELSE IF(K.EQ.K4)THEN
DO I=1,i0
DO J=1,j0
DO L=1,kk
TG(I,J,L)=(T(I,J,L)-TG(I,J,L))/DT
fT(I,J,L)=(q(I,J,L)-fT1(I,J,L))/DT
ENDDO
ENDDO
ENDDO
ELSE
DO I=1,i0
DO J=1,j0
DO L=1,kk
TG(I,J,L)=(TH(I,J,L)-TG(I,J,L))/2.0/DT
fT(I,J,L)=(fT3(I,J,L)-fT1(I,J,L))/2.0/DT
ENDDO
ENDDO
ENDDO
ENDIF
call CALsiq12(U,V,t,q,h,tg,ft,OROG,irec,jrec
& ,I0,j0,kk,df1,df2,f0,pp)
do l=1,kk
do j=1,j0
do i=1,i0
c ampv(i,j,l,k)=ampv1(i,j,l,k)+ampv2(i,j,l,k)
enddo
enddo
enddo
c print*,ampv(12,18,11,kt)
WRITE(6,*)'TIME=',K
c pause 1111
100 continue
CLOSE(20);CLOSE(21);CLOSE(22);CLOSE(23);close(24);close(25)
close(26);close(27)
cc write(filen1,'(a)') 'e:\data\ncep0629\data\siq1.grb'
ccc open(17,file=filen1,
cc & form='unformatted',access='direct',recl=i0*j0*kk*k4)
cc write(17,rec=1) ((((siq1(I,J,K,kt),I=1,i0),J=1,j0)
cc & ,K=1,kk),kt=1,k4)
cc CLOSE(17)
cc write(filen1,'(a)') 'e:\data\ncep0629\data\siq2.grb'
cc open(17,file=filen1,
cc & form='unformatted',access='direct',recl=i0*j0*kk*k4)
cc write(17,rec=1) ((((siq2(I,J,K,kt),I=1,i0),J=1,j0)
cc & ,K=1,kk),kt=1,k4)
cc CLOSE(17)
STOP
END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SUBROUTINE CALsiq12(U,V,T,Q,h,tg,ft,OROG,irec,jrec
& ,I0,j0,k0,df1,df2,f0,pp)
C**************************************************************************
C Compute HEATIGH RATE (SOURCE/SINK) & WATER SOURCE/SINK by liuhz 2009.3 *
C inputU=====> U WIND (m/s) *
C V=====> V WIND (m/s) *
C H ======> GEOPOTENTIAL HEIHGT (GEOPOTENTIAL METER)
C Q ====> SPECIFIC HUMIDITY (G/G) *
C T =====> TEMPERATURE (K) *
C tg =====> local variation of TEMPERATURE *
C ft =====> local variation of SPECIFIC HUMIDITY
C OROG======> OROGRAPHIC HEIHGT(GEOPOTENTIAL METER) *
C output UG(Q1V)=====> ADVECTIVE TRANSFER of HEAT (K/S)
C WG(Q1W)=====> VERTICAL TRANSFEROF HEAT (K/S)
C W(W20) =====> W VELOCITY (HPA/S) ; *
C FV(Q2V)=====> ADVECTIVE TRANSFER of water (G/G/S)
C FW(Q2W)=====> VERTICAL TRANSFEROF water (G/G/S)
C Q1(Q12)=====> HEAT SOURCE/SINK (K/S) ;
C Q2 =====> WATER SOURCE/SINK (G/G/S) ;
C tg (Q1T) =====> local variation of TEMPERATURE (K/S) *
C ft (Q2T) =====> local variation of SPECIFIC HUMIDITY (G/G/S)
C**************************************************************************
C INCLUDE 'PARAMET213'
C PARAMETER (I0=161,J0=113,K0=13,DF1=1.0,DF2=1.0,F0=-10.)
C I0,J0,K0---GRID NUMBER IN X Y Z; DF1-DF2- FIELD RESOLUTION IN X Y;
C F0--FIRST LATITUDE FROM SOUTH TO NORTH
C P(1) --P(K0) EACH LAYEL PRESSURE FROM LOW TO HIHG
C DATA PL/1000.0,925.0,850.0,700.0,600.,500.0,400.0,300.0,250.0,200.
C & ,150.,100.,50./
C ���ò�ֵ��NCEP/NCAR 20��U��V��T��H���ϼ���Q1 & Q2
C ����X��Y��P��T��������
dimension OROG(i0,j0),H(i0,j0,k0),PP(k0),a(i0,j0,k0)
& ,siq1(I0,J0,K0),siq2(I0,J0,K0),AR(I0,J0,K0),P(k0)
DIMENSION U(I0,J0,K0),V(I0,J0,K0),PL(K0),T(I0,J0,K0)
DIMENSION PS(i0,j0),US(i0,j0),VS(i0,j0),WS(i0,j0),ts(i0,j0)
DIMENSION ILAY(i0,j0),q(i0,j0,k0),gtmp(i0,j0,k0)
DIMENSION DIVG(i0,j0,k0),DIVS(i0,j0),W(i0,j0,k0)
dimension wg(i0,j0,k0),ug(i0,j0,k0),tg(i0,j0,k0),TH(i0,j0,k0)
DIMENSION rQ(i0,j0),HT(i0,j0,k0),wt(i0,j0),Q1(i0,j0),Q2(i0,j0)
dimension ft(i0,j0,k0),fv(i0,j0,k0),fw(i0,j0,k0)
dimension ft1(i0,j0,k0),ft3(i0,j0,k0),ff(i0,j0,k0)
COV=3.1415926/180.
c DF=DF2
DLON=DF1*COV
DLAT=DF2*COV
DLOND=DF1
DLATD=DF2
RG=6371229.0
DY=2.0*RG*DLAT
ck=0.286
gama=0.0065
rd=287.04
g=9.80616
dk=g/(rd*gama)
DP=-50.0
WTt=0.0
c DT=6*3600.0
DT=24*3600.0
CP=1005.0
LH=2501000.0
C �����¶�ƽ��---������
do i=2,i0-1
do j=2,j0-1
F1=(F0+FLOAT(J-1)*DLATD)*COV
DX=2.0*RG*COS(F1)*DLON
DO L=1,k0
UG(I,J,L)=U(I,J,L)*(T(I+1,J,L)-T(I-1,J,L))/DX
c & COS((-87.0+1.0*(j-1))*RAD)+
& +V(I,J,L)*(T(I,J+1,L)-T(I,J-1,L))/DY
fv(I,J,L)=U(I,J,L)*(q(I+1,J,L)-q(I-1,J,L))/DX
c & COS((-87.0+1.0*(j-1))*RAD)+
& +V(I,J,L)*(q(I,J+1,L)-q(I,J-1,L))/DY
ENDDO
ENDDO
ENDDO
do j=2,j0-1
F1=(F0+FLOAT(J-1)*DLATD)*COV
DX=2.0*RG*COS(F1)*DLON
DO L=1,k0
UG(1,J,L)=U(1,J,L)*(T(2,J,L)-T(i0-1,J,L))/DX
c & COS((-87.0+1.0*(j-1))*RAD)+
& +V(1,J,L)*(T(1,J+1,L)-T(1,J-1,L))/DY
fv(1,J,L)=U(1,J,L)*(q(2,J,L)-q(i0-1,J,L))/DX
c & COS((-87.0+1.0*(j-1))*RAD)+
& +V(1,J,L)*(q(1,J+1,L)-q(1,J-1,L))/DY
ENDDO
DO L=1,k0
UG(i0,J,L)=U(i0,J,L)*(T(1,J,l)-T(i0-1,J,L))/DX
c & COS((-87.0+1.0*(j-1))*RAD)+
& +V(i0,J,L)*(T(i0,J+1,L)-T(i0,J-1,L))/DY
fv(i0,J,L)=U(i0,J,L)*(q(1,J,l)-q(i0-1,J,L))/DX
c & COS((-87.0+1.0*(j-1))*RAD)+
& +V(i0,J,L)*(q(i0,J+1,L)-q(i0,J-1,L))/DY
enddo
ENDDO
do l=1,k0
CALL BOU1(UG(1,1,L),I0,J0)
CALL BOU1(FV(1,1,L),I0,J0)
enddo
c pause 888
DO I=1,i0
DO J=1,j0
C �����
DO L=1,k0
gtmp(i,j,l)=t(i,j,l)*(1000.0/pp(l))**ck
ENDDO
c ���㶥�㴹ֱ�ٶ�
WT(I,J)=-(TG(I,J,k0)+UG(I,J,k0))/
&((GTMP(I,J,k0)-GTMP(I,J,k0-1))/DP*(PP(k0)/1000.0)**CK)
C ������ڵ��ε���Ͳ�
corog(i,j)=orog(i,j)*g
if(h(i,j,1).gt.orog(i,j))then
ip=1
else
DO L=1,k0-1
IF((H(I,J,L).LE.OROG(I,J)).AND.
& (H(I,J,L+1).GT.OROG(I,J)))THEN
IP=L+1
ELSE
ENDIF
ENDDO
endif
ILAY(I,J)=IP
C ������¡�������ѹ---��Ԫ�������
ts(i,j)=t(i,j,ip)-gama*(orog(i,j)-h(i,j,ip))
C �ڲ������ѹ�������
if(ip.gt.1)then
CC=(H(I,J,IP)-OROG(I,J))/(H(I,J,IP)-H(I,J,IP-1))
CD=(PP(IP-1)/PP(IP))**CC
PS(I,J)=CD*PP(IP)
US(I,J)=U(I,J,IP)+CC*(U(I,J,IP)-U(I,J,IP+1))
VS(I,J)=V(I,J,IP)+CC*(V(I,J,IP)-V(I,J,IP+1))
else
ps(i,j)=pp(ip)*(1-gama*(orog(i,j)-h(i,j,ip))/t(i,j,ip))**dk
us(i,j)=u(i,j,ip)+(U(I,J,IP)-U(I,J,IP+1))*
&(OROG(I,J)-H(I,J,IP))/(H(I,J,IP)-H(I,J,IP+1))
US(I,J)=0.6*US(I,J)
Vs(i,j)=V(i,j,ip)+(V(I,J,IP)-V(I,J,IP+1))*
&(OROG(I,J)-H(I,J,IP))/(H(I,J,IP)-H(I,J,IP+1))
VS(I,J)=0.6*VS(I,J)
endif
enddo
enddo
C ������洹ֱ�ٶ�
do i=2,i0-1
do j=2,j0-1
F1=(F0+FLOAT(J-1)*DLATD)*COV
DX=2.0*RG*COS(F1)*DLON
WS(I,J)=(-1.0)*PS(I,J)*g/Rd/Ts(I,J)*
& (US(I,J)*(OROG(I+1,J)-OROG(I-1,J))/DX
c & COS((-87.0+1.0*(j-1))*RAD)+
& +VS(I,J)*(OROG(I,J+1)-OROG(I,J-1))/DY)
ENDDO
ENDDO
CALL BOU1(WS,I0,J0)
c do j=2,j0-1
cWS(1,J)=(-1.0)*PS(1,J)*g/Rd/Ts(1,J)*
c & US(1,J)*(OROG(2,J)-OROG(i0-1,J))/DX
ccc & COS((-87.0+1.0*(j-1))*RAD)+
c & +VS(I,J)*(OROG(I,J+1)-OROG(I,J-1))/DY
cWS(i0,J)=(-1.0)*PS(i0,J)*g/Rd/Ts(i0,J)*
c & (US(i0,J)*(OROG(1,J)-OROG(i0-1,J))/DX
c & COS((-87.0+1.0*(j-1))*RAD)+
c & +VS(I,J)*(OROG(I,J+1)-OROG(I,J-1))/DY)
c enddo
C �������ɢ��
DO I=2,i0-1
DO J=2,j0-1
F1=(F0+FLOAT(J-1)*DLATD)*COV
DX=2.0*RG*COS(F1)*DLON
TGA=TAN(F1)/RG
DO L=1,k0
DIVG(I,J,L)=(U(I+1,J,L)-U(I-1,J,L))/DX
c & /COS((-87.0+1.0*(j-1))*RAD)/(2*DD)+
& +(V(I,J+1,L)-V(I,J-1,L))/DY-V(I,J,L)*TGA
ENDDO
ENDDO
ENDDO
iI=1
DO J=2,j0-1
F1=(F0+FLOAT(J-1)*DLATD)*COV
DX=2.0*RG*COS(F1)*DLON
TGA=TAN(F1)/RG
DO L=1,k0
DIVG(II,J,L)=(U(II+1,J,L)-U(i0,J,L))/DX
& +(V(II,J+1,L)-V(II,J-1,L))/DY-V(II,J,L)*TGA
ENDDO
ENDDO
iI=i0
DO J=2,j0-1
F1=(F0+FLOAT(J-1)*DLATD)*COV
DX=2.0*RG*COS(F1)*DLON
TGA=TAN(F1)/RG
DO L=1,k0
DIVG(II,J,L)=(U(1,J,L)-U(II-1,J,L))/DX
c & /COS((-87.0+1.0*(j-1))*RAD)/(2*DD)+
& +(V(II,J+1,L)-V(II,J-1,L))/DY-V(II,J,L)*TGA
ENDDO
ENDDO
do l=1,k0
CALL BOU1(DIVG(1,1,l),I0,J0)
enddo
C �������ɢ��
DO I=2,i0-1
DO J=2,j0-1
F1=(F0+FLOAT(J-1)*DLATD)*COV
DX=2.0*RG*COS(F1)*DLON
TGA=TAN(F1)/RG
DIVS(I,J)=(US(I+1,J)-US(I-1,J))/DX
&+(VS(I,J+1)-VS(I,J-1))/DY-VS(I,J)*TGA
ENDDO
ENDDO
iI=1
DO J=2,j0-1
F1=(F0+FLOAT(J-1)*DLATD)*COV
DX=2.0*RG*COS(F1)*DLON
TGA=TAN(F1)/RG
DIVS(II,J)=(US(II+1,J)-US(i0,J))/DX
&+(VS(II,J+1)-VS(II,J-1))/DY-VS(II,J)*TGA
ENDDO
iI=i0
DO J=2,j0-1
F1=(F0+FLOAT(J-1)*DLATD)*COV
DX=2.0*RG*COS(F1)*DLON
TGA=TAN(F1)/RG
DIVS(II,J)=(US(1,J)-US(II-1,J))/DX
c & /COS((-87.0+1.0*(j-1))*RAD)/(2*DD)+
& +(VS(II,J+1)-VS(II,J-1))/DY-VS(II,J)*TGA
ENDDO
CALL BOU1(DIVS,I0,J0)
C ���㴹ֱ�ٶ�
do i=1,i0
do j=2,j0-1
IP=ILAY(I,J)
if(ip.GT.1)THEN
DO L=1,IP-1
W(I,J,L)=-9.99E33
ENDDO
ELSE
ENDIF
W(I,J,IP)=WS(I,J)-(DIVS(I,J)+DIVG(I,J,IP))/2*(PP(IP)-PS(I,J))
do L=IP+1,k0
W(I,J,L)=W(I,J,L-1)-(DIVG(I,J,L-1)+DIVG(I,J,L))/2*DP
ENDDO
MM=(k0-IP+1)*(k0-IP+2)
DO L=IP,k0
LL=L-IP+1
W(I,J,L)=W(I,J,L)-LL*(LL+1)*(W(I,J,k0)-WT(i,j))/REAL(MM)
ENDDO
C ����λ�´�ֱ����
Cif(ip.gt.1)then
Cdo l=1,ip-1
Cwg(i,j,l)=-9.99e33
Cenddo
Celse
Cendif
do l=ip+1,k0-1
wg(i,j,l)=W(I,J,L)*(GTMP(I,J,L+1)-GTMP(I,J,L-1))/2/DP
fw(i,j,l)=W(I,J,L)*(q(I,J,L+1)-q(I,J,L-1))/2/DP
ENDDO
WG(I,J,IP)=W(I,J,IP)*(GTMP(I,J,IP+1)-GTMP(I,J,IP))/DP
WG(I,J,k0)=W(I,J,k0)*(GTMP(I,J,k0)-GTMP(I,J,k0-1))/DP
fW(I,J,IP)=W(I,J,IP)*(q(I,J,IP+1)-q(I,J,IP))/DP
fW(I,J,k0)=W(I,J,k0)*(q(I,J,k0)-q(I,J,k0-1))/DP
ENDDO
ENDDO
do L=ip+1,k0-1
CALL BOU1(WG(1,1,L),I0,J0)
CALL BOU1(fW(1,1,L),I0,J0)
enddo
C ���������
DO I=1,i0
DO J=2,j0-1
IP=ILAY(I,J)
IF(IP.GT.1)THEN
DO L=1,IP-1
HT(I,J,L)=9.99e20
ff(i,j,l)=9.99e20
tg(i,j,l)=9.99e20
ug(i,j,l)=9.99e20
wg(i,j,l)=9.99e20
ft(i,j,l)=9.99e20
fv(i,j,l)=9.99e20
fw(i,j,l)=9.99e20
ENDDO
ELSE
ENDIF
DO L=IP,k0
wg(i,j,l)=WG(I,J,L)*(PP(L)/1000.0)**CK
ft(i,j,l)=-ft(i,j,l)*lh/cp
fv(i,j,l)=-fv(i,j,l)*lh/cp
fw(i,j,l)=-fw(i,j,l)*lh/cp
HT(I,J,L)=TG(I,J,L)+UG(I,J,L)+WG(I,J,L)
ff(I,J,L)=fT(I,J,L)+fv(I,J,L)+fW(I,J,L)
ENDDO
SS=0.0
sx=0.0
do L=IP,k0-1
SI=(HT(I,J,L)+HT(I,J,L+1))*DP/2.0
SS=SS-SI
Sj=(ff(I,J,L)+ff(I,J,L+1))*DP/2.0
Sx=Sx-Sj
ENDDO
Q1(I,J)=SS/G*100.0*CP
q2(i,j)=sx/g*100.0*cp
ENDDO
ENDDO
do L=ip,k0-1
CALL BOU1(HT(1,1,L),I0,J0)
CALL BOU1(ff(1,1,L),I0,J0)
enddo
CALL BOU1(q1,I0,J0)
CALL BOU1(q2,I0,J0)
irec=irec+1
write(20,REC=irec)(((tg(I,J,L),I=1,i0),J=1,j0),l=1,k0-1)
write(21,REC=irec)(((ug(I,J,L),I=1,i0),J=1,j0),l=1,k0-1)
write(22,REC=irec)(((wg(I,J,L),I=1,i0),J=1,j0),l=1,k0-1)
write(23,REC=irec)(((w(I,J,L),I=1,i0),J=1,j0),l=1,k0-1)
write(24,REC=irec)(((ft(I,J,L),I=1,i0),J=1,j0),l=1,k0-1)
write(25,REC=irec)(((fv(I,J,L),I=1,i0),J=1,j0),l=1,k0-1)
write(26,REC=irec)(((fw(I,J,L),I=1,i0),J=1,j0),l=1,k0-1)
jrec=jrec+1
write(27,REC=jrec)((q1(I,J),I=1,i0),J=1,j0)
jrec=jrec+1
write(27,REC=jrec)((q2(I,J),I=1,i0),J=1,j0)
print*,irec,jrec
1999continue
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine smth9(i0,j0,ff,gg)
real ff(i0,j0),gg(i0,j0)
do i=2,i0-1
do j=2,j0-1
s=0.5
gg(i,j)=ff(i,j)+s/2.0*(1-s)*(ff(i+1,j)+ff(i-1,j)+ff(i,j+1)
& +ff(i,j-1)-4*ff(i,j))+s*s/4*(ff(i+1,J+1)+ff(i-1,j-1)+
&ff(i+1,j-1)+ff(i-1,j+1)-4*ff(i,j))
enddo
enddo
i=1
do j=2,j0-1
gg(i,j)=ff(i,j)+s/2.0*(1-s)*(ff(i+1,j)+ff(i0,j)+ff(i,j+1)
& +ff(i,j-1)-4*ff(i,j))+s*s/4*(ff(i+1,J+1)+ff(i0,j-1)+
&ff(i+1,j-1)+ff(i0,j+1)-4*ff(i,j))
enddo
i=i0
do j=2,j0-1
gg(i,j)=ff(i,j)+s/2.0*(1-s)*(ff(1,j)+ff(i-1,j)+ff(i,j+1)
& +ff(i,j-1)-4*ff(i,j))+s*s/4*(ff(1,J+1)+ff(i-1,j-1)+
&ff(i,j-1)+ff(i-1,j+1)-4*ff(i,j))
enddo
j=1
do i=1,i0
gg(i,j)=ff(i,j)
enddo
j=j0
do i=1,i0
gg(i,j)=ff(i,j)
enddo
return
end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SUBROUTINE BOU1(AA,MX,MY)
DIMENSION AA(MX,MY)
SS=.125
M=0
DO 1 I=2,MX-1
AA(I,1)=AA(I,2)
AA(I,MY)=AA(I,MY-1)
1 CONTINUE
DO 2 J=2,MY-1
AA(1,J)=AA(2,J)
AA(MX,J)=AA(MX-1,J)
2 CONTINUE
AA(1,1)=AA(2,1)
AA(1,MY)=AA(2,MY)
AA(MX,1)=AA(MX,2)
AA(MX,MY)=AA(MX,MY-1)
RETURN
END
IDL代码
Q1Q2_Y73: Calculate apparent heat source and moisture sink, after Yanai et al. (1973). Files required:
div_vel_sph.pro
omega_vel.pro
press2alt.pro (by Dominik Brunner)
q1q2_y73.pro地址：http://www.johnny-lin.com/lib.shtml
参考：
http://www.ncl.ucar.edu/Applications/wind.shtml2.
http://read.pudn.com/downloads378/sourcecode/others/1631672/calswa.f__.htm
http://www.johnny-lin.com/lib.shtml
文章来源于微信公众号：气象学家
666
页:
[1]