GrADS　GSM気圧面データ　取得＆描画　バイナリファイル編

GSMデータを用いて、GrADSで描画したいと思います。

GSMデータとは、なんぞやという方は、こちらを参照してください。

2021年10月16日 GrADS　GSMのデータを取得＆描画　NetCDF編

GSMデータの取得

本記事では、京都大学の気象庁データを利用させていただきます。

オリジナルのGSMデータは、国際気象通報式FM92 GRIB 二進形式格子点資料気象通報式(第2版)というデータ形式になっています。

本記事では、そのオリジナルデータ(GRIB2形式)から描画させたいと思います。

Z__C_RJTD_yyyyMMddhhmmss_GSM_GPV_Rjp_Lsurf_FDddhh-ddhh_grib2.binは、地表面データです。

Z__C_RJTD_yyyyMMddhhmmss_GSM_GPV_Rjp_L-pall_FDddhh-ddhh_grib2.binは、気圧面データです。

本記事では、気圧面データに着目します。

yyyyは西暦、MMは月、ddは日、hhは時間、㎜は分、ssは秒です。

wgetコマンドで取得しましょう。

$ wget  http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/2021/07/17/Z__C_RJTD_20210717000000_GSM_GPV_Rjp_L-pall_FD0000-0312_grib2.bin

上では、2021年7月17日の00UTC初期値の気圧面データであるGSMデータを取得しています。予測時間は、84時間後までです。

264時間先までほしいので、合計で3つのファイルを取得します。また、初期値が06,18時の場合、132時間先までですので、合計のファイルは、2つです。

$ wget  http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/2021/07/17/Z__C_RJTD_20210717000000_GSM_GPV_Rjp_L-pall_FD0318-0512_grib2.bin
$ wget  http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/2021/07/17/Z__C_RJTD_20210717000000_GSM_GPV_Rjp_L-pall_FD0518-1100_grib2.bin

wgrib2コマンドでバイナリファイルに変換します。

$ wgrib2 Z__C_RJTD_20210717000000_GSM_GPV_Rjp_L-pall_FD0000-0312_grib2.bin -no_header -bin out00000312.bin
$ wgrib2 Z__C_RJTD_20210717000000_GSM_GPV_Rjp_L-pall_FD0318-0512_grib2.bin -no_header -bin out03180512.bin
$ wgrib2 Z__C_RJTD_20210717000000_GSM_GPV_Rjp_L-pall_FD0518-1100_grib2.bin -no_header -bin out05181100.bin

バイナリファイルの編集

この時のバイナリファイルout????.binは、1000hPaのジオポテンシャル高度、東西風、南北風、気温、鉛直流、相対湿度、975hPaのジオポテンシャル高度、東西風、南北風、….100hPaの上昇流という並び方をしています。

このままだと、GrADSで描画する際、ある高度のジオポテンシャル高度を見たいのに気温が表示されたりとめちゃくちゃになるので、それぞれの変数でまとめるという感じに変更する必要があります。

また、264時間先まで3つのファイルに分かれているので、1つのファイルに統合しようと思います。

Fortranでやります。Pythonもいつかやろうと思います。

バイナリファイルの大きさを確認します。

$ ls -l
-rw-r--r-- 1 tanaka users 194988112 Aug 16 16:46 out00000312.bin
-rw-r--r-- 1 tanaka users  53789824 Aug 16 16:46 out03180512.bin
-rw-r--r-- 1 tanaka users 147922016 Aug 16 16:46 out05121100.bin

out00000312.binのファイルサイズは194988112byte、out03180512.binでは、53789824byte、out05121100.binでは、147922016byteです。

気象業務センターで確認すると、GSMデータは、北緯20~50度と東経120~150度で、気圧面データでは、x方向(経度)格子点間隔0.25の格子点数121、y方向(緯度)格子間隔0.2度の格子点数151です。

