データの時刻の定義に関する情報
https://software.ecmwf.int/wiki/pages/viewpage.action?pageId=56658233
解析値は日に4回 (00:00, 06:00, 12:00, 18:00)
予報値は、各解析値を初期値として、3, 6, 9, 12時間後の値がある。
変数によっては予報値しかないものがあることに注意
使用するデータ（surfaceのデータ）のダウンロードサイト (netCDF)
http://apps.ecmwf.int/datasets/data/interim-full-daily/levtype=sfc/
使用する変数
float longitude(longitude) ;
float latitude(latitude) ;
int time(time) ;
short e(time, latitude, longitude) ;
short tp(time, latitude, longitude)
short p55.162(time, latitude, longitude) ;
p55.162:scale_factor = 0.000815696301596405 ;
p55.162:add_offset = 26.817747768371 ;
p55.162:_FillValue = -32767s ;
p55.162:missing_value = -32767s ;
p55.162:units = "kg m**-2" ;
p55.162:long_name = "Vertical integral of water vapour" ;
short tcwv(time, latitude, longitude) ;
tcwv:scale_factor = 0.00081569645439851 ;
tcwv:add_offset = 26.8177309001802 ;
tcwv:_FillValue = -32767s ;
tcwv:missing_value = -32767s ;
tcwv:units = "kg m**-2" ;
tcwv:long_name = "Total column water vapour" ;
tcwv:standard_name = "lwe_thickness_of_atmosphere_mass_content_of_water_vapor" ;
short p84.162(time, latitude, longitude) ;
p84.162:scale_factor = 1.08452245927098e-07 ;
p84.162:add_offset = 0.000282594957522881 ;
p84.162:_FillValue = -32767s ;
p84.162:missing_value = -32767s ;
p84.162:units = "kg m**-2 s**-1" ;
p84.162:long_name = "Vertical integral of divergence of moisture flux" ;
どの時間ステップを使うかは、変数によって異なる
詳細は、末尾の「ダウンロードしたファイルに関する情報」を参照のこと
スクリプト（NCL)使用
------------------------------
List of the following files:
------------------------------
ERA-I.Evapo.run.sh
ERA-I.Evapo.ncl
ERA-I.IWV_diff.run.sh
ERA-I.IWV_diff.ncl
ERA-I.WVF.run.sh
ERA-I.WVF.ncl
ERA-I.IWV.run.sh
ERA-I.IWV.ncl
ERA-I.Precip.run.sh
ERA-I.Precip.ncl
runncl.sh
------------------------------
Machine info
------------------------------
aofd165.bio.mie-u.ac.jp
/work1/am/2016.PolarLow/ERA-I/Water.Vapor.Budget.ERA-I
======================
ERA-I.Evapo.run.sh
======================
#!/bin/sh
exe=./runncl.sh
ncl=$(basename $0 .run.sh).ncl
if [ ! -f $ncl ]; then
echo Error in $0 : No such file, $ncl
exit 1
fi
indir="."
dirs="$indir"
for dir in $dirs; do
if [ ! -d $dir ]; then
echo Error in $0 : No such directory, $dir
exit 1
fi
done
input="ERA-I_E-P_110101-31.nc"
if [ ! -f $indir/$input ];then
echo Error in $0: No such file, $indir/$input
exit 1
fi
export LANG=C
$exe $ncl "$indir" "$input"
exit 0
----------------------
End of ERA-I.Evapo.run.sh
----------------------
======================
ERA-I.Evapo.ncl
======================
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
begin
month_abbr = (/"","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep", \
"Oct","Nov","Dec"/)
wallClock1 = systemfunc("date") ; retrieve wall clock time
scriptname_i = getenv("NCL_ARG_1")
indir = getenv("NCL_ARG_2")
input = getenv("NCL_ARG_3")
print("")
print(scriptname_i + " is running ...")
scriptname=systemfunc("basename " + scriptname_i + " .ncl")
inname=systemfunc("basename "+ input + " .nc")
outdir="Fig"+"_"+scriptname
system("mkdir -vp "+outdir)
infile=indir+"/"+input
a=addfile(infile,"r")
latmin=40.0
latmax=80.0
lonmin=-30.0
lonmax=90.0
min_lat=71.2
max_lat=76.5
min_lon=-6.0
max_lon=25.5
min_lat2=71.2
max_lat2=78.0
min_lon2=25.501
max_lon2=40.0
m2mm=1000.0
scl=1.E4
hr2s=3600.0
;lon=a->longitude
;lat=a->latitude
time=a->time
time@units="hours since 1900-01-01 00:00:0.0"
mt=dimsizes(time)
;print(mt)
nstart=0
nend=mt-1; 6
nstep=1
e=a->e ;Evaporation, m of water equivalent
lat=a->latitude
lon=a->longitude
e@_FillValue=default_fillvalue(typeof(e))
;print(e@_FillValue)
e2 =(e * e@scale_factor + e@add_offset)*(-1.0)*m2mm ;downward (m->mm)
;;;e2 = e2/hr2s*scl
e2@units="mm"
e2!0="time"
e2!1="lat"
e2!2="lon"
e2&time=time
e2&lat=lat
e2&lon=lon
e2&lat@units="degrees_north"
e2&lon@units="degrees_east"
;print("max(e2)="+max(e2)+" min(e2)="+min(e2))
dim=dimsizes(e2)
ny=dim(1)
nx=dim(2)
e2s=new((/ny,nx/),typeof(e2))
e2s!0="lat"
e2s!1="lon"
e2s&lat=lat
e2s&lon=lon
e2s&lat@units="degrees_north"
e2s&lon@units="degrees_east"
e2s_sub=new((/ny,nx/),typeof(e2))
e2s_sub!0="lat"
e2s_sub!1="lon"
e2s_sub&lat=lat
e2s_sub&lon=lon
e2s_sub&lat@units="degrees_north"
e2s_sub&lon@units="degrees_east"
e2s_sub2=new((/ny,nx/),typeof(e2))
e2s_sub2!0="lat"
e2s_sub2!1="lon"
e2s_sub2&lat=lat
e2s_sub2&lon=lon
e2s_sub2&lat@units="degrees_north"
e2s_sub2&lon@units="degrees_east"
lat2d=new((/ny,nx/),typeof(lat))
lon2d=new((/ny,nx/),typeof(lon))
;printVarSummary(lon)
;printVarSummary(lat)
;printVarSummary(lon2d)
;printVarSummary(lat2d)
do j=0,ny-1
do i=0,nx-1
lat2d(j,i)=lat(j)
lon2d(j,i)=lon(i)
end do
end do
;printVarSummary(e2s)
e2s_ave=new((/mt/),typeof(e2))
;printVarSummary(e2s_ave)
dthr=time(1)-time(0)
dthr2=tofloat(dthr)*2.0
dt=dthr*3600.0
dt2=tofloat(dt*2.0)
;print(dt2)
utc_date = cd_calendar(time, 0)
year = tointeger(utc_date(:,0)) ; Convert to integer for
month = tointeger(utc_date(:,1)) ; use sprinti
day = tointeger(utc_date(:,2))
hour = tointeger(utc_date(:,3))
minute = tointeger(utc_date(:,4))
second = utc_date(:,5)
date_num = sprinti("%0.4i", year)+sprinti("%0.2i", month)+sprinti("%0.2i", day) \
+ sprinti("_%0.2i", hour)
date_str = sprinti("%0.2i", hour)+":00UTC"+sprinti("%0.2i", day)+month_abbr(month)\
+sprinti("%0.4i", year)
time_def=time-6
time_def@units="hours since 1900-01-01 00:00:0.0"
utc_date_def = cd_calendar(time_def, 0)
year = tointeger(utc_date_def(:,0)) ; Convert to integer for
month = tointeger(utc_date_def(:,1)) ; use sprinti
day = tointeger(utc_date_def(:,2))
hour = tointeger(utc_date_def(:,3))
minute = tointeger(utc_date_def(:,4))
second = utc_date_def(:,5)
date_num_def = sprinti("%0.4i", year)+sprinti("%0.2i", month)+sprinti("%0.2i", day) \
+ sprinti("_%0.2i", hour)
date_str_def = sprinti("%0.2i", hour)+":00UTC"+sprinti("%0.2i", day)+month_abbr(month)\
+sprinti("%0.4i", year)
; For Generic mapping tools
datetime_out = sprinti("%0.4i", year)+"-"+sprinti("%0.2i", month)+"-"+\
sprinti("%0.2i", day )+"T" + sprinti("%0.2i", hour)+\
":00:00"
;print(date_str)
figtype="ps" ;"x11"
opts=True
;Coordinates and domain
opts@gsnAddCyclic = False
;Plotting area
opts@vpHeightF = .2 ;1
;opts@vpWidthF = opts@vpHeightF*tofloat(mlon)/tofloat(nlat) ;.4
;Page control
opts@gsnFrame = False
opts@gsnDraw = False
opts@gsnStringFontHeightF = 0.012
opts@lbLabelFontHeightF = 0.012
; Contour line
opts_cn=opts
opts_cn@cnFillOn = False
opts_cn@cnLevelSpacingF = 1
opts_cn@cnLineThicknessF= 1.0
opts_cn@cnInfoLabelOn = False
opts_cn@cnLabelDrawOrder="PostDraw"
opts_cn@cnLineLabelFontHeightF = 0.01 ;8
; Color shade
opts_cs=opts
opts_cs@cnFillOn = True
opts_cs@pmLabelBarOrthogonalPosF = 0.0
opts_cs@pmLabelBarHeightF = 0.04
; Color bar
opts_cs@pmLabelBarOrthogonalPosF = .10 ;move whole thing down
opts_cs@lbTitlePosition = "Right" ;title position
opts_cs@lbTitleFontHeightF= .012 ;make title smaller
opts_cs@lbTitleDirection = "Across" ;title direction
ots_cs_nomap=opts_cs
; MAP
opts_cs@mpProjection = "Stereographic"
opts_cs@mpLimitMode = "Corners"
opts_cs@mpLeftCornerLatF = latmin
opts_cs@mpLeftCornerLonF = lonmin
opts_cs@mpRightCornerLatF = latmax
opts_cs@mpRightCornerLonF = lonmax
opts_cs@mpCenterLonF = 20.
opts_cs@mpCenterLatF = 70. ;This is necessary to fix the wrong value on the WRF file.
opts_cs@mpDataBaseVersion = "HighRes" ;fill land in gray and ocean in transparent
opts_cs@mpAreaMaskingOn = "True"
opts_cs@mpFillOn = False ; True
; VECTOR
vecres = True ;vector only resources
vecres@gsnAddCyclic = False
vecres@gsnDraw = False ;don't draw
vecres@gsnFrame = False ;don't advance frame
vecres@vcGlyphStyle = "LineArrow" ;"FillArrow"
;vecres@vcLineArrowThicknessF = 3 ;change vector thickness
vecres@vcRefLengthF = 0.05 ;define length of vec ref
vecres@vcRefAnnoOrthogonalPosF=-1.0
vecres@vcRefAnnoParallelPosF =1.2
;vecres@vcLabelFontHeightF = 0.08
vecres@gsnRightString = " " ;turn off right string
vecres@gsnLeftString = " " ;turn off left string
vecres@tiXAxisString = " " ;turn off axis label
vecres@vcRefMagnitudeF = 500 ;define vector ref mag
vecres@vcRefAnnoOrthogonalPosF = -0.36 ;move ref vector into plot
xskp=6
yskp=3
; Add a box showing the area of interest
lnres = True
lnres@gsLineColor = "Red"
lnres@gsLineThicknessF = 2.0
res1=opts_cs
res1@cnLinesOn = False
res1@cnLevelSelectionMode = "ExplicitLevels" ;
clev1=(/-0.5,0.0,0.5,1.0,1.5,2.0,2.5,3.0/)
res1@cnLevels = clev1
res2=opts_cs
res2@gsnLeftString = "E [mm/hr]"
res2@gsnCenterString = ""
res2@gsnRightString = ""
res2@cnLinesOn = False
res2@cnLevelSelectionMode = "ExplicitLevels" ;
clev2=(/-0.4,-0.3,-0.2,-0.1, 0.0, 0.1, 0.2, 0.3, 0.4/)
res2@cnLevels = clev2
res22=ots_cs_nomap
res22@cnLinesOn = False
res22@cnLevelSelectionMode = "ExplicitLevels" ;
clev2=(/-0.4,-0.3,-0.2,-0.1, 0.0, 0.1, 0.2, 0.3, 0.4/)
res22@cnLevels = clev2
res22@lbLabelBarOn = False
res22@cnInfoLabelOn = False
;printVarSummary(e)
nt=0
do it=nstart,nend,nstep
print("")
print("it="+it+" "+date_str(it))
templete= scriptname+"_"+inname
figfile = outdir+"/"+ templete + "_" + date_num_def(it)
wks = gsn_open_wks(figtype, figfile)
res1@vpXF = 0.0
res2@vpXF = 0.0
itp=it
res1@gsnLeftString = "E [mm] " +itp
res1@gsnRightString = date_str(itp)
res1@vpYF = 0.85
gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")
plot1=gsn_csm_contour_map(wks,e2(itp,:,:),res1)
print("max(e2)=" + max(e2(itp,:,:)) + " min(e2)=" +min(e2(itp,:,:)))
draw(plot1)
delete(plot1)
res2@gsnLeftString = "E [mm/h]"
res2@gsnRightString = date_str_def(itp)
res2@vpYF = 0.5
gsn_define_colormap(wks,"BlueWhiteOrangeRed")
e2s(:,:)=e2(itp,:,:)/dthr
e2s_sub= where(lat2d.lt.min_lat.or.lat2d.gt.max_lat .or. \
lon2d.lt.min_lon.or.lon2d.gt.max_lon, \
lat2d@_FillValue, e2s)
e2s_sub2=where(lat2d.lt.min_lat2.or.lat2d.gt.max_lat2 .or. \
lon2d.lt.min_lon2.or.lon2d.gt.max_lon2, \
lat2d@_FillValue, e2s)
e2s_sub@_FillValue=e2s@_FillValue
e2s_sub2@_FillValue=e2s@_FillValue
e2s_ave(it) = avg(e2s_sub)+avg(e2s_sub2)
plot2=gsn_csm_contour_map(wks,e2s_sub(:,:),res2)
plot22=gsn_csm_contour(wks,e2s_sub2(:,:),res22)
overlay(plot2,plot22)
dum4 = gsn_add_polyline(wks,plot2,(/min_lon,max_lon,\
max_lon,min_lon,min_lon/),\
(/min_lat,min_lat,max_lat,\
max_lat,min_lat/),lnres)
dum5 = gsn_add_polyline(wks,plot2,(/min_lon2,max_lon2,\
max_lon2,min_lon2,min_lon2/),\
(/min_lat2,min_lat2,max_lat2,\
max_lat2,min_lat2/),lnres)
draw(plot2)
delete(plot2)
;Header
txres = True
txres@txJust="CenterLeft"
txres@txFontHeightF = 0.010
lst = systemfunc("ls -lh "+infile)
strs = str_split(lst, " ")
str1 = "Input: "+strs(8)
str2 = strs(0)+" "+strs(1)+" "+strs(2)+" "+strs(3)+" "+strs(4)+" "+\
strs(5)+" "+strs(6)+" "+strs(7)
gsn_text_ndc(wks,str1,0.05,0.99,txres)
gsn_text_ndc(wks,str2,0.05,0.97,txres)
txres@txJust="CenterCenter"
txres@txFontHeightF = 0.016
text=scriptname
gsn_text_ndc(wks,text,0.5,0.93,txres)
txres@txFontHeightF = 0.016
text=date_str_def(it)
gsn_text_ndc(wks,text,0.5,0.90,txres)
;Footer
txres@txFontHeightF = 0.01
txres@txJust="CenterLeft"
text = systemfunc("export LANC=C; date")
gsn_text_ndc(wks,text, 0.05,0.055,txres)
cwd =systemfunc("pwd")
gsn_text_ndc(wks,"Current dir: "+cwd, 0.05,0.040,txres)
gsn_text_ndc(wks,"Script: "+scriptname, 0.05,0.025,txres)
gsn_text_ndc(wks,"Output dir: "+outdir, 0.05,0.010,txres)
frame(wks)
end do ;it
e2s_ave_kgm2s2=e2s_ave/hr2s*scl
alist= [/datetime_out, e2s_ave, e2s_ave_kgm2s2/]
;print(alist)
;
; TIME SERIES
;
templete= scriptname+"_"+inname
figfile = outdir+"/"+ templete + "_Time.Series"
wks = gsn_open_wks(figtype, figfile)
res = True ;plot mods desired
res@gsnFrame = False ;don't advance frame yet
res@gsnDraw = False ;don't draw plot
res@xyLineThicknessF = 2.0 ;make thicker
res@gsnLeftString=""
res@gsnRightString=""
;printVarSummary(time)
restick = True
restick@ttmFormat = "%h %c %d" ;"%N/%D %H:%M"
time_axis_labels(time,res,restick)
plot = gsn_csm_xy(wks,time,e2s_ave,res)
draw(plot)
delete(res)
;Header
txres = True
txres@txJust="CenterLeft"
txres@txFontHeightF = 0.010
lst = systemfunc("ls -lh "+infile)
strs = str_split(lst, " ")
str1 = "Input: "+strs(8)
str2 = strs(0)+" "+strs(1)+" "+strs(2)+" "+strs(3)+" "+strs(4)+" "+\
strs(5)+" "+strs(6)+" "+strs(7)
gsn_text_ndc(wks,str1,0.05,0.99,txres)
gsn_text_ndc(wks,str2,0.05,0.97,txres)
txres@txJust="CenterCenter"
txres@txFontHeightF = 0.016
text=scriptname
gsn_text_ndc(wks,text,0.5,0.93,txres)
;txres@txFontHeightF = 0.016
;text=chartostring(times0)
;gsn_text_ndc(wks,text,0.5,0.90,txres)
;Footer
txres@txFontHeightF = 0.01
txres@txJust="CenterLeft"
today = systemfunc("export LANC=C; date")
gsn_text_ndc(wks,today, 0.05,0.010,txres)
cwd =systemfunc("pwd")
gsn_text_ndc(wks,"Current dir: "+cwd, 0.05,0.025,txres)
gsn_text_ndc(wks,"Output dir: "+outdir,0.05,0.040,txres)
frame(wks)
;
; ASCII OUT
;
min_lat=71.2
max_lat=76.5
min_lon=-6.0
max_lon=25.5
min_lat2=71.2
max_lat2=78.0
min_lon2=25.501
max_lon2=40.0
header_ofle = (/\
"# "+today,\
"# Input : "+infile,\
"# CWD : "+cwd,\
"# ave area : "+min_lat+" - "+max_lat+" "+min_lon+" - "+max_lon,\
"# ave area2: "+min_lat2+" - "+max_lat2+" "+min_lon2+" - "+max_lon2\
/)
hlist=[/header_ofle/]
ofle=scriptname+"_"+inname+"_Time.Series.txt"
write_table(ofle, "w", hlist, "%s")
write_table(ofle, "a", alist, "%s %11.4f %12.4f")
print("")
print("Output: "+figfile)
print("Output: "+ofle)
print("")
end
----------------------
End of ERA-I.Evapo.ncl
----------------------
======================
ERA-I.IWV_diff.run.sh
======================
#!/bin/sh
exe=./runncl.sh
ncl=$(basename $0 .run.sh).ncl
if [ ! -f $ncl ]; then
echo Error in $0 : No such file, $ncl
exit 1
fi
indir="."
dirs="$indir"
for dir in $dirs; do
if [ ! -d $dir ]; then
echo Error in $0 : No such directory, $dir
exit 1
fi
done
input="ERA-I_IWV_diff_110101-31.nc"
if [ ! -f $indir/$input ];then
echo Error in $0: No such file, $indir/$input
exit 1
fi
export LANG=C
$exe $ncl "$indir" "$input"
exit 0
----------------------
End of ERA-I.IWV_diff.run.sh
----------------------
======================
ERA-I.IWV_diff.ncl
======================
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
begin
month_abbr = (/"","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep", \
"Oct","Nov","Dec"/)
wallClock1 = systemfunc("date") ; retrieve wall clock time
scriptname_i = getenv("NCL_ARG_1")
indir = getenv("NCL_ARG_2")
input = getenv("NCL_ARG_3")
print("")
print(scriptname_i + " is running ...")
scriptname=systemfunc("basename " + scriptname_i + " .ncl")
inname=systemfunc("basename "+ input + " .nc")
outdir="Fig"+"_"+scriptname
system("mkdir -vp "+outdir)
infile=indir+"/"+input
a=addfile(infile,"r")
latmin=40.0
latmax=80.0
lonmin=-30.0
lonmax=90.0
min_lat=71.2
max_lat=76.5
min_lon=-6.0
max_lon=25.5
min_lat2=71.2
max_lat2=78.0
min_lon2=25.501
max_lon2=40.0
m2mm=1000.0
scl=1.E4
hr2s=3600.0
;lon=a->longitude
;lat=a->latitude
time=a->time
time@units="hours since 1900-01-01 00:00:0.0"
mt=dimsizes(time)
;print(mt)
nstart=0
nend=mt-1-1; 6
nstep=1
iwv=a->$"p55.162"$
; p55.162:units = "kg m**-2" ;
; p55.162:long_name = "Vertical integral of water vapour"
tcwv=a->tcwv
; tcwv:units = "kg m**-2" ;
; tcwv:long_name = "Total column water vapour" ;
lat=a->latitude
lon=a->longitude
iwv@_FillValue=default_fillvalue(typeof(iwv))
iwv2 =(iwv * iwv@scale_factor + iwv@add_offset)
tcwv@_FillValue=default_fillvalue(typeof(tcwv))
tcwv2 =(tcwv * tcwv@scale_factor + tcwv@add_offset)
;
; Vertical integral of water vapour, kg m**-2
;
iwv2@units="kg m**-2"
iwv2!0="time"
iwv2!1="lat"
iwv2!2="lon"
iwv2&time=time
iwv2&lat=lat
iwv2&lon=lon
iwv2&lat@units="degrees_north"
iwv2&lon@units="degrees_east"
dim=dimsizes(iwv2)
ny=dim(1)
nx=dim(2)
iwv2s=new((/ny,nx/),typeof(iwv2))
iwv2s!0="lat"
iwv2s!1="lon"
iwv2s&lat=lat
iwv2s&lon=lon
iwv2s&lat@units="degrees_north"
iwv2s&lon@units="degrees_east"
iwv2s_sub=new((/ny,nx/),typeof(iwv2))
iwv2s_sub!0="lat"
iwv2s_sub!1="lon"
iwv2s_sub&lat=lat
iwv2s_sub&lon=lon
iwv2s_sub&lat@units="degrees_north"
iwv2s_sub&lon@units="degrees_east"
iwv2s_sub2=new((/ny,nx/),typeof(iwv2))
iwv2s_sub2!0="lat"
iwv2s_sub2!1="lon"
iwv2s_sub2&lat=lat
iwv2s_sub2&lon=lon
iwv2s_sub2&lat@units="degrees_north"
iwv2s_sub2&lon@units="degrees_east"
;
; Tendency, kg m**-2 s**-1
;
dim=dimsizes(iwv2)
ny=dim(1)
nx=dim(2)
tend2s=new((/ny,nx/),typeof(iwv2s))
tend2s!0="lat"
tend2s!1="lon"
tend2s&lat=lat
tend2s&lon=lon
tend2s&lat@units="degrees_north"
tend2s&lon@units="degrees_east"
tend2s_sub=new((/ny,nx/),typeof(iwv2s))
tend2s_sub!0="lat"
tend2s_sub!1="lon"
tend2s_sub&lat=lat
tend2s_sub&lon=lon
tend2s_sub&lat@units="degrees_north"
tend2s_sub&lon@units="degrees_east"
tend2s_sub2=new((/ny,nx/),typeof(iwv2s))
tend2s_sub2!0="lat"
tend2s_sub2!1="lon"
tend2s_sub2&lat=lat
tend2s_sub2&lon=lon
tend2s_sub2&lat@units="degrees_north"
tend2s_sub2&lon@units="degrees_east"
lat2d=new((/ny,nx/),typeof(lat))
lon2d=new((/ny,nx/),typeof(lon))
do j=0,ny-1
do i=0,nx-1
lat2d(j,i)=lat(j)
lon2d(j,i)=lon(i)
end do
end do
iwv2s_ave=new((/mt/),typeof(iwv2))
tend2s_ave=new((/mt/),typeof(iwv2))
;printVarSummary(cmf2s_ave)
dthr=time(1)-time(0)
dthr2=tofloat(dthr)*2.0
dt=dthr*3600.0
dt2=tofloat(dt*2.0)
;print(dt2)
utc_date = cd_calendar(time, 0)
year = tointeger(utc_date(:,0)) ; Convert to integer for
month = tointeger(utc_date(:,1)) ; use sprinti
day = tointeger(utc_date(:,2))
hour = tointeger(utc_date(:,3))
minute = tointeger(utc_date(:,4))
second = utc_date(:,5)
date_num = sprinti("%0.4i", year)+sprinti("%0.2i", month)+sprinti("%0.2i", day) \
+ sprinti("_%0.2i", hour)
date_str = sprinti("%0.2i", hour)+":00UTC"+sprinti("%0.2i", day)+month_abbr(month)\
+sprinti("%0.4i", year)
time_def=time+6
time_def@units="hours since 1900-01-01 00:00:0.0"
utc_date_def = cd_calendar(time_def, 0)
year = tointeger(utc_date_def(:,0)) ; Convert to integer for
month = tointeger(utc_date_def(:,1)) ; use sprinti
day = tointeger(utc_date_def(:,2))
hour = tointeger(utc_date_def(:,3))
minute = tointeger(utc_date_def(:,4))
second = utc_date_def(:,5)
date_num_def = sprinti("%0.4i", year)+sprinti("%0.2i", month)+sprinti("%0.2i", day) \
+ sprinti("_%0.2i", hour)
date_str_def = sprinti("%0.2i", hour)+":00UTC"+sprinti("%0.2i", day)+month_abbr(month)\
+sprinti("%0.4i", year)
; For Generic mapping tools
datetime_out = sprinti("%0.4i", year)+"-"+sprinti("%0.2i", month)+"-"+\
sprinti("%0.2i", day )+"T" + sprinti("%0.2i", hour)+\
":00:00"
;print(date_str)
figtype="ps" ;"x11"
opts=True
;Coordinates and domain
opts@gsnAddCyclic = False
;Plotting area
opts@vpHeightF = .2 ;1
;opts@vpWidthF = opts@vpHeightF*tofloat(mlon)/tofloat(nlat) ;.4
;Page control
opts@gsnFrame = False
opts@gsnDraw = False
opts@gsnStringFontHeightF = 0.012
opts@lbLabelFontHeightF = 0.012
; Contour line
opts_cn=opts
opts_cn@cnFillOn = False
opts_cn@cnLevelSpacingF = 1
opts_cn@cnLineThicknessF= 1.0
opts_cn@cnInfoLabelOn = False
opts_cn@cnLabelDrawOrder="PostDraw"
opts_cn@cnLineLabelFontHeightF = 0.01 ;8
; Color shade
opts_cs=opts
opts_cs@cnFillOn = True
opts_cs@pmLabelBarOrthogonalPosF = 0.0
opts_cs@pmLabelBarHeightF = 0.04
; Color bar
opts_cs@pmLabelBarOrthogonalPosF = .10 ;move whole thing down
opts_cs@lbTitlePosition = "Right" ;title position
opts_cs@lbTitleFontHeightF= .012 ;make title smaller
opts_cs@lbTitleDirection = "Across" ;title direction
ots_cs_nomap=opts_cs
; MAP
opts_cs@mpProjection = "Stereographic"
opts_cs@mpLimitMode = "Corners"
opts_cs@mpLeftCornerLatF = latmin
opts_cs@mpLeftCornerLonF = lonmin
opts_cs@mpRightCornerLatF = latmax
opts_cs@mpRightCornerLonF = lonmax
opts_cs@mpCenterLonF = 20.
opts_cs@mpCenterLatF = 70. ;This is necessary to fix the wrong value on the WRF file.
opts_cs@mpDataBaseVersion = "HighRes" ;fill land in gray and ocean in transparent
opts_cs@mpAreaMaskingOn = "True"
opts_cs@mpFillOn = False ; True
; VECTOR
vecres = True ;vector only resources
vecres@gsnAddCyclic = False
vecres@gsnDraw = False ;don't draw
vecres@gsnFrame = False ;don't advance frame
vecres@vcGlyphStyle = "LineArrow" ;"FillArrow"
;vecres@vcLineArrowThicknessF = 3 ;change vector thickness
vecres@vcRefLengthF = 0.05 ;define length of vec ref
vecres@vcRefAnnoOrthogonalPosF=-1.0
vecres@vcRefAnnoParallelPosF =1.2
;vecres@vcLabelFontHeightF = 0.08
vecres@gsnRightString = " " ;turn off right string
vecres@gsnLeftString = " " ;turn off left string
vecres@tiXAxisString = " " ;turn off axis label
vecres@vcRefMagnitudeF = 500 ;define vector ref mag
vecres@vcRefAnnoOrthogonalPosF = -0.36 ;move ref vector into plot
xskp=6
yskp=3
; Add a box showing the area of interest
lnres = True
lnres@gsLineColor = "Red"
lnres@gsLineThicknessF = 2.0
res1=opts_cs
res1@cnLinesOn = False
res1@cnLevelSelectionMode = "ExplicitLevels" ;
clev1=(/4, 8, 12, 16, 20, 24, 28, 32/)
res1@cnLevels = clev1
res2=opts_cs
res2@gsnLeftString = "IWV [kgm~S~-2~N]"
res2@gsnCenterString = ""
res2@gsnRightString = ""
res2@cnLinesOn = False
res2@cnLevelSelectionMode = "ExplicitLevels" ;
clev2=(/-0.8,-0.6,-0.4,-0.2, 0.0, 0.2, 0.4, 0.6, 0.8/) ;(/-2.0,-1.6,-1.2,-0.8,-0.4,0.0,0.4,0.8,1.2,1.6,2.0/)
res2@cnLevels = clev2
res22=ots_cs_nomap
res22@cnLinesOn = False
res22@cnLevelSelectionMode = "ExplicitLevels" ;
res22@cnLevels = clev2
res22@lbLabelBarOn = False
res22@cnInfoLabelOn = False
nt=0
do it=nstart,nend,nstep
print("")
print("it="+it+" "+date_str(it))
itp=it
iwv2s(:,:)=iwv2(itp,:,:)
iwv2s@_FillValue=default_fillvalue(typeof(iwv2s))
;printVarSummary(iwv2s)
print("max(iwv2s)=" + max(iwv2s) + " min(iwv2s)=" +min(iwv2s))
iwv2s_sub= where(lat2d.lt.min_lat.or.lat2d.gt.max_lat .or. \
lon2d.lt.min_lon.or.lon2d.gt.max_lon, \
iwv2s@_FillValue, iwv2s)
iwv2s_sub2=where(lat2d.lt.min_lat2.or.lat2d.gt.max_lat2 .or. \
lon2d.lt.min_lon2.or.lon2d.gt.max_lon2, \
iwv2s@_FillValue, iwv2s)
iwv2s_sub@_FillValue=iwv2s@_FillValue
iwv2s_sub2@_FillValue=iwv2s@_FillValue
iwv2s_ave(it) = avg(iwv2s_sub)+avg(iwv2s_sub2)
tend2s(:,:)=(iwv2(itp+1,:,:)-iwv2(itp,:,:))/dt*scl
tend2s@_FillValue=default_fillvalue(typeof(tend2s))
print("max(tend2s)=" + max(tend2s) + " min(tend2s)=" +min(tend2s))
tend2s_sub= where(lat2d.lt.min_lat.or.lat2d.gt.max_lat .or. \
lon2d.lt.min_lon.or.lon2d.gt.max_lon, \
tend2s@_FillValue, tend2s)
tend2s_sub2=where(lat2d.lt.min_lat2.or.lat2d.gt.max_lat2 .or. \
lon2d.lt.min_lon2.or.lon2d.gt.max_lon2, \
tend2s@_FillValue, tend2s)
tend2s_sub@_FillValue=tend2s@_FillValue
tend2s_sub2@_FillValue=tend2s@_FillValue
tend2s_ave(it) = avg(tend2s_sub)+avg(tend2s_sub2)
templete= scriptname+"_"+inname
figfile = outdir+"/"+ templete + "_" + date_num_def(it)
wks = gsn_open_wks(figtype, figfile)
res1@vpXF = 0.0
res2@vpXF = 0.0
res1@gsnLeftString = "IWV [kgm~S~-2~N~]"
res1@gsnRightString = date_str(itp)
res1@vpYF = 0.85
gsn_define_colormap(wks,"WhiteBlueGreenYellowRed");"BlueWhiteOrangeRed")
plot1=gsn_csm_contour_map(wks,iwv2s(:,:),res1)
draw(plot1)
delete(plot1)
res2@gsnLeftString = "Tnd [kgm~S~-2~N~s~S~-2~N~]"
res2@gsnRightString = date_str_def(itp)
res2@vpYF = 0.5
gsn_define_colormap(wks,"BlueWhiteOrangeRed");"WhiteBlueGreenYellowRed")
plot2=gsn_csm_contour_map(wks,tend2s(:,:),res2)
;plot22=gsn_csm_contour(wks,tend2s(:,:),res22)
;overlay(plot2,plot22)
dum4 = gsn_add_polyline(wks,plot2,(/min_lon,max_lon,\
max_lon,min_lon,min_lon/),\
(/min_lat,min_lat,max_lat,\
max_lat,min_lat/),lnres)
dum5 = gsn_add_polyline(wks,plot2,(/min_lon2,max_lon2,\
max_lon2,min_lon2,min_lon2/),\
(/min_lat2,min_lat2,max_lat2,\
max_lat2,min_lat2/),lnres)
draw(plot2)
delete(plot2)
;Header
txres = True
txres@txJust="CenterLeft"
txres@txFontHeightF = 0.010
lst = systemfunc("ls -lh "+infile)
strs = str_split(lst, " ")
str1 = "Input: "+strs(8)
str2 = strs(0)+" "+strs(1)+" "+strs(2)+" "+strs(3)+" "+strs(4)+" "+\
strs(5)+" "+strs(6)+" "+strs(7)
gsn_text_ndc(wks,str1,0.05,0.99,txres)
gsn_text_ndc(wks,str2,0.05,0.97,txres)
txres@txJust="CenterCenter"
txres@txFontHeightF = 0.016
text=scriptname
gsn_text_ndc(wks,text,0.5,0.93,txres)
txres@txFontHeightF = 0.016
text=date_str_def(it)
gsn_text_ndc(wks,text,0.5,0.90,txres)
;Footer
txres@txFontHeightF = 0.01
txres@txJust="CenterLeft"
text = systemfunc("export LANC=C; date")
gsn_text_ndc(wks,text, 0.05,0.055,txres)
cwd =systemfunc("pwd")
gsn_text_ndc(wks,"Current dir: "+cwd, 0.05,0.040,txres)
gsn_text_ndc(wks,"Script: "+scriptname, 0.05,0.025,txres)
gsn_text_ndc(wks,"Output dir: "+outdir, 0.05,0.010,txres)
frame(wks)
end do ;it
alist= [/datetime_out(0:mt-3), iwv2s_ave(0:mt-3), tend2s_ave(0:mt-3)/]
;print(alist)
;
; TIME SERIES
;
templete= scriptname+"_"+inname
figfile = outdir+"/"+ templete + "_Time.Series"
wks = gsn_open_wks(figtype, figfile)
res = True ;plot mods desired
res@gsnFrame = False ;don't advance frame yet
res@gsnDraw = False ;don't draw plot
res@xyLineThicknessF = 2.0 ;make thicker
res@gsnLeftString=""
res@gsnRightString=""
;printVarSummary(time)
restick = True
restick@ttmFormat = "%h %c %d" ;"%N/%D %H:%M"
time_axis_labels(time,res,restick)
plot = gsn_csm_xy(wks,time,iwv2s_ave,res)
draw(plot)
delete(res)
;Header
txres = True
txres@txJust="CenterLeft"
txres@txFontHeightF = 0.010
lst = systemfunc("ls -lh "+infile)
strs = str_split(lst, " ")
str1 = "Input: "+strs(8)
str2 = strs(0)+" "+strs(1)+" "+strs(2)+" "+strs(3)+" "+strs(4)+" "+\
strs(5)+" "+strs(6)+" "+strs(7)
gsn_text_ndc(wks,str1,0.05,0.99,txres)
gsn_text_ndc(wks,str2,0.05,0.97,txres)
txres@txJust="CenterCenter"
txres@txFontHeightF = 0.016
text=scriptname
gsn_text_ndc(wks,text,0.5,0.93,txres)
;txres@txFontHeightF = 0.016
;text=chartostring(times0)
;gsn_text_ndc(wks,text,0.5,0.90,txres)
;Footer
txres@txFontHeightF = 0.01
txres@txJust="CenterLeft"
today = systemfunc("export LANC=C; date")
gsn_text_ndc(wks,today, 0.05,0.010,txres)
cwd =systemfunc("pwd")
gsn_text_ndc(wks,"Current dir: "+cwd, 0.05,0.025,txres)
gsn_text_ndc(wks,"Output dir: "+outdir,0.05,0.040,txres)
frame(wks)
;
; ASCII OUT
;
min_lat=71.2
max_lat=76.5
min_lon=-6.0
max_lon=25.5
min_lat2=71.2
max_lat2=78.0
min_lon2=25.501
max_lon2=40.0
header_ofle = (/\
"# "+today,\
"# Input : "+infile,\
"# CWD : "+cwd,\
"# ave area : "+min_lat+" - "+max_lat+" "+min_lon+" - "+max_lon,\
"# ave area2: "+min_lat2+" - "+max_lat2+" "+min_lon2+" - "+max_lon2\
/)
hlist=[/header_ofle/]
ofle=scriptname+"_"+inname+"_Time.Series.txt"
write_table(ofle, "w", hlist, "%s")
write_table(ofle, "a", alist, "%s %11.4f %12.4f")
print("")
print("Output: "+figfile)
print("Output: "+ofle)
print("")
end
----------------------
End of ERA-I.IWV_diff.ncl
----------------------
======================
ERA-I.WVF.run.sh
======================
#!/bin/sh
exe=./runncl.sh
ncl=$(basename $0 .run.sh).ncl
if [ ! -f $ncl ]; then
echo Error in $0 : No such file, $ncl
exit 1
fi
indir="."
dirs="$indir"
for dir in $dirs; do
if [ ! -d $dir ]; then
echo Error in $0 : No such directory, $dir
exit 1
fi
done
input="ERA-I_WVF_110101-31.nc"
if [ ! -f $indir/$input ];then
echo Error in $0: No such file, $indir/$input
exit 1
fi
export LANG=C
$exe $ncl "$indir" "$input"
exit 0
----------------------
End of ERA-I.WVF.run.sh
----------------------
======================
ERA-I.WVF.ncl
======================
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
begin
month_abbr = (/"","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep", \
"Oct","Nov","Dec"/)
wallClock1 = systemfunc("date") ; retrieve wall clock time
scriptname_i = getenv("NCL_ARG_1")
indir = getenv("NCL_ARG_2")
input = getenv("NCL_ARG_3")
print("")
print(scriptname_i + " is running ...")
scriptname=systemfunc("basename " + scriptname_i + " .ncl")
inname=systemfunc("basename "+ input + " .nc")
outdir="Fig"+"_"+scriptname
system("mkdir -vp "+outdir)
infile=indir+"/"+input
a=addfile(infile,"r")
latmin=40.0
latmax=80.0
lonmin=-30.0
lonmax=90.0
min_lat=71.2
max_lat=76.5
min_lon=-6.0
max_lon=25.5
min_lat2=71.2
max_lat2=78.0
min_lon2=25.501
max_lon2=40.0
m2mm=1000.0
scl=1.E4
hr2s=3600.0
;lon=a->longitude
;lat=a->latitude
time=a->time
time@units="hours since 1900-01-01 00:00:0.0"
mt=dimsizes(time)
;print(mt)
nstart=0
nend=mt-1; 6
nstep=1
divmf=a->$"p84.162"$ ;Vertical integral of DIVERGENCE of moisture flux, kg m**-2 s**-1
lat=a->latitude
lon=a->longitude
divmf@_FillValue=default_fillvalue(typeof(divmf))
divmf2 =(divmf * divmf@scale_factor + divmf@add_offset)
;
; Vertical integral of CONVERGENCE of moisture flux, kg m**-2 s**-1
;
cmf2 = divmf2*(-1.0)*scl ;div->cnv, scl
;printVarSummary(cmf2)
cmf2@units="kg m**-2 s**-1"
cmf2!0="time"
cmf2!1="lat"
cmf2!2="lon"
cmf2&time=time
cmf2&lat=lat
cmf2&lon=lon
cmf2&lat@units="degrees_north"
cmf2&lon@units="degrees_east"
dim=dimsizes(cmf2)
ny=dim(1)
nx=dim(2)
cmf2s=new((/ny,nx/),typeof(cmf2))
cmf2s!0="lat"
cmf2s!1="lon"
cmf2s&lat=lat
cmf2s&lon=lon
cmf2s&lat@units="degrees_north"
cmf2s&lon@units="degrees_east"
cmf2s_sub=new((/ny,nx/),typeof(cmf2))
cmf2s_sub!0="lat"
cmf2s_sub!1="lon"
cmf2s_sub&lat=lat
cmf2s_sub&lon=lon
cmf2s_sub&lat@units="degrees_north"
cmf2s_sub&lon@units="degrees_east"
cmf2s_sub2=new((/ny,nx/),typeof(cmf2))
cmf2s_sub2!0="lat"
cmf2s_sub2!1="lon"
cmf2s_sub2&lat=lat
cmf2s_sub2&lon=lon
cmf2s_sub2&lat@units="degrees_north"
cmf2s_sub2&lon@units="degrees_east"
lat2d=new((/ny,nx/),typeof(lat))
lon2d=new((/ny,nx/),typeof(lon))
do j=0,ny-1
do i=0,nx-1
lat2d(j,i)=lat(j)
lon2d(j,i)=lon(i)
end do
end do
;printVarSummary(cmf2s)
cmf2s_ave=new((/mt/),typeof(cmf2))
;printVarSummary(cmf2s_ave)
dthr=time(1)-time(0)
dthr2=tofloat(dthr)*2.0
dt=dthr*3600.0
dt2=tofloat(dt*2.0)
;print(dt2)
utc_date = cd_calendar(time, 0)
year = tointeger(utc_date(:,0)) ; Convert to integer for
month = tointeger(utc_date(:,1)) ; use sprinti
day = tointeger(utc_date(:,2))
hour = tointeger(utc_date(:,3))
minute = tointeger(utc_date(:,4))
second = utc_date(:,5)
date_num = sprinti("%0.4i", year)+sprinti("%0.2i", month)+sprinti("%0.2i", day) \
+ sprinti("_%0.2i", hour)
date_str = sprinti("%0.2i", hour)+":00UTC"+sprinti("%0.2i", day)+month_abbr(month)\
+sprinti("%0.4i", year)
time_def=time-0
time_def@units="hours since 1900-01-01 00:00:0.0"
utc_date_def = cd_calendar(time_def, 0)
year = tointeger(utc_date_def(:,0)) ; Convert to integer for
month = tointeger(utc_date_def(:,1)) ; use sprinti
day = tointeger(utc_date_def(:,2))
hour = tointeger(utc_date_def(:,3))
minute = tointeger(utc_date_def(:,4))
second = utc_date_def(:,5)
date_num_def = sprinti("%0.4i", year)+sprinti("%0.2i", month)+sprinti("%0.2i", day) \
+ sprinti("_%0.2i", hour)
date_str_def = sprinti("%0.2i", hour)+":00UTC"+sprinti("%0.2i", day)+month_abbr(month)\
+sprinti("%0.4i", year)
; For Generic mapping tools
datetime_out = sprinti("%0.4i", year)+"-"+sprinti("%0.2i", month)+"-"+\
sprinti("%0.2i", day )+"T" + sprinti("%0.2i", hour)+\
":00:00"
;print(date_str)
figtype="ps" ;"x11"
opts=True
;Coordinates and domain
opts@gsnAddCyclic = False
;Plotting area
opts@vpHeightF = .2 ;1
;opts@vpWidthF = opts@vpHeightF*tofloat(mlon)/tofloat(nlat) ;.4
;Page control
opts@gsnFrame = False
opts@gsnDraw = False
opts@gsnStringFontHeightF = 0.012
opts@lbLabelFontHeightF = 0.012
; Contour line
opts_cn=opts
opts_cn@cnFillOn = False
opts_cn@cnLevelSpacingF = 1
opts_cn@cnLineThicknessF= 1.0
opts_cn@cnInfoLabelOn = False
opts_cn@cnLabelDrawOrder="PostDraw"
opts_cn@cnLineLabelFontHeightF = 0.01 ;8
; Color shade
opts_cs=opts
opts_cs@cnFillOn = True
opts_cs@pmLabelBarOrthogonalPosF = 0.0
opts_cs@pmLabelBarHeightF = 0.04
; Color bar
opts_cs@pmLabelBarOrthogonalPosF = .10 ;move whole thing down
opts_cs@lbTitlePosition = "Right" ;title position
opts_cs@lbTitleFontHeightF= .012 ;make title smaller
opts_cs@lbTitleDirection = "Across" ;title direction
ots_cs_nomap=opts_cs
; MAP
opts_cs@mpProjection = "Stereographic"
opts_cs@mpLimitMode = "Corners"
opts_cs@mpLeftCornerLatF = latmin
opts_cs@mpLeftCornerLonF = lonmin
opts_cs@mpRightCornerLatF = latmax
opts_cs@mpRightCornerLonF = lonmax
opts_cs@mpCenterLonF = 20.
opts_cs@mpCenterLatF = 70. ;This is necessary to fix the wrong value on the WRF file.
opts_cs@mpDataBaseVersion = "HighRes" ;fill land in gray and ocean in transparent
opts_cs@mpAreaMaskingOn = "True"
opts_cs@mpFillOn = False ; True
; VECTOR
vecres = True ;vector only resources
vecres@gsnAddCyclic = False
vecres@gsnDraw = False ;don't draw
vecres@gsnFrame = False ;don't advance frame
vecres@vcGlyphStyle = "LineArrow" ;"FillArrow"
;vecres@vcLineArrowThicknessF = 3 ;change vector thickness
vecres@vcRefLengthF = 0.05 ;define length of vec ref
vecres@vcRefAnnoOrthogonalPosF=-1.0
vecres@vcRefAnnoParallelPosF =1.2
;vecres@vcLabelFontHeightF = 0.08
vecres@gsnRightString = " " ;turn off right string
vecres@gsnLeftString = " " ;turn off left string
vecres@tiXAxisString = " " ;turn off axis label
vecres@vcRefMagnitudeF = 500 ;define vector ref mag
vecres@vcRefAnnoOrthogonalPosF = -0.36 ;move ref vector into plot
xskp=6
yskp=3
; Add a box showing the area of interest
lnres = True
lnres@gsLineColor = "Red"
lnres@gsLineThicknessF = 2.0
res1=opts_cs
res1@cnLinesOn = False
res1@cnLevelSelectionMode = "ExplicitLevels" ;
clev1=(/-2.0,-1.6,-1.2,-0.8,-0.4,0.0,0.4,0.8,1.2,1.6,2.0/)
res1@cnLevels = clev1
res2=opts_cs
res2@gsnLeftString = "CNV [10~S~-4~N~kgm~S~-2~N~s~S~-1~N~]"
res2@gsnCenterString = ""
res2@gsnRightString = ""
res2@cnLinesOn = False
res2@cnLevelSelectionMode = "ExplicitLevels" ;
clev2=clev1
res2@cnLevels = clev2
res22=ots_cs_nomap
res22@cnLinesOn = False
res22@cnLevelSelectionMode = "ExplicitLevels" ;
res22@cnLevels = clev2
res22@lbLabelBarOn = False
res22@cnInfoLabelOn = False
nt=0
do it=nstart,nend,nstep
print("")
print("it="+it+" "+date_str(it))
templete= scriptname+"_"+inname
figfile = outdir+"/"+ templete + "_" + date_num_def(it)
wks = gsn_open_wks(figtype, figfile)
res1@vpXF = 0.0
res2@vpXF = 0.0
itp=it
res1@gsnLeftString = "CNV [10~S~-4~N~kgm~S~-2~N~s~S~-1~N~]" +itp
res1@gsnRightString = date_str(itp)
res1@vpYF = 0.85
gsn_define_colormap(wks,"BlueWhiteOrangeRed")
cmf2s(:,:)=cmf2(itp,:,:)
cmf2s@_FillValue=default_fillvalue(typeof(cmf2s))
;printVarSummary(cmf2s)
plot1=gsn_csm_contour_map(wks,cmf2s(:,:),res1)
print("max(cmf2s)=" + max(cmf2s) + " min(cmf2s)=" +min(cmf2s))
draw(plot1)
delete(plot1)
res2@gsnLeftString = "CNV [10~S~-4~N~ kg m~S~-2~N~ s~S~-1~N~]"
res2@gsnRightString = date_str_def(itp)
res2@vpYF = 0.5
gsn_define_colormap(wks,"BlueWhiteOrangeRed")
cmf2s_sub= where(lat2d.lt.min_lat.or.lat2d.gt.max_lat .or. \
lon2d.lt.min_lon.or.lon2d.gt.max_lon, \
cmf2s@_FillValue, cmf2s)
cmf2s_sub2=where(lat2d.lt.min_lat2.or.lat2d.gt.max_lat2 .or. \
lon2d.lt.min_lon2.or.lon2d.gt.max_lon2, \
cmf2s@_FillValue, cmf2s)
cmf2s_sub@_FillValue=cmf2s@_FillValue
cmf2s_sub2@_FillValue=cmf2s@_FillValue
cmf2s_ave(it) = avg(cmf2s_sub)+avg(cmf2s_sub2)
plot2=gsn_csm_contour_map(wks,cmf2s_sub(:,:),res2)
plot22=gsn_csm_contour(wks,cmf2s_sub2(:,:),res22)
overlay(plot2,plot22)
dum4 = gsn_add_polyline(wks,plot2,(/min_lon,max_lon,\
max_lon,min_lon,min_lon/),\
(/min_lat,min_lat,max_lat,\
max_lat,min_lat/),lnres)
dum5 = gsn_add_polyline(wks,plot2,(/min_lon2,max_lon2,\
max_lon2,min_lon2,min_lon2/),\
(/min_lat2,min_lat2,max_lat2,\
max_lat2,min_lat2/),lnres)
draw(plot2)
delete(plot2)
;Header
txres = True
txres@txJust="CenterLeft"
txres@txFontHeightF = 0.010
lst = systemfunc("ls -lh "+infile)
strs = str_split(lst, " ")
str1 = "Input: "+strs(8)
str2 = strs(0)+" "+strs(1)+" "+strs(2)+" "+strs(3)+" "+strs(4)+" "+\
strs(5)+" "+strs(6)+" "+strs(7)
gsn_text_ndc(wks,str1,0.05,0.99,txres)
gsn_text_ndc(wks,str2,0.05,0.97,txres)
txres@txJust="CenterCenter"
txres@txFontHeightF = 0.016
text=scriptname
gsn_text_ndc(wks,text,0.5,0.93,txres)
txres@txFontHeightF = 0.016
text=date_str_def(it)
gsn_text_ndc(wks,text,0.5,0.90,txres)
;Footer
txres@txFontHeightF = 0.01
txres@txJust="CenterLeft"
text = systemfunc("export LANC=C; date")
gsn_text_ndc(wks,text, 0.05,0.055,txres)
cwd =systemfunc("pwd")
gsn_text_ndc(wks,"Current dir: "+cwd, 0.05,0.040,txres)
gsn_text_ndc(wks,"Script: "+scriptname, 0.05,0.025,txres)
gsn_text_ndc(wks,"Output dir: "+outdir, 0.05,0.010,txres)
frame(wks)
end do ;it
cmf2s_ave_kgm2s2=cmf2s_ave
alist= [/datetime_out, cmf2s_ave, cmf2s_ave_kgm2s2/]
;print(alist)
;
; TIME SERIES
;
templete= scriptname+"_"+inname
figfile = outdir+"/"+ templete + "_Time.Series"
wks = gsn_open_wks(figtype, figfile)
res = True ;plot mods desired
res@gsnFrame = False ;don't advance frame yet
res@gsnDraw = False ;don't draw plot
res@xyLineThicknessF = 2.0 ;make thicker
res@gsnLeftString=""
res@gsnRightString=""
;printVarSummary(time)
restick = True
restick@ttmFormat = "%h %c %d" ;"%N/%D %H:%M"
time_axis_labels(time,res,restick)
plot = gsn_csm_xy(wks,time,cmf2s_ave,res)
draw(plot)
delete(res)
;Header
txres = True
txres@txJust="CenterLeft"
txres@txFontHeightF = 0.010
lst = systemfunc("ls -lh "+infile)
strs = str_split(lst, " ")
str1 = "Input: "+strs(8)
str2 = strs(0)+" "+strs(1)+" "+strs(2)+" "+strs(3)+" "+strs(4)+" "+\
strs(5)+" "+strs(6)+" "+strs(7)
gsn_text_ndc(wks,str1,0.05,0.99,txres)
gsn_text_ndc(wks,str2,0.05,0.97,txres)
txres@txJust="CenterCenter"
txres@txFontHeightF = 0.016
text=scriptname
gsn_text_ndc(wks,text,0.5,0.93,txres)
;txres@txFontHeightF = 0.016
;text=chartostring(times0)
;gsn_text_ndc(wks,text,0.5,0.90,txres)
;Footer
txres@txFontHeightF = 0.01
txres@txJust="CenterLeft"
today = systemfunc("export LANC=C; date")
gsn_text_ndc(wks,today, 0.05,0.010,txres)
cwd =systemfunc("pwd")
gsn_text_ndc(wks,"Current dir: "+cwd, 0.05,0.025,txres)
gsn_text_ndc(wks,"Output dir: "+outdir,0.05,0.040,txres)
frame(wks)
;
; ASCII OUT
;
min_lat=71.2
max_lat=76.5
min_lon=-6.0
max_lon=25.5
min_lat2=71.2
max_lat2=78.0
min_lon2=25.501
max_lon2=40.0
header_ofle = (/\
"# "+today,\
"# Input : "+infile,\
"# CWD : "+cwd,\
"# ave area : "+min_lat+" - "+max_lat+" "+min_lon+" - "+max_lon,\
"# ave area2: "+min_lat2+" - "+max_lat2+" "+min_lon2+" - "+max_lon2\
/)
hlist=[/header_ofle/]
ofle=scriptname+"_"+inname+"_Time.Series.txt"
write_table(ofle, "w", hlist, "%s")
write_table(ofle, "a", alist, "%s %11.4f %12.4f")
print("")
print("Output: "+figfile)
print("Output: "+ofle)
print("")
end
----------------------
End of ERA-I.WVF.ncl
----------------------
======================
ERA-I.IWV.run.sh
======================
#!/bin/sh
exe=./runncl.sh
ncl=$(basename $0 .run.sh).ncl
if [ ! -f $ncl ]; then
echo Error in $0 : No such file, $ncl
exit 1
fi
indir="."
dirs="$indir"
for dir in $dirs; do
if [ ! -d $dir ]; then
echo Error in $0 : No such directory, $dir
exit 1
fi
done
input="ERA-I_IWV_110101-31.nc"
if [ ! -f $indir/$input ];then
echo Error in $0: No such file, $indir/$input
exit 1
fi
export LANG=C
$exe $ncl "$indir" "$input"
exit 0
----------------------
End of ERA-I.IWV.run.sh
----------------------
======================
ERA-I.IWV.ncl
======================
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
begin
month_abbr = (/"","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep", \
"Oct","Nov","Dec"/)
wallClock1 = systemfunc("date") ; retrieve wall clock time
scriptname_i = getenv("NCL_ARG_1")
indir = getenv("NCL_ARG_2")
input = getenv("NCL_ARG_3")
print("")
print(scriptname_i + " is running ...")
scriptname=systemfunc("basename " + scriptname_i + " .ncl")
inname=systemfunc("basename "+ input + " .nc")
outdir="Fig"+"_"+scriptname
system("mkdir -vp "+outdir)
infile=indir+"/"+input
a=addfile(infile,"r")
latmin=40.0
latmax=80.0
lonmin=-30.0
lonmax=90.0
min_lat=71.2
max_lat=76.5
min_lon=-6.0
max_lon=25.5
min_lat2=71.2
max_lat2=78.0
min_lon2=25.501
max_lon2=40.0
m2mm=1000.0
scl=1.E4
hr2s=3600.0
;lon=a->longitude
;lat=a->latitude
time=a->time
time@units="hours since 1900-01-01 00:00:0.0"
mt=dimsizes(time)
;print(mt)
nstart=0
nend=mt-1; 6
nstep=1
iwv=a->$"p55.162"$
; p55.162:units = "kg m**-2" ;
; p55.162:long_name = "Vertical integral of water vapour"
tcwv=a->tcwv
; tcwv:units = "kg m**-2" ;
; tcwv:long_name = "Total column water vapour" ;
lat=a->latitude
lon=a->longitude
iwv@_FillValue=default_fillvalue(typeof(iwv))
iwv2 =(iwv * iwv@scale_factor + iwv@add_offset)
tcwv@_FillValue=default_fillvalue(typeof(tcwv))
tcwv2 =(tcwv * tcwv@scale_factor + tcwv@add_offset)
;
; Vertical integral of water vapour, kg m**-2
;
iwv2@units="kg m**-2"
iwv2!0="time"
iwv2!1="lat"
iwv2!2="lon"
iwv2&time=time
iwv2&lat=lat
iwv2&lon=lon
iwv2&lat@units="degrees_north"
iwv2&lon@units="degrees_east"
dim=dimsizes(iwv2)
ny=dim(1)
nx=dim(2)
iwv2s=new((/ny,nx/),typeof(iwv2))
iwv2s!0="lat"
iwv2s!1="lon"
iwv2s&lat=lat
iwv2s&lon=lon
iwv2s&lat@units="degrees_north"
iwv2s&lon@units="degrees_east"
iwv2s_sub=new((/ny,nx/),typeof(iwv2))
iwv2s_sub!0="lat"
iwv2s_sub!1="lon"
iwv2s_sub&lat=lat
iwv2s_sub&lon=lon
iwv2s_sub&lat@units="degrees_north"
iwv2s_sub&lon@units="degrees_east"
iwv2s_sub2=new((/ny,nx/),typeof(iwv2))
iwv2s_sub2!0="lat"
iwv2s_sub2!1="lon"
iwv2s_sub2&lat=lat
iwv2s_sub2&lon=lon
iwv2s_sub2&lat@units="degrees_north"
iwv2s_sub2&lon@units="degrees_east"
;
; Total column water vapour, kg m**-2
;
tcwv2@units="kg m**-2"
tcwv2!0="time"
tcwv2!1="lat"
tcwv2!2="lon"
tcwv2&time=time
tcwv2&lat=lat
tcwv2&lon=lon
tcwv2&lat@units="degrees_north"
tcwv2&lon@units="degrees_east"
dim=dimsizes(tcwv2)
ny=dim(1)
nx=dim(2)
tcwv2s=new((/ny,nx/),typeof(tcwv2))
tcwv2s!0="lat"
tcwv2s!1="lon"
tcwv2s&lat=lat
tcwv2s&lon=lon
tcwv2s&lat@units="degrees_north"
tcwv2s&lon@units="degrees_east"
tcwv2s_sub=new((/ny,nx/),typeof(tcwv2))
tcwv2s_sub!0="lat"
tcwv2s_sub!1="lon"
tcwv2s_sub&lat=lat
tcwv2s_sub&lon=lon
tcwv2s_sub&lat@units="degrees_north"
tcwv2s_sub&lon@units="degrees_east"
tcwv2s_sub2=new((/ny,nx/),typeof(tcwv2))
tcwv2s_sub2!0="lat"
tcwv2s_sub2!1="lon"
tcwv2s_sub2&lat=lat
tcwv2s_sub2&lon=lon
tcwv2s_sub2&lat@units="degrees_north"
tcwv2s_sub2&lon@units="degrees_east"
lat2d=new((/ny,nx/),typeof(lat))
lon2d=new((/ny,nx/),typeof(lon))
do j=0,ny-1
do i=0,nx-1
lat2d(j,i)=lat(j)
lon2d(j,i)=lon(i)
end do
end do
;printVarSummary(cmf2s)
iwv2s_ave=new((/mt/),typeof(iwv2))
tcwv2s_ave=new((/mt/),typeof(tcwv2))
;printVarSummary(cmf2s_ave)
dthr=time(1)-time(0)
dthr2=tofloat(dthr)*2.0
dt=dthr*3600.0
dt2=tofloat(dt*2.0)
;print(dt2)
utc_date = cd_calendar(time, 0)
year = tointeger(utc_date(:,0)) ; Convert to integer for
month = tointeger(utc_date(:,1)) ; use sprinti
day = tointeger(utc_date(:,2))
hour = tointeger(utc_date(:,3))
minute = tointeger(utc_date(:,4))
second = utc_date(:,5)
date_num = sprinti("%0.4i", year)+sprinti("%0.2i", month)+sprinti("%0.2i", day) \
+ sprinti("_%0.2i", hour)
date_str = sprinti("%0.2i", hour)+":00UTC"+sprinti("%0.2i", day)+month_abbr(month)\
+sprinti("%0.4i", year)
time_def=time-0
time_def@units="hours since 1900-01-01 00:00:0.0"
utc_date_def = cd_calendar(time_def, 0)
year = tointeger(utc_date_def(:,0)) ; Convert to integer for
month = tointeger(utc_date_def(:,1)) ; use sprinti
day = tointeger(utc_date_def(:,2))
hour = tointeger(utc_date_def(:,3))
minute = tointeger(utc_date_def(:,4))
second = utc_date_def(:,5)
date_num_def = sprinti("%0.4i", year)+sprinti("%0.2i", month)+sprinti("%0.2i", day) \
+ sprinti("_%0.2i", hour)
date_str_def = sprinti("%0.2i", hour)+":00UTC"+sprinti("%0.2i", day)+month_abbr(month)\
+sprinti("%0.4i", year)
; For Generic mapping tools
datetime_out = sprinti("%0.4i", year)+"-"+sprinti("%0.2i", month)+"-"+\
sprinti("%0.2i", day )+"T" + sprinti("%0.2i", hour)+\
":00:00"
;print(date_str)
figtype="ps" ;"x11"
opts=True
;Coordinates and domain
opts@gsnAddCyclic = False
;Plotting area
opts@vpHeightF = .2 ;1
;opts@vpWidthF = opts@vpHeightF*tofloat(mlon)/tofloat(nlat) ;.4
;Page control
opts@gsnFrame = False
opts@gsnDraw = False
opts@gsnStringFontHeightF = 0.012
opts@lbLabelFontHeightF = 0.012
; Contour line
opts_cn=opts
opts_cn@cnFillOn = False
opts_cn@cnLevelSpacingF = 1
opts_cn@cnLineThicknessF= 1.0
opts_cn@cnInfoLabelOn = False
opts_cn@cnLabelDrawOrder="PostDraw"
opts_cn@cnLineLabelFontHeightF = 0.01 ;8
; Color shade
opts_cs=opts
opts_cs@cnFillOn = True
opts_cs@pmLabelBarOrthogonalPosF = 0.0
opts_cs@pmLabelBarHeightF = 0.04
; Color bar
opts_cs@pmLabelBarOrthogonalPosF = .10 ;move whole thing down
opts_cs@lbTitlePosition = "Right" ;title position
opts_cs@lbTitleFontHeightF= .012 ;make title smaller
opts_cs@lbTitleDirection = "Across" ;title direction
ots_cs_nomap=opts_cs
; MAP
opts_cs@mpProjection = "Stereographic"
opts_cs@mpLimitMode = "Corners"
opts_cs@mpLeftCornerLatF = latmin
opts_cs@mpLeftCornerLonF = lonmin
opts_cs@mpRightCornerLatF = latmax
opts_cs@mpRightCornerLonF = lonmax
opts_cs@mpCenterLonF = 20.
opts_cs@mpCenterLatF = 70. ;This is necessary to fix the wrong value on the WRF file.
opts_cs@mpDataBaseVersion = "HighRes" ;fill land in gray and ocean in transparent
opts_cs@mpAreaMaskingOn = "True"
opts_cs@mpFillOn = False ; True
; VECTOR
vecres = True ;vector only resources
vecres@gsnAddCyclic = False
vecres@gsnDraw = False ;don't draw
vecres@gsnFrame = False ;don't advance frame
vecres@vcGlyphStyle = "LineArrow" ;"FillArrow"
;vecres@vcLineArrowThicknessF = 3 ;change vector thickness
vecres@vcRefLengthF = 0.05 ;define length of vec ref
vecres@vcRefAnnoOrthogonalPosF=-1.0
vecres@vcRefAnnoParallelPosF =1.2
;vecres@vcLabelFontHeightF = 0.08
vecres@gsnRightString = " " ;turn off right string
vecres@gsnLeftString = " " ;turn off left string
vecres@tiXAxisString = " " ;turn off axis label
vecres@vcRefMagnitudeF = 500 ;define vector ref mag
vecres@vcRefAnnoOrthogonalPosF = -0.36 ;move ref vector into plot
xskp=6
yskp=3
; Add a box showing the area of interest
lnres = True
lnres@gsLineColor = "Red"
lnres@gsLineThicknessF = 2.0
res1=opts_cs
res1@cnLinesOn = False
res1@cnLevelSelectionMode = "ExplicitLevels" ;
clev1=(/4, 8, 12, 16, 20, 24, 28, 32/)
res1@cnLevels = clev1
res2=opts_cs
res2@gsnLeftString = "IWV [kgm~S~-2~N]"
res2@gsnCenterString = ""
res2@gsnRightString = ""
res2@cnLinesOn = False
res2@cnLevelSelectionMode = "ExplicitLevels" ;
clev2=clev1
res2@cnLevels = clev2
res22=ots_cs_nomap
res22@cnLinesOn = False
res22@cnLevelSelectionMode = "ExplicitLevels" ;
res22@cnLevels = clev2
res22@lbLabelBarOn = False
res22@cnInfoLabelOn = False
nt=0
do it=nstart,nend,nstep
print("")
print("it="+it+" "+date_str(it))
itp=it
iwv2s(:,:)=iwv2(itp,:,:)
iwv2s@_FillValue=default_fillvalue(typeof(iwv2s))
;printVarSummary(iwv2s)
print("max(iwv2s)=" + max(iwv2s) + " min(iwv2s)=" +min(iwv2s))
iwv2s_sub= where(lat2d.lt.min_lat.or.lat2d.gt.max_lat .or. \
lon2d.lt.min_lon.or.lon2d.gt.max_lon, \
iwv2s@_FillValue, iwv2s)
iwv2s_sub2=where(lat2d.lt.min_lat2.or.lat2d.gt.max_lat2 .or. \
lon2d.lt.min_lon2.or.lon2d.gt.max_lon2, \
iwv2s@_FillValue, iwv2s)
iwv2s_sub@_FillValue=iwv2s@_FillValue
iwv2s_sub2@_FillValue=iwv2s@_FillValue
iwv2s_ave(it) = avg(iwv2s_sub)+avg(iwv2s_sub2)
tcwv2s(:,:)=tcwv2(itp,:,:)
tcwv2s@_FillValue=default_fillvalue(typeof(tcwv2s))
;printVarSummary(tcwv2s)
print("max(tcwv2s)=" + max(tcwv2s) + " min(tcwv2s)=" +min(tcwv2s))
tcwv2s_sub= where(lat2d.lt.min_lat.or.lat2d.gt.max_lat .or. \
lon2d.lt.min_lon.or.lon2d.gt.max_lon, \
tcwv2s@_FillValue, tcwv2s)
tcwv2s_sub2=where(lat2d.lt.min_lat2.or.lat2d.gt.max_lat2 .or. \
lon2d.lt.min_lon2.or.lon2d.gt.max_lon2, \
tcwv2s@_FillValue, tcwv2s)
tcwv2s_sub@_FillValue=tcwv2s@_FillValue
tcwv2s_sub2@_FillValue=tcwv2s@_FillValue
tcwv2s_ave(it) = avg(tcwv2s_sub)+avg(tcwv2s_sub2)
templete= scriptname+"_"+inname
figfile = outdir+"/"+ templete + "_" + date_num_def(it)
wks = gsn_open_wks(figtype, figfile)
res1@vpXF = 0.0
res2@vpXF = 0.0
res1@gsnLeftString = "IWV [kgm~S~-2~N~]" +itp
res1@gsnRightString = date_str(itp)
res1@vpYF = 0.85
gsn_define_colormap(wks,"WhiteBlueGreenYellowRed");"BlueWhiteOrangeRed")
plot1=gsn_csm_contour_map(wks,iwv2s(:,:),res1)
draw(plot1)
delete(plot1)
res2@gsnLeftString = "TCWV [kgm~S~-2~N~]"
res2@gsnRightString = date_str_def(itp)
res2@vpYF = 0.5
gsn_define_colormap(wks,"WhiteBlueGreenYellowRed");"BlueWhiteOrangeRed")
plot2=gsn_csm_contour_map(wks,tcwv2s_sub(:,:),res2)
plot22=gsn_csm_contour(wks,tcwv2s_sub2(:,:),res22)
overlay(plot2,plot22)
dum4 = gsn_add_polyline(wks,plot2,(/min_lon,max_lon,\
max_lon,min_lon,min_lon/),\
(/min_lat,min_lat,max_lat,\
max_lat,min_lat/),lnres)
dum5 = gsn_add_polyline(wks,plot2,(/min_lon2,max_lon2,\
max_lon2,min_lon2,min_lon2/),\
(/min_lat2,min_lat2,max_lat2,\
max_lat2,min_lat2/),lnres)
draw(plot2)
delete(plot2)
;Header
txres = True
txres@txJust="CenterLeft"
txres@txFontHeightF = 0.010
lst = systemfunc("ls -lh "+infile)
strs = str_split(lst, " ")
str1 = "Input: "+strs(8)
str2 = strs(0)+" "+strs(1)+" "+strs(2)+" "+strs(3)+" "+strs(4)+" "+\
strs(5)+" "+strs(6)+" "+strs(7)
gsn_text_ndc(wks,str1,0.05,0.99,txres)
gsn_text_ndc(wks,str2,0.05,0.97,txres)
txres@txJust="CenterCenter"
txres@txFontHeightF = 0.016
text=scriptname
gsn_text_ndc(wks,text,0.5,0.93,txres)
txres@txFontHeightF = 0.016
text=date_str_def(it)
gsn_text_ndc(wks,text,0.5,0.90,txres)
;Footer
txres@txFontHeightF = 0.01
txres@txJust="CenterLeft"
text = systemfunc("export LANC=C; date")
gsn_text_ndc(wks,text, 0.05,0.055,txres)
cwd =systemfunc("pwd")
gsn_text_ndc(wks,"Current dir: "+cwd, 0.05,0.040,txres)
gsn_text_ndc(wks,"Script: "+scriptname, 0.05,0.025,txres)
gsn_text_ndc(wks,"Output dir: "+outdir, 0.05,0.010,txres)
frame(wks)
end do ;it
alist= [/datetime_out, iwv2s_ave, tcwv2s_ave/]
;print(alist)
;
; TIME SERIES
;
templete= scriptname+"_"+inname
figfile = outdir+"/"+ templete + "_Time.Series"
wks = gsn_open_wks(figtype, figfile)
res = True ;plot mods desired
res@gsnFrame = False ;don't advance frame yet
res@gsnDraw = False ;don't draw plot
res@xyLineThicknessF = 2.0 ;make thicker
res@gsnLeftString=""
res@gsnRightString=""
;printVarSummary(time)
restick = True
restick@ttmFormat = "%h %c %d" ;"%N/%D %H:%M"
time_axis_labels(time,res,restick)
plot = gsn_csm_xy(wks,time,iwv2s_ave,res)
draw(plot)
delete(res)
;Header
txres = True
txres@txJust="CenterLeft"
txres@txFontHeightF = 0.010
lst = systemfunc("ls -lh "+infile)
strs = str_split(lst, " ")
str1 = "Input: "+strs(8)
str2 = strs(0)+" "+strs(1)+" "+strs(2)+" "+strs(3)+" "+strs(4)+" "+\
strs(5)+" "+strs(6)+" "+strs(7)
gsn_text_ndc(wks,str1,0.05,0.99,txres)
gsn_text_ndc(wks,str2,0.05,0.97,txres)
txres@txJust="CenterCenter"
txres@txFontHeightF = 0.016
text=scriptname
gsn_text_ndc(wks,text,0.5,0.93,txres)
;txres@txFontHeightF = 0.016
;text=chartostring(times0)
;gsn_text_ndc(wks,text,0.5,0.90,txres)
;Footer
txres@txFontHeightF = 0.01
txres@txJust="CenterLeft"
today = systemfunc("export LANC=C; date")
gsn_text_ndc(wks,today, 0.05,0.010,txres)
cwd =systemfunc("pwd")
gsn_text_ndc(wks,"Current dir: "+cwd, 0.05,0.025,txres)
gsn_text_ndc(wks,"Output dir: "+outdir,0.05,0.040,txres)
frame(wks)
;
; ASCII OUT
;
min_lat=71.2
max_lat=76.5
min_lon=-6.0
max_lon=25.5
min_lat2=71.2
max_lat2=78.0
min_lon2=25.501
max_lon2=40.0
header_ofle = (/\
"# "+today,\
"# Input : "+infile,\
"# CWD : "+cwd,\
"# ave area : "+min_lat+" - "+max_lat+" "+min_lon+" - "+max_lon,\
"# ave area2: "+min_lat2+" - "+max_lat2+" "+min_lon2+" - "+max_lon2\
/)
hlist=[/header_ofle/]
ofle=scriptname+"_"+inname+"_Time.Series.txt"
write_table(ofle, "w", hlist, "%s")
write_table(ofle, "a", alist, "%s %11.4f %12.4f")
print("")
print("Output: "+figfile)
print("Output: "+ofle)
print("")
end
----------------------
End of ERA-I.IWV.ncl
----------------------
======================
ERA-I.Precip.run.sh
======================
#!/bin/sh
exe=./runncl.sh
ncl=$(basename $0 .run.sh).ncl
if [ ! -f $ncl ]; then
echo Error in $0 : No such file, $ncl
exit 1
fi
indir="."
dirs="$indir"
for dir in $dirs; do
if [ ! -d $dir ]; then
echo Error in $0 : No such directory, $dir
exit 1
fi
done
input="ERA-I_E-P_110101-31.nc"
if [ ! -f $indir/$input ];then
echo Error in $0: No such file, $indir/$input
exit 1
fi
export LANG=C
$exe $ncl "$indir" "$input"
exit 0
----------------------
End of ERA-I.Precip.run.sh
----------------------
======================
ERA-I.Precip.ncl
======================
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
begin
month_abbr = (/"","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep", \
"Oct","Nov","Dec"/)
wallClock1 = systemfunc("date") ; retrieve wall clock time
scriptname_i = getenv("NCL_ARG_1")
indir = getenv("NCL_ARG_2")
input = getenv("NCL_ARG_3")
print("")
print(scriptname_i + " is running ...")
scriptname=systemfunc("basename " + scriptname_i + " .ncl")
inname=systemfunc("basename "+ input + " .nc")
outdir="Fig"+"_"+scriptname
system("mkdir -vp "+outdir)
infile=indir+"/"+input
a=addfile(infile,"r")
latmin=40.0
latmax=80.0
lonmin=-30.0
lonmax=90.0
min_lat=71.2
max_lat=76.5
min_lon=-6.0
max_lon=25.5
min_lat2=71.2
max_lat2=78.0
min_lon2=25.501
max_lon2=40.0
m2mm=1000.0
scl=1.E4
hr2s=3600.0
;lon=a->longitude
;lat=a->latitude
time=a->time
time@units="hours since 1900-01-01 00:00:0.0"
mt=dimsizes(time)
;print(mt)
nstart=0
nend=mt-1; 6
nstep=1
tp=a->tp ;Total Precipitation, m
lat=a->latitude
lon=a->longitude
tp@_FillValue=default_fillvalue(typeof(tp))
;print(tp@_FillValue)
tp2 =(tp * tp@scale_factor + tp@add_offset)*m2mm ;downward (m->mm)
;;;e2 = e2/hr2s*scl
tp2@units="mm"
tp2!0="time"
tp2!1="lat"
tp2!2="lon"
tp2&time=time
tp2&lat=lat
tp2&lon=lon
tp2&lat@units="degrees_north"
tp2&lon@units="degrees_east"
;print("max(e2)="+max(e2)+" min(e2)="+min(e2))
dim=dimsizes(tp2)
ny=dim(1)
nx=dim(2)
tp2s=new((/ny,nx/),typeof(tp2))
tp2s!0="lat"
tp2s!1="lon"
tp2s&lat=lat
tp2s&lon=lon
tp2s&lat@units="degrees_north"
tp2s&lon@units="degrees_east"
tp2s_sub=new((/ny,nx/),typeof(tp2))
tp2s_sub!0="lat"
tp2s_sub!1="lon"
tp2s_sub&lat=lat
tp2s_sub&lon=lon
tp2s_sub&lat@units="degrees_north"
tp2s_sub&lon@units="degrees_east"
tp2s_sub2=new((/ny,nx/),typeof(tp2))
tp2s_sub2!0="lat"
tp2s_sub2!1="lon"
tp2s_sub2&lat=lat
tp2s_sub2&lon=lon
tp2s_sub2&lat@units="degrees_north"
tp2s_sub2&lon@units="degrees_east"
lat2d=new((/ny,nx/),typeof(lat))
lon2d=new((/ny,nx/),typeof(lon))
;printVarSummary(lon)
;printVarSummary(lat)
;printVarSummary(lon2d)
;printVarSummary(lat2d)
do j=0,ny-1
do i=0,nx-1
lat2d(j,i)=lat(j)
lon2d(j,i)=lon(i)
end do
end do
;printVarSummary(e2s)
tp2s_ave=new((/mt/),typeof(tp2))
;printVarSummary(tp2s_ave)
dthr=time(1)-time(0)
dthr2=tofloat(dthr)*2.0
dt=dthr*3600.0
dt2=tofloat(dt*2.0)
;print(dt2)
utc_date = cd_calendar(time, 0)
year = tointeger(utc_date(:,0)) ; Convert to integer for
month = tointeger(utc_date(:,1)) ; use sprinti
day = tointeger(utc_date(:,2))
hour = tointeger(utc_date(:,3))
minute = tointeger(utc_date(:,4))
second = utc_date(:,5)
date_num = sprinti("%0.4i", year)+sprinti("%0.2i", month)+sprinti("%0.2i", day) \
+ sprinti("_%0.2i", hour)
date_str = sprinti("%0.2i", hour)+":00UTC"+sprinti("%0.2i", day)+month_abbr(month)\
+sprinti("%0.4i", year)
time_def=time-6
time_def@units="hours since 1900-01-01 00:00:0.0"
utc_date_def = cd_calendar(time_def, 0)
year = tointeger(utc_date_def(:,0)) ; Convert to integer for
month = tointeger(utc_date_def(:,1)) ; use sprinti
day = tointeger(utc_date_def(:,2))
hour = tointeger(utc_date_def(:,3))
minute = tointeger(utc_date_def(:,4))
second = utc_date_def(:,5)
date_num_def = sprinti("%0.4i", year)+sprinti("%0.2i", month)+sprinti("%0.2i", day) \
+ sprinti("_%0.2i", hour)
date_str_def = sprinti("%0.2i", hour)+":00UTC"+sprinti("%0.2i", day)+month_abbr(month)\
+sprinti("%0.4i", year)
; For Generic mapping tools
datetime_out = sprinti("%0.4i", year)+"-"+sprinti("%0.2i", month)+"-"+\
sprinti("%0.2i", day )+"T" + sprinti("%0.2i", hour)+\
":00:00"
;print(date_str)
figtype="ps" ;"x11"
opts=True
;Coordinates and domain
opts@gsnAddCyclic = False
;Plotting area
opts@vpHeightF = .2 ;1
;opts@vpWidthF = opts@vpHeightF*tofloat(mlon)/tofloat(nlat) ;.4
;Page control
opts@gsnFrame = False
opts@gsnDraw = False
opts@gsnStringFontHeightF = 0.012
opts@lbLabelFontHeightF = 0.012
; Contour line
opts_cn=opts
opts_cn@cnFillOn = False
opts_cn@cnLevelSpacingF = 1
opts_cn@cnLineThicknessF= 1.0
opts_cn@cnInfoLabelOn = False
opts_cn@cnLabelDrawOrder="PostDraw"
opts_cn@cnLineLabelFontHeightF = 0.01 ;8
; Color shade
opts_cs=opts
opts_cs@cnFillOn = True
opts_cs@pmLabelBarOrthogonalPosF = 0.0
opts_cs@pmLabelBarHeightF = 0.04
; Color bar
opts_cs@pmLabelBarOrthogonalPosF = .10 ;move whole thing down
opts_cs@lbTitlePosition = "Right" ;title position
opts_cs@lbTitleFontHeightF= .012 ;make title smaller
opts_cs@lbTitleDirection = "Across" ;title direction
ots_cs_nomap=opts_cs
; MAP
opts_cs@mpProjection = "Stereographic"
opts_cs@mpLimitMode = "Corners"
opts_cs@mpLeftCornerLatF = latmin
opts_cs@mpLeftCornerLonF = lonmin
opts_cs@mpRightCornerLatF = latmax
opts_cs@mpRightCornerLonF = lonmax
opts_cs@mpCenterLonF = 20.
opts_cs@mpCenterLatF = 70. ;This is necessary to fix the wrong value on the WRF file.
opts_cs@mpDataBaseVersion = "HighRes" ;fill land in gray and ocean in transparent
opts_cs@mpAreaMaskingOn = "True"
opts_cs@mpFillOn = False ; True
; VECTOR
vecres = True ;vector only resources
vecres@gsnAddCyclic = False
vecres@gsnDraw = False ;don't draw
vecres@gsnFrame = False ;don't advance frame
vecres@vcGlyphStyle = "LineArrow" ;"FillArrow"
;vecres@vcLineArrowThicknessF = 3 ;change vector thickness
vecres@vcRefLengthF = 0.05 ;define length of vec ref
vecres@vcRefAnnoOrthogonalPosF=-1.0
vecres@vcRefAnnoParallelPosF =1.2
;vecres@vcLabelFontHeightF = 0.08
vecres@gsnRightString = " " ;turn off right string
vecres@gsnLeftString = " " ;turn off left string
vecres@tiXAxisString = " " ;turn off axis label
vecres@vcRefMagnitudeF = 500 ;define vector ref mag
vecres@vcRefAnnoOrthogonalPosF = -0.36 ;move ref vector into plot
xskp=6
yskp=3
; Add a box showing the area of interest
lnres = True
lnres@gsLineColor = "Red"
lnres@gsLineThicknessF = 2.0
res1=opts_cs
res1@cnLinesOn = False
res1@cnLevelSelectionMode = "ExplicitLevels" ;
clev1=(/-0.5,0.0,0.5,1.0,1.5,2.0,2.5,3.0/)
res1@cnLevels = clev1
res2=opts_cs
res2@gsnLeftString = "E [mm/hr]"
res2@gsnCenterString = ""
res2@gsnRightString = ""
res2@cnLinesOn = False
res2@cnLevelSelectionMode = "ExplicitLevels" ;
clev2=(/-0.4,-0.3,-0.2,-0.1, 0.0, 0.1, 0.2, 0.3, 0.4/)
res2@cnLevels = clev2
res22=ots_cs_nomap
res22@cnLinesOn = False
res22@cnLevelSelectionMode = "ExplicitLevels" ;
clev2=(/-0.4,-0.3,-0.2,-0.1, 0.0, 0.1, 0.2, 0.3, 0.4/)
res22@cnLevels = clev2
res22@lbLabelBarOn = False
res22@cnInfoLabelOn = False
;printVarSummary(tp)
nt=0
do it=nstart,nend,nstep
print("")
print("it="+it+" "+date_str(it))
templete= scriptname+"_"+inname
figfile = outdir+"/"+ templete + "_" + date_num_def(it)
wks = gsn_open_wks(figtype, figfile)
res1@vpXF = 0.0
res2@vpXF = 0.0
itp=it
res1@gsnLeftString = "P [mm] " +itp
res1@gsnRightString = date_str(itp)
res1@vpYF = 0.85
gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")
plot1=gsn_csm_contour_map(wks,tp2(itp,:,:),res1)
print("max(tp2)=" + max(tp2(itp,:,:)) + " min(tp2)=" +min(tp2(itp,:,:)))
draw(plot1)
delete(plot1)
res2@gsnLeftString = "P [mm/h]"
res2@gsnRightString = date_str_def(itp)
res2@vpYF = 0.5
gsn_define_colormap(wks,"BlueWhiteOrangeRed")
tp2s(:,:)=tp2(itp,:,:)/dthr
tp2s_sub= where(lat2d.lt.min_lat.or.lat2d.gt.max_lat .or. \
lon2d.lt.min_lon.or.lon2d.gt.max_lon, \
lat2d@_FillValue, tp2s)
tp2s_sub2=where(lat2d.lt.min_lat2.or.lat2d.gt.max_lat2 .or. \
lon2d.lt.min_lon2.or.lon2d.gt.max_lon2, \
lat2d@_FillValue, tp2s)
tp2s_sub@_FillValue=tp2s@_FillValue
tp2s_sub2@_FillValue=tp2s@_FillValue
tp2s_ave(it) = avg(tp2s_sub)+avg(tp2s_sub2)
plot2=gsn_csm_contour_map(wks,tp2s_sub(:,:),res2)
plot22=gsn_csm_contour(wks,tp2s_sub2(:,:),res22)
overlay(plot2,plot22)
dum4 = gsn_add_polyline(wks,plot2,(/min_lon,max_lon,\
max_lon,min_lon,min_lon/),\
(/min_lat,min_lat,max_lat,\
max_lat,min_lat/),lnres)
dum5 = gsn_add_polyline(wks,plot2,(/min_lon2,max_lon2,\
max_lon2,min_lon2,min_lon2/),\
(/min_lat2,min_lat2,max_lat2,\
max_lat2,min_lat2/),lnres)
draw(plot2)
delete(plot2)
;Header
txres = True
txres@txJust="CenterLeft"
txres@txFontHeightF = 0.010
lst = systemfunc("ls -lh "+infile)
strs = str_split(lst, " ")
str1 = "Input: "+strs(8)
str2 = strs(0)+" "+strs(1)+" "+strs(2)+" "+strs(3)+" "+strs(4)+" "+\
strs(5)+" "+strs(6)+" "+strs(7)
gsn_text_ndc(wks,str1,0.05,0.99,txres)
gsn_text_ndc(wks,str2,0.05,0.97,txres)
txres@txJust="CenterCenter"
txres@txFontHeightF = 0.016
text=scriptname
gsn_text_ndc(wks,text,0.5,0.93,txres)
txres@txFontHeightF = 0.016
text=date_str_def(it)
gsn_text_ndc(wks,text,0.5,0.90,txres)
;Footer
txres@txFontHeightF = 0.01
txres@txJust="CenterLeft"
text = systemfunc("export LANC=C; date")
gsn_text_ndc(wks,text, 0.05,0.055,txres)
cwd =systemfunc("pwd")
gsn_text_ndc(wks,"Current dir: "+cwd, 0.05,0.040,txres)
gsn_text_ndc(wks,"Script: "+scriptname, 0.05,0.025,txres)
gsn_text_ndc(wks,"Output dir: "+outdir, 0.05,0.010,txres)
frame(wks)
end do ;it
tp2s_ave_kgm2s2=tp2s_ave/hr2s*scl
alist= [/datetime_out, tp2s_ave, tp2s_ave_kgm2s2/]
;print(alist)
;
; TIME SERIES
;
templete= scriptname+"_"+inname
figfile = outdir+"/"+ templete + "_Time.Series"
wks = gsn_open_wks(figtype, figfile)
res = True ;plot mods desired
res@gsnFrame = False ;don't advance frame yet
res@gsnDraw = False ;don't draw plot
res@xyLineThicknessF = 2.0 ;make thicker
res@gsnLeftString=""
res@gsnRightString=""
;printVarSummary(time)
restick = True
restick@ttmFormat = "%h %c %d" ;"%N/%D %H:%M"
time_axis_labels(time,res,restick)
plot = gsn_csm_xy(wks,time,tp2s_ave,res)
draw(plot)
delete(res)
;Header
txres = True
txres@txJust="CenterLeft"
txres@txFontHeightF = 0.010
lst = systemfunc("ls -lh "+infile)
strs = str_split(lst, " ")
str1 = "Input: "+strs(8)
str2 = strs(0)+" "+strs(1)+" "+strs(2)+" "+strs(3)+" "+strs(4)+" "+\
strs(5)+" "+strs(6)+" "+strs(7)
gsn_text_ndc(wks,str1,0.05,0.99,txres)
gsn_text_ndc(wks,str2,0.05,0.97,txres)
txres@txJust="CenterCenter"
txres@txFontHeightF = 0.016
text=scriptname
gsn_text_ndc(wks,text,0.5,0.93,txres)
;txres@txFontHeightF = 0.016
;text=chartostring(times0)
;gsn_text_ndc(wks,text,0.5,0.90,txres)
;Footer
txres@txFontHeightF = 0.01
txres@txJust="CenterLeft"
today = systemfunc("export LANC=C; date")
gsn_text_ndc(wks,today, 0.05,0.010,txres)
cwd =systemfunc("pwd")
gsn_text_ndc(wks,"Current dir: "+cwd, 0.05,0.025,txres)
gsn_text_ndc(wks,"Output dir: "+outdir,0.05,0.040,txres)
frame(wks)
;
; ASCII OUT
;
min_lat=71.2
max_lat=76.5
min_lon=-6.0
max_lon=25.5
min_lat2=71.2
max_lat2=78.0
min_lon2=25.501
max_lon2=40.0
header_ofle = (/\
"# "+today,\
"# Input : "+infile,\
"# CWD : "+cwd,\
"# ave area : "+min_lat+" - "+max_lat+" "+min_lon+" - "+max_lon,\
"# ave area2: "+min_lat2+" - "+max_lat2+" "+min_lon2+" - "+max_lon2\
/)
hlist=[/header_ofle/]
ofle=scriptname+"_"+inname+"_Time.Series.txt"
write_table(ofle, "w", hlist, "%s")
write_table(ofle, "a", alist, "%s %11.4f %12.4f")
print("")
print("Output: "+figfile)
print("Output: "+ofle)
print("")
end
----------------------
End of ERA-I.Precip.ncl
----------------------
======================
runncl.sh
======================
#!/bin/bash
#
# Universal wrapper script for ncl.
# Pass arguments from the command line to environment variables
#
# version 0.1, Thierry Corti, C2SM ETH Zurich
#
E_BADARGS=65
if [ ! -n "$1" ]
then
echo "Usage: `basename $0` script.ncl argument1 argument2 etc."
exit $E_BADARGS
fi
# save number of arguments to environment variable NCL_N_ARG
export NCL_N_ARGS=$#
# save command line arguments to environment variable NCL_ARG_#
for ((index=1; index<=$#; index++))
do
eval export NCL_ARG_$index=\$$index
done
command="ulimit -s unlimited"
echo ==========================================================
echo
echo $command
echo
echo Be careful! The above command will use lots of memory.
echo
echo ==========================================================
$command
# run ncl
ncl -Q -n $1
----------------------
End of runncl.sh
----------------------
残差項を求めるスクリプト
------------------------------
List of the following files:
------------------------------
calc_diff.sh
------------------------------
Machine info
------------------------------
aofd165.bio.mie-u.ac.jp
/work1/am/2016.PolarLow/ERA-I/Water.Vapor.Budget.ERA-I
======================
calc_diff.sh
======================
#!/bin/sh
export LANG=C
prefix=ERA-I
postfix=110101-31_Time.Series
var=WVF
var2=WVF
in1=$prefix"."$var"_"$prefix"_"${var2}_${postfix}.txt
var=Precip
var2=E-P
in2=$prefix"."$var"_"$prefix"_"${var2}_${postfix}.txt
var=Evapo
var2=E-P
in3=$prefix"."$var"_"$prefix"_"${var2}_${postfix}.txt
var=IWV_diff
var2=IWV_diff
in4=$prefix"."$var"_"$prefix"_"${var2}_${postfix}.txt
inlist="$in1 $in2 $in3 $in4"
for input in $inlist; do
if [ ! -f $input ]; then
echo Error in $0 : No such file, $input
exit 1
fi
done
out=ERA-I.Res_ERA-I_Res_110101-31_Time.Series.txt
i=0
while read BUF ; do
ary=(`echo $BUF`) # 配列に格納
if [ ${ary[0]} != "#" ]; then
datetime[$i]=${ary[0]}
term1[$i]=${ary[2]}
i=$(expr $i + 1)
fi
done < $in1
i=0
while read BUF ; do
ary=(`echo $BUF`) # 配列に格納
if [ ${ary[0]} != "#" ]; then
datetime[$i]=${ary[0]}
term2[$i]=${ary[2]}
i=$(expr $i + 1)
fi
done < $in2
i=0
while read BUF ; do
ary=(`echo $BUF`) # 配列に格納
if [ ${ary[0]} != "#" ]; then
datetime[$i]=${ary[0]}
term3[$i]=${ary[2]}
i=$(expr $i + 1)
fi
done < $in3
i=0
while read BUF ; do
ary=(`echo $BUF`) # 配列に格納
if [ ${ary[0]} != "#" ]; then
datetime[$i]=${ary[0]}
term4[$i]=${ary[2]}
i=$(expr $i + 1)
fi
done < $in4
n=$i
echo "# $(pwd)" > $out
echo "# $(hostname)" >> $out
echo "# $(whoami)" >> $out
echo "# $(date)" >> $out
echo "# $in1" >> $out
echo "# $in2" >> $out
echo "# datetime " >> $out
i=0
while [ $i -le $n ]; do
res[$i]=$(echo "scale=4; ${term4[$i]} - ${term1[$i]} + ${term2[$i]} - ${term3[$i]}" |bc)
echo "${datetime[$i]} ${term1[$i]} ${res[$i]}" >> $out
i=$(expr $i + 1)
done
echo
echo $in1
echo $in2
echo $in3
echo $in4
echo
echo $out
echo
exit 0
----------------------
End of calc_diff.sh
----------------------
作図用スクリプト (GMT使用）
------------------------------
List of the following files:
------------------------------
wv.budget.sh
gmtpar.sh
note.sh
------------------------------
Machine info
------------------------------
/work1/am/2016.PolarLow/ERA-I/Water.Vapor.Budget.ERA-I
======================
wv.budget.sh
======================
#!/bin/bash
# Description:
#
# Author: am
#
# Host: aofd165.bio.mie-u.ac.jp
# Directory: /work1/am/2016.PolarLow/ERA-I/Water.Vapor.Budget.ERA-I
#
# Revision history:
# This file is created by /usr/local/mybin/ngmt.sh at 14:59 on 04-19-2017.
. ./gmtpar.sh
echo "Bash script $0 starts."
gmtset PLOT_CLOCK_FORMAT hh
gmtset PLOT_DATE_FORMAT "o y"
gmtset TIME_FORMAT_PRIMARY abbreviated
gmtset ANOT_FONT_SIZE 12
gmtset INPUT_DATE_FORMAT yyyy-mm-dd
gmtset INPUT_CLOCK_FORMAT hh:mm:ss
xrange="2011-01-16T00:00:00/2011-01-26T00:00:00"
yrange1="-1/1.5" #Tnd
yrange2="0/25" #IWV
range1=$xrange/$yrange1
range2=$xrange/$yrange2
size=5/4
prefix=ERA-I
postfix=110101-31_Time.Series
var=IWV
var2=IWV
in0=$prefix"."$var"_"$prefix"_"${var2}_${postfix}.txt
var=WVF
var2=WVF
in1=$prefix"."$var"_"$prefix"_"${var2}_${postfix}.txt
var=Precip
var2=E-P
in2=$prefix"."$var"_"$prefix"_"${var2}_${postfix}.txt
var=Evapo
var2=E-P
in3=$prefix"."$var"_"$prefix"_"${var2}_${postfix}.txt
var=IWV_diff
var2=IWV_diff
in4=$prefix"."$var"_"$prefix"_"${var2}_${postfix}.txt
var=Res
var2=Res
in5=$prefix"."$var"_"$prefix"_"${var2}_${postfix}.txt
inlist="$in0 $in1 $in2 $in3 $in4 $in5"
for input in $inlist; do
if [ ! -f $input ]; then
echo Error in $0 : No such file, $input
exit 1
fi
done
#figdir="../fig/"
#if [ ! -d ${figdir} ];then
# mkdir -p $figdir
#fi
out=$(basename $0 .sh)_$prefix_${postfix}.ps
awk '{if ($1 != "missing" && $1 != "#") print $1, $3}' $in0|\
psxy -R$range2 -JX$size -W5/${black} \
-Bpa1df1d/a5f5:"[kg/m@+2@+]":E \
-X1.5 -Y5 -P -K > $out
awk '{if ($1 != "missing" && $1 != "#") print $1, $3}' $in1|\
psxy -R$range1 -JX$size -W3/${red} \
-O -K >> $out
awk '{if ($1 != "missing" && $1 != "#") print $1, -$3}' $in2|\
psxy -R$range1 -JX$size -W3/${blue} \
-O -K >> $out
awk '{if ($1 != "missing" && $1 != "#") print $1, $3}' $in3|\
psxy -R$range1 -JX$size -W3/${green} \
-Bpa1df1d/a0.5f0.5g100:"[10@+-4 @+kg/m@+2@+/s]":W \
-O -K >> $out
awk '{if ($1 != "missing" && $1 != "#") print $1, $3}' $in4|\
psxy -R$range1 -JX$size -W3/${purple} \
-O -K >> $out
awk '{if ($1 != "missing" && $1 != "#") print $1, $3}' $in5|\
psxy -R$range1 -JX$size -W3/${gray} \
-Bpa1df1d -Bsa1OSn \
-O -K >> $out
xoffset=
yoffset=4
comment=
. ./note.sh
in=${in0}
echo "Done $0"
echo
echo Input : $in0
echo Input : $in1
echo Input : $in2
echo Input : $in3
echo Input : $in4
echo Input : $in5
echo
echo Output : $out
echo
exit 0
----------------------
End of wv.budget.sh
----------------------
======================
gmtpar.sh
======================
#
# GMTのパラメータの設定 (version 4以上対象）
#
# 長さの単位はインチとする
gmtset MEASURE_UNIT INCH
gmtset LABEL_FONT_SIZE 18
gmtset HEADER_FONT_SIZE 20
gmtset ANOT_FONT_SIZE 18
gmtset TICK_PEN 4
# 地図の縦横軸に縞々を使わない
gmtset BASEMAP_TYPE PLAIN
# 紙のサイズはA4
gmtset PAPER_MEDIA A4
# 地図の°の記号
gmtset DEGREE_SYMBOL = degree # Available for ver. 4 or later
# 空白文字の設定
sp="\040" # White space
# =の記号
eq="\075" # equal
# 温度の°の記号
deg="\260" #deg="\312" # degree symbol
# 色のRGB値を設定 (線や点に色をつけるときRGB値を直接指定するより分かりやすい)
red="255/0/0"
orange="255/165/0"
pink="255/192/203"
green="0/255/0"
blue="0/0/255"
midnightblue="25/25/112"
yellow="255/255/0"
gold="255/215/0"
purple="160/32/240"
magenta="255/0/255"
white="255/255/255"
gray="180/180/180"
# 線種を指定
dash="t15_5:15"
dot="t5_5:5"
dotdash="t12_5_5_5:5"
#その他の線種の例
# -W5t20_5:5
# -W5t15_5:5
# -W5t20_5_5_5_5_5:5
# -W5t30_5_5_5:5
# -W5t20_5_10_5:5
# 数字の出力フォーマット(project, grd2xzy等で使う)
gmtset D_FORMAT %lg #デフォルト(規定値)
# 桁数を指定（例：123.45678)
#gmtset D_FORMAT %12.6f
# 状況に応じて
#gmtset D_FORMAT %12.5g
----------------------
End of gmtpar.sh
----------------------
======================
note.sh
======================
# Print time, current working directory and output filename
export LANG=C
currentdir=`pwd`
if [ ${#currentdir} -gt 90 ]; then
curdir1=${currentdir:1:90}
curdir2=${currentdir:91}
else
curdir1=$currentdir
curdir2="\ "
fi
now=`date`
host=`hostname`
#comment=" "
time1=$(ls -l $in | awk '{print $6, $7, $8}')
pstext -JX6/1.2 -R0/1/0/1.2 -N -X${xoffset:-0} -Y${yoffset:-0} << EOF -O >> $out
0 1.1 9 0 1 LM $0 $@
0 0.95 9 0 1 LM ${now}
0 0.80 9 0 1 LM ${host}
0 0.65 9 0 1 LM ${curdir1}
0 0.50 9 0 1 LM ${curdir2}
0 0.35 9 0 1 LM Input: ${in} (${time1})
0 0.1 9 0 1 LM ${comment:-""}
EOF
# Output: ${out}
----------------------
End of note.sh
----------------------
ダウンロードしたファイルに関する情報
------------------------------
List of the following files:
------------------------------
ERA-I_E-P_110101-31.txt
ERA-I_IWV_110101-31.txt
ERA-I_IWV_diff_110101-31.txt
ERA-I_WVF_110101-31.txt
------------------------------
Machine info
------------------------------
/work1/am/2016.PolarLow/ERA-I/Water.Vapor.Budget.ERA-I
======================
ERA-I_E-P_110101-31.txt
======================
netcdf
Final request
Stream:Atmospheric model
Area:90.0°N 180.0°W 30.0°N 180.0°E
Type:Forecast
Dataset:interim_daily
Step:12
Version:1
Type of level:Surface
Time:00:00:00, 12:00:00
Date:20110101 to 20110131
Grid:0.75° x 0.75°
Parameter:Evaporation, Total precipitation
Class:ERA Interim
The status of the request is: complete
Download (9.2MB)
[am@aofd165]
$ ncdump -h ERA-I_E-P_110101-31.nc
netcdf ERA-I_E-P_110101-31 {
dimensions:
longitude = 480 ;
latitude = 81 ;
time = UNLIMITED ; // (62 currently)
variables:
float longitude(longitude) ;
longitude:units = "degrees_east" ;
longitude:long_name = "longitude" ;
float latitude(latitude) ;
latitude:units = "degrees_north" ;
latitude:long_name = "latitude" ;
int time(time) ;
time:units = "hours since 1900-01-01 00:00:0.0" ;
time:long_name = "time" ;
time:calendar = "gregorian" ;
short e(time, latitude, longitude) ;
e:scale_factor = 2.42387284062225e-07 ;
e:add_offset = -0.00505705935012759 ;
e:_FillValue = -32767s ;
e:missing_value = -32767s ;
e:units = "m of water equivalent" ;
e:long_name = "Evaporation" ;
e:standard_name = "lwe_thickness_of_water_evaporation_amount" ;
short tp(time, latitude, longitude) ;
tp:scale_factor = 1.25801245107129e-06 ;
tp:add_offset = 0.0412200359718018 ;
tp:_FillValue = -32767s ;
tp:missing_value = -32767s ;
tp:units = "m" ;
tp:long_name = "Total precipitation" ;
// global attributes:
:Conventions = "CF-1.6" ;
:history = "2017-04-18 23:08:19 GMT by grib_to_netcdf-2.1.0: grib_to_netcdf /data/data04/scratch/_mars-atls19-a562cefde8a29a7288fa0b8b7f9413f7-4VWAct.grib -o /data/data04/scratch/_grib2netcdf-atls09-a562cefde8a29a7288fa0b8b7f9413f7-AFYibY.nc -utime" ;
----------------------
End of ERA-I_E-P_110101-31.txt
----------------------
======================
ERA-I_IWV_110101-31.txt
======================
netcdf
Final request
Stream:Atmospheric model
Area:90.0°N 180.0°W 30.0°N 180.0°E
Type:Analysis
Dataset:interim_daily
Step:0
Version:1
Type of level:Surface
Time:06:00:00, 18:00:00
Date:20110101 to 20110131
Grid:0.75° x 0.75°
Parameter:Total column water vapour, Vertical integral of water vapour
Class:ERA Interim
The status of the request is: complete
Download (9.2MB)
[~/2016.PolarLow/ERA-I/Water.Vapor.Budget.ERA-I]
[am@aofd165]
$ ncdump -h ERA-I_IWV_110101-31.nc netcdf ERA-I_IWV_110101-31 {
dimensions:
longitude = 480 ;
latitude = 81 ;
time = UNLIMITED ; // (62 currently)
variables:
float longitude(longitude) ;
longitude:units = "degrees_east" ;
longitude:long_name = "longitude" ;
float latitude(latitude) ;
latitude:units = "degrees_north" ;
latitude:long_name = "latitude" ;
int time(time) ;
time:units = "hours since 1900-01-01 00:00:0.0" ;
time:long_name = "time" ;
time:calendar = "gregorian" ;
short p55.162(time, latitude, longitude) ;
p55.162:scale_factor = 0.000815696301596405 ;
p55.162:add_offset = 26.817747768371 ;
p55.162:_FillValue = -32767s ;
p55.162:missing_value = -32767s ;
p55.162:units = "kg m**-2" ;
p55.162:long_name = "Vertical integral of water vapour" ;
short tcwv(time, latitude, longitude) ;
tcwv:scale_factor = 0.00081569645439851 ;
tcwv:add_offset = 26.8177309001802 ;
tcwv:_FillValue = -32767s ;
tcwv:missing_value = -32767s ;
tcwv:units = "kg m**-2" ;
tcwv:long_name = "Total column water vapour" ;
tcwv:standard_name = "lwe_thickness_of_atmosphere_mass_content_of_water_vapor" ;
// global attributes:
:Conventions = "CF-1.6" ;
:history = "2017-04-19 04:58:38 GMT by grib_to_netcdf-2.1.0: grib_to_netcdf /data/data04/scratch/_mars-atls19-a562cefde8a29a7288fa0b8b7f9413f7-RFvLIm.grib -o /data/data04/scratch/_grib2netcdf-atls02-a562cefde8a29a7288fa0b8b7f9413f7-EmTCCN.nc -utime" ;
}
----------------------
End of ERA-I_IWV_110101-31.txt
----------------------
======================
ERA-I_IWV_diff_110101-31.txt
======================
netcdf
Final request
Stream:Atmospheric model
Area:90.0°N 180.0°W 30.0°N 180.0°E
Type:Analysis
Dataset:interim_daily
Step:0
Version:1
Type of level:Surface
Time:00:00:00, 12:00:00
Date:20110101 to 20110201
Grid:0.75° x 0.75°
Parameter:Total column water vapour, Vertical integral of water vapour
Class:ERA Interim
The status of the request is: complete
Download (9.5MB)
Request output:
$ ncdump -h ERA-I_IWV_diff.nc netcdf ERA-I_IWV_diff {
dimensions:
longitude = 480 ;
latitude = 81 ;
time = UNLIMITED ; // (64 currently)
variables:
float longitude(longitude) ;
longitude:units = "degrees_east" ;
longitude:long_name = "longitude" ;
float latitude(latitude) ;
latitude:units = "degrees_north" ;
latitude:long_name = "latitude" ;
int time(time) ;
time:units = "hours since 1900-01-01 00:00:0.0" ;
time:long_name = "time" ;
time:calendar = "gregorian" ;
short p55.162(time, latitude, longitude) ;
p55.162:scale_factor = 0.000828587749480852 ;
p55.162:add_offset = 27.2483292700157 ;
p55.162:_FillValue = -32767s ;
p55.162:missing_value = -32767s ;
p55.162:units = "kg m**-2" ;
p55.162:long_name = "Vertical integral of water vapour" ;
short tcwv(time, latitude, longitude) ;
tcwv:scale_factor = 0.000828585734857865 ;
tcwv:add_offset = 27.2483453344748 ;
tcwv:_FillValue = -32767s ;
tcwv:missing_value = -32767s ;
tcwv:units = "kg m**-2" ;
tcwv:long_name = "Total column water vapour" ;
tcwv:standard_name = "lwe_thickness_of_atmosphere_mass_content_of_water_vapor" ;
// global attributes:
:Conventions = "CF-1.6" ;
:history = "2017-04-19 07:58:23 GMT by grib_to_netcdf-2.1.0: grib_to_netcdf /data/data04/scratch/_mars-atls17-a562cefde8a29a7288fa0b8b7f9413f7-34RMEK.grib -o /data/data04/scratch/_grib2netcdf-atls18-a562cefde8a29a7288fa0b8b7f9413f7-g__qJ1.nc -utime" ;
}
----------------------
End of ERA-I_IWV_diff_110101-31.txt
----------------------
======================
ERA-I_WVF_110101-31.txt
======================
netcdf
Final request
Stream:Atmospheric model
Area:90.0°N 180.0°W 30.0°N 180.0°E
Type:Analysis
Dataset:interim_daily
Step:0
Version:1
Type of level:Surface
Time:06:00:00, 18:00:00
Date:20110101 to 20110131
Grid:0.75° x 0.75°
Parameter:Vertical integral of divergence of moisture flux, Vertical integral of eastward water vapour flux, Vertical integral of northward water vapour flux
Class:ERA Interim
The status of the request is: complete
Download (13.8MB)
[~/2016.PolarLow/ERA-I/Water.Vapor.Budget.ERA-I]
[am@aofd165]
$ ncdump -h ERA-I_WVF_110101-31.nc
netcdf ERA-I_WVF_110101-31 {
dimensions:
longitude = 480 ;
latitude = 81 ;
time = UNLIMITED ; // (62 currently)
variables:
float longitude(longitude) ;
longitude:units = "degrees_east" ;
longitude:long_name = "longitude" ;
float latitude(latitude) ;
latitude:units = "degrees_north" ;
latitude:long_name = "latitude" ;
int time(time) ;
time:units = "hours since 1900-01-01 00:00:0.0" ;
time:long_name = "time" ;
time:calendar = "gregorian" ;
short p71.162(time, latitude, longitude) ;
p71.162:scale_factor = 0.039770833045231 ;
p71.162:add_offset = 634.89771712254 ;
p71.162:_FillValue = -32767s ;
p71.162:missing_value = -32767s ;
p71.162:units = "kg m**-1 s**-1" ;
p71.162:long_name = "Vertical integral of eastward water vapour flux" ;
short p72.162(time, latitude, longitude) ;
p72.162:scale_factor = 0.0310094879738548 ;
p72.162:add_offset = 556.210935685701 ;
p72.162:_FillValue = -32767s ;
p72.162:missing_value = -32767s ;
p72.162:units = "kg m**-1 s**-1" ;
p72.162:long_name = "Vertical integral of northward water vapour flux" ;
short p84.162(time, latitude, longitude) ;
p84.162:scale_factor = 1.08452245927098e-07 ;
p84.162:add_offset = 0.000282594957522881 ;
p84.162:_FillValue = -32767s ;
p84.162:missing_value = -32767s ;
p84.162:units = "kg m**-2 s**-1" ;
p84.162:long_name = "Vertical integral of divergence of moisture flux" ;
// global attributes:
:Conventions = "CF-1.6" ;
:history = "2017-04-19 04:06:16 GMT by grib_to_netcdf-2.1.0: grib_to_netcdf /data/data04/scratch/_mars-atls17-a562cefde8a29a7288fa0b8b7f9413f7-_iqnAb.grib -o /data/data04/scratch/_grib2netcdf-atls20-a562cefde8a29a7288fa0b8b7f9413f7-kVW9gd.nc -utime" ;
}
----------------------
End of ERA-I_WVF_110101-31.txt
----------------------