使い方
N.MWDH観測点の2-5Hz,15時〜一時間の中央値の場合
./allmedian.sh _N.MWDH_U_bp2-5_sqr 15
allmedian.sh
#!/bin/bash
#./allmedian.sh _N.MMOH_U_bp2-5_sqr 15
for yr in {19..20}
do
for month in {1..9}
do
for day in {1..9}
do
./medi2 ${yr}0${month}/${yr}0${month}0${day}${1}.s ${2} >> allmedi${1}_${2}.txt
done
for day in {10..31}
do
./medi2 ${yr}0${month}/${yr}0${month}${day}${1}.s ${2} >> allmedi${1}_${2}.txt
done
done
for month in {10..12}
do
for day in {1..9}
do
./medi2 ${yr}${month}/${yr}${month}0${day}${1}.s ${2} >> allmedi${1}_${2}.txt
done
for day in {10..31}
do
./medi2 ${yr}${month}/${yr}${month}${day}${1}.s ${2} >> allmedi${1}_${2}.txt
done
done
done
########
yr=21
for month in {1..5}
do
for day in {1..9}
do
./medi2 ${yr}0${month}/${yr}0${month}0${day}${1}.s ${2} >> allmedi${1}_${2}.txt
done
for day in {10..31}
do
./medi2 ${yr}0${month}/${yr}0${month}${day}${1}.s ${2} >> allmedi${1}_${2}.txt
done
done
# for month in {10..12}
#
# for day in {1..9}
# ./medi2 ${yr}${month}/${yr}${month}0${day}${1}.s ${2} >> allmedi${1}_${2}.txt
# done
#
# for day in {10..31}
# ./medi2 ${yr}${month}/${yr}${month}${day}${1}.s ${2} >> allmedi${1}_${2}.txt
# done
# done
スクリプト中に含まれるmedi2は以下のプログラムから
sac_median_2.f95
program main
implicit none
integer i, a(40), max, npts, hour
real b(70),c(2), dt, med
real, allocatable :: x(:), y(:)
character :: file1*80, arg*10, khd(48)*4, name*10
if(iargc()>0) then
call getarg (1, file1) !コマンドラインで引数!
call getarg (2, arg) !コマンドラインで引数所得
read(arg,*) hour
endif
max=10000000
allocate (x(max),y(max))
call arsac(file1,b,a,khd,x,max)
npts = a(10)
dt = b(1)
do i = hour*360000+1,(hour+1)*360000+1
y(i-hour*360000) = x(i)
enddo
call median(y, med, 360000)
med = med*1e-9 !nm/sからm/sへ変換
! write(*,*)med
name = '20'//file1(6:11)
write(*,*)name, med
end program main
!########################################################
subroutine arsac(fn,fhd,nhd,khd,y,maxy)
real fhd(70),y(maxy)
integer nhd(40),i,i1,i2,i3,npts
character khd(48)*4,fn*80
open(1,file=fn,access='direct',recl=632,status='old')
read(1,rec=1)(fhd(i1),i1=1,70),(nhd(i2),i2=1,40),(khd(i3),i3=1,48)
npts=nhd(10)
close(1)
open(1,file=fn,access='direct',recl=4,status='old')
do i=1,npts
! read(1,rec=i+158,end=100)y(i) Mar.1 H.Aoyama
read(1,rec=i+158)y(i)
enddo
! 100 continue
close(1)
return
end subroutine arsac
!#####################################################
subroutine median(a, med, n)
implicit none
integer,intent(in):: n
real,intent(inout):: a(n)
real,intent(inout):: med
integer i, m
real am
!バブルソート(遅い)
! do i = 1, n-1 !sort a
! am = minval(a(i+1:n))
! m = minloc(a(i+1:n), 1) + i
! if (a(i) > am) then
! a(m) = a(i)
! a(i) = am
! endif
! enddo
call quicksort(a,1,360001)
if(mod(n, 2) /= 0) then !median of a(:)
med = a((n-1)/2+1)
else
med = 0.5*(a(n/2)+a(n/2+1))
endif
end subroutine median
!##############################################
recursive subroutine quicksort(a, first, last)
implicit none
real a(*)
integer first, last
integer i, j
real x, t
x = a( (first+last) / 2 )
i = first
j = last
do
do while (a(i) < x)
i=i+1
end do
do while (x < a(j))
j=j-1
end do
if (i >= j) exit
t = a(i); a(i) = a(j); a(j) = t
i=i+1
j=j-1
end do
if (first < i - 1) call quicksort(a, first, i - 1)
if (j + 1 < last) call quicksort(a, j + 1, last)
end subroutine quicksort