3年分の振幅計算（1時間の中央値）

#使い方
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
```