z方向（1000hPa・975hPa・950hPa・925hPa・900hPa・850hPa・800hPa・700hPa・600hPa・500hPa・400hPa・300hPa・250hPa・200hPa・150hPa・100hPa）の16個あります。

物理量として、ジオポテンシャル高度、東西風、南北風、気温、鉛直流、相対湿度(~300hPaまで)があります。

out00000312.binは、時間間隔3時間、84時間先までのファイルですので、t方向に29です。

1データが4byte、物理量が6つありますが、相対湿度は300hPaまでしかありません。

4×5×151×121×16×29=169,554,880

4×1×151×121×12×29=25,433,232

169,554,880+25,433,232=194,988,112

同様にout03180512.binでは、時間間隔6時間、87~132時間先までですので、t方向に8です。

4×5×151×121×16×8=46,773,760

4×1×151×121×12×8=7,016,064

46,773,760+ 7,016,064 =53,789,824

同様にout05181100.binでは、時間間隔6時間、135~264時間先までですので、t方向に22です。

4×5×151×121×16×22=128,627,840

4×1×151×121×12×22=19,294,176

128,627,840+19,294,176=147,922,016

想定したファイルサイズになったので、プログラムを書きます。

implicit none
real,allocatable::out0312(:),out0512(:),out1100(:)
real,allocatable::gz0312(:,:,:,:),uu0312(:,:,:,:),vv0312(:,:,:,:),temp0312(:,:,:,:),ww0312(:,:,:,:),rh0312(:,:,:,:)
real,allocatable::gz0512(:,:,:,:),uu0512(:,:,:,:),vv0512(:,:,:,:),temp0512(:,:,:,:),ww0512(:,:,:,:),rh0512(:,:,:,:)
real,allocatable::gz1100(:,:,:,:),uu1100(:,:,:,:),vv1100(:,:,:,:),temp1100(:,:,:,:),ww1100(:,:,:,:),rh1100(:,:,:,:)
real,allocatable::gz(:,:,:,:),uu(:,:,:,:),vv(:,:,:,:),temp(:,:,:,:),ww(:,:,:,:),rh(:,:,:,:)

integer::xx,yy,zz,tt,num,tdata,rnum
integer,parameter::xxdim=121,yydim=151,zzdim=16,humdim=12
integer,parameter::out0312ttdim=29,out0512ttdim=8,out1100ttdim=22,ttdim=59
integer::out0312recl,out0512recl,out1100recl
out0312recl = 4 * 5 * xxdim * yydim * zzdim * out0312ttdim + 4 * 1 * xxdim * yydim * humdim * out0312ttdim
out0512recl = 4 * 5 * xxdim * yydim * zzdim * out0512ttdim + 4 * 1 * xxdim * yydim * humdim * out0512ttdim
out1100recl = 4 * 5 * xxdim * yydim * zzdim * out1100ttdim + 4 * 1 * xxdim * yydim * humdim * out1100ttdim

