awk
- 行単位でデータファイルを様々に加工する awk(おーく)は、入力データ１行ごとに各列の様々な処理を行うコマンドである。 ここで、行とはファイルデータの横のライン、列とは縦のラインである。 awkを用いると、例えば「ファイルの２列目を２倍する」というような操作もコマンドライン一行で可能になる。
基本形は、
% awk '条件文{実行文}' ファイル名である。'は「Shiftキー＋数字の７」である。括弧は「中括弧」{}でなくてはならない。 また、条件文は省略可能である。 例えば、３列あるデータファイルのうち、１列目はそのまま、２列目は２倍、３列目は３で割りたいというときは、% awk '{print $1,$2*2.0,$3/3.0}' datafileとする。ここで、print文中の$1は１列目を表す。各列は,で区切る。print文中ではC言語の主な演算記号が使える。 以下に、awkによる処理の例を挙げる。% awk '$1 > 4.0{print}' datafile ある条件に一致した行のみを出力する。この例は一列目が4.0よりも大きければ出力。 print文で列を指定しなければ全列が出力される。$0と明示的に指定することも可能。% awk '$1 <= 5.0 {print $3} $1 > 10.0 {print $2}' datafile 複数の条件によって出力パターンを変える。 例では、１列目が5以上なら3列目を出力し、1列目が10より大きければ2列目を出力。% awk `$1 == "S"{print}` datafile 条件文の中で比較対象が文字の場合は、""で括る。% awk '$1 !~ "#"{print}' datafile 読み飛ばしたい行には行頭に#をつけ、#をコメントアウト行と認識させる。 #で始まる行はスキップする。正確には、1列目に#が含まれていなければ出力する。 ~はその列にある文字が含まれている時という条件を表し、 !は直後の条件を否定（not）する。% awk '{printf("%d %f %s\n",$1,$2,$3)}' datafile 出力書式を指定したいときに用いる。 printfの用法はC言語のprintfに準ずる。% awk 'NF >0 {print $1}' datafile 空行をスキップさせる。 NFはawkの変数で、現在読み込んでいる行の列(フィールド)の数を表す。% awk '{if ($1 == "B" || $1 == "N") print "B"; else print "C"}' datafile if/else文。複雑な条件文を組み込む。 if文以下の||は「または」(or)を表す。「かつ」(and)は&&。% awk '{if ($1 !~ "#") if ($1 == 0 && $2 == 0) print $0,0.5,"r"; else print $0,$1/($1+$2); else print "# is escaped."}' datafile 同じくif/else文。二重に入れ子になっているとき。 上記の例は二行に分けてあるが、実際は一行につなげて書く。% awk 'NR==1{s0=$1 } NR!=1{r=$1-s0; s0=$1; print r}' datafile 前後の行での処理。例は、１列目において、１行前の値を今の行の値から引いて出力する。 １行目では今の値を保持させるのみとし、２行目以降で引き算を行っている。 NRはNFと同様にawkの変数で、処理をした行の数を表す。 s0というのは、自分で設定した変数。ここに前の行の一列目の値を保持させている。awk内の変数を挙げる。- NF ：列の数
- NR ：行の数。
- FILENAME ：今読み込んでいるファイルの名前
2011年11月10日木曜日
使いこなせ！シェルスクリプト！！データ解析
http://www.e.ics.nara-wu.ac.jp/~nogu/tips/unix_command.html
2011年10月26日水曜日
chage_axcel6
#! /bin/tcsh -f
#####################################
# cloud change
# chage of axcel
#####################################
#####################################
filename1=change
filename2=read_data6
filename3=output_data
ffilename1=psea
ffilename2=u
ffilename3=v
ffilename4=srf_u
ffilename5=srf_v
ofname1=output_axcel
ofname2=output_axcel_c
####################################
####################################
#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'
*psea
'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'
*u
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename2}.grd'
*'set t ${t}'
'set x 1 418'
'set y 1 418'
k=1
while (k<=${zgrid})
'set z 'k''
'd u'
*say 'yoshiyoshi'
*say k
k=k+1
endwhile
'disable fwrite'
'c'
*v
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename3}.grd'
*'set t ${t}'
'set x 1 418'
'set y 1 418'
k=1
while (k<=${zgrid})
'set z 'k''
'd v'
*say 'yoshiyoshi'
*say k
k=k+1
endwhile
'disable fwrite'
'c'
*srf_u
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename4}.grd'
*'set t ${t}'
'set x 1 418'
'set y 1 418'
'd usrf'
'disable fwrite'
'c'
*srf_v
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename5}.grd'
*'set t ${t}'
'set x 1 418'
'set y 1 418'
'd vsrf'
'disable fwrite'
'c'
'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'
write(*,*) 'time_step pn plonn platn v1n'
write(32,*) '${ofilename2}'
write(32,*) 'time_step pn plonn platn output_srfc output_c1000 output_c800 output_c600 output_c400 output_c200'
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
!filename1=change
!filename2=read_data
!filename3=output_data
!ffilename1=psea
!ffilename2=u
!ffilename3=v
!ffilename4=srf_u
!ffilename5=srf_v
real(kind=4), dimension(xgrid,ygrid) :: ${ffilename1}
real(kind=4), dimension(xgrid,ygrid,zgrid) :: ${ffilename2}
real(kind=4), dimension(xgrid,ygrid,zgrid) :: ${ffilename3}
real(kind=4), dimension(xgrid,ygrid) :: ${ffilename4}
real(kind=4), dimension(xgrid,ygrid) :: ${ffilename5}
real(kind=4), dimension(zgrid) :: output_c
real(kind=4), dimension(ygrid) :: corr_f
real(kind=4), dimension(xgrid,ygrid,zgrid) :: zcurl
real(kind=4), dimension(xgrid,ygrid) :: zcurls
real(kind=4) :: p, axcel_avcwc, l_axcel
real(kind=4) :: pn, ppn, lonn, latn, plonn, platn, v1n, v2n, vlonn, vlatn, vmaxn, x
real(kind=4) :: zcurl_sum, s_sum, cx, cy, output_srfc
real(kind=4) :: f, rad, pai, omega
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=ngrid3*4)
read (unit=12,rec=1) ${ffilename2}
open (unit=13,file='${ffilename3}.grd',form='unformatted',access='direct',recl=ngrid3*4)
read (unit=13,rec=1) ${ffilename3}
open (unit=14,file='${ffilename4}.grd',form='unformatted',access='direct',recl=ngrid2*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 u
!do k = 1, ${zgrid}
!write(*,*) 'u(300,300,k) ', u(300,300,k)
!enddo
!!!!!!!check v
!do k = 1, ${zgrid}
!write(*,*) 'v(100,100,k) ', v(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}
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! cal circulation
!!!! centered difference
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!set coriolis
!rad=30*3.14/180
!f=sin(rad)
pai=3.141592
omega=2*pai/60/60/24
write(*,*) '!!!set coriolis'
do j = 1, ${ygrid}
rad=j*${gsize2}
rad=rad${f_lat}
rad=rad*pai/180
corr_f(j)=2*omega*sin(rad)
!write(*,*) 'corr_f', corr_f(j)
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!cal zcurl
cx=${xgrid}-1
cy=${ygrid}-1
do k = 1, ${zgrid}
do i = 2, cx
do j = 2, cy
zcurl(i,j,k)=(v(i+1,j,k)-v(i-1,j,k))-(u(i,j+1,k)-u(i,j-1,k))
zcurl(i,j,k)=zcurl(i,j,k)/(2*${gsize1})
zcurl(i,j,k)= zcurl(i,j,k) + corr_f(j)
enddo
enddo
enddo
!do k = 1, 16
!write(*,*) 'zcurl(300,300,k) 300,300,k', zcurl(300,300,k), 300, 300, k
!enddo
zcurl_sum=0
do k = 1, ${zgrid}
do i = 2, ${xgrid}
do j = 2, ${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
zcurl_sum = zcurl_sum + zcurl(i,j,k)
step = step + 1
! write(*,*) 'zcurl_sum zcurl(i,j,k) step i j k', zcurl_sum, zcurl(i,j,k), step, i, j, k
! write(*,*) pkxn, pkyn, 'step', 'cwc_sum', step, cwc_sum, l_axcel, i, j
!write(*,*) 'yoshiyoshi'
endif
l_axcel=0
enddo
enddo
output_c(k)=zcurl_sum*(step*${gsize1}*${gsize1})
!write(*,*) 'step output_c(k)', step, output_c(k)
zcurl_sum=0
step=0
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!cal srf_zcurl
cx=${xgrid}-1
cy=${ygrid}-1
do i = 2, cx
do j = 2, cy
zcurls(i,j)=(srf_v(i+1,j)-srf_v(i-1,j))-(srf_u(i,j+1)-srf_u(i,j-1))
zcurls(i,j)=zcurls(i,j)/(2*${gsize1})
zcurls(i,j)= zcurls(i,j) + corr_f(j)
enddo
enddo
!do k = 1, 16
!write(*,*) 'zcurl(300,300,k) 300,300,k', zcurl(300,300,k), 300, 300, k
!enddo
zcurl_sum=0
do i = 2, ${xgrid}
do j = 2, ${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
zcurl_sum = zcurl_sum + zcurls(i,j)
step = step + 1
! write(*,*) 'zcurl_sum zcurls(i,j) step i j k', zcurl_sum, zcurls(i,j), step, i, j, k
! write(*,*) pkxn, pkyn, 'step', 'cwc_sum', step, cwc_sum, l_axcel, i, j
!write(*,*) 'yoshiyoshi'
endif
l_axcel=0
enddo
enddo
output_srfc=zcurl_sum*(step*${gsize1}*${gsize1})
!write(*,*) 'step output_srfc(k)', step, output_srfc
zcurl_sum=0
step=0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! output data
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(*,*) ${time_step}, pn, plonn, platn, output_srfc, output_c(1), output_c(7), output_c(9), output_c(11), output_c(14)
write(32,*) ${time_step}, pn, plonn, platn, output_srfc, output_c(1), output_c(7), output_c(9), output_c(11), output_c(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 ${ofname2}.txt ${ofname2}.dat
#####################################
# cloud change
# chage of axcel
#####################################
#####################################
filename1=change
filename2=read_data6
filename3=output_data
ffilename1=psea
ffilename2=u
ffilename3=v
ffilename4=srf_u
ffilename5=srf_v
ofname1=output_axcel
ofname2=output_axcel_c
####################################
####################################
#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'
*psea
'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'
*u
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename2}.grd'
*'set t ${t}'
'set x 1 418'
'set y 1 418'
k=1
while (k<=${zgrid})
'set z 'k''
'd u'
*say 'yoshiyoshi'
*say k
k=k+1
endwhile
'disable fwrite'
'c'
*v
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename3}.grd'
*'set t ${t}'
'set x 1 418'
'set y 1 418'
k=1
while (k<=${zgrid})
'set z 'k''
'd v'
*say 'yoshiyoshi'
*say k
k=k+1
endwhile
'disable fwrite'
'c'
*srf_u
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename4}.grd'
*'set t ${t}'
'set x 1 418'
'set y 1 418'
'd usrf'
'disable fwrite'
'c'
*srf_v
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename5}.grd'
*'set t ${t}'
'set x 1 418'
'set y 1 418'
'd vsrf'
'disable fwrite'
'c'
'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'
write(*,*) 'time_step pn plonn platn v1n'
write(32,*) '${ofilename2}'
write(32,*) 'time_step pn plonn platn output_srfc output_c1000 output_c800 output_c600 output_c400 output_c200'
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
!filename1=change
!filename2=read_data
!filename3=output_data
!ffilename1=psea
!ffilename2=u
!ffilename3=v
!ffilename4=srf_u
!ffilename5=srf_v
real(kind=4), dimension(xgrid,ygrid) :: ${ffilename1}
real(kind=4), dimension(xgrid,ygrid,zgrid) :: ${ffilename2}
real(kind=4), dimension(xgrid,ygrid,zgrid) :: ${ffilename3}
real(kind=4), dimension(xgrid,ygrid) :: ${ffilename4}
real(kind=4), dimension(xgrid,ygrid) :: ${ffilename5}
real(kind=4), dimension(zgrid) :: output_c
real(kind=4), dimension(ygrid) :: corr_f
real(kind=4), dimension(xgrid,ygrid,zgrid) :: zcurl
real(kind=4), dimension(xgrid,ygrid) :: zcurls
real(kind=4) :: p, axcel_avcwc, l_axcel
real(kind=4) :: pn, ppn, lonn, latn, plonn, platn, v1n, v2n, vlonn, vlatn, vmaxn, x
real(kind=4) :: zcurl_sum, s_sum, cx, cy, output_srfc
real(kind=4) :: f, rad, pai, omega
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=ngrid3*4)
read (unit=12,rec=1) ${ffilename2}
open (unit=13,file='${ffilename3}.grd',form='unformatted',access='direct',recl=ngrid3*4)
read (unit=13,rec=1) ${ffilename3}
open (unit=14,file='${ffilename4}.grd',form='unformatted',access='direct',recl=ngrid2*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 u
!do k = 1, ${zgrid}
!write(*,*) 'u(300,300,k) ', u(300,300,k)
!enddo
!!!!!!!check v
!do k = 1, ${zgrid}
!write(*,*) 'v(100,100,k) ', v(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}
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! cal circulation
!!!! centered difference
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!set coriolis
!rad=30*3.14/180
!f=sin(rad)
pai=3.141592
omega=2*pai/60/60/24
write(*,*) '!!!set coriolis'
do j = 1, ${ygrid}
rad=j*${gsize2}
rad=rad${f_lat}
rad=rad*pai/180
corr_f(j)=2*omega*sin(rad)
!write(*,*) 'corr_f', corr_f(j)
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!cal zcurl
cx=${xgrid}-1
cy=${ygrid}-1
do k = 1, ${zgrid}
do i = 2, cx
do j = 2, cy
zcurl(i,j,k)=(v(i+1,j,k)-v(i-1,j,k))-(u(i,j+1,k)-u(i,j-1,k))
zcurl(i,j,k)=zcurl(i,j,k)/(2*${gsize1})
zcurl(i,j,k)= zcurl(i,j,k) + corr_f(j)
enddo
enddo
enddo
!do k = 1, 16
!write(*,*) 'zcurl(300,300,k) 300,300,k', zcurl(300,300,k), 300, 300, k
!enddo
zcurl_sum=0
do k = 1, ${zgrid}
do i = 2, ${xgrid}
do j = 2, ${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
zcurl_sum = zcurl_sum + zcurl(i,j,k)
step = step + 1
! write(*,*) 'zcurl_sum zcurl(i,j,k) step i j k', zcurl_sum, zcurl(i,j,k), step, i, j, k
! write(*,*) pkxn, pkyn, 'step', 'cwc_sum', step, cwc_sum, l_axcel, i, j
!write(*,*) 'yoshiyoshi'
endif
l_axcel=0
enddo
enddo
output_c(k)=zcurl_sum*(step*${gsize1}*${gsize1})
!write(*,*) 'step output_c(k)', step, output_c(k)
zcurl_sum=0
step=0
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!cal srf_zcurl
cx=${xgrid}-1
cy=${ygrid}-1
do i = 2, cx
do j = 2, cy
zcurls(i,j)=(srf_v(i+1,j)-srf_v(i-1,j))-(srf_u(i,j+1)-srf_u(i,j-1))
zcurls(i,j)=zcurls(i,j)/(2*${gsize1})
zcurls(i,j)= zcurls(i,j) + corr_f(j)
enddo
enddo
!do k = 1, 16
!write(*,*) 'zcurl(300,300,k) 300,300,k', zcurl(300,300,k), 300, 300, k
!enddo
zcurl_sum=0
do i = 2, ${xgrid}
do j = 2, ${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
zcurl_sum = zcurl_sum + zcurls(i,j)
step = step + 1
! write(*,*) 'zcurl_sum zcurls(i,j) step i j k', zcurl_sum, zcurls(i,j), step, i, j, k
! write(*,*) pkxn, pkyn, 'step', 'cwc_sum', step, cwc_sum, l_axcel, i, j
!write(*,*) 'yoshiyoshi'
endif
l_axcel=0
enddo
enddo
output_srfc=zcurl_sum*(step*${gsize1}*${gsize1})
!write(*,*) 'step output_srfc(k)', step, output_srfc
zcurl_sum=0
step=0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! output data
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(*,*) ${time_step}, pn, plonn, platn, output_srfc, output_c(1), output_c(7), output_c(9), output_c(11), output_c(14)
write(32,*) ${time_step}, pn, plonn, platn, output_srfc, output_c(1), output_c(7), output_c(9), output_c(11), output_c(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 ${ofname2}.txt ${ofname2}.dat
change_axcel5
#! /bin/tcsh -f
#####################################
# cwc change
# chage of axcel
#####################################
#####################################
filename1=change
filename2=read_data5
filename3=output_data
ffilename1=psea
ffilename2=cwc #[g/kg]
ofname1=output_axcel
ofname2=output_axcel_cwc
####################################
#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'
*psea
'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 x 1 418'
'set y 1 418'
k=1
while (k<=${zgrid})
'set z 'k''
'd cwc*1000'
*say 'yoshiyoshi'
*say k
k=k+1
endwhile
'disable fwrite'
'c'
'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'
write(*,*) 'time_step pn plonn platn'
write(32,*) '${ofilename2}'
write(32,*) 'time_step pn plonn platn cwc1000 cwc800 cwc600 cwc400 cwc200 sumcwc1000 sumcwc800 sumcwc600 sumcwc400 sumcwc200 allcwc'
write(32,*) 'time_step pn plonn platn cwc1000 cwc800 cwc600 cwc400 cwc200 sumcwc1000 sumcwc800 sumcwc600 sumcwc400 sumcwc200 allcwc'
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
!filename1=change
!filename2=read_data
!filename3=output_data
!ffilename1=psea
!ffilename2=cwc
!ofname1=output_axcel
!ofname2=output_axcel_cwc
real(kind=4), dimension(xgrid,ygrid) :: ${ffilename1}
real(kind=4), dimension(xgrid,ygrid,zgrid) :: ${ffilename2}
real(kind=4), dimension(zgrid) :: output_cwc
real(kind=4), dimension(zgrid) :: output_cwc_av
real(kind=4), dimension(zgrid) :: output_cwc_zintegral
real(kind=4), dimension(zgrid) :: output_cwc_av_zintegral
real(kind=4) :: pn, ppn, lonn, latn, plonn, platn, v1n, v2n, vlonn, vlatn, vmaxn, x
real(kind=4) :: p, axcel_avcwc, l_axcel
real(kind=4) :: cwc_sum, av_cwc, allcwc
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=ngrid3*4)
read (unit=12,rec=1) ${ffilename2}
open (unit=31, file='${ofname1}.txt', position='append', form='FORMATTED')
open (unit=32, file='${ofname2}.txt', position='append', form='FORMATTED')
!!!!!!!check cwc
!do k = 1, ${zgrid}
!write(*,*) 'cwc(100,300,k) ', cwc(100,300,k)
!enddo
!!!!!!!check cwc
!do k = 1, ${zgrid}
!do i = 1, ${xgrid}
!do j = 1, ${ygrid}
! if (cwc(i,j,k) > 0) then
! write(*,*) 'i, j, k, cwc(i,j,k) ', i, j, k, cwc(i,j,k)
! endif
!enddo
!enddo
!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}
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!cwc 1~zgrid
step=0
cwc_sum=0
if (${time_step} == 0) then
do i = 1, ${xgrid}
do j = 1, ${ygrid}
do k = 1, ${zgrid}
cwc(i,j,k)=0
enddo
enddo
enddo
endif
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
cwc_sum = cwc_sum + cwc(i,j,k)
step = step + 1
! write(*,*) 'cwc_sum cwc(i,j,k) step i j k', cwc_sum, cwc(i,j,k), step, i, j, k
! write(*,*) pkxn, pkyn, 'step', 'cwc_sum', step, cwc_sum, l_axcel, i, j
endif
l_axcel=0
enddo
enddo
output_cwc(k) = cwc_sum
output_cwc_av(k) = cwc_sum / step
output_cwc_zintegral(k) = output_cwc_zintegral(k) + output_cwc(k)
output_cwc_av_zintegral(k) = output_cwc_av_zintegral(k) + output_cwc_av(k)
! write(*,*) 'k output_cwc(k) output_cwc_av(k)', k, output_cwc(k), output_cwc_av(k)
step=0
cwc_sum=0
enddo
allcwc=0
do k =1 ,${zgrid}
allcwc = allcwc + output_cwc_zintegral(k)
!write(*,*) 'k, allcwc, output_cwc_zintegral(k)', k, allcwc, output_cwc_zintegral(k)
enddo
do k =1 ,${zgrid}
!write(*,*) 'k output_cwc(k) output_cwc_av(k), output_cwc_zintegral, output_cwc_av_zintegral', k, output_cwc(k), output_cwc_av(k), output_cwc_zintegral(k), output_cwc_zintegral(k)
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! output data
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(31,*) ${time_step}, pn, plonn, platn
write(*,*) ${time_step}, pn, plonn, platn
!write(32,*) ${time_step}, pn, plonn, platn, avcwc1000, avcwc800, avcwc600, avcwc400, avcwc200, sumcwc1000, sumcwc800, sumcwc600, sumcwc400, sumcwc200, allcwc
write(32,*) ${time_step}, pn, plonn, platn, output_cwc(1), output_cwc(7), output_cwc(9), output_cwc(11), output_cwc(14), output_cwc_zintegral(1), output_cwc_zintegral(7), output_cwc_zintegral(9), output_cwc_zintegral(11), output_cwc_zintegral(14), allcwc
write(*,*) ${time_step}, pn, plonn, platn, output_cwc(1), output_cwc(7), output_cwc(9), output_cwc(11), output_cwc(14), output_cwc_zintegral(1), output_cwc_zintegral(7), output_cwc_zintegral(9), output_cwc_zintegral(11), output_cwc_zintegral(14), allcwc
close(11)
close(12)
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`
our=`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 ${ofname2}.txt ${ofname2}.dat
cp ${ofname3}.txt ${ofname3}.dat
#####################################
# cwc change
# chage of axcel
#####################################
#####################################
filename1=change
filename2=read_data5
filename3=output_data
ffilename1=psea
ffilename2=cwc #[g/kg]
ofname1=output_axcel
ofname2=output_axcel_cwc
####################################
#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'
*psea
'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 x 1 418'
'set y 1 418'
k=1
while (k<=${zgrid})
'set z 'k''
'd cwc*1000'
*say 'yoshiyoshi'
*say k
k=k+1
endwhile
'disable fwrite'
'c'
'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'
write(*,*) 'time_step pn plonn platn'
write(32,*) '${ofilename2}'
write(32,*) 'time_step pn plonn platn cwc1000 cwc800 cwc600 cwc400 cwc200 sumcwc1000 sumcwc800 sumcwc600 sumcwc400 sumcwc200 allcwc'
write(32,*) 'time_step pn plonn platn cwc1000 cwc800 cwc600 cwc400 cwc200 sumcwc1000 sumcwc800 sumcwc600 sumcwc400 sumcwc200 allcwc'
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
!filename1=change
!filename2=read_data
!filename3=output_data
!ffilename1=psea
!ffilename2=cwc
!ofname1=output_axcel
!ofname2=output_axcel_cwc
real(kind=4), dimension(xgrid,ygrid) :: ${ffilename1}
real(kind=4), dimension(xgrid,ygrid,zgrid) :: ${ffilename2}
real(kind=4), dimension(zgrid) :: output_cwc
real(kind=4), dimension(zgrid) :: output_cwc_av
real(kind=4), dimension(zgrid) :: output_cwc_zintegral
real(kind=4), dimension(zgrid) :: output_cwc_av_zintegral
real(kind=4) :: pn, ppn, lonn, latn, plonn, platn, v1n, v2n, vlonn, vlatn, vmaxn, x
real(kind=4) :: p, axcel_avcwc, l_axcel
real(kind=4) :: cwc_sum, av_cwc, allcwc
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=ngrid3*4)
read (unit=12,rec=1) ${ffilename2}
open (unit=31, file='${ofname1}.txt', position='append', form='FORMATTED')
open (unit=32, file='${ofname2}.txt', position='append', form='FORMATTED')
!!!!!!!check cwc
!do k = 1, ${zgrid}
!write(*,*) 'cwc(100,300,k) ', cwc(100,300,k)
!enddo
!!!!!!!check cwc
!do k = 1, ${zgrid}
!do i = 1, ${xgrid}
!do j = 1, ${ygrid}
! if (cwc(i,j,k) > 0) then
! write(*,*) 'i, j, k, cwc(i,j,k) ', i, j, k, cwc(i,j,k)
! endif
!enddo
!enddo
!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}
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!cwc 1~zgrid
step=0
cwc_sum=0
if (${time_step} == 0) then
do i = 1, ${xgrid}
do j = 1, ${ygrid}
do k = 1, ${zgrid}
cwc(i,j,k)=0
enddo
enddo
enddo
endif
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
cwc_sum = cwc_sum + cwc(i,j,k)
step = step + 1
! write(*,*) 'cwc_sum cwc(i,j,k) step i j k', cwc_sum, cwc(i,j,k), step, i, j, k
! write(*,*) pkxn, pkyn, 'step', 'cwc_sum', step, cwc_sum, l_axcel, i, j
endif
l_axcel=0
enddo
enddo
output_cwc(k) = cwc_sum
output_cwc_av(k) = cwc_sum / step
output_cwc_zintegral(k) = output_cwc_zintegral(k) + output_cwc(k)
output_cwc_av_zintegral(k) = output_cwc_av_zintegral(k) + output_cwc_av(k)
! write(*,*) 'k output_cwc(k) output_cwc_av(k)', k, output_cwc(k), output_cwc_av(k)
step=0
cwc_sum=0
enddo
allcwc=0
do k =1 ,${zgrid}
allcwc = allcwc + output_cwc_zintegral(k)
!write(*,*) 'k, allcwc, output_cwc_zintegral(k)', k, allcwc, output_cwc_zintegral(k)
enddo
do k =1 ,${zgrid}
!write(*,*) 'k output_cwc(k) output_cwc_av(k), output_cwc_zintegral, output_cwc_av_zintegral', k, output_cwc(k), output_cwc_av(k), output_cwc_zintegral(k), output_cwc_zintegral(k)
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! output data
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(31,*) ${time_step}, pn, plonn, platn
write(*,*) ${time_step}, pn, plonn, platn
!write(32,*) ${time_step}, pn, plonn, platn, avcwc1000, avcwc800, avcwc600, avcwc400, avcwc200, sumcwc1000, sumcwc800, sumcwc600, sumcwc400, sumcwc200, allcwc
write(32,*) ${time_step}, pn, plonn, platn, output_cwc(1), output_cwc(7), output_cwc(9), output_cwc(11), output_cwc(14), output_cwc_zintegral(1), output_cwc_zintegral(7), output_cwc_zintegral(9), output_cwc_zintegral(11), output_cwc_zintegral(14), allcwc
write(*,*) ${time_step}, pn, plonn, platn, output_cwc(1), output_cwc(7), output_cwc(9), output_cwc(11), output_cwc(14), output_cwc_zintegral(1), output_cwc_zintegral(7), output_cwc_zintegral(9), output_cwc_zintegral(11), output_cwc_zintegral(14), allcwc
close(11)
close(12)
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`
our=`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 ${ofname2}.txt ${ofname2}.dat
cp ${ofname3}.txt ${ofname3}.dat
change_axcel4
! /bin/tcsh -f
#####################################
# cloud change
# chage of axcel
# rain before diff hour ~now time2
#####################################
#####################################
filename1=change
filename2=read_data4
filename3=output_data
ffilename1=psea
ffilename2=rain1
ffilename3=rain2
ofname1=output_axcel
ofname2=output_axcel_rain
####################################
####################################
#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
#: ${start1:=199201070000}
: ${start1:=199201041200}
: ${end1:=199201082100}
#: ${end1:=199201050000}
year1=`echo ${start1} | cut -c 1-4`
month1=`echo ${start1} | cut -c 5-6`
day1=`echo ${start1} | cut -c 7-8`
hour1=`echo ${start1} |cut -c 9-10`
h1=00
echo start time ${year1} ${month1} ${day1} ${hour1} ${h1}
echo end time ${end1}
: ${time1:=${start1}}
: ${start1:=${start1}}
#timestep
timestep=1
time_step=0
diff=1
#now time2
year2=`echo ${start1} | cut -c 1-4`
month2=`echo ${start1} | cut -c 5-6`
day2=`echo ${start1} | cut -c 7-8`
hour2=`echo ${start1} |cut -c 9-10`
h2=00
#### next time step
#hour=`expr $hour + 1`
hour2=`expr ${hour1} + ${diff}`
hour2=`echo ${hour2} | awk '{printf("%02d",$1)}'`
if [ $hour2 -ge 24 ]; then
hour2=`expr ${hour2} - 24`
# day=`expr $day + 1`
day2=`echo ${day2} | awk '{printf("%02d",$1+1)}'`
if [ $day2 -ge 32 ];then
day2=01
# month=`expr $month + 1`
month2=`echo ${month2} | awk '{printf("%02d",$1+1)}'`
if [ $month2 -ge 13 ];then
month2=01
year2=`expr $year2 + 1`
fi
fi
fi
time2=${year2}${month2}${day2}${hour2}${h2}
echo start before time1 ${year1} ${month1} ${day1} ${hour1} ${h1}
echo start now time2 ${year2} ${month2} ${day2} ${hour2} ${h2}
#########################
#########################
#########################
#########################
i=0
while [ ${time1} -le ${end1} ]
do
time_step=`expr ${timestep} \* ${i}`
echo time_step ${time_step}
echo start before time1 ${year1} ${month1} ${day1} ${hour1} ${h1}
echo start now time2 ${year2} ${month2} ${day2} ${hour2} ${h2}
##decide m
if [ ${month2} -eq 01 ]; then
m2=JAN
elif [ ${month2} -eq 02 ]; then
m2=FEB
elif [ ${month2} -eq 03 ]; then
m2=MAR
elif [ ${month2} -eq 04 ]; then
m2=APR
elif [ ${month2} -eq 05 ]; then
m2=MAY
elif [ ${month2} -eq 06 ]; then
m2=JUN
elif [ ${month2} -eq 07 ]; then
m2=JUL
elif [ ${month2} -eq 08 ]; then
m2=AUG
elif [ ${month2} -eq 09 ]; then
m2=SEP
elif [ ${month2} -eq 10 ]; then
m2=OCT
elif [ ${month2} -eq 11 ]; then
m2=NOV
elif [ ${month2} -eq 12 ]; then
m2=DEM
fi
##############################################
#read grads data
: ${prefix:=anl_p} #want data
echo ${time1}
echo ${time2}
### time2 rain2 now time
cat > ${filename2}.gs <<EOF
'reinit'
'open _NHMMRPPFCSVSTD1_${time2}.ctl'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename1}.grd'
*'set t ${t}'
'set x 1 ${xgrid}'
'set y 1 ${ygrid}'
!'set x 1 418'
!'set y 1 418'
'd psea'
'disable fwrite'
'c'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename3}.grd'
'set x 1 ${xgrid}'
'set y 1 ${ygrid}'
'd rain'
'disable fwrite'
'c'
'quit'
EOF
grads -blc ${filename2}.gs
### time1 rain1
cat > ${filename2}.gs <<EOF
'reinit'
'open _NHMMRPPFCSVSTD1_${time1}.ctl'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename2}.grd'
'set x 1 ${xgrid}'
'set y 1 ${ygrid}'
'd rain'
'disable fwrite'
'c'
'quit'
EOF
grads -blc ${filename2}.gs
########################################
####program read_data###################
########################################
if [ ${time1} -eq ${start1} ]; 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'
write(*,*) 'time_step pn plonn platn'
write(32,*) '${ofilename2}'
write(32,*) 'time_step pn plonn platn area_rain av_rain'
write(*,*) 'time_step pn plonn platn area_rain av_rain'
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
! integer, parameter :: time1=${time1}, time2=${time2},
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
!filename1=change_cloud
!filename2=read_data
!filename3=output_data
!ffilename1=psea
!ffilename2=rain1
!ffilename3=rain2
!ofname1=output_axcel
!ofname2=output_axcel_rain
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) :: rain_d
real(kind=4) :: pn, ppn, lonn, latn, plonn, platn, v1n, v2n, vlonn, vlatn, vmaxn, x
real(kind=4) :: p, axcel_avcwc, l_axcel
real(kind=4) :: rain_sum, av_rain, area_rain
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=31, file='${ofname1}.txt', position='append', form='FORMATTED')
open (unit=32, file='${ofname2}.txt', position='append', form='FORMATTED')
!!!!!!!check rain1,2, d
!do i = 100, 300
! do j = 100, 300
!if (rain1(i,j) < rain2(i,j)) then
!write(*,*) 'rain1(i,j) rain2(i,j)', rain1(i,j), rain2(i,j)
!endif
!enddo
!enddo
!do i = 200, ${xgrid}
!do j = 200, ${ygrid}
!psea(i,j)=psea(i,j)/100
!write(*,*) 'p ', psea(150,300)
!enddo
!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
!write(*,*) 'yoshiyoshi', step
!write(*,*) 'kxn, kyn', kxn, kyn
!step=step+1
endif
endif
n=n+1
enddo
enddo
pkxn=kxn
pkyn=kyn
!pmin location
plonn=kxn*${gsize2}${f_lon}
platn=kyn*${gsize2}${f_lat}
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!cal rain
step=0
rain_sum=0
if (${time_step} == 0) then
do i = 1, ${xgrid}
do j = 1, ${ygrid}
rain1(i,j)=0
rain2(i,j)=0
enddo
enddo
endif
!ck rain>0
do i = 1, ${xgrid}
do j = 1, ${ygrid}
if (rain1(i,j) < 0) then
rain1(i,j)=0
endif
enddo
enddo
do i = 1, ${xgrid}
do j = 1, ${ygrid}
if (rain2(i,j) < 0) then
rain2(i,j)=0
endif
enddo
enddo
do i = 1, ${xgrid}
do j = 1, ${ygrid}
rain_d(i,j)=rain2(i,j)-rain1(i,j)
!if(rain_d(i,j) > 0) then
!write(*,*) 'rain_d(i,j) rain2(i,j) rain1(i,j)', rain_d(i,j), rain2(i,j), rain1(i,j)
!endif
enddo
enddo
do i = 1, ${xgrid}
do j = 1, ${ygrid}
if (rain_d(i,j) < 0) then
rain_d(i,j)=0
endif
enddo
enddo
rain_sum=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
rain_sum = rain_sum + rain_d(i,j)
step = step + 1
!write(*,*) 'rain_sum rain_d(i,j) step i j', rain_sum, rain_d(i,j), step, i, j
!write(*,*) pkxn, pkyn, 'step', 'rain_sum', step, rain_sum, l_axcel, i, j
endif
l_axcel=0
enddo
enddo
area_rain = rain_sum
av_rain = rain_sum / step
write(*,*) 'area_rain av_rain step', area_rain, av_rain, step
! real(kind=4) :: rain_sum, av_rain, area_rain
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! output data
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(31,*) ${time_step}, pn, plonn, platn
write(*,*) ${time_step}, pn, plonn, platn
write(32,*) ${time_step}, pn, plonn, platn, area_rain, av_rain
write(*,*) ${time_step}, pn, plonn, platn, area_rain, av_rain
!write(32,*) 'start before time1', ${time1}
!write(32,*) 'start now time2', ${time2}
close(11)
close(12)
close(13)
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
#### next time1 step
#hour=`expr $hour + 1`
hour1=`expr $hour1 + $timestep`
hour11=$hour1
hour1=`echo ${hour1} | awk '{printf("%02d",$1)}'`
if [ $hour1 -ge 24 ]; then
hour1=00
# day=`expr $day + 1`
day1=`echo ${day1} | awk '{printf("%02d",$1+1)}'`
if [ $day1 -ge 32 ];then
day1=01
# month=`expr $month + 1`
month1=`echo ${month1} | awk '{printf("%02d",$1+1)}'`
if [ $month1 -ge 13 ];then
month1=01
year1=`expr $year1 + 1`
fi
fi
fi
time1=${year1}${month1}${day1}${hour1}${h1}
#### next time2 step
#hour=`expr $hour + 1`
#"hour2=`expr $hour1 + $diff`
#day2=`expr $day1 + $diff`
hour2=`expr $hour2 + $timestep`
echo hour1 $hour1
echo hour2 $hour2
hour2=`echo ${hour2} | awk '{printf("%02d",$1)}'`
if [ $hour2 -ge 24 ]; then
hour2=`expr ${hour2} - 24 | awk '{printf("%02d",$1)}'`
echo hour2 $hour2
# day=`expr $day + 1`
day2=`echo ${day2} | awk '{printf("%02d",$1+1)}'`
if [ $day2 -ge 32 ];then
day2=01
# month=`expr $month + 1`
month2=`echo ${month2} | awk '{printf("%02d",$1+1)}'`
if [ $month2 -ge 13 ];then
month2=01
year2=`expr $year2 + 1`
fi
fi
fi
time2=${year2}${month2}${day2}${hour2}${h2}
done
cp ${ofname1}.txt ${ofname1}.dat
cp ${ofname2}.txt ${ofname2}.dat
#####################################
# cloud change
# chage of axcel
# rain before diff hour ~now time2
#####################################
#####################################
filename1=change
filename2=read_data4
filename3=output_data
ffilename1=psea
ffilename2=rain1
ffilename3=rain2
ofname1=output_axcel
ofname2=output_axcel_rain
####################################
####################################
#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
#: ${start1:=199201070000}
: ${start1:=199201041200}
: ${end1:=199201082100}
#: ${end1:=199201050000}
year1=`echo ${start1} | cut -c 1-4`
month1=`echo ${start1} | cut -c 5-6`
day1=`echo ${start1} | cut -c 7-8`
hour1=`echo ${start1} |cut -c 9-10`
h1=00
echo start time ${year1} ${month1} ${day1} ${hour1} ${h1}
echo end time ${end1}
: ${time1:=${start1}}
: ${start1:=${start1}}
#timestep
timestep=1
time_step=0
diff=1
#now time2
year2=`echo ${start1} | cut -c 1-4`
month2=`echo ${start1} | cut -c 5-6`
day2=`echo ${start1} | cut -c 7-8`
hour2=`echo ${start1} |cut -c 9-10`
h2=00
#### next time step
#hour=`expr $hour + 1`
hour2=`expr ${hour1} + ${diff}`
hour2=`echo ${hour2} | awk '{printf("%02d",$1)}'`
if [ $hour2 -ge 24 ]; then
hour2=`expr ${hour2} - 24`
# day=`expr $day + 1`
day2=`echo ${day2} | awk '{printf("%02d",$1+1)}'`
if [ $day2 -ge 32 ];then
day2=01
# month=`expr $month + 1`
month2=`echo ${month2} | awk '{printf("%02d",$1+1)}'`
if [ $month2 -ge 13 ];then
month2=01
year2=`expr $year2 + 1`
fi
fi
fi
time2=${year2}${month2}${day2}${hour2}${h2}
echo start before time1 ${year1} ${month1} ${day1} ${hour1} ${h1}
echo start now time2 ${year2} ${month2} ${day2} ${hour2} ${h2}
#########################
#########################
#########################
#########################
i=0
while [ ${time1} -le ${end1} ]
do
time_step=`expr ${timestep} \* ${i}`
echo time_step ${time_step}
echo start before time1 ${year1} ${month1} ${day1} ${hour1} ${h1}
echo start now time2 ${year2} ${month2} ${day2} ${hour2} ${h2}
##decide m
if [ ${month2} -eq 01 ]; then
m2=JAN
elif [ ${month2} -eq 02 ]; then
m2=FEB
elif [ ${month2} -eq 03 ]; then
m2=MAR
elif [ ${month2} -eq 04 ]; then
m2=APR
elif [ ${month2} -eq 05 ]; then
m2=MAY
elif [ ${month2} -eq 06 ]; then
m2=JUN
elif [ ${month2} -eq 07 ]; then
m2=JUL
elif [ ${month2} -eq 08 ]; then
m2=AUG
elif [ ${month2} -eq 09 ]; then
m2=SEP
elif [ ${month2} -eq 10 ]; then
m2=OCT
elif [ ${month2} -eq 11 ]; then
m2=NOV
elif [ ${month2} -eq 12 ]; then
m2=DEM
fi
##############################################
#read grads data
: ${prefix:=anl_p} #want data
echo ${time1}
echo ${time2}
### time2 rain2 now time
cat > ${filename2}.gs <<EOF
'reinit'
'open _NHMMRPPFCSVSTD1_${time2}.ctl'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename1}.grd'
*'set t ${t}'
'set x 1 ${xgrid}'
'set y 1 ${ygrid}'
!'set x 1 418'
!'set y 1 418'
'd psea'
'disable fwrite'
'c'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename3}.grd'
'set x 1 ${xgrid}'
'set y 1 ${ygrid}'
'd rain'
'disable fwrite'
'c'
'quit'
EOF
grads -blc ${filename2}.gs
### time1 rain1
cat > ${filename2}.gs <<EOF
'reinit'
'open _NHMMRPPFCSVSTD1_${time1}.ctl'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename2}.grd'
'set x 1 ${xgrid}'
'set y 1 ${ygrid}'
'd rain'
'disable fwrite'
'c'
'quit'
EOF
grads -blc ${filename2}.gs
########################################
####program read_data###################
########################################
if [ ${time1} -eq ${start1} ]; 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'
write(*,*) 'time_step pn plonn platn'
write(32,*) '${ofilename2}'
write(32,*) 'time_step pn plonn platn area_rain av_rain'
write(*,*) 'time_step pn plonn platn area_rain av_rain'
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
! integer, parameter :: time1=${time1}, time2=${time2},
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
!filename1=change_cloud
!filename2=read_data
!filename3=output_data
!ffilename1=psea
!ffilename2=rain1
!ffilename3=rain2
!ofname1=output_axcel
!ofname2=output_axcel_rain
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) :: rain_d
real(kind=4) :: pn, ppn, lonn, latn, plonn, platn, v1n, v2n, vlonn, vlatn, vmaxn, x
real(kind=4) :: p, axcel_avcwc, l_axcel
real(kind=4) :: rain_sum, av_rain, area_rain
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=31, file='${ofname1}.txt', position='append', form='FORMATTED')
open (unit=32, file='${ofname2}.txt', position='append', form='FORMATTED')
!!!!!!!check rain1,2, d
!do i = 100, 300
! do j = 100, 300
!if (rain1(i,j) < rain2(i,j)) then
!write(*,*) 'rain1(i,j) rain2(i,j)', rain1(i,j), rain2(i,j)
!endif
!enddo
!enddo
!do i = 200, ${xgrid}
!do j = 200, ${ygrid}
!psea(i,j)=psea(i,j)/100
!write(*,*) 'p ', psea(150,300)
!enddo
!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
!write(*,*) 'yoshiyoshi', step
!write(*,*) 'kxn, kyn', kxn, kyn
!step=step+1
endif
endif
n=n+1
enddo
enddo
pkxn=kxn
pkyn=kyn
!pmin location
plonn=kxn*${gsize2}${f_lon}
platn=kyn*${gsize2}${f_lat}
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!cal rain
step=0
rain_sum=0
if (${time_step} == 0) then
do i = 1, ${xgrid}
do j = 1, ${ygrid}
rain1(i,j)=0
rain2(i,j)=0
enddo
enddo
endif
!ck rain>0
do i = 1, ${xgrid}
do j = 1, ${ygrid}
if (rain1(i,j) < 0) then
rain1(i,j)=0
endif
enddo
enddo
do i = 1, ${xgrid}
do j = 1, ${ygrid}
if (rain2(i,j) < 0) then
rain2(i,j)=0
endif
enddo
enddo
do i = 1, ${xgrid}
do j = 1, ${ygrid}
rain_d(i,j)=rain2(i,j)-rain1(i,j)
!if(rain_d(i,j) > 0) then
!write(*,*) 'rain_d(i,j) rain2(i,j) rain1(i,j)', rain_d(i,j), rain2(i,j), rain1(i,j)
!endif
enddo
enddo
do i = 1, ${xgrid}
do j = 1, ${ygrid}
if (rain_d(i,j) < 0) then
rain_d(i,j)=0
endif
enddo
enddo
rain_sum=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
rain_sum = rain_sum + rain_d(i,j)
step = step + 1
!write(*,*) 'rain_sum rain_d(i,j) step i j', rain_sum, rain_d(i,j), step, i, j
!write(*,*) pkxn, pkyn, 'step', 'rain_sum', step, rain_sum, l_axcel, i, j
endif
l_axcel=0
enddo
enddo
area_rain = rain_sum
av_rain = rain_sum / step
write(*,*) 'area_rain av_rain step', area_rain, av_rain, step
! real(kind=4) :: rain_sum, av_rain, area_rain
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! output data
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(31,*) ${time_step}, pn, plonn, platn
write(*,*) ${time_step}, pn, plonn, platn
write(32,*) ${time_step}, pn, plonn, platn, area_rain, av_rain
write(*,*) ${time_step}, pn, plonn, platn, area_rain, av_rain
!write(32,*) 'start before time1', ${time1}
!write(32,*) 'start now time2', ${time2}
close(11)
close(12)
close(13)
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
#### next time1 step
#hour=`expr $hour + 1`
hour1=`expr $hour1 + $timestep`
hour11=$hour1
hour1=`echo ${hour1} | awk '{printf("%02d",$1)}'`
if [ $hour1 -ge 24 ]; then
hour1=00
# day=`expr $day + 1`
day1=`echo ${day1} | awk '{printf("%02d",$1+1)}'`
if [ $day1 -ge 32 ];then
day1=01
# month=`expr $month + 1`
month1=`echo ${month1} | awk '{printf("%02d",$1+1)}'`
if [ $month1 -ge 13 ];then
month1=01
year1=`expr $year1 + 1`
fi
fi
fi
time1=${year1}${month1}${day1}${hour1}${h1}
#### next time2 step
#hour=`expr $hour + 1`
#"hour2=`expr $hour1 + $diff`
#day2=`expr $day1 + $diff`
hour2=`expr $hour2 + $timestep`
echo hour1 $hour1
echo hour2 $hour2
hour2=`echo ${hour2} | awk '{printf("%02d",$1)}'`
if [ $hour2 -ge 24 ]; then
hour2=`expr ${hour2} - 24 | awk '{printf("%02d",$1)}'`
echo hour2 $hour2
# day=`expr $day + 1`
day2=`echo ${day2} | awk '{printf("%02d",$1+1)}'`
if [ $day2 -ge 32 ];then
day2=01
# month=`expr $month + 1`
month2=`echo ${month2} | awk '{printf("%02d",$1+1)}'`
if [ $month2 -ge 13 ];then
month2=01
year2=`expr $year2 + 1`
fi
fi
fi
time2=${year2}${month2}${day2}${hour2}${h2}
done
cp ${ofname1}.txt ${ofname1}.dat
cp ${ofname2}.txt ${ofname2}.dat
change_axcel3.sh
#! /bin/tcsh -f
#####################################
# qv change
# chage of axcel
#####################################
#####################################
filename1=change
filename2=read_data3
filename3=output_data
ffilename1=psea
ffilename2=qv
ofname1=output_axcel
ofname2=output_axcel_qv
####################################
#model grid size
xgrid=418
ygrid=418
zgrid=24
#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 ${xgrid}'
'set y 1 ${ygrid}'
!'set x 1 418'
!'set y 1 418'
'd psea'
'disable fwrite'
'c'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename2}.grd'
'set x 1 418'
'set y 1 418'
k=1
while (k<=${zgrid})
'set z 'k''
'd qv*1000'
*say 'yoshiyoshi'
*say k
k=k+1
endwhile
'disable fwrite'
'c'
'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'
write(*,*) 'time_step pn plonn platn'
write(32,*) '${ofilename2}'
write(32,*) 'time_step pn plonn platn qv20m qv1180m qv5330m qv7830m'
write(*,*) 'time_step pn plonn platn qv20m qv1180m qv5330m qv7830m'
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
!filename1=change
!filename2=read_data
!filename3=output_data
!ffilename1=psea
!ffilename2=qv
!ofname1=output_axcel
!ofname2=output_axcel_qv
real(kind=4), dimension(xgrid,ygrid) :: ${ffilename1}
real(kind=4), dimension(xgrid,ygrid,zgrid) :: ${ffilename2}
real(kind=4), dimension(zgrid) :: output_qv
real(kind=4) :: pn, ppn, lonn, latn, plonn, platn, v1n, v2n, vlonn, vlatn, vmaxn, x
real(kind=4) :: p, axcel_avcwc, l_axcel
real(kind=4) :: qv_sum, avqv
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=ngrid3*4)
read (unit=12,rec=1) ${ffilename2}
open (unit=31, file='${ofname1}.txt', position='append', form='FORMATTED')
open (unit=32, file='${ofname2}.txt', position='append', form='FORMATTED')
!!!!!!!check qv
!do k = 1, ${zgrid}
!write(*,*) 'qv(100,100,k) ', qv(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
!write(*,*) 'yoshiyoshi', psea(i,j)
!write(*,*) 'kxn, kyn', kxn, kyn
!step=step+1
endif
endif
n=n+1
enddo
enddo
pkxn=kxn
pkyn=kyn
!pmin location
plonn=kxn*${gsize2}${f_lon}
platn=kyn*${gsize2}${f_lat}
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
step=0
qv_sum=0
if (${time_step} == 0) then
do i = 1, ${xgrid}
do j = 1, ${ygrid}
do k = 1, ${zgrid}
qv(i,j,k)=0
enddo
enddo
enddo
endif
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
if(qv(i,j,k)<100) then
qv_sum = qv_sum + qv(i,j,k)
step = step + 1
!write(*,*) 'qv_sum qv(i,j,k) step i j k', qv_sum, qv(i,j,k), step, i, j, k
!write(*,*) pkxn, pkyn, 'step', 'qv_sum', step, qv_sum, l_axcel, i, j
endif
endif
l_axcel=0
enddo
enddo
output_qv(k) = qv_sum / step
!write(*,*) 'output_qv(k) k', output_qv(k), step, k
step=0
qv_sum=0
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! output data
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(31,*) ${time_step}, pn, plonn, platn
write(*,*) ${time_step}, pn, plonn, platn
write(32,*) ${time_step}, pn, plonn, platn, output_qv(2), output_qv(10), output_qv(20), output_qv(24)
write(*,*) ${time_step}, pn, plonn, platn, output_qv(2), output_qv(10), output_qv(20), output_qv(24)
close(11)
close(12)
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`
our=`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 ${ofname2}.txt ${ofname2}.dat
cp ${ofname3}.txt ${ofname3}.dat
#####################################
# qv change
# chage of axcel
#####################################
#####################################
filename1=change
filename2=read_data3
filename3=output_data
ffilename1=psea
ffilename2=qv
ofname1=output_axcel
ofname2=output_axcel_qv
####################################
#model grid size
xgrid=418
ygrid=418
zgrid=24
#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 ${xgrid}'
'set y 1 ${ygrid}'
!'set x 1 418'
!'set y 1 418'
'd psea'
'disable fwrite'
'c'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename2}.grd'
'set x 1 418'
'set y 1 418'
k=1
while (k<=${zgrid})
'set z 'k''
'd qv*1000'
*say 'yoshiyoshi'
*say k
k=k+1
endwhile
'disable fwrite'
'c'
'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'
write(*,*) 'time_step pn plonn platn'
write(32,*) '${ofilename2}'
write(32,*) 'time_step pn plonn platn qv20m qv1180m qv5330m qv7830m'
write(*,*) 'time_step pn plonn platn qv20m qv1180m qv5330m qv7830m'
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
!filename1=change
!filename2=read_data
!filename3=output_data
!ffilename1=psea
!ffilename2=qv
!ofname1=output_axcel
!ofname2=output_axcel_qv
real(kind=4), dimension(xgrid,ygrid) :: ${ffilename1}
real(kind=4), dimension(xgrid,ygrid,zgrid) :: ${ffilename2}
real(kind=4), dimension(zgrid) :: output_qv
real(kind=4) :: pn, ppn, lonn, latn, plonn, platn, v1n, v2n, vlonn, vlatn, vmaxn, x
real(kind=4) :: p, axcel_avcwc, l_axcel
real(kind=4) :: qv_sum, avqv
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=ngrid3*4)
read (unit=12,rec=1) ${ffilename2}
open (unit=31, file='${ofname1}.txt', position='append', form='FORMATTED')
open (unit=32, file='${ofname2}.txt', position='append', form='FORMATTED')
!!!!!!!check qv
!do k = 1, ${zgrid}
!write(*,*) 'qv(100,100,k) ', qv(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
!write(*,*) 'yoshiyoshi', psea(i,j)
!write(*,*) 'kxn, kyn', kxn, kyn
!step=step+1
endif
endif
n=n+1
enddo
enddo
pkxn=kxn
pkyn=kyn
!pmin location
plonn=kxn*${gsize2}${f_lon}
platn=kyn*${gsize2}${f_lat}
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
step=0
qv_sum=0
if (${time_step} == 0) then
do i = 1, ${xgrid}
do j = 1, ${ygrid}
do k = 1, ${zgrid}
qv(i,j,k)=0
enddo
enddo
enddo
endif
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
if(qv(i,j,k)<100) then
qv_sum = qv_sum + qv(i,j,k)
step = step + 1
!write(*,*) 'qv_sum qv(i,j,k) step i j k', qv_sum, qv(i,j,k), step, i, j, k
!write(*,*) pkxn, pkyn, 'step', 'qv_sum', step, qv_sum, l_axcel, i, j
endif
endif
l_axcel=0
enddo
enddo
output_qv(k) = qv_sum / step
!write(*,*) 'output_qv(k) k', output_qv(k), step, k
step=0
qv_sum=0
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! output data
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(31,*) ${time_step}, pn, plonn, platn
write(*,*) ${time_step}, pn, plonn, platn
write(32,*) ${time_step}, pn, plonn, platn, output_qv(2), output_qv(10), output_qv(20), output_qv(24)
write(*,*) ${time_step}, pn, plonn, platn, output_qv(2), output_qv(10), output_qv(20), output_qv(24)
close(11)
close(12)
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`
our=`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 ${ofname2}.txt ${ofname2}.dat
cp ${ofname3}.txt ${ofname3}.dat
change_axcel2.sh
#! /bin/tcsh -f
#####################################
# cloud change
# chage of axcel
#####################################
#####################################
filename1=change_cloud
filename2=read_data2
filename3=output_data
ffilename1=psea
ffilename2=v # 1000hPa
ffilename3=allcwc
ffilename4=rv #1e+4
ffilename5=rvsrf #1e+4
ffilename6=t #[K]
ffilename7=tsrf #[K]
ofname1=output_axcel
#ofname2=output_axcel_rv
ofname3=output_axcel_ept
####################################
####################################
#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'
*psea
'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'
*v (u*u+v*v)**0.5 1000hPa
'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'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename6}.grd'
'set x 1 418'
'set y 1 418'
k=1
while (k<=${zgrid})
'set z 'k''
'd t'
*say 'yoshiyoshi'
*say k
k=k+1
endwhile
'disable fwrite'
'c'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename7}.grd'
*'set t ${t}'
'set x 1 418'
'set y 1 418'
'd tsrf'
'disable fwrite'
'c'
'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')
open (unit=33, file='${ofname3}.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'
write(33,*) '${ofilename3}'
write(33,*) 'time_step pn plonn platn v1n aveptsrf output_ept1000 output_ept800 output_ept600 output_ept400'
write(*,*) 'time_step pn plonn platn v1n aveptsrf output_ept1000 output_ept800 output_ept600 output_ept400'
close(31)
!#close(32)
close(33)
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 # 1000hPa
!ffilename3=allcwc
!ffilename4=rv #1e+4
!ffilename5=rvsrf #1e+4
!ffilename6=t #[K]
!ffilename7=tsrf #[K]
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(xgrid,ygrid,zgrid) :: ${ffilename6}
real(kind=4), dimension(xgrid,ygrid) :: ${ffilename7}
real(kind=4), dimension(zgrid) :: output_rv
real(kind=4), dimension(zgrid) :: output_ept
real(kind=4), dimension(zgrid) :: cal_p
real(kind=4), dimension(xgrid,ygrid,zgrid) :: pt
real(kind=4), dimension(xgrid,ygrid,zgrid) :: es
real(kind=4), dimension(xgrid,ygrid,zgrid) :: ws
real(kind=4), dimension(xgrid,ygrid,zgrid) :: ept
real(kind=4), dimension(xgrid,ygrid) :: pts
real(kind=4), dimension(xgrid,ygrid) :: ess
real(kind=4), dimension(xgrid,ygrid) :: wss
real(kind=4), dimension(xgrid,ygrid) :: epts
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
real(kind=4) :: ept_sum, epts_sum, avept, aveptsrf, rde, cpe, le, c1, c2, c3, c4
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=16,file='${ffilename6}.grd',form='unformatted',access='direct',recl=ngrid3*4)
read (unit=16,rec=1) ${ffilename6}
open (unit=17,file='${ffilename7}.grd',form='unformatted',access='direct',recl=ngrid2*4)
read (unit=17,rec=1) ${ffilename7}
open (unit=31, file='${ofname1}.txt', position='append', form='FORMATTED')
!open (unit=32, file='${ofname2}.txt', position='append', form='FORMATTED')
open (unit=33, file='${ofname3}.txt', position='append', form='FORMATTED')
!!!!!!!check rv
!do k = 1, ${zgrid}
!write(*,*) 'rv(100,100,k) ', rv(100,100,k)
!enddo
!!!!!!!check t
!do k = 1, ${zgrid}
!write(*,*) 't(100,100,k) ', t(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}
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! cal ept
!ept 1~zgrid
!equipment
cal_p(1)=1000
cal_p(2)=975
cal_p(3)=950
cal_p(4)=925
cal_p(5)=900
cal_p(6)=850
cal_p(7)=800
cal_p(8)=700
cal_p(9)=600
cal_p(10)=500
cal_p(11)=400
cal_p(12)=300
cal_p(13)=250
cal_p(14)=200
cal_p(15)=150
cal_p(16)=100
rde=287
cpe=1004
le=2500000
c2=rde/cpe
c3=le/cpe
!ept_sum, avept, aveptsrf, rde, cpe, le, pt, c1, c2, c3, c4
!rde, cpe, le, pt, c1, c2, c3, c4
pt(i,j,k)=1
es(i,j,k)=1
ws(i,j,k)=1
ept(i,j,k)=1
!cal ept z 1~zgrid
do k = 1, ${zgrid}
do i = 1, ${xgrid}
do j = 1, ${ygrid}
c1=1000/cal_p(k)
pt(i,j,k)=c1**c2
pt(i,j,k)=pt(i,j,k)*t(i,j,k)
es(i,j,k)=17.27*(t(i,j,k)-273.15)/(t(i,j,k)-35.86) !#tetens
es(i,j,k)=exp(es(i,j,k))
ws(i,j,k)=0.622*es(i,j,k)/cal_p(k)
ept(i,j,k)=le*ws(i,j,k)/cpe
ept(i,j,k)=ept(i,j,k)/t(i,j,k)
ept(i,j,k)=pt(i,j,k)*exp(ept(i,j,k))
enddo
enddo
enddo
!!!!!!!check ept
!do k = 1, ${zgrid}
!write(*,*) 'ept(100,100,k) ', ept(100,100,k)
!enddo
!ept 1~zgrid
step=0
ept_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!
ept_sum = ept_sum + ept(i,j,k)
step = step + 1
!write(*,*) 'ept_sum ept(i,j,k) step i j k', ept_sum, ept(i,j,k), step, i, j, k
!write(*,*) pkxn, pkyn, 'step', 'ept_sum', step, ept_sum, l_axcel, i, j
endif
l_axcel=0
enddo
enddo
output_ept(k) = ept_sum / step
!write(*,*) 'output_ept(k) k', output_ept(k), step, k
step=0
ept_sum=0
enddo
! real(kind=4) :: ept_sum, avept, aveptsrf, rde, cpe, le, c1, c2, c3, c4, avept, aveptsrf
!ept srf
!cal ept srf use kousiki
!cal ept z 1~zgrid
do i = 1, ${xgrid}
do j = 1, ${ygrid}
c1=1000/1000
pts(i,j)=c1**c2
pts(i,j)=pts(i,j)*tsrf(i,j)
ess(i,j)=17.27*(tsrf(i,j)-273.15)/(tsrf(i,j)-35.86) !#tetens
ess(i,j)=exp(ess(i,j))
wss(i,j)=0.622*ess(i,j)/1000
epts(i,j)=le*wss(i,j)/cpe
epts(i,j)=epts(i,j)/tsrf(i,j)
epts(i,j)=pts(i,j)*exp(epts(i,j))
enddo
enddo
!do i = 150, 200
!write(*,*) 'pts, ess, wss epts(272,294)', pts(272,294), ess(272,294), wss(272,294), epts(272,294)
!enddo
step=0
epts_sum=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
if (epts(i,j)<1000) then
epts_sum = epts_sum + epts(i,j)
step = step + 1
!write(*,*) 'epts_sum epts(i,j) step i j', epts_sum, epts(i,j), step, i, j
!write(*,*) pkxn, pkyn, 'step', 'epts_sum', step, ept_sum, l_axcel, i, j
endif
endif
l_axcel=0
enddo
enddo
aveptsrf = epts_sum / step
!write(*,*) 'aveptsrf ', aveptsrf, step
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! 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)
!write(33,*) 'time_step pn plonn platn v1n aveptsrf output_ept1000 output_ept800 output_ept600 output_ept400'
write(33,*) ${time_step}, pn, plonn, platn, v1n, aveptsrf, output_ept(1), output_ept(7), output_ept(9), output_ept(11)
write(*,*) ${time_step}, pn, plonn, platn, v1n, aveptsrf, output_ept(1), output_ept(7), output_ept(9), output_ept(11)
close(11)
close(12)
close(13)
close(14)
close(15)
close(16)
close(17)
close(31)
!close(32)
close(33)
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 ${ofname2}.txt ${ofname2}.dat
cp ${ofname3}.txt ${ofname3}.dat
#####################################
# cloud change
# chage of axcel
#####################################
#####################################
filename1=change_cloud
filename2=read_data2
filename3=output_data
ffilename1=psea
ffilename2=v # 1000hPa
ffilename3=allcwc
ffilename4=rv #1e+4
ffilename5=rvsrf #1e+4
ffilename6=t #[K]
ffilename7=tsrf #[K]
ofname1=output_axcel
#ofname2=output_axcel_rv
ofname3=output_axcel_ept
####################################
####################################
#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'
*psea
'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'
*v (u*u+v*v)**0.5 1000hPa
'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'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename6}.grd'
'set x 1 418'
'set y 1 418'
k=1
while (k<=${zgrid})
'set z 'k''
'd t'
*say 'yoshiyoshi'
*say k
k=k+1
endwhile
'disable fwrite'
'c'
'set gxout fwrite'
'set undef dfile'
'set fwrite ${ffilename7}.grd'
*'set t ${t}'
'set x 1 418'
'set y 1 418'
'd tsrf'
'disable fwrite'
'c'
'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')
open (unit=33, file='${ofname3}.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'
write(33,*) '${ofilename3}'
write(33,*) 'time_step pn plonn platn v1n aveptsrf output_ept1000 output_ept800 output_ept600 output_ept400'
write(*,*) 'time_step pn plonn platn v1n aveptsrf output_ept1000 output_ept800 output_ept600 output_ept400'
close(31)
!#close(32)
close(33)
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 # 1000hPa
!ffilename3=allcwc
!ffilename4=rv #1e+4
!ffilename5=rvsrf #1e+4
!ffilename6=t #[K]
!ffilename7=tsrf #[K]
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(xgrid,ygrid,zgrid) :: ${ffilename6}
real(kind=4), dimension(xgrid,ygrid) :: ${ffilename7}
real(kind=4), dimension(zgrid) :: output_rv
real(kind=4), dimension(zgrid) :: output_ept
real(kind=4), dimension(zgrid) :: cal_p
real(kind=4), dimension(xgrid,ygrid,zgrid) :: pt
real(kind=4), dimension(xgrid,ygrid,zgrid) :: es
real(kind=4), dimension(xgrid,ygrid,zgrid) :: ws
real(kind=4), dimension(xgrid,ygrid,zgrid) :: ept
real(kind=4), dimension(xgrid,ygrid) :: pts
real(kind=4), dimension(xgrid,ygrid) :: ess
real(kind=4), dimension(xgrid,ygrid) :: wss
real(kind=4), dimension(xgrid,ygrid) :: epts
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
real(kind=4) :: ept_sum, epts_sum, avept, aveptsrf, rde, cpe, le, c1, c2, c3, c4
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=16,file='${ffilename6}.grd',form='unformatted',access='direct',recl=ngrid3*4)
read (unit=16,rec=1) ${ffilename6}
open (unit=17,file='${ffilename7}.grd',form='unformatted',access='direct',recl=ngrid2*4)
read (unit=17,rec=1) ${ffilename7}
open (unit=31, file='${ofname1}.txt', position='append', form='FORMATTED')
!open (unit=32, file='${ofname2}.txt', position='append', form='FORMATTED')
open (unit=33, file='${ofname3}.txt', position='append', form='FORMATTED')
!!!!!!!check rv
!do k = 1, ${zgrid}
!write(*,*) 'rv(100,100,k) ', rv(100,100,k)
!enddo
!!!!!!!check t
!do k = 1, ${zgrid}
!write(*,*) 't(100,100,k) ', t(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}
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! cal ept
!ept 1~zgrid
!equipment
cal_p(1)=1000
cal_p(2)=975
cal_p(3)=950
cal_p(4)=925
cal_p(5)=900
cal_p(6)=850
cal_p(7)=800
cal_p(8)=700
cal_p(9)=600
cal_p(10)=500
cal_p(11)=400
cal_p(12)=300
cal_p(13)=250
cal_p(14)=200
cal_p(15)=150
cal_p(16)=100
rde=287
cpe=1004
le=2500000
c2=rde/cpe
c3=le/cpe
!ept_sum, avept, aveptsrf, rde, cpe, le, pt, c1, c2, c3, c4
!rde, cpe, le, pt, c1, c2, c3, c4
pt(i,j,k)=1
es(i,j,k)=1
ws(i,j,k)=1
ept(i,j,k)=1
!cal ept z 1~zgrid
do k = 1, ${zgrid}
do i = 1, ${xgrid}
do j = 1, ${ygrid}
c1=1000/cal_p(k)
pt(i,j,k)=c1**c2
pt(i,j,k)=pt(i,j,k)*t(i,j,k)
es(i,j,k)=17.27*(t(i,j,k)-273.15)/(t(i,j,k)-35.86) !#tetens
es(i,j,k)=exp(es(i,j,k))
ws(i,j,k)=0.622*es(i,j,k)/cal_p(k)
ept(i,j,k)=le*ws(i,j,k)/cpe
ept(i,j,k)=ept(i,j,k)/t(i,j,k)
ept(i,j,k)=pt(i,j,k)*exp(ept(i,j,k))
enddo
enddo
enddo
!!!!!!!check ept
!do k = 1, ${zgrid}
!write(*,*) 'ept(100,100,k) ', ept(100,100,k)
!enddo
!ept 1~zgrid
step=0
ept_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!
ept_sum = ept_sum + ept(i,j,k)
step = step + 1
!write(*,*) 'ept_sum ept(i,j,k) step i j k', ept_sum, ept(i,j,k), step, i, j, k
!write(*,*) pkxn, pkyn, 'step', 'ept_sum', step, ept_sum, l_axcel, i, j
endif
l_axcel=0
enddo
enddo
output_ept(k) = ept_sum / step
!write(*,*) 'output_ept(k) k', output_ept(k), step, k
step=0
ept_sum=0
enddo
! real(kind=4) :: ept_sum, avept, aveptsrf, rde, cpe, le, c1, c2, c3, c4, avept, aveptsrf
!ept srf
!cal ept srf use kousiki
!cal ept z 1~zgrid
do i = 1, ${xgrid}
do j = 1, ${ygrid}
c1=1000/1000
pts(i,j)=c1**c2
pts(i,j)=pts(i,j)*tsrf(i,j)
ess(i,j)=17.27*(tsrf(i,j)-273.15)/(tsrf(i,j)-35.86) !#tetens
ess(i,j)=exp(ess(i,j))
wss(i,j)=0.622*ess(i,j)/1000
epts(i,j)=le*wss(i,j)/cpe
epts(i,j)=epts(i,j)/tsrf(i,j)
epts(i,j)=pts(i,j)*exp(epts(i,j))
enddo
enddo
!do i = 150, 200
!write(*,*) 'pts, ess, wss epts(272,294)', pts(272,294), ess(272,294), wss(272,294), epts(272,294)
!enddo
step=0
epts_sum=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
if (epts(i,j)<1000) then
epts_sum = epts_sum + epts(i,j)
step = step + 1
!write(*,*) 'epts_sum epts(i,j) step i j', epts_sum, epts(i,j), step, i, j
!write(*,*) pkxn, pkyn, 'step', 'epts_sum', step, ept_sum, l_axcel, i, j
endif
endif
l_axcel=0
enddo
enddo
aveptsrf = epts_sum / step
!write(*,*) 'aveptsrf ', aveptsrf, step
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! 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)
!write(33,*) 'time_step pn plonn platn v1n aveptsrf output_ept1000 output_ept800 output_ept600 output_ept400'
write(33,*) ${time_step}, pn, plonn, platn, v1n, aveptsrf, output_ept(1), output_ept(7), output_ept(9), output_ept(11)
write(*,*) ${time_step}, pn, plonn, platn, v1n, aveptsrf, output_ept(1), output_ept(7), output_ept(9), output_ept(11)
close(11)
close(12)
close(13)
close(14)
close(15)
close(16)
close(17)
close(31)
!close(32)
close(33)
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 ${ofname2}.txt ${ofname2}.dat
cp ${ofname3}.txt ${ofname3}.dat
chage_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_axcel1
ofname2=output_axcel_rv
####################################
####################################
#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 ${ofname2}.txt ${ofname2}.dat
#####################################
# 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_axcel1
ofname2=output_axcel_rv
####################################
####################################
#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 ${ofname2}.txt ${ofname2}.dat
2011年10月19日水曜日
gnuplot 2つの要素をプロット
#! /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_vamx
plot_file12=axcel_pmin_cwc
plot_file13=axcel_vamx_cwc
################################################
################################################
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 [990: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 y2tics 0, 5
set ytics nomirror
set xlabel "time step [h]"
set ylabel "[hPa]"
set y2label "[m/s]"
set term postscript
set output "${plot_file11}.ps"
replot
################################################
set yrange [990:1010]
set y2range [0:0.1]
plot "${read_file}" using 1:2 axis x1y1 with line title "${plot_file1}" linetype 3 linecolor rgbcolor "red" linewidth 3,\
"${read_file}" using 1:6 axis x1y2 with line title "${plot_file2}" linetype 3 linecolor rgbcolor "blue" linewidth 3
set y2tics 0, 0.02
set ytics nomirror
set xlabel "time step [h]"
set ylabel "[hPa]"
set y2label "[g/m2]"
set term postscript
set output "${plot_file12}.ps"
replot
################################################
set yrange [0:30]
set y2range [0:0.1]
plot "${read_file}" using 1:5 axis x1y1 with line title "${plot_file3}" linetype 3 linecolor rgbcolor "green" linewidth 3,\
"${read_file}" using 1:6 axis x1y2 with line title "${plot_file2}" linetype 3 linecolor rgbcolor "blue" linewidth 3
set y2tics 0, 0.02
set ytics nomirror
set xlabel "time step [h]"
set ylabel "[m/s]"
set y2label "[g/m2]"
set term postscript
set output "${plot_file13}.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
ps2pdf ${plot_file12}.ps 12_${plot_file12}.pdf
rm ${plot_file12}.ps
ps2pdf ${plot_file13}.ps 13_${plot_file13}.pdf
rm ${plot_file13}.ps
################################################
#plot_axcel1
################################################
read_file=output_axcelrv.dat
plot_file1=axcel_pmin
plot_file2=axcel_avcwc
plot_file3=axcel_vmax
################################################
plot_file11=axcel_pmin_vamx
plot_file12=axcel_pmin_cwc
plot_file13=axcel_vamx_cwc
################################################
################################################
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 [990: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 y2tics 0, 5
set ytics nomirror
set xlabel "time step [h]"
set ylabel "[hPa]"
set y2label "[m/s]"
set term postscript
set output "${plot_file11}.ps"
replot
################################################
set yrange [990:1010]
set y2range [0:0.1]
plot "${read_file}" using 1:2 axis x1y1 with line title "${plot_file1}" linetype 3 linecolor rgbcolor "red" linewidth 3,\
"${read_file}" using 1:6 axis x1y2 with line title "${plot_file2}" linetype 3 linecolor rgbcolor "blue" linewidth 3
set y2tics 0, 0.02
set ytics nomirror
set xlabel "time step [h]"
set ylabel "[hPa]"
set y2label "[g/m2]"
set term postscript
set output "${plot_file12}.ps"
replot
################################################
set yrange [0:30]
set y2range [0:0.1]
plot "${read_file}" using 1:5 axis x1y1 with line title "${plot_file3}" linetype 3 linecolor rgbcolor "green" linewidth 3,\
"${read_file}" using 1:6 axis x1y2 with line title "${plot_file2}" linetype 3 linecolor rgbcolor "blue" linewidth 3
set y2tics 0, 0.02
set ytics nomirror
set xlabel "time step [h]"
set ylabel "[m/s]"
set y2label "[g/m2]"
set term postscript
set output "${plot_file13}.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
ps2pdf ${plot_file12}.ps 12_${plot_file12}.pdf
rm ${plot_file12}.ps
ps2pdf ${plot_file13}.ps 13_${plot_file13}.pdf
rm ${plot_file13}.ps
登録:
投稿 (Atom)