NHMの計算結果からpminとvmaxそして、それぞれの位置を求めるシェルスクリプト

NHMの計算結果からpminとvmaxそして、それぞれの位置を求めるシェルスクリプト

find_n_v_pmin.sh



#! /bin/tcsh -f

#####################################
filename=find_n_v_pmin


#type yyyymmddhhhh start and end time
: ${start:=199201020000}
: ${end:=199201081000}
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=6

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


: ${prefix:=anl_p}    #want data

    cat > data.gs <<EOF
'reinit'
'open _NHMMRPPFCSVSTD1_${time}.ctl'

'set gxout fwrite'
'set undef dfile'
'set fwrite ${filename}.grd'
*'set t ${t}'
'set x 1 400'
'set y 1 300'
'd psea'
'disable fwrite'
'c'
'set gxout fwrite'
'set undef dfile'
'set fwrite v_${filename}.grd'
*'set t ${t}'
'set x 1 400'
'set y 1 300'
'd sqrt(u*u+v*v)'
'disable fwrite'
'c'

'quit'

#157
EOF
    grads -blc data.gs
########################################
####program read_data###################
########################################

if [ ${time} -eq ${start} ]; then
##start file content
cat <<EOF > start.f90
program start
 implicit none
open (unit=7, file='${filename}.txt', form='FORMATTED')
write(7,*) '${filename} NH'
write(7,*) 'time_step mini_psea xgrid ygrid lon lat vmax vlon vlat'
write(*,*) 'time_step mini_psea xgrid ygrid lon lat vmax vlon vlat'
close(7)
end program start
EOF
##output data.txt by f90
f90 -o start.out start.f90
./start.out
##start
fi


####find min psea
cat <<EOF > ${filename}.f90
program find_pmin
 implicit none
 integer, parameter :: xgrid=400, ygrid=300, rgrid=10
 integer :: n, i, j, kx, ky, is, ie, js, je
 real(kind=4), dimension(xgrid,ygrid) :: p_d
 real(kind=4), dimension(xgrid,ygrid) :: v_d
 real(kind=4) :: p, pp, lon, lat, v, vv, vlon ,vlat
open (unit=17,file='${filename}.grd',form='unformatted',access='direct',recl=xgrid*ygrid*4)
read (unit=17,rec=1) p_d
open (unit=18,file='v_${filename}.grd',form='unformatted',access='direct',recl=xgrid*ygrid*4)
read (unit=18,rec=1) v_d

open (unit=7, file='${filename}.txt', position='append', form='FORMATTED')

p=2000
do i = 139, 277 !155~180
   do j = 150, ygrid !150 @NH
    pp= p_d(i-1,j) + p_d(i+1,j) + p_d(i,j-1) + p_d(i,j-1)
    pp= pp / 4
      if (p_d(i,j) < pp) then  
    if (p_d(i,j) < p) then
        kx=i
        ky=j
        p=p_d(i,j)
!         write(*,*) n, ${time_step}, i, j, p
!         write(7,*) n, ${time_step}, i, j, p
    endif
      endif  
   n=n+1 
   enddo
enddo

lon=kx*0.18+130
lat=ky*0.18-30


is=kx-10
js=ky-10
ie=kx+10
je=ky+10

do  i = is, ie
    do j = js, je
    vv= v_d(i-1,j) + v_d(i+1,j) + v_d(i,j-1) + v_d(i,j-1)
    vv= vv / 4
      if (v_d(i,j) > vv) then  
    if (v_d(i,j) > v) then
        kx=i
        ky=j
        v=v_d(i,j)  
    endif
      endif  
    enddo
enddo

vlon=kx*0.18+130
vlat=ky*0.18-30


write(*,*)  ${time_step}, p, kx, ky, lon, lat, v, vlon, vlat
write(7,*)  ${time_step}, p, kx, ky, lon, lat, v, vlon, vlat 

close(7)
close(17)
close(18)

end program  find_pmin
EOF


##output data.txt by f90
f90 -o ${filename}.out ${filename}.f90
./${filename}.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