allocate(gz0312(xxdim,yydim,zzdim,out0312ttdim),uu0312(xxdim,yydim,zzdim,out0312ttdim),vv0312(xxdim,yydim,zzdim,out0312ttdim),&
temp0312(xxdim,yydim,zzdim,out0312ttdim),ww0312(xxdim,yydim,zzdim,out0312ttdim),rh0312(xxdim,yydim,zzdim,out0312ttdim))
allocate(gz0512(xxdim,yydim,zzdim,out0512ttdim),uu0512(xxdim,yydim,zzdim,out0512ttdim),vv0512(xxdim,yydim,zzdim,out0512ttdim),&
temp0512(xxdim,yydim,zzdim,out0512ttdim),ww0512(xxdim,yydim,zzdim,out0512ttdim),rh0512(xxdim,yydim,zzdim,out0512ttdim))
allocate(gz1100(xxdim,yydim,zzdim,out1100ttdim),uu1100(xxdim,yydim,zzdim,out1100ttdim),vv1100(xxdim,yydim,zzdim,out1100ttdim),&
temp1100(xxdim,yydim,zzdim,out1100ttdim),ww1100(xxdim,yydim,zzdim,out1100ttdim),rh1100(xxdim,yydim,zzdim,out1100ttdim))
allocate(gz(xxdim,yydim,zzdim,ttdim),uu(xxdim,yydim,zzdim,ttdim),vv(xxdim,yydim,zzdim,ttdim),temp(xxdim,yydim,zzdim,ttdim),&
ww(xxdim,yydim,zzdim,ttdim),rh(xxdim,yydim,zzdim,ttdim))
allocate(out0312(out0312recl/4),out0512(out0512recl/4),out1100(out1100recl/4))
open(10,file="out00000312.bin",access="direct",recl=out0312recl)
open(11,file="out03180512.bin",access="direct",recl=out0512recl)
open(12,file="out05181100.bin",access="direct",recl=out1100recl)
read(10,rec=1)out0312(:)
read(11,rec=1)out0512(:)
read(12,rec=1)out1100(:)
close(10)
close(11)
close(12)
!read
num=1
do tt=1,out0312ttdim
do zz=1,humdim
     call sub0312(gz0312)
     call sub0312(uu0312)
     call sub0312(vv0312)
     call sub0312(temp0312)
     call sub0312(ww0312)
     call sub0312(rh0312)
enddo
do zz=humdim+1,zzdim
     call sub0312(gz0312)
     call sub0312(uu0312)
     call sub0312(vv0312)
     call sub0312(temp0312)
     call sub0312(ww0312)
     do yy=1,yydim
     do xx=1,xxdim
        rh0312(xx,yy,zz,tt)=0.0
     enddo
     enddo
enddo
enddo
num=1
do tt=1,out0512ttdim
do zz=1,humdim
     call sub0512(gz0512)
     call sub0512(uu0512)
     call sub0512(vv0512)
     call sub0512(temp0512)
     call sub0512(ww0512)
     call sub0512(rh0512)
enddo
do zz=humdim+1,zzdim
     call sub0512(gz0512)
     call sub0512(uu0512)
     call sub0512(vv0512)
     call sub0512(temp0512)
     call sub0512(ww0512)
     do yy=1,yydim
     do xx=1,xxdim
        rh0512(xx,yy,zz,tt)=0.0
     enddo
     enddo
enddo
enddo

num=1
do tt=1,out1100ttdim
do zz=1,humdim
     call sub1100(gz1100)
     call sub1100(uu1100)
     call sub1100(vv1100)
     call sub1100(temp1100)
     call sub1100(ww1100)
     call sub1100(rh1100)
enddo
do zz=humdim+1,zzdim
     call sub1100(gz1100)
     call sub1100(uu1100)
     call sub1100(vv1100)
     call sub1100(temp1100)
     call sub1100(ww1100)
     do yy=1,yydim
     do xx=1,xxdim
        rh1100(xx,yy,zz,tt)=0.0
     enddo
     enddo
enddo
enddo

!write
!write(*,*)gz1100(xxdim,yydim,zzdim,:)

