2014.08.02
ขั้นตอนการอ่านค่าข้อมูลลมจาก
ECMWF เพื่อนำมาประมาณค่าให้ได้ความละเอียดเท่า DEM เพื่อใช้ประมวลผลต่อไปในงานฝนที่เกี่ยวเนื่องกับภูเขา
ขั้นตอน
1. โหลดข้อมูลลมจากเวป ECMWF ณเวลาที่ต้องการ
ผมใช้วันที่ 2010-08-10 ณ เวลา 12UTC เลือกพิกัด lon100E to 140E,lat 10S to 40N ความละเอียดข้อมูลที่ 1.5 องศา เลือกฟอร์แมทแบบ GRIB
wgrib wind_surface_2010081018utc_EPR_ICP.grib |
wgrib -i wind_surface_2010081012utc_EPR_ICP.grib -o
wind_surface_2010081012utc_EPR_ICP.bin
3. อ่านข้อมูลในโปรแกรมฟอร์แทรนตามโค้ดข้างล่างนี้
program read_write2grads
!------------------------------------------------
! 2014.08.02 by Nattapon Mahavik
! Read wind data from ECMWF at one timestep
! first use WGRIB to decode GRIB format to
binary
! then, restructure and export to binary for
GRADS
! note that shifting of u and v for x axis
occurs
! I experiment but not yet found reason
!------------------------------------------------
implicit none
integer,parameter ::
ncol=29,nrow=35,ntim=4 !-x and y
dimension
real,dimension(ncol,nrow) :: u,v,un,vn !-wind orographic and moisture
convergence flux
integer,parameter :: sizeofreal=4 !-record size
integer :: i,j,m,n
integer, parameter :: shft_xu=1, shft_xv=3
!-shifting in x-axis of u,v by comparing to origial data
open (10,file='wind_surface_2010081012_EPR_ICP.bin',form='unformatted',access='direct',recl=ncol*nrow*sizeofreal)
read(10,rec=3)u
read(10,rec=4)v
close(10)
do i=1,ncol
do
j=1,nrow
m=ncol-i+1 !-col inverse
n=nrow-j+1 !-row inverse
un(i- shft_xu,n)=u(i,j)
vn(i- shft_xv,n)=v(i,j)
enddo
enddo
open
(20,file='wind_surface_2010081012_EPR_ICP.grd',form='unformatted',access='direct',recl=ncol*nrow*sizeofreal)
write(20, rec=1)un
write(20, rec=2)vn
close(20)
print*,'done..'
endprogram read_write2grads
4. เปิดในโปรแกรม
GRADS ซึ่งต้องสร้าง Control ไฟล์ ก่อนให้ชื่อว่า “ctl_wind_surface_2010081012_EPR_ICP_ext.ctl”
dset ^wind_surface_2010081012_EPR_ICP.grd
title wind surface 10m ECMWF 1.5 deg
undef 9.999E+20
xdef 29 linear 99.000000 1.500000
ydef 35 linear -10.500000 1.5
zdef 1 linear 1 1
tdef 1 linear 12Z10aug2010 6hr
vars 2
u 0
99 100 zonal wind [m/s]
v 0
99 100 meridional wind [m/s]
endvars
ใช้คำสั่งในการเปิดดังนี้
open
ctl_wind_surface_2010081012_EPR_ICP_ext.ctl
set gxout shaded
d u
d u;v
cbarn
ผลที่ได้ดังภาพข้างล่าง ทำการตรวจสอบแล้วพบว่าไฟล์ไบนารี่ที่ผมสร้างนั้นตรงกับข้อมูลตั้งต้นจากฟอร์แมท
GRIB ที่ดาวน์โหลดมาจาก เวป ECMWF (โปรแกรม grib2ctl เพื่อแปลงเป็นไบนารี่ที่ GRADS อ่านได้)แสดงว่าโปรแกรมถูกต้อง แต่ข้อสงสัยคือ
หากผมไม่ได้ใส่ค่า Shift
ของแกน X ในลม U,V เข้าไปจะทำให้ภาพที่ได้จากไบนารี่ที่ผมสร้างไม่ตรงกับข้อมูลตั้งต้น
ผมตรวจสอบหากข้อมูลแล้วพบว่าแกน X เลื่อนไป 1 และ 3 จุดภาพใน U,V ตามลำดับ ยังหาสาเหตุไม่ได้
แต่อย่างไรก็ตามหากใช้ข้อมูลแค่หนึ่งช่วงเวลา ก็สามารถนำขั้นตอนนี้ไปใช้ได้
ตัวอย่างภาพเลื่อนตามที่บอกดังข้างล่าง
ไม่มีความคิดเห็น:
แสดงความคิดเห็น