NCL计算Plumb波活动通量

是冉冉升起的冉

结果跟摸鱼的对比了一下，但感觉在大西洋地区有一点差异，大体趋势差不多，20°低纬地区有一点奇怪，然后垂直方向没有对比；有问题联系我ing；

1. ;///////////////////////////////////////////////////////////////
   
2. ;  Calculate 3-D wave-activity flux:plumb[1985]
   
3. ;///////////////////////////////////////////////////////////////
   
4. ;  Written by Ran,20240606.
   
5. ;  E-mail: 602023280016@smail.nju.edu.cn
   
6. 
   
7. ;If there is a problem with the code calculation, please contact me by email.
   
8. 
   
9. ;Input data：hgt t u v
   
10. 
    
11. ;Note that the data downloaded by ERA5 is the gravitational potential (unites:m^2/s^2).
    
12. ;The ncep data is the gravity potential height z, which needs to be multiplied by g.
    
13. 
    
14. ;Input data must include the entire longitude！！！ Because to take the zonal deviation!!
    
15. 
    
16. ;Return value -  Four multi-dimentional arrays containing:
    
17. ;                    x-component:  Fx
    
18. ;                    y-component:  Fy
    
19. ;                    z-component:  Fz
    
20. ;                    z_ano
    
21. 
    
22. ;Usage: [/Fx,Fy,Fz,z_ano/] = plumb_pLLL(hgt,t,u,v)
    
23. 
    
24. 
    
25. undef("plumb_pLLL")
    
26. function plumb_pLLL(hgt[*][*][*][*]:numeric,t[*][*][*][*]:numeric,u[*][*][*][*]:numeric,v[*][*][*][*]:numeric)
    
27. local re,sclhgt,k,pi,dName,level,lat,lon,f,wmg,ga,z,coslat,sin2lat,leveltmp,ftmp,coslattmp,sin2lattmp,t_ano,u_ano,v_ano,t_mean,dt_meandz,\
    
28.         S,Stmp,dvzdlon,duzdlon,dtzdlon
    
29. begin
    
30. ;基本量
    
31. 
    
32. ; Radius of the earth
    
33. re      = 6.37*10^6;6378388
    
34. 
    
35. ; scale height
    
36. sclhgt  = 8000.
    
37. 
    
38. k=0.286
    
39. 
    
40. ; pi
    
41. pi      = atan(1.0)*4.
    
42. printVarSummary(hgt)
    
43. dName   = getvardims(hgt)
    
44. if (any(ismissing(dName(1:)))) then
    
45.    print("fatal: plumb_pLLL: requires that all the rightmost dimensions be named")
    
46.    exit
    
47. end if
    
48.     level   = hgt&$dName(1)$
    
49.     lat     = hgt&$dName(2)$
    
50.     lon     = hgt&$dName(3)$
    
51. 
    
52. ; Coriolis parameter
    
53. f       = lat(:)
    
54. f       = (/2.*2.*pi/(60.*60.*24.)*sin(pi/180.*f)/)
    
55. f@_FillValue = hgt@_FillValue
    
56. 
    
57. ;Rotational angular velocity of the earth
    
58. wmg=2.*pi/(60.*60.*24.)
    
59. 
    
60. ; Gravitational acceleration
    
61. ga      = 9.80665
    
62. z=hgt*ga
    
63. copy_VarMeta(hgt, z)
    
64. ; cosine
    
65. coslat  = cos(lat(:)*pi/180.)
    
66. coslat@_FillValue = default_fillvalue(typeof(coslat))
    
67. coslat  = (/where(coslat.le.0.,coslat@_FillValue,coslat)/)
    
68. 
    
69. ; sin2lat
    
70. sin2lat  = sin(2*(lat(:)*pi/180.))
    
71. sin2lat@_FillValue = default_fillvalue(typeof(sin2lat))
    
72. sin2lat  = (/where(sin2lat.le.0.,sin2lat@_FillValue,sin2lat)/)
    
73. 
    
74. 
    
75. ; 1-D -> 4-D
    
76. leveltmp   = conform_dims(dimsizes(hgt),level,1)
    
77. ftmp       = conform_dims(dimsizes(hgt),f,2)
    
78. coslattmp  = conform_dims(dimsizes(hgt),coslat,2)
    
79. sin2lattmp = conform_dims(dimsizes(hgt),sin2lat,2)
    
80. ;Zonal deviation
    