do tdata=1,ttdim
do zz=1,zzdim
do yy=1,yydim
do xx=1,xxdim
   if (tdata <= out0312ttdim)then
      gz(xx,yy,zz,tdata)=gz0312(xx,yy,zz,tdata)
      uu(xx,yy,zz,tdata)=uu0312(xx,yy,zz,tdata)
      vv(xx,yy,zz,tdata)=vv0312(xx,yy,zz,tdata)
      temp(xx,yy,zz,tdata)=temp0312(xx,yy,zz,tdata)
      ww(xx,yy,zz,tdata)=ww0312(xx,yy,zz,tdata)
      rh(xx,yy,zz,tdata)=rh0312(xx,yy,zz,tdata)

   elseif(tdata <= out0312ttdim+out0512ttdim)then
      tt=tdata - out0312ttdim
      gz(xx,yy,zz,tdata)=gz0512(xx,yy,zz,tt)
      uu(xx,yy,zz,tdata)=uu0512(xx,yy,zz,tt)
      vv(xx,yy,zz,tdata)=vv0512(xx,yy,zz,tt)
      temp(xx,yy,zz,tdata)=temp0512(xx,yy,zz,tt)
      ww(xx,yy,zz,tdata)=ww0512(xx,yy,zz,tt)
      rh(xx,yy,zz,tdata)=rh0512(xx,yy,zz,tt)
   else
      tt=tdata - out0312ttdim - out0512ttdim
      gz(xx,yy,zz,tdata)=gz1100(xx,yy,zz,tt)
      uu(xx,yy,zz,tdata)=uu1100(xx,yy,zz,tt)
      vv(xx,yy,zz,tdata)=vv1100(xx,yy,zz,tt)
      temp(xx,yy,zz,tdata)=temp1100(xx,yy,zz,tt)
      ww(xx,yy,zz,tdata)=ww1100(xx,yy,zz,tt)
      rh(xx,yy,zz,tdata)=rh1100(xx,yy,zz,tt)
   endif
enddo
enddo
enddo
enddo
!write(*,*)ww(:,:,:,:)
open(20,file="gsm.bin",status="unknown",access="direct",form="unformatted",recl=4*xxdim*yydim*zzdim)
rnum=1
do tt=1,ttdim
   write(20,rec=rnum)gz(:,:,:,tt)
   rnum=rnum+1
   write(20,rec=rnum)uu(:,:,:,tt)
   rnum=rnum+1
   write(20,rec=rnum)vv(:,:,:,tt)
   rnum=rnum+1
   write(20,rec=rnum)temp(:,:,:,tt)
   rnum=rnum+1
   write(20,rec=rnum)ww(:,:,:,tt)
   rnum=rnum+1
   write(20,rec=rnum)rh(:,:,:,tt)
   rnum=rnum+1
enddo
close(20)

stop

contains

subroutine sub0312(var)
  implicit none
  real::var(xxdim,yydim,zzdim,out0312ttdim)
  do yy=1,yydim
  do xx=1,xxdim
        var(xx,yy,zz,tt)=out0312(num)
        num=num+1
  enddo
  enddo
end subroutine 
subroutine sub0512(var)
  implicit none
  real::var(xxdim,yydim,zzdim,out0512ttdim)
  do yy=1,yydim
  do xx=1,xxdim
        var(xx,yy,zz,tt)=out0512(num)
        num=num+1
  enddo
  enddo
end subroutine 
subroutine sub1100(var)
  implicit none
  real::var(xxdim,yydim,zzdim,out1100ttdim)
  do yy=1,yydim
  do xx=1,xxdim
        var(xx,yy,zz,tt)=out1100(num)
        num=num+1
  enddo
  enddo
end subroutine 
  
end

コントロールファイルの作成

GrADSで描画するために、コントロールファイルを作る必要があります。

気象業務センターで確認すると、日本域のGSMデータは、日本とその近海(北緯20度～50度、東経120度〜150度)で、x方向(経度)格子間隔0.25度の格子点数121個、y方向(緯度)格子間隔0.2度の格子点数151個です。

また、気圧面データは、z方向（1000hPa・975hPa・950hPa・925hPa・900hPa・850hPa・800hPa・700hPa・600hPa・500hPa・400hPa・300hPa・250hPa・200hPa・150hPa・100hPa）の16個あります。

この情報を下にコントロールファイルを作成します。

GSM(日本域)

DSET gsm.bin
UNDEF 999999999E+20
XDEF 121 LINEAR 120 0.25
YDEF 151 LINEAR 20 0.2
ZDEF 16 LEVELS 1000 975 950 925 900 850 800 700 600 500 400 300 250 200 150 100
TDEF 59 LINEAR 00z01jan0000 3hr
VARS 6
gz 16 99 gz
u 16 99 u
v 16 99 v
temp 16 99 temp
w 16 99 w
rh 16 99 rh
ENDVARS