!---------------------------------------------------------------------
!  Fortran program to calculate local time comoposite of Tbb (by Kawashima)
!  1時間毎のTbbデータ[ 配列：東西×南北×時間 = im x jm x lm]から、
!  local time composite (各時刻（00, 01, 02,....23LT) の平均値）を計算する. 
!  各点の値および、空間方向に平滑化した値を出力する. 
!  lm: 入力データの時間方向の数 
!  nday: compositeを取る日数 
!  ism,jsm : 東西x南北（2*ism+1)x(2*jsm+1)点について平均
!-------------------------------------------------------------------
      PARAMETER (im=1752,jm=1502,lm=72)  
      PARAMETER (ld=24,lini=10,nday=3)  
      PARAMETER (ism=10,jsm=10)  
      INTEGER i,j,k,l,m
      REAL TBB(im,jm,lm)    !入力データ
      REAL TBBC(im,jm,ld)   !compositeデータ
      REAL TBBCS(im,jm,ld)    !平滑化compositeデータ
      REAL TBBA(im,jm)      ! compositeデータの24時間平均値
      REAL TBBAS(im,jm)      !平滑化compositeデータの24時間平均値
      REAL TBBC2(im,jm,2*ld+1)   !compositeデータの出力用配列
      REAL TBBCS2(im,jm,2*ld+1)   !平滑化compositeデータの出力用配列

! 入力データを開く
       open (10,file='tbb_combined.dat',
     *      form='unformatted',access='direct',recl=im)

! 出力データを開く
       open (20,file='tbb_composite.dat',
     *      form='unformatted',access='direct',recl=im)


      irecin=0
      irecout=0

      do 100 l=1,lm    ! 時間のループ
      write (6,*)'input L=', L

! 入力データの読込み
      do 50 j=1,jm
      irecin=irecin+1   
      read(10,rec=irecin)  (TBB(i,j,l),i=1,im)
 50   continue

      write (6,*)'TBB(500,500)=', TBB(500,500,l)

 100  continue         !時間のループ終わり

!!!計算確認用ダミーデータの作成!!!!!!!!!
!!      do 60 l=1,lm
!!      do 60 j=1,jm
!!      do 60 i=1,im
!!      TBB(i,j,l)=l
!! 60   continue
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! compositeデータの初期化
      do 200 l=1,ld     
      do 200 j=1,jm
      do 200 i=1,im
      TBBC(i,j,l)=0.0
      TBBCS(i,j,l)=0.0 
 200  CONTINUE

! compositeデータの平均値の初期化
      do 210 j=1,jm
      do 210 i=1,im
      TBBA(i,j)=0.0
      TBBAS(i,j)=0.0
 210  CONTINUE

!以下、入力データtbb_combined.datが 00 LT (=09UTC) から始まっている場合を想定している.
!この時刻を起点に、1時間ごとにずらしつつ、24時間おきのデータを加算していく

      do 300 l=1,ld                  !時刻(00,01,02,...,23LT)についてのループ
      do 300 m=l,l+ld*(nday-1),ld    !ld=24 時間おき
      do 300 j=1,jm
      do 300 i=1,im
      TBBC(i,j,l)=TBBC(i,j,l)+TBB(i,j,m)/nday
 300  continue


! TBBCを東西南北方向に平均（2*ism+1)x(2*jsm+1) をとることで平滑化
      do 320 l=1,ld                 
      do 320 j=1+jsm,jm-jsm
      do 320 i=1+ism,im-ism
      do 320 m=j-jsm,j+jsm
      do 320 k=i-ism,i+ism
      TBBCS(i,j,l)=TBBCS(i,j,l)+TBBC(k,m,l)/(2*ism+1)/(2*jsm+1)
 320  continue


! 日平均からの偏差を見るため、日平均値を計算する
      do 350 l=1,ld     
      do 350 j=1,jm
      do 350 i=1,im
      TBBA(i,j)=TBBA(i,j)+TBBC(i,j,l)/ld 
      TBBAS(i,j)=TBBAS(i,j)+TBBCS(i,j,l)/ld 
 350  CONTINUE

!作図等の便宜上、00～23LTを2サイクル分＋00LT (2*24+1=49時刻分）を出力するため、
!TBBC2（im,jm,2*ld+1)にTBBCを入れる

      do 400 l=1,ld     
      do 400 j=1,jm
      do 400 i=1,im
      TBBC2(i,j,l)=TBBC(i,j,l)
      TBBCS2(i,j,l)=TBBCS(i,j,l)
 400  CONTINUE

      do 410 l=ld+1,2*ld
      do 410 j=1,jm
      do 410 i=1,im
      TBBC2(i,j,l)=TBBC(i,j,l-ld)
      TBBCS2(i,j,l)=TBBCS(i,j,l-ld)
 410  CONTINUE

      do 420 j=1,jm
      do 420 i=1,im
      TBBC2(i,j,2*ld+1)=TBBC(i,j,1)
      TBBCS2(i,j,2*ld+1)=TBBCS(i,j,1)
 420  CONTINUE

! 以下、データの出力
      write (6,*) 'L,TBBC2(500,500,l),TBBCS2(500,500,l)'
      do 600 l=1,2*ld+1

!  動作確認のため、ある点(ここでは i,j = 500,500) の値を画面上に表示する
      write (6,*) L,TBBC2(500,500,l),TBBCS2(500,500,l)

! compositeデータの出力
      do 500 j=1,jm
      irecout=irecout+1
 500  write(20,rec=irecout)  (TBBC2(i,j,l),i=1,im)

! compositeの日平均からの偏差の出力
      do 520 j=1,jm
      irecout=irecout+1
 520  write(20,rec=irecout)  (TBBC2(i,j,l)-TBBA(i,j),i=1,im)

! 平滑化compositeの日平均からの偏差の出力
      do 540 j=1,jm
      irecout=irecout+1
 540  write(20,rec=irecout)  (TBBCS2(i,j,l),i=1,im)

! 平滑化compositeの日平均からの偏差の出力
      do 560 j=1,jm
      irecout=irecout+1
 560  write(20,rec=irecout)  (TBBCS2(i,j,l)-TBBAS(i,j),i=1,im)


 600  continue

      stop
      end 
