：從零搭建近岸海域模擬與避坑指南)
簡介POM普林斯頓海洋模型是一套面向海洋數(shù)值模擬初學者與科研人員的開源建模工具可在f77、f90、g77、pgf77、CVF等多種編譯器下直接運行適合用于河道水流、河口及近岸海域的水動力與溫鹽場模擬研究。資源包共21個文件約658KB包含f與c源碼、n與dat輸入數(shù)據(jù)、out運行輸出、runpom2k執(zhí)行腳本、change版本記錄、pdf中文說明書及zip/tar工具包等覆蓋從編譯、算例運行到結果查看的完整流程。POM2K提供兩個可直接上手的實例河道模型參數(shù)已內(nèi)置于代碼中河口模型則附帶原始網(wǎng)格文件其輸入文件由GRID-DATA/GRID.f生成便于用戶理解網(wǎng)格與初始場的構造方式。借助GRID-DATA中的插值代碼可將地形、溫鹽、表面風場及熱通量等資料插值到網(wǎng)格點生成初始條件與強迫條件為自建模型提供有效工具。目前已有997人學習下載適合希望快速搭建POM建模環(huán)境、掌握算例配置與數(shù)據(jù)前處理思路的讀者參考。1. POM 普林斯頓海洋模型從零建立一套可跑通的海域模擬如果你手頭只有一份岸線數(shù)據(jù)、幾組風場和開邊界水位卻要在一周內(nèi)搭出一套能出流場和溫鹽結構的近岸模型POMPrinceton Ocean Model大概率是你繞不開的選擇。它屬于三維原始方程海洋模型采用 sigma 坐標貼合地形垂向用二階湍流閉合時間積分走的是外模式與內(nèi)模式分離的顯式差分。這套結構決定了它對近岸、河口、陸架海域的適應性很強也決定了它對時間步長和網(wǎng)格質(zhì)量相當敏感。我見過不少人拿到源碼包后卡在編譯和參數(shù)配置上其實真正難的不是代碼本身而是搞清楚每個輸入文件對應哪塊物理過程。這篇筆記按“資源是什么、怎么配、坑在哪”的順序拆一遍適合做近岸水動力、溫鹽輸運、工程海域論證的從業(yè)者也適合剛接觸海洋數(shù)值模擬的研究生照著復現(xiàn)。2. POM 的坐標系與差分框架為什么它適合近岸而不是深海2.1 sigma 坐標與內(nèi)外模式分離的取舍POM 最核心的設計是垂向 sigma 坐標也就是把水深歸一化到 0 到 1 之間自由表面和底部地形都能被貼合。這樣做的好處是近岸淺水區(qū)的垂向分辨率不會被水深變化吃掉河口和陸架坡折帶都能保持同樣的垂向?qū)訑?shù)。代價是斜壓梯度項在陡坡處容易產(chǎn)生誤差所以實際配置時地形坡度一般控制在 0.05 以內(nèi)超過這個值要么平滑地形要么加密水平網(wǎng)格。時間積分上POM 把正壓模態(tài)和斜壓模態(tài)拆開。外模式負責自由表面和垂向平均流時間步長受 CFL 條件限制通常取幾秒到十幾秒內(nèi)模式負責溫鹽和三維流場步長可以放大到外模式的幾十倍。這種分離讓計算量降下來但也帶來一個隱患內(nèi)外模式之間的耦合頻率如果設得不合理會出現(xiàn)質(zhì)量不守恒或者表面重力波失真。常見做法是讓內(nèi)外模式步長比控制在 20 到 50 之間具體看網(wǎng)格分辨率和最大水深。2.2 控制方程與湍流閉合的落地含義POM 求解的是 Boussinesq 近似下的三維原始方程包括動量方程、連續(xù)方程、溫度鹽度輸運方程和狀態(tài)方程。動量方程里保留了科氏力、壓強梯度力、湍流黏性項和底摩擦。湍流閉合用的是 Mellor-Yamada 二階閉合方案也就是常說的 2.5 階模型它通過求解湍動能和湍流宏觀尺度兩個預報量來給出垂向混合系數(shù)。對使用者來說這套閉合方案意味著你不需要自己調(diào)混合層深度但需要關注兩個參數(shù)一個是湍動能產(chǎn)生項里的常數(shù)另一個是壁面函數(shù)里的粗糙度。前者在源碼里通常以常量形式出現(xiàn)不建議隨意改后者和底摩擦系數(shù)直接相關近岸模擬里一般取 0.002 到 0.005 之間。如果模擬出來的底層流速偏大先查底摩擦再查垂向分層是否被過度混合。2.3 網(wǎng)格與地形文件的準備步驟POM 的輸入網(wǎng)格是正交曲線網(wǎng)格水平方向用 Arakawa C 網(wǎng)格離散也就是速度分量分別位于網(wǎng)格邊中點標量位于網(wǎng)格中心。準備網(wǎng)格時第一步是確定研究區(qū)域和開邊界位置第二步是生成正交曲線網(wǎng)格第三步是把水深插值到網(wǎng)格點上并做平滑。常見做法是用正交曲線網(wǎng)格生成工具先出網(wǎng)格文件再導出 POM 需要的 fortran 格式。水深數(shù)據(jù)一般來自實測或公開海圖插值后要做三次平滑每次平滑系數(shù)不超過 0.5否則會抹掉真實地形特征。下面是一個把水深數(shù)據(jù)插值到 POM 網(wǎng)格并平滑的 Python 片段用的是常見科學計算庫。import numpy as np from scipy.interpolate import griddata # 讀取原始水深散點lon, lat, depth raw np.loadtxt(depth_raw.txt) lon, lat, depth raw[:, 0], raw[:, 1], raw[:, 2] # 讀取 POM 網(wǎng)格點坐標 grid np.loadtxt(pom_grid.txt) gx, gy grid[:, 0], grid[:, 1] # 插值到網(wǎng)格點近岸用線性外?？蓳Q cubic depth_grid griddata((lon, lat), depth, (gx, gy), methodlinear) # 填補插值空洞避免陸地邊界出現(xiàn) NaN mask np.isnan(depth_grid) if mask.any(): depth_grid[mask] griddata( (lon, lat), depth, (gx[mask], gy[mask]), methodnearest ) # 三次平滑系數(shù) 0.3保留地形趨勢 smoothed depth_grid.copy() for _ in range(3): smoothed[1:-1, 1:-1] ( smoothed[1:-1, 1:-1] * 0.4 (smoothed[:-2, 1:-1] smoothed[2:, 1:-1] smoothed[1:-1, :-2] smoothed[1:-1, 2:]) * 0.15 ) np.savetxt(depth_pom.txt, smoothed, fmt%.3f)這段代碼的邏輯是先做散點插值再處理邊界空洞最后做三次加權平滑。參數(shù)上平滑系數(shù) 0.3 是經(jīng)驗值地形復雜的河口可以降到 0.2開闊陸架可以升到 0.4。插值方法選 linear 是為了避免 cubic 在數(shù)據(jù)稀疏區(qū)產(chǎn)生過沖如果原始水深點足夠密可以換成 cubic 提升光滑度。輸出格式按 POM 讀取要求寫成純文本每行一個網(wǎng)格點的水深值。3. 編譯與輸入文件配置把源碼跑起來的關鍵幾步3.1 編譯環(huán)境與 Makefile 調(diào)整POM 源碼通常是 Fortran 77 和 Fortran 90 混合編譯前先確認編譯器支持固定格式和自由格式混編。常見做法是用 gfortran在 Makefile 里指定-ffixed-line-length-132和-fno-second-underscore避免長行被截斷和符號名不匹配。如果用的是 ifort對應選項是-extend-source和-assume nounderscore。編譯順序一般是先編公共塊和參數(shù)模塊再編主程序。源碼包里通常有一個makefile或Makefile里面已經(jīng)寫好了依賴關系但路徑需要按本機環(huán)境改。我一般會先跑一遍make clean再make看第一個報錯出現(xiàn)在哪個文件。如果是找不到netcdf.inc說明輸出用了 NetCDF需要先裝好庫并在 Makefile 里加-I和-L路徑。# 清理舊目標文件 make clean # 指定編譯器和 NetCDF 路徑后編譯 make FCgfortran \ FFLAGS-O2 -ffixed-line-length-132 -fno-second-underscore \ NETCDF_INC/usr/local/include \ NETCDF_LIB/usr/local/lib -lnetcdf這里-O2是常規(guī)優(yōu)化調(diào)試階段可以換成-g -O0方便定位數(shù)組越界。NETCDF_INC和NETCDF_LIB按實際安裝路徑改如果不用 NetCDF 輸出可以把相關源文件從 Makefile 里注釋掉。編譯成功后會在當前目錄生成可執(zhí)行文件名字通常是pom.exe或pom2k。3.2 輸入文件清單與參數(shù)對應關系POM 的輸入文件一般包括網(wǎng)格文件、水深文件、初始溫鹽文件、開邊界文件、風場和熱通量文件。每個文件對應一個 Fortran 讀取單元單元號在源碼的params或infile里定義。下面這張表列出常見文件和關鍵參數(shù)方便對照檢查。文件名內(nèi)容關鍵參數(shù)常見取值grid.dat網(wǎng)格坐標網(wǎng)格點數(shù)、正交性由生成工具決定depth.dat水深最小水深不小于 2 minit_ts.dat初始溫鹽表層溫度、鹽度按季節(jié)實測obc.dat開邊界水位、流速、溫鹽潮汐或?qū)崪ywind.dat風場風速、風向10 m 高度heat.dat熱通量短波、長波、感熱按月份配置時最容易出錯的是開邊界。POM 支持輻射邊界和指定邊界兩種輻射邊界適合外海開邊界指定邊界適合有實測潮位的斷面。如果開邊界給的是潮汐調(diào)和常數(shù)需要在源碼里加一個潮汐預報模塊或者提前把水位時間序列算好寫進obc.dat。風場和熱通量的時間步長要和模型時間步長對齊否則會插值出錯。3.3 時間步長與 CFL 條件的實際估算時間步長的選取直接決定模型能不能穩(wěn)定跑完。外模式步長受 CFL 條件限制公式是dt dx / sqrt(g * H_max)其中dx是最小網(wǎng)格間距H_max是最大水深g是重力加速度。舉個例子最小網(wǎng)格 500 米最大水深 100 米算出來dt要小于 15 秒左右。實際配置時取 10 秒留余量。內(nèi)模式步長通常是外模式的 20 到 50 倍但不要超過 60 倍否則斜壓調(diào)整跟不上。如果模擬區(qū)域有強潮汐內(nèi)外步長比取小一點比如 20如果是陸架環(huán)流可以取到 40。跑之前先用理想地形做一次冷啟動看自由表面是否發(fā)散再換真實地形。# 冷啟動檢查先用 10 秒外步長跑 1 天 ./pom.exe run_control.txt # 查看輸出日志里的 CFL 數(shù)和表面水位極值 grep -E CFL|SSH pom.log | tail -20如果日志里 CFL 數(shù)持續(xù)大于 1說明步長偏大或者地形有尖峰需要減小步長或平滑地形。表面水位極值如果出現(xiàn)異常大值先查開邊界是否給了不合理的水位再查底摩擦是否過小。4. 避坑與排查跑 POM 時最容易翻車的五個地方4.1 編譯通過但運行報數(shù)組越界現(xiàn)象是編譯沒問題一運行就Segmentation fault或者Fortran runtime error: Index out of bounds。原因通常是網(wǎng)格文件里的點數(shù)比源碼里聲明的數(shù)組維度大或者水深文件的行數(shù)和網(wǎng)格點數(shù)不一致。解決方法是先檢查params文件里的IM和JM是否和網(wǎng)格文件頭一致再逐行核對水深文件的行數(shù)。我一般會寫一個簡單的 Python 腳本統(tǒng)計行數(shù)和網(wǎng)格文件對比。4.2 自由表面發(fā)散或水位異?，F(xiàn)象是跑了幾十步后表面水位指數(shù)增長日志里 CFL 數(shù)越來越大。原因可能是外模式步長超過 CFL 限制或者開邊界水位給了不合理的突變。解決方法是先把外步長減半試跑如果穩(wěn)定了說明是步長問題如果還發(fā)散檢查開邊界文件的時間序列是否有跳變必要時做線性插值平滑。4.3 溫鹽出現(xiàn)負值或超過物理范圍現(xiàn)象是溫鹽場跑著跑著出現(xiàn)負鹽度或者溫度超過 40 度。原因通常是平流項用了中心差分在強梯度區(qū)產(chǎn)生數(shù)值振蕩。POM 里可以選擇迎風差分或者加水平擴散。解決方法是把平流方案改成迎風或者把水平擴散系數(shù)調(diào)大但不要超過 100 m2/s否則會過度抹平鋒面。4.4 底摩擦系數(shù)取錯導致底層流速失真現(xiàn)象是表層流速合理底層流速偏大或偏小。原因是底摩擦系數(shù)沒有按底質(zhì)類型調(diào)整。沙質(zhì)海底取 0.002 到 0.003泥質(zhì)海底取 0.004 到 0.005。如果模擬區(qū)域底質(zhì)不均可以按區(qū)域給不同系數(shù)但要在源碼里加判斷。我一般會先用一個中間值跑一遍再和實測流速對比微調(diào)。4.5 輸出文件太大導致磁盤寫滿現(xiàn)象是跑了一半程序卡住日志沒有報錯但輸出目錄不再更新。原因是 NetCDF 或二進制輸出沒有做時間降采樣每個時間步都寫。解決方法是在輸出模塊里加一個輸出間隔比如每 24 個內(nèi)模式步寫一次或者只輸出表層和底層。如果已經(jīng)寫滿先清理舊文件再改輸出頻率重跑。5. 進階技巧用理想實驗驗證配置再上真實海域5.1 用理想地形做冷啟動驗證真實海域配置復雜直接上容易把配置錯誤和物理問題混在一起。我習慣先用一個理想化的矩形海域做冷啟動地形是平底或者簡單斜坡開邊界給一個正弦潮汐初始溫鹽是常數(shù)。跑一天后檢查表面水位是否和開邊界同相位流速是否對稱。如果理想實驗都跑不通真實海域更不用談。import numpy as np import matplotlib.pyplot as plt # 讀取理想實驗的表面水位輸出 ssh np.loadtxt(ideal_ssh.txt) time np.arange(len(ssh)) * 3600 # 假設每小時輸出一次 # 畫水位時間序列檢查是否和開邊界同相位 plt.plot(time / 3600, ssh, labelmodel SSH) plt.axhline(0, colorgray, linestyle--) plt.xlabel(Time (hour)) plt.ylabel(SSH (m)) plt.legend() plt.savefig(ssh_check.png, dpi150)這段代碼用來快速檢查水位相位。如果模型水位和開邊界相位差超過半個周期說明開邊界輻射條件有問題需要檢查輻射速度參數(shù)。參數(shù)上輻射速度一般取當?shù)刂亓Σㄋ贉\水區(qū)約等于sqrt(g * H)。5.2 用實測數(shù)據(jù)做溫鹽剖面校驗理想實驗通過后換真實地形和實測開邊界。跑完一個潮周期后把模擬的溫鹽剖面和實測 CTD 數(shù)據(jù)對比。常見做法是取同一位置同一時刻的剖面畫在同一張圖上。如果表層溫度偏差超過 1 度先查熱通量文件的時間插值如果底層鹽度偏差超過 0.5查底摩擦和垂向混合系數(shù)。校驗項可接受偏差優(yōu)先排查表面水位0.1 m開邊界、底摩擦表層溫度1.0 度熱通量、短波穿透底層鹽度0.5垂向混合、底摩擦表層流速0.1 m/s風場、開邊界5.3 長期運行的穩(wěn)定性習慣跑長期模擬時我一般會先跑 10 天檢查能量是否守恒再決定是否延長。能量檢查的方法是看動能和勢能的日變化是否合理如果總能量持續(xù)增長說明有數(shù)值不穩(wěn)定。另外每次改配置后都從冷啟動重跑不要熱啟動疊加否則錯誤會累積。從那以后我每次改完參數(shù)都強制走一遍理想實驗確認沒問題再上真實海域。希望幫到你。本文還有配套的精品資源點擊獲取 POM普林斯頓海洋模型實戰(zhàn)：從零搭建近岸海域模擬與避坑指南
簡介POM普林斯頓海洋模型是一套面向海洋數(shù)值模擬初學者與科研人員的開源建模工具可在f77、f90、g77、pgf77、CVF等多種編譯器下直接運行適合用于河道水流、河口及近岸海域的水動力與溫鹽場模擬研究。資源包共21個文件約658KB包含f與c源碼、n與dat輸入數(shù)據(jù)、out運行輸出、runpom2k執(zhí)行腳本、change版本記錄、pdf中文說明書及zip/tar工具包等覆蓋從編譯、算例運行到結果查看的完整流程。POM2K提供兩個可直接上手的實例河道模型參數(shù)已內(nèi)置于代碼中河口模型則附帶原始網(wǎng)格文件其輸入文件由GRID-DATA/GRID.f生成便于用戶理解網(wǎng)格與初始場的構造方式。借助GRID-DATA中的插值代碼可將地形、溫鹽、表面風場及熱通量等資料插值到網(wǎng)格點生成初始條件與強迫條件為自建模型提供有效工具。目前已有997人學習下載適合希望快速搭建POM建模環(huán)境、掌握算例配置與數(shù)據(jù)前處理思路的讀者參考。1. POM 普林斯頓海洋模型從零建立一套可跑通的海域模擬如果你手頭只有一份岸線數(shù)據(jù)、幾組風場和開邊界水位卻要在一周內(nèi)搭出一套能出流場和溫鹽結構的近岸模型POMPrinceton Ocean Model大概率是你繞不開的選擇。它屬于三維原始方程海洋模型采用 sigma 坐標貼合地形垂向用二階湍流閉合時間積分走的是外模式與內(nèi)模式分離的顯式差分。這套結構決定了它對近岸、河口、陸架海域的適應性很強也決定了它對時間步長和網(wǎng)格質(zhì)量相當敏感。我見過不少人拿到源碼包后卡在編譯和參數(shù)配置上其實真正難的不是代碼本身而是搞清楚每個輸入文件對應哪塊物理過程。這篇筆記按“資源是什么、怎么配、坑在哪”的順序拆一遍適合做近岸水動力、溫鹽輸運、工程海域論證的從業(yè)者也適合剛接觸海洋數(shù)值模擬的研究生照著復現(xiàn)。2. POM 的坐標系與差分框架為什么它適合近岸而不是深海2.1 sigma 坐標與內(nèi)外模式分離的取舍POM 最核心的設計是垂向 sigma 坐標也就是把水深歸一化到 0 到 1 之間自由表面和底部地形都能被貼合。這樣做的好處是近岸淺水區(qū)的垂向分辨率不會被水深變化吃掉河口和陸架坡折帶都能保持同樣的垂向?qū)訑?shù)。代價是斜壓梯度項在陡坡處容易產(chǎn)生誤差所以實際配置時地形坡度一般控制在 0.05 以內(nèi)超過這個值要么平滑地形要么加密水平網(wǎng)格。時間積分上POM 把正壓模態(tài)和斜壓模態(tài)拆開。外模式負責自由表面和垂向平均流時間步長受 CFL 條件限制通常取幾秒到十幾秒內(nèi)模式負責溫鹽和三維流場步長可以放大到外模式的幾十倍。這種分離讓計算量降下來但也帶來一個隱患內(nèi)外模式之間的耦合頻率如果設得不合理會出現(xiàn)質(zhì)量不守恒或者表面重力波失真。常見做法是讓內(nèi)外模式步長比控制在 20 到 50 之間具體看網(wǎng)格分辨率和最大水深。2.2 控制方程與湍流閉合的落地含義POM 求解的是 Boussinesq 近似下的三維原始方程包括動量方程、連續(xù)方程、溫度鹽度輸運方程和狀態(tài)方程。動量方程里保留了科氏力、壓強梯度力、湍流黏性項和底摩擦。湍流閉合用的是 Mellor-Yamada 二階閉合方案也就是常說的 2.5 階模型它通過求解湍動能和湍流宏觀尺度兩個預報量來給出垂向混合系數(shù)。對使用者來說這套閉合方案意味著你不需要自己調(diào)混合層深度但需要關注兩個參數(shù)一個是湍動能產(chǎn)生項里的常數(shù)另一個是壁面函數(shù)里的粗糙度。前者在源碼里通常以常量形式出現(xiàn)不建議隨意改后者和底摩擦系數(shù)直接相關近岸模擬里一般取 0.002 到 0.005 之間。如果模擬出來的底層流速偏大先查底摩擦再查垂向分層是否被過度混合。2.3 網(wǎng)格與地形文件的準備步驟POM 的輸入網(wǎng)格是正交曲線網(wǎng)格水平方向用 Arakawa C 網(wǎng)格離散也就是速度分量分別位于網(wǎng)格邊中點標量位于網(wǎng)格中心。準備網(wǎng)格時第一步是確定研究區(qū)域和開邊界位置第二步是生成正交曲線網(wǎng)格第三步是把水深插值到網(wǎng)格點上并做平滑。常見做法是用正交曲線網(wǎng)格生成工具先出網(wǎng)格文件再導出 POM 需要的 fortran 格式。水深數(shù)據(jù)一般來自實測或公開海圖插值后要做三次平滑每次平滑系數(shù)不超過 0.5否則會抹掉真實地形特征。下面是一個把水深數(shù)據(jù)插值到 POM 網(wǎng)格并平滑的 Python 片段用的是常見科學計算庫。import numpy as np from scipy.interpolate import griddata # 讀取原始水深散點lon, lat, depth raw np.loadtxt(depth_raw.txt) lon, lat, depth raw[:, 0], raw[:, 1], raw[:, 2] # 讀取 POM 網(wǎng)格點坐標 grid np.loadtxt(pom_grid.txt) gx, gy grid[:, 0], grid[:, 1] # 插值到網(wǎng)格點近岸用線性外?？蓳Q cubic depth_grid griddata((lon, lat), depth, (gx, gy), methodlinear) # 填補插值空洞避免陸地邊界出現(xiàn) NaN mask np.isnan(depth_grid) if mask.any(): depth_grid[mask] griddata( (lon, lat), depth, (gx[mask], gy[mask]), methodnearest ) # 三次平滑系數(shù) 0.3保留地形趨勢 smoothed depth_grid.copy() for _ in range(3): smoothed[1:-1, 1:-1] ( smoothed[1:-1, 1:-1] * 0.4 (smoothed[:-2, 1:-1] smoothed[2:, 1:-1] smoothed[1:-1, :-2] smoothed[1:-1, 2:]) * 0.15 ) np.savetxt(depth_pom.txt, smoothed, fmt%.3f)這段代碼的邏輯是先做散點插值再處理邊界空洞最后做三次加權平滑。參數(shù)上平滑系數(shù) 0.3 是經(jīng)驗值地形復雜的河口可以降到 0.2開闊陸架可以升到 0.4。插值方法選 linear 是為了避免 cubic 在數(shù)據(jù)稀疏區(qū)產(chǎn)生過沖如果原始水深點足夠密可以換成 cubic 提升光滑度。輸出格式按 POM 讀取要求寫成純文本每行一個網(wǎng)格點的水深值。3. 編譯與輸入文件配置把源碼跑起來的關鍵幾步3.1 編譯環(huán)境與 Makefile 調(diào)整POM 源碼通常是 Fortran 77 和 Fortran 90 混合編譯前先確認編譯器支持固定格式和自由格式混編。常見做法是用 gfortran在 Makefile 里指定-ffixed-line-length-132和-fno-second-underscore避免長行被截斷和符號名不匹配。如果用的是 ifort對應選項是-extend-source和-assume nounderscore。編譯順序一般是先編公共塊和參數(shù)模塊再編主程序。源碼包里通常有一個makefile或Makefile里面已經(jīng)寫好了依賴關系但路徑需要按本機環(huán)境改。我一般會先跑一遍make clean再make看第一個報錯出現(xiàn)在哪個文件。如果是找不到netcdf.inc說明輸出用了 NetCDF需要先裝好庫并在 Makefile 里加-I和-L路徑。# 清理舊目標文件 make clean # 指定編譯器和 NetCDF 路徑后編譯 make FCgfortran \ FFLAGS-O2 -ffixed-line-length-132 -fno-second-underscore \ NETCDF_INC/usr/local/include \ NETCDF_LIB/usr/local/lib -lnetcdf這里-O2是常規(guī)優(yōu)化調(diào)試階段可以換成-g -O0方便定位數(shù)組越界。NETCDF_INC和NETCDF_LIB按實際安裝路徑改如果不用 NetCDF 輸出可以把相關源文件從 Makefile 里注釋掉。編譯成功后會在當前目錄生成可執(zhí)行文件名字通常是pom.exe或pom2k。3.2 輸入文件清單與參數(shù)對應關系POM 的輸入文件一般包括網(wǎng)格文件、水深文件、初始溫鹽文件、開邊界文件、風場和熱通量文件。每個文件對應一個 Fortran 讀取單元單元號在源碼的params或infile里定義。下面這張表列出常見文件和關鍵參數(shù)方便對照檢查。文件名內(nèi)容關鍵參數(shù)常見取值grid.dat網(wǎng)格坐標網(wǎng)格點數(shù)、正交性由生成工具決定depth.dat水深最小水深不小于 2 minit_ts.dat初始溫鹽表層溫度、鹽度按季節(jié)實測obc.dat開邊界水位、流速、溫鹽潮汐或?qū)崪ywind.dat風場風速、風向10 m 高度heat.dat熱通量短波、長波、感熱按月份配置時最容易出錯的是開邊界。POM 支持輻射邊界和指定邊界兩種輻射邊界適合外海開邊界指定邊界適合有實測潮位的斷面。如果開邊界給的是潮汐調(diào)和常數(shù)需要在源碼里加一個潮汐預報模塊或者提前把水位時間序列算好寫進obc.dat。風場和熱通量的時間步長要和模型時間步長對齊否則會插值出錯。3.3 時間步長與 CFL 條件的實際估算時間步長的選取直接決定模型能不能穩(wěn)定跑完。外模式步長受 CFL 條件限制公式是dt dx / sqrt(g * H_max)其中dx是最小網(wǎng)格間距H_max是最大水深g是重力加速度。舉個例子最小網(wǎng)格 500 米最大水深 100 米算出來dt要小于 15 秒左右。實際配置時取 10 秒留余量。內(nèi)模式步長通常是外模式的 20 到 50 倍但不要超過 60 倍否則斜壓調(diào)整跟不上。如果模擬區(qū)域有強潮汐內(nèi)外步長比取小一點比如 20如果是陸架環(huán)流可以取到 40。跑之前先用理想地形做一次冷啟動看自由表面是否發(fā)散再換真實地形。# 冷啟動檢查先用 10 秒外步長跑 1 天 ./pom.exe run_control.txt # 查看輸出日志里的 CFL 數(shù)和表面水位極值 grep -E CFL|SSH pom.log | tail -20如果日志里 CFL 數(shù)持續(xù)大于 1說明步長偏大或者地形有尖峰需要減小步長或平滑地形。表面水位極值如果出現(xiàn)異常大值先查開邊界是否給了不合理的水位再查底摩擦是否過小。4. 避坑與排查跑 POM 時最容易翻車的五個地方4.1 編譯通過但運行報數(shù)組越界現(xiàn)象是編譯沒問題一運行就Segmentation fault或者Fortran runtime error: Index out of bounds。原因通常是網(wǎng)格文件里的點數(shù)比源碼里聲明的數(shù)組維度大或者水深文件的行數(shù)和網(wǎng)格點數(shù)不一致。解決方法是先檢查params文件里的IM和JM是否和網(wǎng)格文件頭一致再逐行核對水深文件的行數(shù)。我一般會寫一個簡單的 Python 腳本統(tǒng)計行數(shù)和網(wǎng)格文件對比。4.2 自由表面發(fā)散或水位異?，F(xiàn)象是跑了幾十步后表面水位指數(shù)增長日志里 CFL 數(shù)越來越大。原因可能是外模式步長超過 CFL 限制或者開邊界水位給了不合理的突變。解決方法是先把外步長減半試跑如果穩(wěn)定了說明是步長問題如果還發(fā)散檢查開邊界文件的時間序列是否有跳變必要時做線性插值平滑。4.3 溫鹽出現(xiàn)負值或超過物理范圍現(xiàn)象是溫鹽場跑著跑著出現(xiàn)負鹽度或者溫度超過 40 度。原因通常是平流項用了中心差分在強梯度區(qū)產(chǎn)生數(shù)值振蕩。POM 里可以選擇迎風差分或者加水平擴散。解決方法是把平流方案改成迎風或者把水平擴散系數(shù)調(diào)大但不要超過 100 m2/s否則會過度抹平鋒面。4.4 底摩擦系數(shù)取錯導致底層流速失真現(xiàn)象是表層流速合理底層流速偏大或偏小。原因是底摩擦系數(shù)沒有按底質(zhì)類型調(diào)整。沙質(zhì)海底取 0.002 到 0.003泥質(zhì)海底取 0.004 到 0.005。如果模擬區(qū)域底質(zhì)不均可以按區(qū)域給不同系數(shù)但要在源碼里加判斷。我一般會先用一個中間值跑一遍再和實測流速對比微調(diào)。4.5 輸出文件太大導致磁盤寫滿現(xiàn)象是跑了一半程序卡住日志沒有報錯但輸出目錄不再更新。原因是 NetCDF 或二進制輸出沒有做時間降采樣每個時間步都寫。解決方法是在輸出模塊里加一個輸出間隔比如每 24 個內(nèi)模式步寫一次或者只輸出表層和底層。如果已經(jīng)寫滿先清理舊文件再改輸出頻率重跑。5. 進階技巧用理想實驗驗證配置再上真實海域5.1 用理想地形做冷啟動驗證真實海域配置復雜直接上容易把配置錯誤和物理問題混在一起。我習慣先用一個理想化的矩形海域做冷啟動地形是平底或者簡單斜坡開邊界給一個正弦潮汐初始溫鹽是常數(shù)。跑一天后檢查表面水位是否和開邊界同相位流速是否對稱。如果理想實驗都跑不通真實海域更不用談。import numpy as np import matplotlib.pyplot as plt # 讀取理想實驗的表面水位輸出 ssh np.loadtxt(ideal_ssh.txt) time np.arange(len(ssh)) * 3600 # 假設每小時輸出一次 # 畫水位時間序列檢查是否和開邊界同相位 plt.plot(time / 3600, ssh, labelmodel SSH) plt.axhline(0, colorgray, linestyle--) plt.xlabel(Time (hour)) plt.ylabel(SSH (m)) plt.legend() plt.savefig(ssh_check.png, dpi150)這段代碼用來快速檢查水位相位。如果模型水位和開邊界相位差超過半個周期說明開邊界輻射條件有問題需要檢查輻射速度參數(shù)。參數(shù)上輻射速度一般取當?shù)刂亓Σㄋ贉\水區(qū)約等于sqrt(g * H)。5.2 用實測數(shù)據(jù)做溫鹽剖面校驗理想實驗通過后換真實地形和實測開邊界。跑完一個潮周期后把模擬的溫鹽剖面和實測 CTD 數(shù)據(jù)對比。常見做法是取同一位置同一時刻的剖面畫在同一張圖上。如果表層溫度偏差超過 1 度先查熱通量文件的時間插值如果底層鹽度偏差超過 0.5查底摩擦和垂向混合系數(shù)。校驗項可接受偏差優(yōu)先排查表面水位0.1 m開邊界、底摩擦表層溫度1.0 度熱通量、短波穿透底層鹽度0.5垂向混合、底摩擦表層流速0.1 m/s風場、開邊界5.3 長期運行的穩(wěn)定性習慣跑長期模擬時我一般會先跑 10 天檢查能量是否守恒再決定是否延長。能量檢查的方法是看動能和勢能的日變化是否合理如果總能量持續(xù)增長說明有數(shù)值不穩(wěn)定。另外每次改配置后都從冷啟動重跑不要熱啟動疊加否則錯誤會累積。從那以后我每次改完參數(shù)都強制走一遍理想實驗確認沒問題再上真實海域。希望幫到你。本文還有配套的精品資源點擊獲取