LoginSignup
1
1

More than 5 years have passed since last update.

WRFの出力から定点における時間—高度断面を描く

Last updated at Posted at 2019-04-05

2019-04-06_07-36
/work05/manda/WRF.POST/NCL/H30.7.GOUU/SOUNDING
$ WRF_TIME-HEIGHT.sh
script name =WRF_TIME-HEIGHT.ncl
namelist =WRF_TIME-HEIGHT.nml


List of the following files:

WRF_TIME-HEIGHT.sh
WRF_TIME-HEIGHT.ncl
WRF_PL_TIME-HEIGHT.sh
WRF_OBS_PL_TIME-HEIGHT.sh

WRF_TIME-HEIGHT.sh
#!/bin/bash

RUNNAME=H30.07.ERA.MYNN.NoSF01
INDIR_ROOT=/work05/manda/WRF.Result
INDIR=${INDIR_ROOT}/${RUNNAME}
WRFOUT_TIME=2018-07-04_00:00:00
DOMAIN=d03
WRFOUT=wrfout_${DOMAIN}_${WRFOUT_TIME}

STNNM=MATSUE
YYYY=2018
PLAT=35.45
PLON=133.07

OFLE=$(basename $0 .sh)_${STNNM}_${RUNNAME}.${DOMAIN}.txt

if [ ! -f ${INDIR}/${WRFOUT} ]; then
echo ERROR in $0 : NO SUCH FILE, ${INDIR}/$WRFOUT
fi

NCL=$(basename $0 .sh).ncl
NML=$(basename $0 .sh).nml
RUN=runncl.sh

cat <<EOF>$NML
STNNM   =${STNNM}
YYYY    =${YYYY}
PLAT    =${PLAT}
PLON    =${PLON}
RUNNAME =${RUNNAME}
INDIR   =${INDIR}
WRFOUT  =${WRFOUT}
OFLE    =${OFLE}
EOF

$RUN $NCL "${NML}"

exit 0

#----------------------
# End of WRF_TIME-HEIGHT.sh
#----------------------
WRF_TIME-HEIGHT.ncl
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin

MO=7
DS=5
DE=8
HS=0
DH=3
INIT_ETIME=24
interp_levels    = (/300,400,500,600,700,800,850,900,925,950,975,1000/)
INIT_OUT_TIME="2018-07-05_00:00:00"
LAST_OUT_TIME="2018-07-08_00:00:00"


scriptname_in = getenv("NCL_ARG_1")
scriptname=systemfunc("basename " + scriptname_in + " .ncl")
nml    = getenv("NCL_ARG_2")

print("script name ="+scriptname_in)
print("namelist    ="+nml)
print("")



print("READ NAMELIST FILE")

STNNM=systemfunc("grep STNNM "+nml+ "|cut -f2 -d'='|cut -f1 -d','")
YYYY=toint(systemfunc("grep YYYY "+nml+ "|cut -f2 -d'='|cut -f1 -d','"))
PLAT=tofloat(systemfunc("grep PLAT "+nml+ "|cut -f2 -d'='|cut -f1 -d','"))
PLON=tofloat(systemfunc("grep PLON "+nml+ "|cut -f2 -d'='|cut -f1 -d','"))
RUNNAME=systemfunc("grep RUNNAME "+nml+ "|cut -f2 -d'='|cut -f1 -d','")
INDIR=systemfunc("grep INDIR "+nml+ "|cut -f2 -d'='|cut -f1 -d','")
WRFOUT=systemfunc("grep WRFOUT "+nml+ "|cut -f2 -d'='|cut -f1 -d','")
OFLE=systemfunc("grep OFLE "+nml+ "|cut -f2 -d'='|cut -f1 -d','")
print("DONE READ NAMELIST FILE.")
print("")

print("CHECK NAMELIST FILE")

