因應intel MPI轉換wrfout格式

Table of contents

背景
程式設計
呼叫與執行結果
- 呼叫
- 結果

背景

作為總結國網上wrfout與mcip之間的衝突，這篇的做法並不是非常光彩(但真的是沒辦法中的辦法了)。就是：
1. 找一個國網mcip程式可以接受的wrfout模版(一般工作站所產生的netcdf)，
2. 將國網上跑的wrfout結果(pnetcdf)讀出來
3. 將內容填進模板
詳細的問題補充於平行與序列化libnetcdf.settings的差異。
憑藉的理由(python處理的強項)：
- python對於檔案是否是一般netcdf還是pnetcdf所產生並不敏感（相對fortran程式）。
- python處理上GB的wrfout檔案，並沒有延遲或龜速的情況。
- python程式設計要比準備fortran程式庫、修改平行運作的程式錯誤，要容易的多了。
流程方案

程式	版本	編譯	說明
wrf.exe	4.0/4.2.1等等不拘	國網上的intel mpi	必須是pnetcdf版本、國網上不接受序列版本netcdf的IO方式
TransWrfout	3.11	python	將pnetcdf之wrfout結果填入指定序列netcdf模版(本程式)
mcip	5.0/5.1	ifort	無mpi之序列版本。不能使用平行作業，詳見mcip平行運作方案

目前環保署公版模式並未提供mcip程式與編譯建議。國網官方測試也只限制在wrf/cmaq的主程式，都未有詳細交代mcip的執行細節。
可將序列處理、序列mcip的過程放在背景執行，前景可以繼續執行其他domain的real/wrf，因此即使是序列執行，也不會影響預報的進度或造成速度瓶頸。

程式設計

模版的準備

去除Time維度，以減少檔案容量
分別準備grid45/grid09/grid03等等3個範圍
模板和真正執行的wrf之間，必須有相同的版次。
- 檢視：ncdump -h $nc|grep WRF

陣列變數

先恢復時間維度的長度
按照不同階層數的維度進行數據的移轉
如果陣列階層數相同，但維度的長度不同，則不處理。

nc = netCDF4.Dataset(fname,'r+')
nc0= netCDF4.Dataset(fnam0,'r')
V =[list(filter(lambda x:nc.variables[x].ndim==j, [i for i in nc.variables])) for j in [1,2,3,4]]
V0=[list(filter(lambda x:nc0.variables[x].ndim==j, [i for i in nc0.variables])) for j in [1,2,3,4]]

v='V10'
nt,nrow,ncol=nc0[v].shape
print(nt)
for t in range(nt):
  nc[v][t,:,:]=nc0[v][t,:,:]
for v in V[3]:
  nc[v][:,:,:,:]=nc0[v][:,:,:,:]
for v in V[2]:
  if v not in V0[2]:continue
  nc[v][:,:,:]=nc0[v][:,:,:]
for v in V[1]:
  if nc[v].shape!=nc0[v].shape:continue
  nc[v][:,:]=nc0[v][:,:]
for v in V[0]:
  nc[v][:]=nc0[v][:]

屬性轉移

事先使用ncdump讀出檔案的屬性變數項目備用

sinotec2@lgn303 /work/sinotec2/cmaqruns/forecast/grid45/wrfout
$ head *.txt
                :TITLE = " OUTPUT FROM WRF V4.2 MODEL" ;
                :START_DATE = "2023-03-08_00:00:00" ;
                :SIMULATION_START_DATE = "2023-03-08_00:00:00" ;
                :WEST-EAST_GRID_DIMENSION = 221 ;
                :SOUTH-NORTH_GRID_DIMENSION = 129 ;
                :BOTTOM-TOP_GRID_DIMENSION = 25 ;
                :DX = 45000.f ;
                :DY = 45000.f ;
                :AERCU_OPT = 0 ;
                :AERCU_FCT = 1.f ;
...

有減號的變數要規避，會造成執行錯誤
新版有restart的內容，無關mcip之應用

with open('att.txt', 'r') as f:
  var=[i.split(':')[1].split('=')[0].replace(' ','') for i in f if len(i)>0]
for v in var:
  if '-' in v:continue
  try:
    exec('nc.'+v+'=nc0.'+v)
  except:
    continue

檔案管理

引數：$BEGD
取代之前從$gfs目錄連結到$fcst之作法
讀取檔案目錄名稱
- 目錄：即執行wrf的目錄(pwd)
- 檔名：逐日標記(dates)
輸出目標targ
- targ=fcst+'/grid45/wrfout'
- mcip使用的$InMetFiles仍為固定檔名之暫存檔
  - nam0=pwd+'/wrfout_d0'+ad0+'_'+dates[jj]+'_00:00:00'
以下範例為台灣杉3號

tdy=sys.argv[1]
bdate=datetime.datetime.strptime(tdy,"%Y-%m-%d")
nd=12
dates=[(bdate+datetime.timedelta(days=i)).strftime("%Y-%m-%d") for i in range(nd)]
...
pwd=subprocess.check_output("pwd" ,shell=True).decode('utf8').strip('\n')
fcst='/work/sinotec2/cmaqruns/forecast'
targ=fcst+'/grid45/wrfout'

ads={1:['3'],2:['1','2']}
gds={1:['03'],2:['45','09']}
ndms=1
if pwd[-3:]=='45k':
  ndms=2
with open(targ+'/att.txt', 'r') as f:
  var=[i.split(':')[1].split('=')[0].replace(' ','') for i in f if len(i)>0]
for ii in range(ndms):
  ad=ads[ndms][ii]
  ad0=ad
  if ndms==1:ad0='1'
  ftemp=targ+'/wrfout_d0'+ad+'_template'
  for jj in range(nd):
    fnam0=pwd+'/wrfout_d0'+ad0+'_'+dates[jj]+'_00:00:00'
    if not os.path.isfile(fnam0):continue
    nc0= netCDF4.Dataset(fnam0,'r')
    fname=targ+'/wrfout_d0'+ad+'_'+str(jj)
    os.system('test -e '+fname+' && rm -f '+fname+';cp '+ftemp+' '+fname)
    nc = netCDF4.Dataset(fname,'r+')

序列版mcip的執行

for ii in range(ndms):
...    
  gd=gds[ndms][ii]
  os.system('/usr/bin/sbatch '+fcst+'/csh_mcip.sh '+gd)

呼叫與執行結果

呼叫

在$gfs/fcst.cs中，執行wrf.exe後隨即在背景中執行

...
  # wrf
    mpiexec.hydra -bootstrap slurm -n 120 ../run/wrf.exe
    $BIN/sub $BIN/trans_wrfout.py $BEGD
...

結果

用最簡便的方式解決。
wrfout檔案會增加至少1/4

Download: 轉換wrfout格式：trans_wrfout.py