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