print("STNNM:  "+STNNM)
print("YYYY:   "+YYYY)
print("PLAT:   "+PLAT)
print("PLON:   "+PLON)
print("RUNNAME:"+RUNNAME)
print("INDIR:  "+INDIR)
print("WRFOUT: "+WRFOUT)
print("OFLE  : "+OFLE)

print("DONE CHECK NAMELIST FILE.")
print("")



a = addfile(INDIR+"/"+WRFOUT,"r")
;print(a)

times=wrf_user_getvar(a,"times",-1)
;printVarSummary(times)
dim=dimsizes(times)
nt=dim(0)
do n=0,nt-1
if(times(n) .eq. INIT_OUT_TIME)then
nout_s=n
end if
if(times(n) .eq. LAST_OUT_TIME)then
nout_e=n
end if
end do ;n

delete(dim)

print("OUTPUT START TIME="+times(nout_s))
print("OUTPUT END TIME  ="+times(nout_e))


print("OUTPUT: "+OFLE)

PWD=systemfunc("pwd")
NOW=systemfunc("date -R")

header_ofle = (/\
"# PWD: "+ PWD, \
"# NOW: "+ NOW, \
"# SCRIPT: "+scriptname_in, \
"# STNNM: "+ STNNM, \
"# PLAT  : "+ PLAT, \
"# PLON  : "+ PLON, \
"# INDIR: "+ INDIR, \
"#    hr             datetime     P      U        V     EPT", \
"#                              hPa     m/s     m/s       K"/)
hlist=[/header_ofle/]
write_table(OFLE, "w", hlist, "%s")

ETIME=tofloat(INIT_ETIME)

do n=nout_s,nout_e; TIME LOOP

print("ETIME="+ETIME+" TIME="+times(n))

lat2d = wrf_user_getvar(a,"XLAT",n)    ; size = (nlat,nlon)
lon2d = wrf_user_getvar(a,"XLONG",n)    ; size = (nlat,nlon)
;printVarSummary(lat2d)
;printVarSummary(lon2d)

p3  = wrf_user_getvar(a,"pressure",n)
pi = rcm2points (lat2d, lon2d, p3, PLAT, PLON, 2)
dim=dimsizes(pi)
km=dim(0)
p1=new((/km/),typeof(pi),pi@_FillValue)

; EQUIVALENT POTENTIAL TEMPERATURE (K)
pfull=wrf_user_getvar(a,"pres",n)
qv   =wrf_user_getvar(a,"QVAPOR",n)
tk   =wrf_user_getvar(a,"tk",n)

pfi = rcm2points (lat2d, lon2d, pfull, PLAT, PLON, 2)
qvi = rcm2points (lat2d, lon2d, qv   , PLAT, PLON, 2)
tki = rcm2points (lat2d, lon2d, tk   , PLAT, PLON, 2)

pf3=new((/km,1,1/),typeof(pfi),pfi@_FillValue)
qv3=new((/km,1,1/),typeof(qvi),qvi@_FillValue)
tk3=new((/km,1,1/),typeof(tki),tki@_FillValue)

pf3(:,0,0)=pfi(:,0)
qv3(:,0,0)=qvi(:,0)
tk3(:,0,0)=tki(:,0)

ept=wrf_eth(qv3,tk3,pf3)
;printVarSummary(ept)
delete([/qv3,tk3,pf3/])


dim=dimsizes(interp_levels)
NLEV=dim(0)
po=new((/NLEV/),float,default_fillvalue("float"))

uvm = wrf_user_getvar(a,"uvmet",n)    ; size = (nlat,nlon)
um=uvm(0,:,:,:)
vm=uvm(1,:,:,:)
ui = rcm2points (lat2d, lon2d, um, PLAT, PLON, 2)
vi = rcm2points (lat2d, lon2d, vm, PLAT, PLON, 2)
dim=dimsizes(ui)
km=dim(0)
u1  =new((/km/),typeof(ui), ui@_FillValue)
v1  =new((/km/),typeof(vi), vi@_FillValue)
ept1=new((/km/),typeof(ept),default_fillvalue(typeof(ept)))