81. t_ano   =dim_rmvmean_n_Wrap(t, 3)
    
82. u_ano   =dim_rmvmean_n_Wrap(u, 3)
    
83. v_ano   =dim_rmvmean_n_Wrap(v, 3)
    
84. z_ano   =dim_rmvmean_n_Wrap(z, 3)
    
85. ;printVarSummary(t_ano)
    
86. ;printVarSummary(z_ano)
    
87. 
    
88. ;Zonal mean
    
89. t_mean=dim_avg_n_Wrap(t, 3)
    
90. dt_meandz = center_finite_diff_n(t_mean,-sclhgt*log(level/1000.),False,0,1)
    
91. 
    
92. 
    
93. S=dt_meandz+k*t_mean/sclhgt;
    
94. copy_VarCoords_1(hgt, S)
    
95. ;printVarSummary(S)
    
96. ; 3-D -> 4-D
    
97. Stmp = conform_dims(dimsizes(hgt),S,(/0,1,2/))
    
98. ;printVarSummary(Stmp)
    
99. 
    
100. ; dvz/dlon
     
101. dvzdlon      =  center_finite_diff_n(v_ano*z_ano,lon*pi/180.,True,0,3)
     
102. ; duz/dlon
     
103. duzdlon      =  center_finite_diff_n(u_ano*z_ano,lon*pi/180.,True,0,3)
     
104. ; dtz/dlon
     
105. dtzdlon      =  center_finite_diff_n(t_ano*z_ano,lon*pi/180.,True,0,3)
     
106. 
     
107. ;printVarSummary(sin2lattmp)
     
108. ;printVarSummary(dvzdlon)
     
109. ;printVarSummary(wmg)
     
110. ;printVarSummary(coslattmp)
     
111. 
     
112. ; x-component
     
113. Fx      = leveltmp/1000.*coslattmp*v_ano*v_ano-1./(2.*wmg*re*sin2lattmp)*dvzdlon
     
114. ;printVarSummary(Fx)
     
115. 
     
116. ; y-component
     
117. Fy      = leveltmp/1000.*coslattmp*(-u_ano*v_ano+1./(2.*wmg*re*sin2lattmp)*duzdlon)
     
118. 
     
119. ; z-component
     
120. Fz      = leveltmp/1000.*coslattmp*(ftmp/Stmp*(v_ano*t_ano-1./(2.*wmg*re*sin2lattmp)*dtzdlon))
     
121. 
     
122. copy_VarCoords(hgt,Fx)
     
123. copy_VarCoords(Fx,Fy)
     
124. copy_VarCoords(Fx,Fz)
     
125. 
     
126. Fx@units = "m^2/s^2"
     
127. Fx@var_desc = "wave-activity flux"
     
128. Fx@long_name = "x-component of wave-activity flux (Plumb1985)"
     
129. 
     
130. Fy@units = "m^2/s^2"
     
131. Fy@var_desc = "wave-activity flux"
     
132. Fy@long_name = "y-component of wave-activity flux (Plumb1985)"
     
133. 
     
134. Fz@units = "Pa m/s^2"
     
135. Fz@var_desc = "wave-activity flux"
     
136. Fz@long_name = "z-component of wave-activity flux (Plumb1985)"
     
137. 
     
138. return [/Fx,Fy,Fz,z_ano/]
     
139. 
     
140. end

复制代码

是冉冉升起的冉

Fx计算的时候发现少打了一个括号ing (v_ano*v_ano-1./(2.*wmg*re*sin2lattmp)*dvzdlon)

隽今、

您好，我想问一下我下载的NCEP数据里纬度有负的（-90，90），在计算 f 的时候要把负的地方手动变成正值吗？还是不管？

是冉冉升起的冉

隽今、 发表于 2025-3-16 11:22
您好，我想问一下我下载的NCEP数据里纬度有负的（-90，90），在计算 f 的时候要把负的地方手动变成正值吗？ ...

确实没想过南半球的问题哈哈哈哈，但我觉得计算的时候，南半球直接f=-2Ωsinφ这样吧。

Ausiore

是冉冉升起的冉 发表于 2025-3-17 10:55
确实没想过南半球的问题哈哈哈哈，但我觉得计算的时候，南半球直接f=-2Ωsinφ这样吧。

应该不需要，南半球的纬度本身就是负的，f=2*Ω*sin(lat)，南半球的lat是负的，所以f<0，这很合理