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