;Array Reversal https://www.ncl.ucar.edu/Document/Language/array_reverse.shtml
p1(:)  = pi(::-1,0)
u1(:)  = ui(::-1,0)
v1(:)  = vi(::-1,0)
ept1(:)= ept(::-1,0,0)

;do k=0,km-1
;print(p1(k)+" "+u1(k)+" "+v1(k))
;end do ;k
;print("")

po=tofloat(interp_levels)
uo  = linint1 (p1,u1, False, po, 0)
vo  = linint1 (p1,v1, False, po, 0)
epto= linint1 (p1,ept1, False, po, 0)
delete([/p1,u1,v1,ept1/])

;SURFACE DATA
psfc   =wrf_user_getvar(a,"PSFC",n)
t2     =wrf_user_getvar(a,"T2",n)
q2     =wrf_user_getvar(a,"Q2",n)
uvm10=wrf_user_getvar(a,"uvmet10",n)
um10=uvm10(0,:,:)
vm10=uvm10(1,:,:)

u10i = rcm2points (lat2d, lon2d, um10, PLAT, PLON, 2)
v10i = rcm2points (lat2d, lon2d, vm10, PLAT, PLON, 2)
psi = rcm2points (lat2d, lon2d, psfc, PLAT, PLON, 2)
t2i = rcm2points (lat2d, lon2d, t2  , PLAT, PLON, 2)
q2i = rcm2points (lat2d, lon2d, q2  , PLAT, PLON, 2)
;printVarSummary(psi)

pf3=new((/1,1,1/),typeof(psi),pfi@_FillValue)
qv3=new((/1,1,1/),typeof(q2i),qvi@_FillValue)
tk3=new((/1,1,1/),typeof(t2i),tki@_FillValue)
pf3(:,0,0)=psi(:)
qv3(:,0,0)=q2i(:)
tk3(:,0,0)=t2i(:)

epts=wrf_eth(qv3,tk3,pf3)
;printVarSummary(epts)
delete([/qv3,tk3,pf3/])

do k=0,NLEV-1
if(ismissing(uo(k)) .eq. False)then

alist= [/ETIME, times(n), po(k), uo(k), vo(k), epto(k) /]
write_table(OFLE, "a", alist, "%7.1f %s %6.0f %7.3f %7.3f %7.3f %5.1f %8.5f %7.3f")

;print(po(k)+" "+uo(k)+" "+vo(k)+" "+epto(k))
end if
end do ;k

psihpa=psi(0)/100.
alist= [/ETIME, times(n), psihpa, u10i(0), v10i(0), epts(0,0,0) /]
write_table(OFLE, "a", alist, "%7.1f %s %6.0f %7.3f %7.3f %7.3f %5.1f %8.5f %7.3f")

;print(psihpa+" "+u10i(0)+" "+v10i(0)+" "+epts(0,0,0))



ETIME=ETIME+toint(DH)
end do ;n TIME-LOOP

end

#----------------------
# End of WRF_TIME-HEIGHT.ncl
#----------------------
WRF_PL_TIME-HEIGHT.sh
#!/bin/bash
# Description:
#
# Author: manda
#
# Host: calypso.bosai.go.jp
# Directory: /work05/manda/18.H30.GOUU/SOUNDING
#
# Revision history:
#  This file is created by /home/manda/mybin/ngmt.sh at 17:44 on 04-04-2019.

. ./gmtpar.sh
echo "Bash script $0 starts."
gmtset HEADER_OFFSET -0.5c
gmtset HEADER_FONT_SIZE 16

RUNNAME=H30.07.ERA.MYNN.NoSF01
RNSHORT=CNTL
DOMAIN=d03
DATASET=$RNSHORT

#STNNO=47600
#STNNM=WAJIMA
STNNO=47741
STNNM=MATSUE

INIT_TIME=00Z04JUL2018

VSCL=0.02 #1m/s=>0.02 inch ~ 0.025*2 cm ~ 0.25*2 mm

