#! /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
0 件のコメント:
コメントを投稿