3時間雨量の計算

[2017年 2月 22日 水曜日 15:54:18 JST]

[~/2017.Heavy.Rain/Radar.Precip.3hr]

[am@aofd165]

$ radar.precip.3hr.run.sh 20160830 09

............

./radar.precip.3hr.run.sh

Input:

intm.data.1hr/R1h.20160830_07.bin

intm.data.1hr/R1h.20160830_08.bin

intm.data.1hr/R1h.20160830_09.bin

Output: output/R3h.20160830_09.bin

$ grads -bcp "r3h.gs"

.............

Output: R3h.20160830_09.eps

[2017年 2月 22日 水曜日 15:56:51 JST]

[~/2017.Heavy.Rain/Radar.Precip.3hr]

[am@aofd165]

$ srcdump.sh radar.precip.3hr.run.sh 1hr.to.3hr/*.f90 radar.precip.1hr.run.sh 10min.to.1hr/*.f90 read_jmaradar/*.c read_jmaradar/Makefile r3h.gs>tmp.txt

------------------------------

List of the following files:

------------------------------

radar.precip.3hr.run.sh

1hr.to.3hr/1hr_to_3hr.f90

radar.precip.1hr.run.sh

10min.to.1hr/10min_to_1hour.f90

r3h.gs

------------------------------

Machine info

------------------------------

aofd165.bio.mie-u.ac.jp

/work1/am/2017.Heavy.Rain/Radar.Precip.3hr

Wed Feb 22 15:54:18 JST 2017

======================

radar.precip.3hr.run.sh

======================

#!/bin/bash

rdir=/mnt/raid112MANDA/KakenA.H28/JMA-Radar

nx=2560

ny=3360

x0=118.006250

y0=20.004167

dx=0.012500

dy=0.008333

itmdir1hr="intm.data.1hr"

outdir="output"

mkdir -vp $outdir

exe1=radar.precip.1hr.run.sh

if [ ! -f $exe1 ]; then

echo

echo Error in $0 : No such file, $exe1

echo

exit 1

fi

if [ $# -ne 2 ]; then

echo ERROR in $0: Wrong argument.

echo

echo Usage: $0 yyyymmdd hh

echo

exit 1

fi

yyyymmdd_i=$1

hh_i=$2

opt1="-r $rdir -m $nx -n $ny -x $x0 -y $y0 -d $dx -e $dy -i $itmdir1hr"

if [ $hh_i = "00" ]; then

yyyymmdda[0]=$(date -d"${yyyymmdd_i}"-1days '+%Y%m%d')

hha[0]="22"

yyyymmdda[1]=$(date -d"${yyyymmdd_i}"-1days '+%Y%m%d')

hha[1]="23"

yyyymmdda[2]=${yyyymmdd_i}

hha[2]="00"

elif [ $hh_i = "01" ]; then

yyyymmdda[0]=$(date -d"${yyyymmdd_i}"-1days '+%Y%m%d')

hha[0]="23"

yyyymmdda[2]=${yyyymmdd_i}

hha[2]="01"

yyyymmdda[1]=${yyyymmdd_i}

hha[1]="00"

else

yyyymmdda[0]=${yyyymmdd_i}

hha[0]=$(printf %02d $(expr ${hh_i} - 2))

yyyymmdda[1]=${yyyymmdd_i}

hha[1]=$(printf %02d $(expr ${hh_i} - 1))

yyyymmdda[2]=${yyyymmdd_i}

hha[2]=${hh_i}

fi

i=0

while [ $i -le 2 ]; do

yyyymmdd=${yyyymmdda[$i]}

hh=${hha[$i]}

echo $yyyymmdd $hh

itm1hr[$i]=R1h.${yyyymmdd}_${hh}.bin

$exe1 $opt1 $yyyymmdd $hh ${itm1hr[$i]}

if [ $? -ne 0 ]; then

echo

echo "Error in $0 while running $exe1"

echo

echo $yyyymmdd $opt1 $hh ${itm1hr[$i]}

echo

exit 1

fi

i=$(expr $i + 1)

done

echo

echo Hourly averaged data:

i=0

while [ $i -le 2 ]; do

echo "${itmdir1hr}/${itm1hr[$i]}"

i=$(expr $i + 1)

done

echo

exe=1hr_to_3hr

if [ ! -f $exe ]; then

echo

echo ERROR in $0: No such file, $exe

echo

exit 1

fi

output=${outdir}/R3h.${yyyymmdd_i}_${hh_i}.bin

namelist="namelist.${exe}.txt"

cat <<EOF >$namelist

&para

im=${nx},

jm=${ny},

input0="${itmdir1hr}/${itm1hr[0]}",

input1="${itmdir1hr}/${itm1hr[1]}",

input2="${itmdir1hr}/${itm1hr[2]}",

output="${output}",

&end

EOF

command="$exe < $namelist"

echo

echo $command

echo

$exe < $namelist

if [ $? -ne 0 ]; then

echo

echo "Error in $0 while running the following:"

echo " $command"

echo

exit 1

fi

echo

echo $0

echo

echo Input:

i=0

while [ $i -le 2 ]; do

echo "${itmdir1hr}/${itm1hr[$i]}"

i=$(expr $i + 1)

done

echo

echo "Output: $output"

echo

yyyy=${yyyymmdd_i:0:4}

mm=${yyyymmdd_i:4:2}

dd=${yyyymmdd_i:6:2}

hh=${hh_i}

if [ $mm = "01" ]; then mmm="Jan"; fi

if [ $mm = "02" ]; then mmm="Feb"; fi

if [ $mm = "03" ]; then mmm="Mar"; fi

if [ $mm = "04" ]; then mmm="Apr"; fi

if [ $mm = "05" ]; then mmm="May"; fi

if [ $mm = "06" ]; then mmm="Jun"; fi

if [ $mm = "07" ]; then mmm="Jul"; fi

if [ $mm = "08" ]; then mmm="Aug"; fi

if [ $mm = "09" ]; then mmm="Sep"; fi

if [ $mm = "10" ]; then mmm="Oct"; fi

if [ $mm = "11" ]; then mmm="Nov"; fi

if [ $mm = "12" ]; then mmm="Dec"; fi

ctl=r3h.ctl

cat <<EOF>$ctl

dset ${output}

undef -999.9

options template

xdef ${nx} LINEAR ${x0} ${dx}

ydef ${ny} LINEAR ${y0} ${dy}

zdef 1 LEVELS 1000

tdef 1 LINEAR ${hh}:00Z${dd}${mmm}${yyyy} 1hr

vars 1

rr 0 0 rainfall

endvars

EOF

echo

echo GrADS control file, $ctl

echo

cat $ctl

echo

exit 0

----------------------

End of radar.precip.3hr.run.sh

----------------------

======================

1hr.to.3hr/1hr_to_3hr.f90

======================

program one_to_three

! Description:

!

! Host: aofd165.fish.nagasaki-u.ac.jp

! Directory: /work2/kunoki/to_manda_sensei/Tools/JMA_radar/10min.to.1hr

!

! ifortの場合以下のようにオプション-assume bytereclを付けてコンパイルすること。

! $ ifort -assume byterecl hoge.f90

! こうしないと、"recl"がbyte単位ではなく、word単位になってしまう

!

! use

! implicit none

real,allocatable:: rr(:,:),rr1(:,:)

integer im,jm

character(len=500)::input0, input1, input2

character(len=500)::output

namelist /para/im,jm,input0, input1, input2, output

read(*,nml=para)

print '(A,A)','Input: ',trim(input0)

print '(A,A)','Input: ',trim(input1)

print '(A,A)','Input: ',trim(input2)

print '(A,A)','Output: ',trim(output)

print *

allocate(rr(im,jm),rr1(im,jm))

rr(:,:)=0.0

open(10,file=input0,form="unformatted",access="direct",recl=4*im*jm,action="read")

read(10,rec=1)rr1

rr(:,:)=rr(:,:)+rr1(:,:)

close(10)

open(10,file=input1,form="unformatted",access="direct",recl=4*im*jm,action="read")

read(10,rec=1)rr1

rr(:,:)=rr(:,:)+rr1(:,:)

close(10)

open(10,file=input2,form="unformatted",access="direct",recl=4*im*jm,action="read")

read(10,rec=1)rr1

rr(:,:)=rr(:,:)+rr1(:,:)

close(10)

!rr(:,:)=rr(:,:)/6.0

! We do NOT need averaging since we want 3-hr precipitation (mm/3hr).

!

! Unit of rr1 is mm/hr.

!

open(10,file=output,form="unformatted",access="direct",recl=4*im*jm)

write(10,rec=1)rr

close(10)

end program one_to_three

----------------------

End of 1hr.to.3hr/1hr_to_3hr.f90

----------------------

======================

radar.precip.1hr.run.sh

======================

#!/bin/sh

# Default values

rdir=/mnt/raid112MANDA/KakenA.H28/JMA-Radar

nx=2560

ny=3360

x0=118.006250

y0=20.004167

dx=0.012500

dy=0.008333

itmdir1h="intm.data.1hr"

while getopts r:m:n:x:y:d:e:i: OPT; do

case $OPT in

"r" ) flagr="true" ; rdir="$OPTARG" ;;

"m" ) flagm="true" ; nx="$OPTARG" ;;

"n" ) flagn="true" ; ny="$OPTARG" ;;

"x" ) flagx="true" ; x0="$OPTARG" ;;

"y" ) flagy="true" ; y0="$OPTARG" ;;

"d" ) flagd="true" ; dx="$OPTARG" ;;

"e" ) flage="true" ; dy="$OPTARG" ;;

"i" ) flago="true" ; itmdir1hr="$OPTARG" ;;

* ) echo "Usage: $0 [-r rdir] [-m nx] [-n ny] [-x x0] [-y y0] [-d dx] [-e dy] [-i itmdir1h] yyyymmdd hh output" 1>&2

esac

done

shift $(expr $OPTIND - 1) # 不要となったオプション部分を shift コマンドで切り捨てる

#echo

#echo "rdir: $rdir"

#echo "nx = $nx"

#echo "ny = $ny"

#echo "x0 = $x0"

#echo "y0 = $y0"

#echo "dx = $dx"

#echo "dy = $dy"

#echo "itmdir1h: $itmdir1h"

#echo

if [ $# -ne 3 ]; then

echo ERROR in $0: Wrong argument.

echo

echo Usage: $0 yyyymmdd hh output

echo

exit 1

fi

yyyymmdd_i=$1

hh_i=$2

output=$3

ss="00"

if [ $hh_i -lt 0 -o $hh_i -gt 23 ]; then

echo

echo ERROR in $0

echo "hh must be 0 <= hh <= 23."

echo

exit 1

fi

if [ $hh_i = "00" ]; then

yyyymmdda[0]=$(date -d"${yyyymmdd_i}"-1days '+%Y%m%d')

hha[0]="23"

mia[0]="10"

yyyymmdda[1]=${yyyymmdda[0]}

hha[1]="23"

mia[1]="20"

yyyymmdda[2]=${yyyymmdda[0]}

hha[2]="23"

mia[2]="30"

yyyymmdda[3]=${yyyymmdda[0]}

hha[3]="23"

mia[3]="40"

yyyymmdda[4]=${yyyymmdda[0]}

hha[4]="23"

mia[4]="50"

yyyymmdda[5]=${yyyymmdd_i}

hha[5]="00"

mia[5]="00"

else

yyyymmdda[0]=${yyyymmdd_i}

hha[0]=$(printf %02d $(expr ${hh_i} - 1))

mia[0]="10"

yyyymmdda[1]=${yyyymmdd_i}

hha[1]=${hha[0]}

mia[1]="20"

yyyymmdda[2]=${yyyymmdd_i}

hha[2]=${hha[0]}

mia[2]="30"

yyyymmdda[3]=${yyyymmdd_i}

hha[3]=${hha[0]}

mia[3]="40"

yyyymmdda[4]=${yyyymmdd_i}

hha[4]=${hha[0]}

mia[4]="50"

yyyymmdda[5]=${yyyymmdd_i}

hha[5]=${hh_i}

mia[5]="00"

fi

mkdir -vp $itmdir1hr

itmdir="intm.data"

mkdir -vp $itmdir

dirs="$rdir $itmdir $itmdir1hr"

for dir in $dirs; do

if [ ! -d $dir ]; then

echo

echo Error in $0 : No such directory, $dir

echo while checking dirs

exit 1

fi

done

prefix="Z__C_RJTD_"

postfix="_RDR_JMAGPV_Ggis1km_Prr10lv_ANAL_grib2.bin"

cwd=$(pwd)

i=0

while [ $i -le 5 ]; do

yyyymmdd=${yyyymmdda[$i]}

yyyy=${yyyymmdd:0:4}

mm=${yyyymmdd:4:2}

dd=${yyyymmdd:6:2}

hh=${hha[$i]}

mi=${mia[$i]}

echo "$yyyy $mm $dd $hh $mi"

indir=$rdir"/"$yyyy"/"$mm"/"$dd

if [ ! -d $indir ]; then

echo

echo Error in $0 : No such directory, $dir

echo

exit 1

fi

yyyymmddhhmiss=${yyyy}${mm}${dd}${hh}${mi}${ss}

input[$i]=${indir}/${prefix}${yyyymmddhhmiss}${postfix}

if [ ! -f ${input[$i]} ]; then

for f in ${indir}/${prefix}${yyyymmddhhmiss}_RDR_JMAGPV__grib2.tar; do

echo

echo "Going to ${indir} ..."

echo

cd ${indir}

tar -xvf ${f}

if [ $? -ne 0 ]; then

echo

echo ERROR in $0 while running tar command.

echo tar file: ${f}

echo

echo "Going back to ${cwd} ..."

echo

cd ${cwd}

echo "ABNORMAL END"

echo

exit 1

fi

done

echo

echo "Going back to ${cwd} ..."

echo

cd ${cwd}

fi

if [ ! -f ${input[$i]} ]; then

echo

echo Error in $0 : No such file, ${input[$i]}

echo

exit 1

fi

intmfile[$i]=${itmdir}/${yyyy}${mm}${dd}_${hh}${mi}${sc}_radar10m.bin

command="./jmaradar2bin ${input[$i]} ${intmfile[$i]}"

echo

echo $command

echo

$command

if [ $? -ne 0 ]; then

echo

echo "Error in $0 while running the following:"

echo " $command"

echo

exit 1

fi

i=$(expr $i + 1)

done #i

i=0

echo

while [ $i -le 5 ]; do

echo "Input file: ${input[$i]}"

i=$(expr $i + 1)

done #i

i=0

echo

while [ $i -le 5 ]; do

echo "Intermidiate file: ${intmfile[$i]}"

i=$(expr $i + 1)

done #i

echo

exe=10min_to_1hour

if [ ! -f $exe ]; then

echo

echo ERROR in $0: No such file, $exe

echo

exit 1

fi

namelist="namelist.${exe}.txt"

cat <<EOF >$namelist

&para

im=${nx},

jm=${ny},

input0="${intmfile[0]}",

input1="${intmfile[1]}"

input2="${intmfile[2]}"

input3="${intmfile[3]}"

input4="${intmfile[4]}"

input5="${intmfile[5]}"

output="${itmdir1h}/${output}"

&end

EOF

command="$exe < $namelist"

echo

echo $command

echo

$exe < $namelist

if [ $? -ne 0 ]; then

echo

echo "Error in $0 while running the following:"

echo " $command"

echo

exit 1

fi

yyyy=${yyyymmdd_i:0:4}

mm=${yyyymmdd_i:4:2}

dd=${yyyymmdd_i:6:2}

hh=${hh_i}

if [ $mm = "01" ]; then mmm="Jan"; fi

if [ $mm = "02" ]; then mmm="Feb"; fi

if [ $mm = "03" ]; then mmm="Mar"; fi

if [ $mm = "04" ]; then mmm="Apr"; fi

if [ $mm = "05" ]; then mmm="May"; fi

if [ $mm = "06" ]; then mmm="Jun"; fi

if [ $mm = "07" ]; then mmm="Jul"; fi

if [ $mm = "08" ]; then mmm="Aug"; fi

if [ $mm = "09" ]; then mmm="Sep"; fi

if [ $mm = "10" ]; then mmm="Oct"; fi

if [ $mm = "11" ]; then mmm="Nov"; fi

if [ $mm = "12" ]; then mmm="Dec"; fi

ctl=r1h.ctl

cat <<EOF>$ctl

dset ${output}

undef -999.9

options template

xdef ${nx} LINEAR ${x0} ${dx}

ydef ${ny} LINEAR ${y0} ${dy}

zdef 1 LEVELS 1000

tdef 1 LINEAR ${hh}:00Z${dd}${mmm}${yyyy} 1hr

vars 1

rr 0 0 rainfall

endvars

EOF

echo

echo "Output time: ${yyyymmdd_i} ${hh_i}"

echo

#echo

#echo GrADS control file, $ctl

#echo

#cat $ctl

#echo

exit 0

----------------------

End of radar.precip.1hr.run.sh

----------------------

======================

10min.to.1hr/10min_to_1hour.f90

======================

program ten2one

! Description:

!

! Host: aofd165.fish.nagasaki-u.ac.jp

! Directory: /work2/kunoki/to_manda_sensei/Tools/JMA_radar/10min.to.1hr

!

! ifortの場合以下のようにオプション-assume bytereclを付けてコンパイルすること。

! $ ifort -assume byterecl hoge.f90

! こうしないと、"recl"がbyte単位ではなく、word単位になってしまう

!

! use

! implicit none

real,allocatable:: rr(:,:),rr1(:,:)

integer im,jm

character(len=500)::input0, input1, input2, input3, input4, input5

character(len=500)::output

namelist /para/im,jm,input0, input1, input2, input3, input4, input5, output

read(*,nml=para)

print '(A,A)','Input: ',trim(input0)

print '(A,A)','Input: ',trim(input1)

print '(A,A)','Input: ',trim(input2)

print '(A,A)','Input: ',trim(input3)

print '(A,A)','Input: ',trim(input4)

print '(A,A)','Input: ',trim(input5)

print '(A,A)','Output: ',trim(output)

print *

allocate(rr(im,jm),rr1(im,jm))

rr(:,:)=0.0

open(10,file=input0,form="unformatted",access="direct",recl=4*im*jm,action="read")

read(10,rec=1)rr1

rr(:,:)=rr(:,:)+rr1(:,:)

close(10)

open(10,file=input1,form="unformatted",access="direct",recl=4*im*jm,action="read")

read(10,rec=1)rr1

rr(:,:)=rr(:,:)+rr1(:,:)

close(10)

open(10,file=input2,form="unformatted",access="direct",recl=4*im*jm,action="read")

read(10,rec=1)rr1

rr(:,:)=rr(:,:)+rr1(:,:)

close(10)

open(10,file=input3,form="unformatted",access="direct",recl=4*im*jm,action="read")

read(10,rec=1)rr1

rr(:,:)=rr(:,:)+rr1(:,:)

close(10)

open(10,file=input4,form="unformatted",access="direct",recl=4*im*jm,action="read")

read(10,rec=1)rr1

rr(:,:)=rr(:,:)+rr1(:,:)

close(10)

open(10,file=input5,form="unformatted",access="direct",recl=4*im*jm,action="read")

read(10,rec=1)rr1

rr(:,:)=rr(:,:)+rr1(:,:)

close(10)

rr(:,:)=rr(:,:)/6.0

! NOTE:

! Since unit of rr1 is mm/hr, not mm/10min, we need to

! divide rr by six.

open(10,file=output,form="unformatted",access="direct",recl=4*im*jm)

write(10,rec=1)rr

close(10)

end program ten2one

----------------------

End of 10min.to.1hr/10min_to_1hour.f90

----------------------

======================

r3h.gs

======================

ctl='r3h.ctl'

'open 'ctl

'cc'

'q dims'

line=sublin(result,5)

datetime=subwrd(line,6)

hh=substr(datetime,1,2)

dd=substr(datetime,4,2)

mmm=substr(datetime,6,3)

yyyy=substr(datetime,9,4)

timeutc=hh'00UTC' dd mmm yyyy

say

say timeutc

say

if(mmm='JAN');mm='01';endif

if(mmm='FEB');mm='02';endif

if(mmm='MAR');mm='03';endif

if(mmm='APR');mm='04';endif

if(mmm='MAY');mm='05';endif

if(mmm='JUN');mm='06';endif

if(mmm='JUL');mm='07';endif

if(mmm='AUG');mm='08';endif

if(mmm='SEP');mm='09';endif

if(mmm='OCT');mm='10';endif

if(mmm='NOV');mm='11';endif

if(mmm='DEC');mm='12';endif

'color 10 100 10 -kind white->white->rainbow'

'set gxout shade2'

'd rr'

'cbarn'

'draw title 'timeutc' R`b3h`n [mm/3hr]'

fig='R3h.'yyyy''mm''dd'_'hh'.eps'

'gxprint 'fig

say

say 'Output: 'fig

say

'close 1'

'quit'

----------------------

End of r3h.gs

----------------------