VEC_LEG_MAG=10
VSCL_MAG=$(echo "scale=4; ${VEC_LEG_MAG}*${VSCL}"|bc)
PI=3.14159265359

range=18/102/760/1020
size=X5/-5
xanot=a24f12
yanot=a50f50

xlab="Elapsed${sp}time${sp}from${sp}${INIT_TIME}${sp}[hr]"
ylab="hPa"
title="${DATASET}${sp}${STNNM}"
anot=${xanot}:${xlab}:/${yanot}:${ylab}::.${title}:WSne

#indir=${STNNO}${STNNM}
in=WRF_TIME-HEIGHT_${STNNM}_${RUNNAME}.${DOMAIN}.txt
if [ ! -f $in ]; then
  echo Error in $0 : No such file, $in
  exit 1
fi

out=$(basename $0 .sh)_${STNNM}_${RUNNAME}.${DOMAIN}.ps

# EPT

GRD=$(basename $0 .sh).EPT.grd
XYZ=$(basename $0 .sh).EPT.xyz

awk '{if($1!="#" && $3!=-999.9 && $6!=-999.9) print $1,$3,$6}' $in |\
blockmean -R$range -I12/20 >$XYZ

surface $XYZ -R$range -T1 -I6/10 -G$GRD


grdcontour $GRD -J$size -R$range -A4f12 -C2 -X1.5 -Y2 -K -P >$out

#awk '{ if($1!="#" && $3!=-999.9) print $1,$3}' $in |\
#psxy -R -J$size -Sc0.1 -G0 -O -K >>$out

# WIND
awk -v s=$VSCL -v pi=${PI} \
'{if($1!="#"&&$9!=-999.9&&$10!=-999.9) \
printf "%12.6f %12.6f %12.6f %12.6f\n", \
$1,$3,(180.0/pi)*atan2(-$5,$4),sqrt($4*$4+$5*$5)*s}' $in| \
psxy -R -J$size -Sv0.01/0.1/0.05n0.2 -G0 -O -K >>$out
#===========================================================
# NOTE: NEGATIVE OF $5 IS USED IN THE ABOVE AWK SCRIPT, 
# SINCE Y-AXIS IS REVERSED.
#===========================================================

psbasemap -J$size -R$range -B$anot -O -K >>$out


yoffset=5

echo "0.85 0.2 270 ${VEC_LEG_MAG}" |\
awk -v s=$VSCL '{print $1,$2, 270-$3,$4*s}' |\
psxy -R0/1/0/1 -JX5/0.5 -Sv0.01/0.1/0.05n0.2 -G0 -O -K -Y${yoffset} >>$out

echo "0.9 0.2 10 0 1 CM E" |\
pstext -R -JX -O -K >>$out

echo "0.85 0.2 180 ${VEC_LEG_MAG}" |\
awk -v s=$VSCL '{print $1,$2, 270-$3,$4*s}' |\
psxy -R0/1/0/1 -JX5/0.5 -Sv0.01/0.1/0.05n0.2 -G0 -O -K >>$out

echo "0.85 0.7 10 0 1 CM N" |\
pstext -R -JX -O -K >>$out

echo "0.77 0.3 14 0 1 CM ${VEC_LEG_MAG}m/s" |\
pstext -R -JX -O -K >>$out


xoffset=
yoffset=0

export LANG=C

curdir1=$(pwd)
now=$(date)
host=$(hostname)

time=$(ls -l ${in} | awk '{print $6, $7, $8}')
time1=$(ls -l ${out} | awk '{print $6, $7, $8}')

pstext <<EOF -JX6/1.5 -R0/1/0/1.5 -N -X${xoffset:-0} -Y${yoffset:-0} -O >> $out
0 1.50  9 0 1 LM $0 $@
0 1.35  9 0 1 LM ${now}
0 1.20  9 0 1 LM ${host}
0 1.05  9 0 1 LM ${curdir1}
0 0.90  9 0 1 LM Input:  ${in}  (${time})
0 0.75  9 0 1 LM Output: ${out} (${time1})
EOF
#0 0.60  9 0 1 LM Input:
#0 0.45  9 0 1 LM Input:

