(2022/11/13發表於個人臉書)
博士後心得Part 6.2，來寫我把一個邊界層/對流參數法，MYNN-EDMF ( Mellor-Yamada-Nakanishi-Niino Eddy-Diffusivity/Mass-Flux)，放進GFDL氣候模式AM4的過程和心得。MYNN-EDMF是一個在WRF裡頭的參數法，JPL的合作者是發展者之一，我們會想把這個參數法放進AM4，是想看這個參數法能不能改進AM4在海洋層積雲的模擬結果。
我大概花了一年半左右做這工作。首先把MYNN-EDMF的程式放進AM4裡，包括寫介面、新增變數之類，這大概花了一兩個月。剩下的一年多就是不停的測試，發現各種問題，然後修改。雖然最後模擬結果不甚理想，但整個過程我學到非常多東西，不僅是學習如何處理研究上遇到的困難，還親眼見識和參與老闆Leo的思考過程，收穫良多。
1. Cloud scheme incompatibility
MYNN-EDMF處理的變數是總水含量 (total water mixing ratio; q_t) 和 冰-水位溫 (ice-liquid water potential temperature; θ_li)。這兩個變數在沒有降水的濕過程是守恆的，所以在計算紊流的混合過程時，不需要考慮來源項和消散項 (source and sink terms)，方程式會比較簡單。很多邊界層的研究和模式常常使用q_t和θ_li這兩個變數。
要把q_t和θ_l轉成常用的溫度、水氣、雲水含量、和雲量等等，MYNN-EDMF使用一個概率密度函數的雲參數方案 (Probability Density Function cloud scheme; PDF cloud scheme)，將一組 (q_t, θ_li) 對應到一組 (T, q, q_l, q_a)。這可以來計算紊流混合之前和之後的溫度、水氣、雲水含量、和雲量，並計算紊流對這些變數造成的改變。
但是，AM4已經有自己的雲參數化方案，Tiedtke方案，想讓MYNN-EDMF的PDF方案和Tiedtke方案共存，會遇到非常大的問題。比方說想計算紊流造成雲量的變化量，可以用兩種方式計算：
1: dφ/dt = (φEDMFt+ᇫt - φEDMFt ) / ᇫt
2: dφ/dt = (φEDMFt+ᇫt - φAM4t ) / ᇫt
φ = {T, q, q_l, q_a}，t是時間，ᇫt是模式的時間步長 (time step; AM4是30分鐘)，下標EDMF代表由MYNN-EDMF自己的PDF雲方案算出，下標AM4代表由AM4的雲方案算出的。
方式一是MYNN-EDMF原本的方式，在內部是一致的，但要回傳給AM4時會遇到問題。比方說假設MYNN-EDMF計算出紊流使得雲量減少30%，但AM4可能根本沒有雲，雲量不可能因為紊流而減少。
方式二的問題是dφ/dt同時包含紊流和雲參數方案的影響。比方說紊流很弱對雲量沒有影響，但是B方式得到的dφ/dt可能不是零，因為φEDMFt+ᇫt 不等於 φAM4t。
此外，要怎麼將MYNN-EDMF算出的雲量和雲水變化傳給AM4的Tiedtke方案，又是另一個很難處理的問題。Tiedtke方案是預報型方案 (prognostic scheme)，考慮了大尺度抬升、雲內濕空氣和雲外乾空氣的混合、輻射和紊流冷卻對雲的生成和消散等過程。MYNN-EDMF算出的雲量和雲水變化，是要把其當作Tiedtke方案的額外過程，還是取代已有的過程，並不是很明確。
我做了一些測試，得到的結論是在AM4同時使用MYNN-EDMF的PDF雲方案和Tiedtke雲方案是條死路，走不下去。老闆Leo、我、還有合作者在想有沒有辦法不使用MYNN-EDMF的PDF雲方案，我們最後想出了兩個方法：
讓MYNN-EDMF去處理紊流對溫度和水氣的影響，把雲量和雲水視為普通物質 (passive tracer)，紊流會移動雲量和雲水，但不牽扯任何生成或消散過程，這個方式類似於AM4原本的邊界層參數化方案。不過，這個方式需要大幅修改MYNN-EDMF程式，也失去使用q_t和θ_li的好處。
還是讓MYNN-EDMF去處理q_t和θ_li，但去診斷出紊流ED和MF對雲水和雲量的改變量。知道雲水的改變量後，就可以回推紊流對溫度和水氣的改變量。這個方式不需要使用PDF雲方案，還保有使用q_t和θ_l的好處。這方法需要多寫個程式根據MYNN-EDMF的資訊，診斷出紊流ED和MF對雲水和雲量的改變量，實務上也是可行的。
我們最後決定試B方法，不過遇到了MF逸出率可能是負值的問題 (negative detrainment rate)，這也讓我們和合作者有些激烈的討論。下個小節來簡述這個逸出率負值的問題，以及兩邊激烈的討論。
2. Negative detrainment rate in MF
B方法要診斷出紊流ED和MF對雲水和雲量的改變量。
ED的話，我用AM4原本就有的紊流擴散程式 (vertical diffusion solver) 配上MYNN-EDMF提供的紊流擴散係數 (eddy-diffusion coefficients)，就可以診斷出來。因為AM4的紊流擴散程式和MYNN-EDMF在不同的模組 (一個在physics_up，另一個在physics_down)，我花了一點時間處理了模組間變數傳遞的技術問題。
MF的話，需要質量通量 (mass flux)和逸出量 (detrainment rate) 來診斷出MF對網格平均的雲水和雲量的改變量。我重新推導了次網格的質量通量如何影響網格平均量的方程式，在假設上升氣流只佔很小一部分的網格區域，Eddy-flux/Source和Detrainment/Subsidence的方程式是同等的。久違的推導方程式。
MYNN-EDMF有提供質量通量，但沒有逸出量，所以我自己在裡頭寫了程式去算，用這個方程式：
(1/M_i) dM_i/dz = ε_i - δ_i
M是質量通量，z是高度，ε是逸入率 (fractional entrainment rate)，δ是逸出率 (fractional detrainment rate)，i是第i個上升氣流塊 (plume)。
在MYNN-EDMF中，假設上升氣流的範圍不隨高度改變，逸入率ε是隨機的，由此計算上升氣流的垂直速度以及質量通量。因為M_i和ε_i是已知的，就可以算出δ_i。
在算δ_i時，我發現在某些高度上δ_i是負值，這是不合理的。從定義上來說，δ_i必須是正值，要不然δ_i就會是逸入率，而不是逸出率了。δ_i是負值也代表MYNN-EDMF違反質量守恆定律，因為逸入量 (M_i * ε_i) 不足以彌補質量通量的變化 (dM_i/dz)，有些質量必須無中生有來滿足質量守恆定律。
我老闆Leo覺得δ_i是負值，MYNN-EDMF違反質量守恆定律是非常嚴重的問題，除非能解決，要不然GFDL這邊不可能接受把這個MYNN-EDMF方案放進模式裡頭。JPL合作者Joao和Kay承認δ_i是負值的確是個問題，但他們不覺得這必須是最優先解決的事項。他們的論點是模式裡頭用很簡化的上升氣流塊 (plume) 來表示複雜的積雲對流，是因為數學上比較好處理，物理上，沒有人知道積雲對流適不適合用上升氣流塊來表示。既然都能接受用上升氣流塊來表示複雜的積雲對流這個很大的假設，或許把力氣花在改進上升氣流塊的其他部分，而不是把逸出率負值當成最優先的事項。
這個討論持續了好幾週，有時候氣氛蠻火爆的。我覺得兩邊的說法都有道理，也都有問題，比較像是做研究的信念不同。最後JPL那邊同意要優先解決逸出率負值的問題，我們討論出幾個方案：
調整MYNN-EDMF的上升氣流範圍，來滿足質量守恆定律。
MYNN-EDMF還是假設上升氣流範圍不隨高度改變，用迭代法找出逸入率來同時滿足上升速度方程式和質量通量方程式。
發展新的逸入率ε計算方式
我們最後採用方案1，因為方案2 & 3的工程太浩大。
3. Tests, tests, and tests
我們最後對MYNN-EDMF做了兩個主要修改 (1) 不使用PDF雲方案; (2) 調整MYNN-EDMF的上升氣流範圍，來滿足質量守恆定律。程式寫完以後，就是不斷地測試、測試、測試。測試程式有沒有bugs，測試模擬會不會爆掉，看模擬結果的表現，然後做很多敏感度測試。
我主要用AM4的一維氣柱模式 (single-column model; SCM) 來測試新增的程式有沒有問題，有沒有bugs。也用SCM跑了海洋層積雲的觀測實驗 DYCOMS-II RF01，把結果和原本的AM4、大渦流模擬結果 (large-eddy simulation)、以及觀測結果做些比較。
SCM測試完了以後，就跑8-30天的全球大氣模擬。第一次跑模式就爆了，因為MYNN-EDMF在某些大氣條件會計算出很不合理的數值，像是2,000 K/day的加熱率，於是就爆掉了。模式爆掉的話，就要去找出問題在哪，是在哪個地方爆掉，然後做些修改，重複測試。
順利找到問題並解決的話，會跑5-10年左右的全球大氣模擬，跟原本的AM4和觀測資料比較。模式結果看起來沒問題的話，會跑30年左右的全球大氣模擬，然後做些比較仔細的結果分析。
我們放MYNN-EDMF進AM4的主要目的是想能不能減少AM4在海洋層積雲區域的誤差。沒想到 (可能也不意外)，MYNN-EDMF讓誤差更大，基本上雲都不見了@@
我做了些檢查和分析，初步結論是問題可能出在MYNN-EDMF跟AM4 Tiedtke雲參數方案的耦合 (coupling)。Tiedtke方案是計算非對流區的雲量，非對流區的垂直速度是用網格平均的垂直速度和對流區的質量通量來計算。因為MYNN-EDMF的質量通量很大，在相同的網格平均的垂直速度下，非對流區的上升速度被抑制，甚至會變成下沉運動，不利雲的生成。另外，Tiedtke方案有一項是雲跟乾空氣的混合效應，稱為 erosion。MYNN-EDMF逸出的雲都被erosion消滅了，不是很清楚為什麼。
看起來由於MYNN-EDMF上升氣流太強，把雲都弄不見了。所以我們調整MYNN-EDMF裡頭的參數，做了一連串的敏感度測試，看能不能把雲帶回來。這些測試包括調整MF的逸入率、邊界層高度太低或紊流動能太小的話把MF關掉、在海洋層積雲地區把MF關掉等等，不過結果都不是很理想。
因為MYNN-EDMF的結果不好，加上JPL的主要合作者Kay離開去私人公司，我們就沒有繼續做更多的模擬和測試。基本上，這個MYNN-EDMF的研究計畫就停止了。
4. Epilogue
我花了一年多的時間 (2020年底到2022年初) 試著把MYNN-EDMF放進AM4裡，整個過程不是很順，一直不斷地的碰壁、碰壁、碰壁。我一開始完全沒有意識到MYNN-EDMF的PDF雲參數方案會造成問題; 在診斷MF的逸出率時，才發現有負值的問題。
在這些不斷碰壁的過程中，我親眼見識和參與老闆Leo的思考過程，我們一同想辦法來面對這些困難。我想，如果是我自己一個人負責這計劃，我應該是很快就投降了，想不到這麼多ideas。偉大的科學家，不是順遂的做研究，找到好題目，做出好結果。而是在研究碰壁後要怎麼處理，怎麼想辦法繞過去，或是轉向其他方向，這才是科學精神。
我覺得我在Leo身上學到最重要的兩個東西是：
1. 科學上的不容妥協之物
MYNN-EDMF違反質量守恆定律是非常嚴重的問題，Leo很堅持這點。雖然用上升氣流塊來表示積雲對流有自己的問題，但還是不應該有逸出率是負值的問題。
我自己可以感覺到JPL合作者Joao和Kay一開始不是很想處理逸出率是負值的問題，是Leo堅持要改。雖然這是雙邊的合作計畫，但科學上的信念還是有不容妥協之物。
2. 嘗試的機會
雖然MYNN-EDMF的初步結果很差，但我們還是做了很多後續測試，因為Leo覺得還是要給它一些機會。
我覺得這種心態在研究上很重要。結果不如預期，但還是給它一些機會，看能不能改進。測試過程中會學到一些經驗，如果只因為結果很差，就直接放棄了，老實說會蠻可惜的。
最後寫一些雜記。
雖然我花了很多時間和心力成功的把MYNN-EDMF放進AM4裡，也學到非常寶貴的經驗，但不是很能把成果發表在期刊上。一是因為AM4 + MYNN-EDMF的模擬結果太差，負面的結果不是很好發表，二是這個工作並沒有很重要的科學發現，主要是處理技術議題。除了真正在做參數法發展和模式開發的人，一般的大氣研究學者對這個MYNN-EDMF的工作可能沒啥興趣。
我是有打算把MYNN-EDMF工作寫一份詳細的報告，留給自己和GFDL。這一篇心得文可算是初稿吧。
Figure. Cloud scheme incompatibility between the MYNN-EDMF and AM4.
Figure. AM4 single column model (SCM) results for the DYCOMS RF01 case. (left) cloud fraction averaged over the simulated hour 3-4. (right) cloud liquid specific humidity averaged over the simulated hour 3-4. Colors line indicate different SCM configurations: the standard AM4 (red), the standard AM4 with its radiation scheme (pink), the AM4 with MYNN ED scheme (orange), the AM4 with the MYNN-EDMF scheme using different entrainment length scale Lo from 100 meter (green; small entrainment) to 25 (purple; large entrainment). Black dash line indicates the large-eddy simulation results provide by one EDMF CPT member, Professor Georgios Matheou. Black dots on the right panel are aircraft measurement taken from Zhu et al. (2005).
Figure: 1980-1984 AMIP, annual mean TOA net↓ solar flux bias w.r.t CERES-EBAF