!= flux 計算のための大気層分割プログラム
!
!= module of dividing vertical layer for flux calculation
!
! Authors::   Masanori Onishi
! Version::   $Id: main_dcrtm.f90,v 1.00 onishi$
! Tag Name::  $Name: dcrtm-20170202 $
! Copyright:: Copyright (C) GFD Dennou Club, 2017. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]

!2016/12/22
!2017/01/05
!2017/02/02

module mod_optdepitr

  ! 格子点設定
  ! grid point setting
  !
  use gridset, only: GridsetSet2, &
    &                imax2,       &
    &                jmax2,       &
    &                kmax2

  ! 2方向近似放射モデル
  ! two stream approximation
  !
  use rad_rte_two_stream_app_itg2,                &
    &          only: RadRTETwoStreamAppInit_itg2, &
    &                RadRTETwoStreamAppLW_itg2


  ! 宣言文 ; Declaration statements
  !
  implicit none

  private

  public :: DivLayCount

contains

  subroutine DivLayCount (  &
    & eps,                  &
    & DelOptDep, SSA, AF,   &
    & PFsrf, PFtop, Albedo, &
    & NLay                  )

    real(8), intent(in ) :: eps
    real(8), intent(in ) :: DelOptDep
    real(8), intent(in ) :: SSA
    real(8), intent(in ) :: AF
    real(8), intent(in ) :: PFsrf
    real(8), intent(in ) :: PFtop
    real(8), intent(in ) :: Albedo
    integer, intent(out) :: NLay

    real(8), allocatable :: xyz_SSA          (:, :, :)
    real(8), allocatable :: xyz_AF           (:, :, :)
    real(8), allocatable :: xyr_OptDep       (:, :, :)
    real(8), allocatable :: xy_SurfAlbedo    (:, :)
    real(8), allocatable :: xyr_PFInted      (:, :, :)
    real(8), allocatable :: xy_SurfPFInted   (:, :)
    real(8), allocatable :: xy_SurfDPFDTInted(:, :)
    real(8), allocatable :: xyr_RadUwFlux    (:, :, :)
    real(8), allocatable :: xyr_RadDwFlux    (:, :, :)
    real(8), allocatable :: xyra_DelRadUwFlux(:, :, :, :)
    real(8), allocatable :: xyra_DelRadDwFlux(:, :, :, :)

    integer :: k_itr, k_r
    real(8) :: FXUP1, FXUP2, FXDN1, FXDN2
    real(8) :: TRUP, TRDN, preTRUP, preTRDN
    real(8) :: ERUP, ERDN
    integer :: MAXNUM

    MAXNUM = 20 ! 2017/02/02

    NLay = 1
    
    preTRUP = -1.0_8
    preTRDN = -1.0_8

    do k_itr = 0, MAXNUM

      kmax2 = 2**k_itr

      call GridsetSet2(imax2, jmax2, kmax2)

      allocate (xyz_SSA          (0:imax2-1, 1:jmax2, 1:kmax2) )
      allocate (xyz_AF           (0:imax2-1, 1:jmax2, 1:kmax2) )
      allocate (xyr_OptDep       (0:imax2-1, 1:jmax2, 0:kmax2) )
      allocate (xy_SurfAlbedo    (0:imax2-1, 1:jmax2         ) )
      allocate (xyr_PFInted      (0:imax2-1, 1:jmax2, 0:kmax2) )
      allocate (xy_SurfPFInted   (0:imax2-1, 1:jmax2         ) )
      allocate (xy_SurfDPFDTInted(0:imax2-1, 1:jmax2         ) )
      allocate (xyr_RadUwFlux    (0:imax2-1, 1:jmax2, 0:kmax2) )
      allocate (xyr_RadDwFlux    (0:imax2-1, 1:jmax2, 0:kmax2) )
      allocate (xyra_DelRadUwFlux(0:imax2-1, 1:jmax2, 0:kmax2, 0:1) )
      allocate (xyra_DelRadDwFlux(0:imax2-1, 1:jmax2, 0:kmax2, 0:1) )
 
      ! フラックスの計算
      ! calculate flux  
      !

      ! dcpam5
      xy_SurfAlbedo    (:, :)    = 0.0_8 !Albedo
      xy_SurfDPFDTInted(:, :)    = 0.0_8
      xyz_SSA = SSA
      xyz_AF  = AF

      do k_r = 0, kmax2
        xyr_OptDep(:,:,k_r) = DelOptDep/kmax2 * (kmax2 - k_r)
        !print *, xyr_OptDep(:,:,k_r)
      end do

!    if(.true.) then
      do k_r = 0, kmax2
        xyr_PFInted(:,:,k_r) = PFsrf + (PFtop - PFsrf)*k_r/kmax2
        !print *, xyr_PFInted(:,:,k_r)
      end do
      xy_SurfPFInted   (:, :)    = xyr_PFInted (:,:,0)
!    else
!      xyr_PFInted = 0.0_8
!      xy_SurfPFInted = 1.0_8
!    end if

      call RadRTETwoStreamAppLW_itg2(                      &
        & xyz_SSA, xyz_AF,                                 & ! (in)
        & xyr_OptDep,                                      & ! (in)
        & xy_SurfAlbedo,                                   & ! (in)
        & xyr_PFInted, xy_SurfPFInted, xy_SurfDPFDTInted,  & ! (in)
        & xyr_RadUwFlux, xyr_RadDwFlux,                    & ! (out)
        & xyra_DelRadUwFlux, xyra_DelRadDwFlux             & ! (out)
        & )

      FXUP1 = xyr_RadUwFlux(0,1,0)
      FXUP2 = xyr_RadUwFlux(0,1,kmax2)
      FXDN1 = xyr_RadDwFlux(0,1,0)
      FXDN2 = xyr_RadDwFlux(0,1,kmax2)

      if(FXUP1 == 0.0_8) then
        FXUP1 = 1.0_8
      end if

      if(FXDN2 == 0.0_8) then
        FXDN2 = 1.0_8
      end if

      TRUP  = FXUP2/FXUP1
      TRDN  = FXDN1/FXDN2
      !TRUP  = FXUP2!/FXUP1
      !TRDN  = FXDN1!/FXDN2

      !print *, kmax2, TRUP, TRDN
      ERUP = (TRUP - preTRUP)/preTRUP
      ERDN = (TRDN - preTRDN)/preTRDN
      !print *, 'ERR', ERUP, ERDN

      if( (abs(ERUP) < eps) .and. (abs(ERDN) < eps) ) then
        !print *, 'Flux calculation is enough'
        NLay = 2**(k_itr-1)
        exit
      else
        preTRUP = TRUP
        preTRDN = TRDN
      end if

      deallocate (xyz_SSA          )
      deallocate (xyz_AF           )
      deallocate (xyr_OptDep       )
      deallocate (xy_SurfAlbedo    )
      deallocate (xyr_PFInted      )
      deallocate (xy_SurfPFInted   )
      deallocate (xy_SurfDPFDTInted)
      deallocate (xyr_RadUwFlux    )
      deallocate (xyr_RadDwFlux    )
      deallocate (xyra_DelRadUwFlux)
      deallocate (xyra_DelRadDwFlux)

    end do ! k_itr

    if( (kmax2 == 2**MAXNUM).and.(NLay == 1) ) then
      print *, 'Flux calculation is not enough: divided layer number=', kmax2
      stop
    end if

  end subroutine DivLayCount

end module mod_optdepitr
