gnuplot ２つの数値をプロット
#! /bin/tcsh -f
################################################
#plot_axcel1
################################################
read_file=output_axcelrv.dat
plot_file1=axcel_pmin
plot_file2=axcel_avcwc
plot_file3=axcel_vmax
plot_file11=axcel_pmin_avcwc
gnuplot <<EOF
################################################
################################################
plot "${read_file}" using 1:2 with line title "${plot_file1}" linetype 3 linecolor rgbcolor "red" linewidth 3
set title "${plot_file1}"
set xlabel "time step [h]"
set ylabel "[hPa]"
set term postscript
set output "${plot_file1}.ps"
replot
################################################
plot "${read_file}" using 1:6 with line title "${plot_file2}" linetype 3 linecolor rgbcolor "blue" linewidth 3
set title "${plot_file2}"
set xlabel "time step [h]"
set ylabel "[g/m2]"
set term postscript
set output "${plot_file2}.ps"
replot
################################################
plot "${read_file}" using 1:5 with line title "${plot_file3}" linetype 3 linecolor rgbcolor "green" linewidth 3
set title "${plot_file3}"
set xlabel "time step [h]"
set ylabel "[m/s]"
set term postscript
set output "${plot_file3}.ps"
replot
################################################
### double plot#################################
################################################
set yrange [980:1010]
set y2range [0:30]
plot "${read_file}" using 1:2 axis x1y1 with line title "${plot_file1}" linetype 3 linecolor rgbcolor "red" linewidth 3,\
"${read_file}" using 1:5 axis x1y2 with line title "${plot_file3}" linetype 3 linecolor rgbcolor "green" linewidth 3
set xlabel "time step [h]"
set ylabel "[hPa]"
set y2label "[m/s]"
set term postscript
set output "${plot_file11}.ps"
replot
################################################
EOF
################################################
################################################
################################################
################################################
################################################
### make pdf ###################################
################################################
################################################
################################################
################################################
################################################
ps2pdf ${plot_file1}.ps 1_${plot_file1}.pdf
rm ${plot_file1}.ps
ps2pdf ${plot_file2}.ps 2_${plot_file2}.pdf
rm ${plot_file2}.ps
ps2pdf ${plot_file3}.ps 3_${plot_file3}.pdf
rm ${plot_file3}.ps
ps2pdf ${plot_file11}.ps 11_${plot_file11}.pdf
rm ${plot_file11}.ps
2011年10月19日水曜日
axcel 解析 時系列変化
change_axcel1.sh
#! /bin/tcsh -f
#####################################
# cloud change
# chage of axcel
#####################################
#####################################
filename1=change_cloud
filename2=read_data
filename3=output_data
ffilename1=psea
ffilename2=v
ffilename3=allcwc
ffilename4=rv #1e+4
ffilename5=rvsrf #1e+4
ofname1=output_axcel
ofname2=output_axcelrv
####################################
####################################
#model grid size
xgrid=418
ygrid=418
zgrid=16
#include +,=
f_lon=+129.88
f_lat=-45.12
####################################
# set grid for cal
#20km=1grid is 0.18
xgrid_a=56 #10=55.5 grid
ygrid_a=56 #10=55.5 grid
gsize1=20 # 20km or 5km
gsize2=0.18 #grid size 20km=0.18, 5km=0.045
#definition region
lon_l=10.00
lat_l=10.00
#for cal
lon_l2=15.00
lat_l2=15.00
#length from center
ll_km=1000 #[km]
ll_m=`expr $ll_km \* 1000` #[m]
####################################
####################################
#type yyyymmddhhhh start and end time
: ${start:=199201041200}
: ${end:=199201082100}
#: ${end:=199201041500}
year=`echo ${start} | cut -c 1-4`
month=`echo ${start} | cut -c 5-6`
day=`echo ${start} | cut -c 7-8`
hour=`echo ${start} |cut -c 9-10`
h=00
echo start time ${year} ${month} ${day} ${hour} ${h}
echo end time ${end}
: ${time:=${start}}
#timestep
timestep=1
time_step=0
i=0
while [ ${time} -le ${end} ]
do
time_step=`expr ${timestep} \* ${i}`
echo time_step ${time_step}
##decide m
if [ ${month} -eq 01 ]; then
m=JAN
elif [ ${month} -eq 02 ]; then
m=FEB
elif [ ${month} -eq 03 ]; then
m=MAR
elif [ ${month} -eq 04 ]; then
m=APR
elif [ ${month} -eq 05 ]; then
m=MAY
elif [ ${month} -eq 06 ]; then
m=JUN
elif [ ${month} -eq 07 ]; then
m=JUL
elif [ ${month} -eq 08 ]; then
m=AUG
elif [ ${month} -eq 09 ]; then
m=SEP
elif [ ${month} -eq 10 ]; then
m=OCT
elif [ ${month} -eq 11 ]; then
m=NOV
elif [ ${month} -eq 12 ]; then
m=DEC
fi
##decide m
#echo $hour
time=${year}${month}${day}${hour}${h}
#echo yoshiyoshi
echo time $time
echo month $m
t=`expr ${day} \* 4 - 3`
t=`expr ${hour} / 6 + ${t}`
##############################################
#read grads data
: ${prefix:=anl_p} #want data
cat > ${filename2}.gs <<EOF
'reinit'
'open _NHMMRPPFCSVSTD1_${time}.ctl'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename1}.grd'
*'set t ${t}'
'set x 1 418'
'set y 1 418'
'd psea'
'disable fwrite'
'c'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename2}.grd'
*'set t ${t}'
'set x 1 418'
'set y 1 418'
'd sqrt(u*u+v*v)'
'disable fwrite'
'c'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename3}.grd'
*'set t ${t}'
'set x 1 418'
'set y 1 418'
'c=vint(lev(z=1), cwc, 100)'
'd c'
'disable fwrite'
'c'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename4}.grd'
'set x 1 418'
'set y 1 418'
k=1
while (k<=${zgrid})
'set z 'k''
'a=hcurl(u,v)'
'd a*10000'
*say 'yoshiyoshi'
*say k
k=k+1
endwhile
'disable fwrite'
'c'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename5}.grd'
*'set t ${t}'
'set x 1 418'
'set y 1 418'
'a=hcurl(usrf,vsrf)'
'd a*10000'
'disable fwrite'
'c'
*ffilename1=psea
*ffilename2=v
*ffilename3=allcwc
*ffilename4=rv
*ffilename5=rv_srf
'quit'
EOF
grads -blc ${filename2}.gs
########################################
####program read_data###################
########################################
if [ ${time} -eq ${start} ]; then
##start : make output file
cat <<EOF > start.f90
program start
implicit none
open (unit=31, file='${ofname1}.txt', form='FORMATTED')
open (unit=32, file='${ofname2}.txt', form='FORMATTED')
write(31,*) '${ofilename1}'
write(31,*) 'time_step pn plonn platn v1n allarea_avcwc axcel_avcwc'
write(*,*) 'time_step pn plonn platn v1n allarea_avcwc axcel_avcwc'
write(32,*) '${ofilename2}'
write(32,*) 'time_step pn plonn platn v1n axcel_avcwc avrvsrf output_rv1000 output_rv800 output_rv600 output_rv400 output_rv200'
write(*,*) 'time_step pn plonn platn v1n axcel_avcwc avrvsrf output_rv1000 output_rv800 output_rv600 output_rv400 output_rv200'
close(31)
close(32)
end program start
EOF
##output data.txt by f90
f90 -o start.out start.f90
./start.out
##start
fi
####find min psea
cat <<EOF > ${filename2}.f90
program ${filename2}
implicit none
integer, parameter :: xgrid=${xgrid}, ygrid=${ygrid}, zgrid=${zgrid}, rgrid=10
real, parameter :: gsize1=${gsize1}, gsize2=${gsize2}
real, parameter :: ll_km=${ll_km}, ll_m=${ll_m}
integer :: n, i, j, k, step, is, ie, js, je
integer :: kxn, kyn, pkxn, pkyn
integer :: ngrid2, ngrid3, size
!ffilename1=psea
!ffilename2=v
!ffilename3=allcwc
!ffilename4=rv #1e+4
!ffilename5=rv_srf #1e+4
real(kind=4), dimension(xgrid,ygrid) :: ${ffilename1}
real(kind=4), dimension(xgrid,ygrid) :: ${ffilename2}
real(kind=4), dimension(xgrid,ygrid) :: ${ffilename3}
real(kind=4), dimension(xgrid,ygrid,zgrid) :: ${ffilename4}
real(kind=4), dimension(xgrid,ygrid) :: ${ffilename5}
real(kind=4), dimension(zgrid) :: output_rv
real(kind=4) :: pn, ppn, lonn, latn, plonn, platn, v1n, v2n, vlonn, vlatn, vmaxn, x
real(kind=4) :: cwc_sum, avcwc, allarea_avcwc
real(kind=4) :: p, axcel_avcwc, l_axcel
real(kind=4) :: rv_sum, avrv, avrvsrf
ngrid2=xgrid*ygrid
ngrid3=xgrid*ygrid*zgrid
open (unit=11,file='${ffilename1}.grd',form='unformatted',access='direct',recl=ngrid2*4)
read (unit=11,rec=1) ${ffilename1}
open (unit=12,file='${ffilename2}.grd',form='unformatted',access='direct',recl=ngrid2*4)
read (unit=12,rec=1) ${ffilename2}
open (unit=13,file='${ffilename3}.grd',form='unformatted',access='direct',recl=ngrid2*4)
read (unit=13,rec=1) ${ffilename3}
open (unit=14,file='${ffilename4}.grd',form='unformatted',access='direct',recl=ngrid3*4)
read (unit=14,rec=1) ${ffilename4}
open (unit=15,file='${ffilename5}.grd',form='unformatted',access='direct',recl=ngrid2*4)
read (unit=15,rec=1) ${ffilename5}
open (unit=31, file='${ofname1}.txt', position='append', form='FORMATTED')
open (unit=32, file='${ofname2}.txt', position='append', form='FORMATTED')
!!!!!!!check
!do k = 1, ${zgrid}
!write(*,*) 'rv(100,100,k) ', rv(100,100,k)
!enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!find axcel pmin lon,lat
pn=2000
do i = 83, 277 !145~180
do j = 255, 388 !1~25 @NH
ppn= psea(i-1,j) + psea(i+1,j) + psea(i,j-1) + psea(i,j-1)
ppn= ppn / 4
if (psea(i,j) < ppn) then
if (psea(i,j) < pn) then
pn=psea(i,j)
kxn=i
kyn=j
endif
endif
n=n+1
enddo
enddo
pkxn=kxn
pkyn=kyn
!pmin location
plonn=kxn*${gsize2}${f_lon}
platn=kyn*${gsize2}${f_lat}
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!find axcel vmax vlon,vlat
v1n=1
is=kxn-10
js=kyn-10
ie=kxn+10
je=kyn+10
do i = is, ie
do j = js, je
v2n= v(i-1,j) + v(i+1,j) + v(i,j-1) + v(i,j-1)
v2n= v2n / 4
if (v(i,j) > v2n) then
if (v(i,j) > v1n) then
kxn=i
kyn=j
v1n=v(i,j)
endif
endif
enddo
enddo
vlonn=kxn*${gsize2}${f_lon}
vlatn=kyn*${gsize2}${f_lat}
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! cloud water content
allcwc(i,j)=allcwc(i,j)*1000
!change [kg/kg] to [g/kg]
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!allarea_avcwc area1:-25~20,150~200(160E)
cwc_sum=0
step=0
do i = 111, 388
do j= 111, 361
!do i = 1, ${xgrid}
! do j= 1, ${ygrid}
cwc_sum = cwc_sum + allcwc(i,j)
step= step + 1
enddo
enddo
allarea_avcwc = cwc_sum / step
!write(*,*) 'step allarea_avcwc ', step, allarea_avcwc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! a axcel_avcwc
cwc_sum=0
step=0
do i = 1 , ${xgrid}
do j = 1 , ${ygrid}
l_axcel=(i-pkxn)**2+(j-pkyn)**2
l_axcel=(l_axcel)**0.5
l_axcel=l_axcel*${gsize1} !between axcel(pmin) and point
if (l_axcel < ${ll_km}) then
cwc_sum=cwc_sum+ allcwc(i,j)
step = step + 1
!write(*,*) pkxn, pkyn, 'step', 'cwc_sum', step, cwc_sum, l_axcel, i, j
endif
l_axcel=0
enddo
enddo
axcel_avcwc = cwc_sum / step
write(*,*) 'step axcel_avcwc ', step, axcel_avcwc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! cal rv
!rv 1~zgrid
step=0
rv_sum=0
do k = 1, ${zgrid}
do i = 1, ${xgrid}
do j = 1, ${ygrid}
l_axcel=(i-pkxn)**2+(j-pkyn)**2
l_axcel=(l_axcel)**0.5
l_axcel=l_axcel*${gsize1} !between axcel(pmin) and point
if (l_axcel < ${ll_km}) then
rv_sum = rv_sum + rv(i,j,k)
step = step + 1
!write(*,*) 'output_rv(k) rv(i,j,k) step i j k', output_rv(k), rv(i,j,k), step, i, j, k
!write(*,*) pkxn, pkyn, 'step', 'rv_sum', step, rv_sum, l_axcel, i, j
endif
l_axcel=0
enddo
enddo
output_rv(k) = rv_sum / step
!write(*,*) 'output_rv(k) k', output_rv(k), step, k
step=0
rv_sum=0
enddo
!!!!!!
do k = 1, ${zgrid}
write(*,*) 'output_rv(k) step k', output_rv(k), step, k
enddo
!rvsrf
step=0
rv_sum=0
avrv=0
do i = 1, ${xgrid}
do j = 1, ${ygrid}
l_axcel=(i-pkxn)**2+(j-pkyn)**2
l_axcel=(l_axcel)**0.5
l_axcel=l_axcel*${gsize1} !between axcel(pmin) and point
if (l_axcel < ${ll_km}) then
rv_sum=rv_sum + rvsrf(i,j)
step = step + 1
!write(*,*) pkxn, pkyn, 'step', 'rv_sum', step, rv_sum, l_axcel, i, j
endif
l_axcel=0
enddo
enddo
avrvsrf = rv_sum / step
!write(*,*) 'step avrvsrf', rv_sum, avrvsrf
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! output data
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(31,*) ${time_step}, pn, plonn, platn, v1n, allarea_avcwc, axcel_avcwc
write(*,*) ${time_step}, pn, plonn, platn, v1n, allarea_avcwc, axcel_avcwc
write(32,*) ${time_step}, pn, plonn, platn, v1n, axcel_avcwc, avrvsrf, output_rv(1), output_rv(7), output_rv(9), output_rv(11), output_rv(14)
write(*,*) ${time_step}, pn, plonn, platn, v1n, axcel_avcwc, avrvsrf, output_rv(1), output_rv(7), output_rv(9), output_rv(11), output_rv(14)
close(11)
close(12)
close(13)
close(14)
close(15)
close(31)
close(32)
end program
EOF
##output data.txt by f90
f90 -o ${filename2}.out ${filename2}.f90
./${filename2}.out
#echo ${day}
#echo "yoshiyoshi"
i=`expr $i + 1`
#echo $i
#hour=`expr $hour + 1`
hour=`expr $hour + $timestep`
hour=`echo ${hour} | awk '{printf("%02d",$1)}'`
if [ $hour -eq 24 ]; then
hour=00
# day=`expr $day + 1`
day=`echo ${day} | awk '{printf("%02d",$1+1)}'`
if [ $day -ge 32 ];then
day=01
# month=`expr $month + 1`
month=`echo ${month} | awk '{printf("%02d",$1+1)}'`
if [ $month -ge 13 ];then
month=01
year=`expr $year + 1`
fi
fi
fi
time=${year}${month}${day}${hour}${h}
done
cp ${ofname1}.txt ${ofname1}.dat
cp ${ofname1}.txt ${ofname2}.dat
#! /bin/tcsh -f
#####################################
# cloud change
# chage of axcel
#####################################
#####################################
filename1=change_cloud
filename2=read_data
filename3=output_data
ffilename1=psea
ffilename2=v
ffilename3=allcwc
ffilename4=rv #1e+4
ffilename5=rvsrf #1e+4
ofname1=output_axcel
ofname2=output_axcelrv
####################################
####################################
#model grid size
xgrid=418
ygrid=418
zgrid=16
#include +,=
f_lon=+129.88
f_lat=-45.12
####################################
# set grid for cal
#20km=1grid is 0.18
xgrid_a=56 #10=55.5 grid
ygrid_a=56 #10=55.5 grid
gsize1=20 # 20km or 5km
gsize2=0.18 #grid size 20km=0.18, 5km=0.045
#definition region
lon_l=10.00
lat_l=10.00
#for cal
lon_l2=15.00
lat_l2=15.00
#length from center
ll_km=1000 #[km]
ll_m=`expr $ll_km \* 1000` #[m]
####################################
####################################
#type yyyymmddhhhh start and end time
: ${start:=199201041200}
: ${end:=199201082100}
#: ${end:=199201041500}
year=`echo ${start} | cut -c 1-4`
month=`echo ${start} | cut -c 5-6`
day=`echo ${start} | cut -c 7-8`
hour=`echo ${start} |cut -c 9-10`
h=00
echo start time ${year} ${month} ${day} ${hour} ${h}
echo end time ${end}
: ${time:=${start}}
#timestep
timestep=1
time_step=0
i=0
while [ ${time} -le ${end} ]
do
time_step=`expr ${timestep} \* ${i}`
echo time_step ${time_step}
##decide m
if [ ${month} -eq 01 ]; then
m=JAN
elif [ ${month} -eq 02 ]; then
m=FEB
elif [ ${month} -eq 03 ]; then
m=MAR
elif [ ${month} -eq 04 ]; then
m=APR
elif [ ${month} -eq 05 ]; then
m=MAY
elif [ ${month} -eq 06 ]; then
m=JUN
elif [ ${month} -eq 07 ]; then
m=JUL
elif [ ${month} -eq 08 ]; then
m=AUG
elif [ ${month} -eq 09 ]; then
m=SEP
elif [ ${month} -eq 10 ]; then
m=OCT
elif [ ${month} -eq 11 ]; then
m=NOV
elif [ ${month} -eq 12 ]; then
m=DEC
fi
##decide m
#echo $hour
time=${year}${month}${day}${hour}${h}
#echo yoshiyoshi
echo time $time
echo month $m
t=`expr ${day} \* 4 - 3`
t=`expr ${hour} / 6 + ${t}`
##############################################
#read grads data
: ${prefix:=anl_p} #want data
cat > ${filename2}.gs <<EOF
'reinit'
'open _NHMMRPPFCSVSTD1_${time}.ctl'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename1}.grd'
*'set t ${t}'
'set x 1 418'
'set y 1 418'
'd psea'
'disable fwrite'
'c'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename2}.grd'
*'set t ${t}'
'set x 1 418'
'set y 1 418'
'd sqrt(u*u+v*v)'
'disable fwrite'
'c'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename3}.grd'
*'set t ${t}'
'set x 1 418'
'set y 1 418'
'c=vint(lev(z=1), cwc, 100)'
'd c'
'disable fwrite'
'c'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename4}.grd'
'set x 1 418'
'set y 1 418'
k=1
while (k<=${zgrid})
'set z 'k''
'a=hcurl(u,v)'
'd a*10000'
*say 'yoshiyoshi'
*say k
k=k+1
endwhile
'disable fwrite'
'c'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename5}.grd'
*'set t ${t}'
'set x 1 418'
'set y 1 418'
'a=hcurl(usrf,vsrf)'
'd a*10000'
'disable fwrite'
'c'
*ffilename1=psea
*ffilename2=v
*ffilename3=allcwc
*ffilename4=rv
*ffilename5=rv_srf
'quit'
EOF
grads -blc ${filename2}.gs
########################################
####program read_data###################
########################################
if [ ${time} -eq ${start} ]; then
##start : make output file
cat <<EOF > start.f90
program start
implicit none
open (unit=31, file='${ofname1}.txt', form='FORMATTED')
open (unit=32, file='${ofname2}.txt', form='FORMATTED')
write(31,*) '${ofilename1}'
write(31,*) 'time_step pn plonn platn v1n allarea_avcwc axcel_avcwc'
write(*,*) 'time_step pn plonn platn v1n allarea_avcwc axcel_avcwc'
write(32,*) '${ofilename2}'
write(32,*) 'time_step pn plonn platn v1n axcel_avcwc avrvsrf output_rv1000 output_rv800 output_rv600 output_rv400 output_rv200'
write(*,*) 'time_step pn plonn platn v1n axcel_avcwc avrvsrf output_rv1000 output_rv800 output_rv600 output_rv400 output_rv200'
close(31)
close(32)
end program start
EOF
##output data.txt by f90
f90 -o start.out start.f90
./start.out
##start
fi
####find min psea
cat <<EOF > ${filename2}.f90
program ${filename2}
implicit none
integer, parameter :: xgrid=${xgrid}, ygrid=${ygrid}, zgrid=${zgrid}, rgrid=10
real, parameter :: gsize1=${gsize1}, gsize2=${gsize2}
real, parameter :: ll_km=${ll_km}, ll_m=${ll_m}
integer :: n, i, j, k, step, is, ie, js, je
integer :: kxn, kyn, pkxn, pkyn
integer :: ngrid2, ngrid3, size
!ffilename1=psea
!ffilename2=v
!ffilename3=allcwc
!ffilename4=rv #1e+4
!ffilename5=rv_srf #1e+4
real(kind=4), dimension(xgrid,ygrid) :: ${ffilename1}
real(kind=4), dimension(xgrid,ygrid) :: ${ffilename2}
real(kind=4), dimension(xgrid,ygrid) :: ${ffilename3}
real(kind=4), dimension(xgrid,ygrid,zgrid) :: ${ffilename4}
real(kind=4), dimension(xgrid,ygrid) :: ${ffilename5}
real(kind=4), dimension(zgrid) :: output_rv
real(kind=4) :: pn, ppn, lonn, latn, plonn, platn, v1n, v2n, vlonn, vlatn, vmaxn, x
real(kind=4) :: cwc_sum, avcwc, allarea_avcwc
real(kind=4) :: p, axcel_avcwc, l_axcel
real(kind=4) :: rv_sum, avrv, avrvsrf
ngrid2=xgrid*ygrid
ngrid3=xgrid*ygrid*zgrid
open (unit=11,file='${ffilename1}.grd',form='unformatted',access='direct',recl=ngrid2*4)
read (unit=11,rec=1) ${ffilename1}
open (unit=12,file='${ffilename2}.grd',form='unformatted',access='direct',recl=ngrid2*4)
read (unit=12,rec=1) ${ffilename2}
open (unit=13,file='${ffilename3}.grd',form='unformatted',access='direct',recl=ngrid2*4)
read (unit=13,rec=1) ${ffilename3}
open (unit=14,file='${ffilename4}.grd',form='unformatted',access='direct',recl=ngrid3*4)
read (unit=14,rec=1) ${ffilename4}
open (unit=15,file='${ffilename5}.grd',form='unformatted',access='direct',recl=ngrid2*4)
read (unit=15,rec=1) ${ffilename5}
open (unit=31, file='${ofname1}.txt', position='append', form='FORMATTED')
open (unit=32, file='${ofname2}.txt', position='append', form='FORMATTED')
!!!!!!!check
!do k = 1, ${zgrid}
!write(*,*) 'rv(100,100,k) ', rv(100,100,k)
!enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!find axcel pmin lon,lat
pn=2000
do i = 83, 277 !145~180
do j = 255, 388 !1~25 @NH
ppn= psea(i-1,j) + psea(i+1,j) + psea(i,j-1) + psea(i,j-1)
ppn= ppn / 4
if (psea(i,j) < ppn) then
if (psea(i,j) < pn) then
pn=psea(i,j)
kxn=i
kyn=j
endif
endif
n=n+1
enddo
enddo
pkxn=kxn
pkyn=kyn
!pmin location
plonn=kxn*${gsize2}${f_lon}
platn=kyn*${gsize2}${f_lat}
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!find axcel vmax vlon,vlat
v1n=1
is=kxn-10
js=kyn-10
ie=kxn+10
je=kyn+10
do i = is, ie
do j = js, je
v2n= v(i-1,j) + v(i+1,j) + v(i,j-1) + v(i,j-1)
v2n= v2n / 4
if (v(i,j) > v2n) then
if (v(i,j) > v1n) then
kxn=i
kyn=j
v1n=v(i,j)
endif
endif
enddo
enddo
vlonn=kxn*${gsize2}${f_lon}
vlatn=kyn*${gsize2}${f_lat}
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! cloud water content
allcwc(i,j)=allcwc(i,j)*1000
!change [kg/kg] to [g/kg]
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!allarea_avcwc area1:-25~20,150~200(160E)
cwc_sum=0
step=0
do i = 111, 388
do j= 111, 361
!do i = 1, ${xgrid}
! do j= 1, ${ygrid}
cwc_sum = cwc_sum + allcwc(i,j)
step= step + 1
enddo
enddo
allarea_avcwc = cwc_sum / step
!write(*,*) 'step allarea_avcwc ', step, allarea_avcwc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! a axcel_avcwc
cwc_sum=0
step=0
do i = 1 , ${xgrid}
do j = 1 , ${ygrid}
l_axcel=(i-pkxn)**2+(j-pkyn)**2
l_axcel=(l_axcel)**0.5
l_axcel=l_axcel*${gsize1} !between axcel(pmin) and point
if (l_axcel < ${ll_km}) then
cwc_sum=cwc_sum+ allcwc(i,j)
step = step + 1
!write(*,*) pkxn, pkyn, 'step', 'cwc_sum', step, cwc_sum, l_axcel, i, j
endif
l_axcel=0
enddo
enddo
axcel_avcwc = cwc_sum / step
write(*,*) 'step axcel_avcwc ', step, axcel_avcwc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! cal rv
!rv 1~zgrid
step=0
rv_sum=0
do k = 1, ${zgrid}
do i = 1, ${xgrid}
do j = 1, ${ygrid}
l_axcel=(i-pkxn)**2+(j-pkyn)**2
l_axcel=(l_axcel)**0.5
l_axcel=l_axcel*${gsize1} !between axcel(pmin) and point
if (l_axcel < ${ll_km}) then
rv_sum = rv_sum + rv(i,j,k)
step = step + 1
!write(*,*) 'output_rv(k) rv(i,j,k) step i j k', output_rv(k), rv(i,j,k), step, i, j, k
!write(*,*) pkxn, pkyn, 'step', 'rv_sum', step, rv_sum, l_axcel, i, j
endif
l_axcel=0
enddo
enddo
output_rv(k) = rv_sum / step
!write(*,*) 'output_rv(k) k', output_rv(k), step, k
step=0
rv_sum=0
enddo
!!!!!!
do k = 1, ${zgrid}
write(*,*) 'output_rv(k) step k', output_rv(k), step, k
enddo
!rvsrf
step=0
rv_sum=0
avrv=0
do i = 1, ${xgrid}
do j = 1, ${ygrid}
l_axcel=(i-pkxn)**2+(j-pkyn)**2
l_axcel=(l_axcel)**0.5
l_axcel=l_axcel*${gsize1} !between axcel(pmin) and point
if (l_axcel < ${ll_km}) then
rv_sum=rv_sum + rvsrf(i,j)
step = step + 1
!write(*,*) pkxn, pkyn, 'step', 'rv_sum', step, rv_sum, l_axcel, i, j
endif
l_axcel=0
enddo
enddo
avrvsrf = rv_sum / step
!write(*,*) 'step avrvsrf', rv_sum, avrvsrf
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! output data
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(31,*) ${time_step}, pn, plonn, platn, v1n, allarea_avcwc, axcel_avcwc
write(*,*) ${time_step}, pn, plonn, platn, v1n, allarea_avcwc, axcel_avcwc
write(32,*) ${time_step}, pn, plonn, platn, v1n, axcel_avcwc, avrvsrf, output_rv(1), output_rv(7), output_rv(9), output_rv(11), output_rv(14)
write(*,*) ${time_step}, pn, plonn, platn, v1n, axcel_avcwc, avrvsrf, output_rv(1), output_rv(7), output_rv(9), output_rv(11), output_rv(14)
close(11)
close(12)
close(13)
close(14)
close(15)
close(31)
close(32)
end program
EOF
##output data.txt by f90
f90 -o ${filename2}.out ${filename2}.f90
./${filename2}.out
#echo ${day}
#echo "yoshiyoshi"
i=`expr $i + 1`
#echo $i
#hour=`expr $hour + 1`
hour=`expr $hour + $timestep`
hour=`echo ${hour} | awk '{printf("%02d",$1)}'`
if [ $hour -eq 24 ]; then
hour=00
# day=`expr $day + 1`
day=`echo ${day} | awk '{printf("%02d",$1+1)}'`
if [ $day -ge 32 ];then
day=01
# month=`expr $month + 1`
month=`echo ${month} | awk '{printf("%02d",$1+1)}'`
if [ $month -ge 13 ];then
month=01
year=`expr $year + 1`
fi
fi
fi
time=${year}${month}${day}${hour}${h}
done
cp ${ofname1}.txt ${ofname1}.dat
cp ${ofname1}.txt ${ofname2}.dat
2011年10月10日月曜日
シェルスクリプト 四則演算
四則演算
シェルスクリプトで数値の演算を行いたい場合は、「expr」コマンドを使用します。例えは、変数に5と3を足した数値を格納したい場合、「変 数=5+3」と書くことはできません。変数には「5+3」という文字列が格納されることになります。シェルスクリプトでは、「変数=`expr 5 + 3`」と書きます。「expr」コマンドで使用できる算術演算子は以下のとおりです。
| 算術演算子 | 意味 |
|---|---|
| a + b | aとbの和 |
| a - b | aとbの差 |
| a \* b | aとbの積 |
| a / b | aとbの商 |
| a % b | aとbの剰余 |
expr_test.sh
#!/bin/sh
a=`expr 5 + 3`
b=`expr 5 - 3`
ab=`expr $a + $b`
c=`expr 5 \* 3`
d=`expr 5 / 3`
e=`expr 5 % 3`
echo "a=$a"
echo "b=$b"
echo "ab=$ab"
echo "c=$c"
echo "d=$d"
echo "e=$e"
a=`expr 5 + 3`
b=`expr 5 - 3`
ab=`expr $a + $b`
c=`expr 5 \* 3`
d=`expr 5 / 3`
e=`expr 5 % 3`
echo "a=$a"
echo "b=$b"
echo "ab=$ab"
echo "c=$c"
echo "d=$d"
echo "e=$e"
gnuplotをシェルスクリプトで動かす
plot_cloudwater1.sh
#! /bin/tcsh -f
################################################
#plot cloudwater
################################################
read_file=output_data.dat
plot_file1=allarea_avcwc
gnuplot <<EOF
plot "${read_file}" using 1:10 with line title "${plot_file1}" linetype 2 linecolor rgbcolor "blue" linewidth 3
set title "${plot_file1}"
set xlabel "time step [h]"
set ylabel "[g/m2]"
set term postscript
set output "${plot_file1}.ps"
replot
! is comentout
EOF
ps2pdf ${plot_file1}.ps 1_${plot_file1}.pdf
rm ${plot_file1}.ps
#! /bin/tcsh -f
################################################
#plot cloudwater
################################################
read_file=output_data.dat
plot_file1=allarea_avcwc
gnuplot <<EOF
plot "${read_file}" using 1:10 with line title "${plot_file1}" linetype 2 linecolor rgbcolor "blue" linewidth 3
set title "${plot_file1}"
set xlabel "time step [h]"
set ylabel "[g/m2]"
set term postscript
set output "${plot_file1}.ps"
replot
! is comentout
EOF
ps2pdf ${plot_file1}.ps 1_${plot_file1}.pdf
rm ${plot_file1}.ps
gnuplot スペシャル
基本的なgnuplotのサイト
http://t16web.lanl.gov/Kawano/gnuplot/intro/basic.html
gnuplotでplotする色、太さの変更
http://d.hatena.ne.jp/yrgnah_yats/20090802/1249188972
より
gnuplotで線の色を変える
epsにすると線幅が細くなるので、太くしたければlw 10と大きめにする。lc rgb "[color name]"で用いる色の名前[color name]については
List of known color names:
■white #ffffff = 255 255 255
■black #000000 = 0 0 0
■grey #bebebe = 190 190 190
■red #ff0000 = 255 0 0
■yellow #ffff00 = 255 255 0
■green #00ff00 = 0 255 0
■light-green #90ee90 = 144 238 144
■dark-green #006400 = 0 100 0
■blue #0000ff = 0 0 255
■cyan #00ffff = 0 255 255
■magenta #ff00ff = 255 0 255
■turquoise #40e0d0 = 64 224 208
■pink #ffc0cb = 255 192 203
■salmon #fa8072 = 250 128 114
■khaki #f0e68c = 240 230 140
■orange #ffa500 = 255 165 0
■purple #a020f0 = 160 32 240
さらに線種の一覧も出してみた。
同様に
では線種の比較をやっている。pngやgifで出力すればすべて実線になるが、epsだと破線などになる。
また
http://t16web.lanl.gov/Kawano/gnuplot/intro/basic.html
gnuplotでplotする色、太さの変更
http://d.hatena.ne.jp/yrgnah_yats/20090802/1249188972
より
gnuplotで線の色を変える
線種の変更
例えば赤い破線で書きたければ次のようにすればよい。plot sin(x) linetype 2 linecolor rgbcolor "red" linewidth 10 plot sin(x) lt 2 lc rgb "red" lw 10 plot sin(x) lt 2 lc rgb "#ff0000" lw 10 plot sin(x) lt 2 lc 1ltはlinetype(線種：実線・破線・色・太さが決められている)、lcはlinecolor、lwはlinewidthを表す。rgbはrgbcolorの省略形。ltは-1から10まで指定できる。lwは太さで1-20の範囲で自分は使っている。もう一つ例を示す。
plot sin(x) lt 3 lc rgb "yellow" lw 10, cos(x) lt 4 lc rgb "green" lw 15, \ 0.5 lt 2 lc rgb "blue" lw 2 , -0.5 lt 1 lc rgb "salmon" lw 5 set term post eps enhanced color set output "test.eps" replot
epsにすると線幅が細くなるので、太くしたければlw 10と大きめにする。lc rgb "[color name]"で用いる色の名前[color name]については
show palette colornamesとすれば一覧が表示される。一部紹介。
List of known color names:
■white #ffffff = 255 255 255
■black #000000 = 0 0 0
■grey #bebebe = 190 190 190
■red #ff0000 = 255 0 0
■yellow #ffff00 = 255 255 0
■green #00ff00 = 0 255 0
■light-green #90ee90 = 144 238 144
■dark-green #006400 = 0 100 0
■blue #0000ff = 0 0 255
■cyan #00ffff = 0 255 255
■magenta #ff00ff = 255 0 255
■turquoise #40e0d0 = 64 224 208
■pink #ffc0cb = 255 192 203
■salmon #fa8072 = 250 128 114
■khaki #f0e68c = 240 230 140
■orange #ffa500 = 255 165 0
■purple #a020f0 = 160 32 240
さらに線種の一覧も出してみた。
plot -1 lt -1 lw 10,\ 0 lt 0 lw 10,\ 1 lt 1 lw 10,\ 2 lt 2 lw 10,\ 3 lt 3 lw 10,\ 4 lt 4 lw 10,\ 5 lt 5 lw 10,\ 6 lt 6 lw 10,\ 7 lt 7 lw 10,\ 8 lt 8 lw 10,\ 9 lt 9 lw 10,\ 10 lt 10 lw 10 set term post eps enhanced color set output "test.eps" replot
同様に
では線種の比較をやっている。pngやgifで出力すればすべて実線になるが、epsだと破線などになる。
また
testと打ち込めば線種の比較ができる。
2011年7月31日日曜日
番号一覧 マーク、色
番号一覧
色
| 番号 | 名前 | 色 |
|---|---|---|
| 0 | background | ■■■ |
| 1 | foreground | ■■■ |
| 2 | red | ■■■ |
| 3 | green | ■■■ |
| 4 | dark blue | ■■■ |
| 5 | light blue | ■■■ |
| 6 | magenta | ■■■ |
| 7 | yellow | ■■■ |
| 8 | orange | ■■■ |
| 9 | purple | ■■■ |
| 10 | yellow/green | ■■■ |
| 11 | medium blue | ■■■ |
| 12 | dark yellow | ■■■ |
| 13 | aqua | ■■■ |
| 14 | dark purple | ■■■ |
| 15 | gray | ■■■ |
マークの種類
| 番号 | マーク | |||||||||||||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | (なし) | |||||||||||||||||||||||||||||||||||||||
| 1 | ＋ | |||||||||||||||||||||||||||||||||||||||
| 2 | ○ | |||||||||||||||||||||||||||||||||||||||
| 3 | ● | |||||||||||||||||||||||||||||||||||||||
| 4 | □ | |||||||||||||||||||||||||||||||||||||||
| 5 | ■ | |||||||||||||||||||||||||||||||||||||||
| 6 | × | |||||||||||||||||||||||||||||||||||||||
| 7 | ◇ | |||||||||||||||||||||||||||||||||||||||
| 8 | △ | |||||||||||||||||||||||||||||||||||||||
| 9 | ▲ | |||||||||||||||||||||||||||||||||||||||
| 10 | ○ | |||||||||||||||||||||||||||||||||||||||
| 11 | ● |
登録:
投稿 (Atom)