rm -vf $GRD $XYZ

echo OUTPUT:
ls -lh --time-style=long-iso $out

exit

#----------------------
# End of WRF_PL_TIME-HEIGHT.sh
#----------------------
WRF_OBS_PL_TIME-HEIGHT.sh
#!/bin/bash
# Description:
#
# Author: manda
#
# Host: calypso.bosai.go.jp
# Directory: /work05/manda/18.H30.GOUU/SOUNDING
#
# Revision history:
#  This file is created by /home/manda/mybin/ngmt.sh at 17:44 on 04-04-2019.

. ./gmtpar.sh
echo "Bash script $0 starts."
gmtset HEADER_OFFSET -0.5c
gmtset HEADER_FONT_SIZE 16

RUNNAME=H30.07.ERA.MYNN.NoSF01
RNSHORT=CNTL
DOMAIN=d03
DATASET=$RNSHORT

#STNNO=47600
#STNNM=WAJIMA
STNNO=47741
STNNM=MATSUE

INIT_TIME=00Z04JUL2018

VSCL=0.02 #1m/s=>0.02 inch ~ 0.025*2 cm ~ 0.25*2 mm

VEC_LEG_MAG=10
VSCL_MAG=$(echo "scale=4; ${VEC_LEG_MAG}*${VSCL}"|bc)
PI=3.14159265359

range=18/102/760/1020
size=X5/-5
xanot=a24f12
yanot=a50f50

xlab="Elapsed${sp}time${sp}from${sp}${INIT_TIME}${sp}[hr]"
ylab="hPa"
anot=${xanot}:${xlab}:/${yanot}:${ylab}:WSne

#indir=${STNNO}${STNNM}
in=WRF_TIME-HEIGHT_${STNNM}_${RUNNAME}.${DOMAIN}.txt
if [ ! -f $in ]; then
  echo Error in $0 : No such file, $in
  exit 1
fi

out=$(basename $0 .sh)_${STNNM}_${RUNNAME}.${DOMAIN}.ps

# EPT

GRD=$(basename $0 .sh).EPT.grd
XYZ=$(basename $0 .sh).EPT.xyz

awk '{if($1!="#" && $3!=-999.9 && $6!=-999.9) print $1,$3,$6}' $in |\
blockmean -R$range -I12/20 >$XYZ

surface $XYZ -R$range -T1 -I6/10 -G$GRD


grdcontour $GRD -J$size -R$range -A4f12 -C2 -X1.5 -Y2 -K -P >$out

#awk '{ if($1!="#" && $3!=-999.9) print $1,$3}' $in |\
#psxy -R -J$size -Sc0.1 -G0 -O -K >>$out

# WIND
awk -v s=$VSCL -v pi=${PI} \
'{if($1!="#"&&$9!=-999.9&&$10!=-999.9) \
printf "%12.6f %12.6f %12.6f %12.6f\n", \
$1,$3,(180.0/pi)*atan2(-$5,$4),sqrt($4*$4+$5*$5)*s}' $in| \
psxy -R -J$size -Sv0.01/0.1/0.05n0.2 -G0 -O -K >>$out
#===========================================================
# NOTE: NEGATIVE OF $5 IS USED IN THE ABOVE AWK SCRIPT, 
# SINCE Y-AXIS IS REVERSED.
#===========================================================


#==============================================================
# OBSERVATION
#
indir=/work05/manda/18.H30.GOUU/SOUNDING/${STNNO}${STNNM}
in1=${indir}/${STNNO}${STNNM}_OUT.txt
if [ ! -f $in1 ]; then
  echo Error in $0 : No such file, $in1
  exit 1
fi


# EPT

GRD=$(basename $0 .sh).EPT.grd
XYZ=$(basename $0 .sh).EPT.xyz

