!------------------------------------------------------------------
!  Fortran program to calculate fractional area from HIMAWARI Tbb data 
!  (by Kawashima, 2023/07/27)
!  Tbbデータ[配列：IM x JM = 6000 x 6000] から、
!  特定の長方形領域 (I= is to ie, J=js to je) における
!  各温度帯 (T=TH0(K)からDTH(K)おき)の格子数および占有率(in %)を求めるプログラム
!  さらに、格子数および占有率(in %)の時刻毎のコンポジットも作成する。
!  LM: 読み込むべきHIMAWARIデータの数
!  NM: 温度帯の数
!  MM: compositeの時刻の数(=24) 
!  NDAY:読み込むHIMAWARIデータの日数 (=LM/MM)
!   ifort -convert big_endian calc-fraction-1K.f 
!   表示用 ctlファイルは tbb-fraction-1K.ctl と tbb-fraction-composite-1K.ctl
!-------------------------------------------------------------------
      PARAMETER (IM=6000,JM=6000,LM=480)
      PARAMETER (NM=100,MM=24,NDAY=20)  
      CHARACTER filename_in*100
      REAL TBB(IM,JM)      ! 元のHIMAWARIデータ
      REAL TBBR(IM,JM)     ! 元のHIMAWARIデータの緯度方向の順番を逆にしたもの(南→北の順)
      REAL TH(NM+1)          ! 各温度帯の閾値 (NM+1個)
      INTEGER NUM(NM,LM)   ! 個数 
      INTEGER NUMC(NM,MM)  ! 個数 (composite)
      REAL FRAC(NM,LM)     ! 占有率  (in %)
      REAL FRACC(NM,MM)    ! 占有率 (composite, in %)

!     温度の閾値(DTHおきに、NM 個）の設定      
      TH0=220.0
      DTH=1.0

      do N=1,NM+1
      TH(N)=TH0+DTH*(N-1)
      end do 

! count 数の初期化      
      do L=1,LM
         do N=1,NM
           NUM(N,L)=0
         end do 
      end do 

      do M=1,MM
         do N=1,NM
           NUMC(N,M)=0
         end do
      end do 

!    対象領域(長方形領域)の設定
!    HIMAWARIデータの格子点は0.02度おき
!    1 degree =50 grid (0.02 degree interval)
      
      is=2500    ! lon=134.99E
      ie=2650    ! lon=137.99E
      js=4900    ! lat=37.99N
      je=5000    ! lat=39.99N

      
! open input data file name (YYYYMMDDHHMM_TIR_01.dat)
! 入力データ名のリストを開く
      open (10,file='filename_in.txt',status='old')

! 出力データを開く(各データ)  
      open (40,file='tbb-fraction-1K.dat',
     *          form='unformatted',access='direct',recl=nm)

! 出力データを開く(local time compositeデータ, compositeを作らないなら不要)  
      open (50,file='tbb-fraction-composite-1K.dat',
     *          form='unformatted',access='direct',recl=nm)


      irecout=0
      irecout2=0      
      
      do 1000 L=1,LM            ! 時間のループ

      irecin=0

      write (6,*)'L=', L

! 入力データ名の読み出し
      read (10,'(a)') filename_in

! 入力データ名を画面上に書き出す
      write(6,*) filename_in

! 上で読み出した名前の入力データを開く
       open (30,file=filename_in,convert="little_endian",
     *          form='unformatted',access='direct',recl=im)

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

!  元のHIMAWARIデータはy方向の並びが高緯度->低緯度となっていて分かりにくいので、
!  低緯度->高緯度の並びのデータ(TBBR)に直す

      do 200 J=1,JM
      do 200 I=1,IM
      TBBR(I,J)=TBB(I,JM+1-J)
 200  CONTINUE

!    I=is～ie, J=js～je の長方形の領域について
!    TH(N) <= TBBR < TH(N+1) の格子数を計算する

      do 220 N=1,NM

      do 210 J=js,je
      do 210 I=is,ie
      IF (TBBR(I,J) .ge. TH(N) .and. TBBR(I,J) .lt. TH(N+1)) THEN
           NUM(N,L)=NUM(N,L)+1
        ELSE
           NUM(N,L)=NUM(N,L)
      END IF
 210  CONTINUE

!     格子点数から占有面積率(in %)に変換する
      FRAC(N,L)=REAL(NUM(N,L))/(ie-is+1)/(je-js+1)*100.      
 220  CONTINUE
      
!     check       
         write(6,*)'L,N,Range,NUMBER,FRACION'
      do 230 N=1,NM
         write(6,*) L,N,TH(N),TH(N)+1.0,NUM(N,L),FRAC(N,L)
 230  CONTINUE

!    格子点数、占有率データの書きだし
      irecout=irecout+1
      write(40,rec=irecout) (REAL(NUM(N,L)),N=1,NM) 

      irecout=irecout+1
      write(40,rec=irecout) (FRAC(N,L),N=1,NM) 
       
 1000 continue                            !時間のループ終了


!    local time コンポジットを作らないなら、これ以下は削除か読み飛ばすようにする (GO TO 700)
!    以降、各時刻のコンポジットの作成. 
!    L=MからL=M+(NDAY-1)*MMまで、MM=24個おきに加算していく作業をM=1(=00LT)からM=MM(=23LT)まで行う

      do 400 N=1,NM                  ! 温度帯の数についてのループ
      do 400 M=1,MM                  ! 時刻のループ (00 LT to 23LT)
      do 400 L=M,M+(NDAY-1)*MM,MM    ! L=MからL=M+(NDAY-1)*MMまで、MM=24個おきに加算する
      NUMC(N,M)=NUMC(N,M)+NUM(N,L)/NDAY
 400  CONTINUE

      DO 500 M=1,MM      
      DO 500 N=1,NM
      FRACC(N,M)=REAL(NUMC(N,M))/(ie-is+1)/(je-js+1)*100.      
 500  CONTINUE

! check 
       write(6,*) 'M,N,NUMC(N,M),FRACC(N,M)'            
      DO 550 M=1,MM      
      DO 550 N=1,NM      
       write(6,*) M,N,NUMC(N,M),FRACC(N,M)      
 550  CONTINUE
      
      do 600 M=1,MM

      irecout2=irecout2+1
      write(50,rec=irecout2) (REAL(NUMC(N,M)),N=1,NM) 

      irecout2=irecout2+1
      write(50,rec=irecout2) (FRACC(N,M),N=1,NM) 
      
 600  CONTINUE

 700  CONTINUE
      
      stop
      end 

      
