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