awk '{if($1!="#" && $3!=-999.9 && $14!=-999.9) print $1,$3,$14}' $in1 |\
blockmean -R$range -I12/20 >$XYZ

surface $XYZ -R$range -T1 -I6/10 -G$GRD


grdcontour $GRD -J$size -R$range -W3/255/0/0 -A4f12+kred -C2 -O -K>>$out

#awk '{ if($1!="#" && $3!=-999.9) print $1,$3}' $in |\
#psxy -R -J$size -Sc0.1 -G0 -O -K >>$out

# WIND
awk -v s=$VSCL -v pi=${PI} \
'{if($1!="#"&&$9!=-999.9&&$10!=-999.9) \
printf "%12.6f %12.6f %12.6f %12.6f\n", \
$1,$3,(180.0/pi)*atan2(-$10,$9),sqrt($9*$9+$10*$10)*s}' $in1| \
psxy -R -J$size -Sv0.01/0.1/0.05n0.2 -Gred -O -K >>$out
#===========================================================
# NOTE: NEGATIVE OF $10 IS USED IN THE ABOVE AWK SCRIPT, 
# SINCE Y-AXIS IS REVERSED.
#===========================================================


#
# END OF OBSERVATION
#==============================================================


psbasemap -J$size -R$range -B$anot -O -K >>$out


yoffset=5

echo "0.85 0.2 270 ${VEC_LEG_MAG}" |\
awk -v s=$VSCL '{print $1,$2, 270-$3,$4*s}' |\
psxy -R0/1/0/1 -JX5/0.5 -Sv0.01/0.1/0.05n0.2 -G0 -O -K -Y${yoffset} >>$out

echo "0.9 0.2 10 0 1 CM E" |\
pstext -R -JX -O -K >>$out

echo "0.85 0.2 180 ${VEC_LEG_MAG}" |\
awk -v s=$VSCL '{print $1,$2, 270-$3,$4*s}' |\
psxy -R0/1/0/1 -JX5/0.5 -Sv0.01/0.1/0.05n0.2 -G0 -O -K >>$out

echo "0.85 0.7 10 0 1 CM N" |\
pstext -R -JX -O -K >>$out

echo "0.77 0.35 14 0 1 CM ${VEC_LEG_MAG}m/s" |\
pstext -R -JX -O -K >>$out

echo "0 0.35 14 0 1 LM ${STNNM}" |\
pstext -R -JX -O -K >>$out

echo "0.25 0.35 14 0 1 LM SONDE" |\
pstext -R -JX -Gred -O -K >>$out

echo "0.4 0.35 14 0 1 LM ${RNSHORT}" |\
pstext -R -JX -O -K >>$out


xoffset=
yoffset=0

export LANG=C

curdir1=$(pwd)
now=$(date)
host=$(hostname)

time=$(ls -l ${in} | awk '{print $6, $7, $8}')
time1=$(ls -l ${in1} | awk '{print $6, $7, $8}')
time2=$(ls -l ${out} | awk '{print $6, $7, $8}')

pstext <<EOF -JX6/1.5 -R0/1/0/1.5 -N -X${xoffset:-0} -Y${yoffset:-0} -O >> $out
0 1.50  9 0 1 LM $0 $@
0 1.35  9 0 1 LM ${now}
0 1.20  9 0 1 LM ${host}
0 1.05  9 0 1 LM ${curdir1}
0 0.90  9 0 1 LM Input:  ${in}  (${time})
0 0.75  9 0 1 LM Input : ${in1} (${time1})
0 0.60  9 0 1 LM Output: ${out} (${time2})
EOF
#0 0.60  9 0 1 LM Input:
#0 0.45  9 0 1 LM Input:

rm -vf $GRD $XYZ

echo OUTPUT:
ls -lh --time-style=long-iso $out

exit

#----------------------
# End of WRF_OBS_PL_TIME-HEIGHT.sh
#----------------------
1
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
1
1