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



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
#----------------------
```