!= 物理量変換 プログラム
!
!= convert physical value module
!
! Authors::   Masanori Onishi
! Version::   $Id: convert.f90,v 1.00 onishi$
! Tag Name::  $Name: dcrtm-20161018 $
! Copyright:: Copyright (C) GFD Dennou Club, 2016. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]


module convert
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !

  ! モジュール引用 ; USE statements
  !

  ! 物理定数設定1
  ! physical constants setteing 1
  !
!  use constants0, only: PI

  ! 物理定数設定2
  ! physical constants setteing 2
  !
  use planetconst, only: Grav


  implicit none

  private

  public :: KtoTauSSA

contains 

  subroutine KtoTauSSA (                        & ! 2016/08/05
    & kmax, mmax, wmax,                         &
    & r_Press, r_Temp, rm_MixRatio, m_MolWt,    &
    & rmw_RayScatK, rmw_LAbsorpK, rmw_CAbsorpK, &
    & rw_OptDepTOA, zw_SingleScatA              )
    !
    ! convert absorption and scattering coefficient 
    ! to optial depth and single scattering albedo
    !
    integer, intent(in ) :: kmax
    integer, intent(in ) :: mmax
    integer, intent(in ) :: wmax
    real(8), intent(in ) :: r_Press       (0:kmax)
    real(8), intent(in ) :: r_Temp        (0:kmax)
    real(8), intent(in ) :: rm_MixRatio   (0:kmax, 0:mmax)
    real(8), intent(in ) :: m_MolWt       (        0:mmax)
    real(8), intent(in ) :: rmw_RayScatK  (0:kmax, 0:mmax, 0:wmax)
    real(8), intent(in ) :: rmw_LAbsorpK  (0:kmax, 0:mmax, 0:wmax)
    real(8), intent(in ) :: rmw_CAbsorpK  (0:kmax, 0:mmax, 0:wmax)
    real(8), intent(out) :: rw_OptDepTOA  (0:kmax,         0:wmax)
    real(8), intent(out) :: zw_SingleScatA(1:kmax,         0:wmax)
    integer :: k_r, k_m, k_w
    real(8) :: r_MolAve                   (0:kmax                )
    real(8) :: z_MolAve                   (1:kmax                )
    real(8) :: zm_MixRatio                (0:kmax, 0:mmax        )
    real(8) :: zw_OptDepS                 (1:kmax,         0:wmax)
    real(8) :: zw_OptDepA                 (1:kmax,         0:wmax)


    do k_r = 0, kmax
      r_MolAve(k_r) = 0.0_8
      do k_m = 0, mmax
        r_MolAve(k_r) = r_MolAve(k_r) + m_MolWt(k_m)*rm_MixRatio(k_r,k_m)
      end do
    end do

    do k_r = 1, kmax

      z_MolAve(k_r) = ( &
        & r_MolAve(k_r-1) * r_Press(k_r-1)/r_Temp(k_r-1) + &
        & r_MolAve(k_r  ) * r_Press(k_r  )/r_Temp(k_r  ) ) &
        & /(r_Press(k_r-1)/r_Temp(k_r-1) + r_Press(k_r)/r_Temp(k_r))

      do k_m = 0, mmax
        zm_MixRatio(k_r,k_m) = ( &
          & rm_MixRatio(k_r-1,k_m) * r_Press(k_r-1)/r_Temp(k_r-1) + &
          & rm_MixRatio(k_r  ,k_m) * r_Press(k_r  )/r_Temp(k_r  ) ) &
          & /(r_Press(k_r-1)/r_Temp(k_r-1) + r_Press(k_r)/r_Temp(k_r))
      end do

    end do 

    zw_OptDepS = 0.0_8
    do k_r = 1, kmax
      do k_m = 0, mmax
        if(zm_MixRatio(k_r,k_m) .ne. 0.0_8) then
          zw_OptDepS(k_r,:) = zw_OptDepS(k_r,:) + &
            &                  ( &
            &                    (rmw_RayScatK(k_r-1,k_m,:)                            )* &
            &                    m_MolWt(k_m)*r_Press(k_r-1)*rm_MixRatio(k_r-1,k_m)/r_Temp(k_r-1) + &
            &                    (rmw_RayScatK(k_r  ,k_m,:)                            )* &
            &                    m_MolWt(k_m)*r_Press(k_r  )*rm_MixRatio(k_r  ,k_m)/r_Temp(k_r  ) )/&
            &                    ( r_Press(k_r-1)*rm_MixRatio(k_r-1,k_m)/r_Temp(k_r-1) +    &
            &                      r_Press(k_r  )*rm_MixRatio(k_r  ,k_m)/r_Temp(k_r  ) )*   &
            &                    zm_MixRatio(k_r,k_m)
        end if
      end do !k_m in k_r
    end do ! k_r

    ! test
    !do k_m = 1, mmax
    !  print *, 'k_m', k_m
    !  print *, 'zvalue'
    !  do k_r = 1, kmax
    !    if(zm_MixRatio(k_r,k_m) .ne. 0.0_8) then
    !      print *,  + &
    !        &                  ( &
    !        &                    (rmw_LAbsorpK(k_r-1,k_m,:)                            )* &
    !        &                    m_MolWt(k_m)*r_Press(k_r-1)*rm_MixRatio(k_r-1,k_m)/r_Temp(k_r-1) + &
    !        &                    (rmw_LAbsorpK(k_r  ,k_m,:)                            )* &
    !        &                    m_MolWt(k_m)*r_Press(k_r  )*rm_MixRatio(k_r  ,k_m)/r_Temp(k_r  ) )/&
    !        &                    ( r_Press(k_r-1)*rm_MixRatio(k_r-1,k_m)/r_Temp(k_r-1) +    &
    !        &                    r_Press(k_r  )*rm_MixRatio(k_r  ,k_m)/r_Temp(k_r  ) )*   &
    !        &                    zm_MixRatio(k_r,k_m)
    !    end if
    !  end do !k_m in k_r
    !end do ! k_r
    ! test

    zw_OptDepA = 0.0_8
    do k_r = 1, kmax
      do k_m = 1, mmax
        if(zm_MixRatio(k_r,k_m) .ne. 0.0_8) then
          zw_OptDepA(k_r,:) = zw_OptDepA(k_r,:) + &
            &                  ( &
            &                    (rmw_LAbsorpK(k_r-1,k_m,:)                            )* &
            &                    m_MolWt(k_m)*r_Press(k_r-1)*rm_MixRatio(k_r-1,k_m)/r_Temp(k_r-1) + &
            &                    (rmw_LAbsorpK(k_r  ,k_m,:)                            )* &
            &                    m_MolWt(k_m)*r_Press(k_r  )*rm_MixRatio(k_r  ,k_m)/r_Temp(k_r  ) )/&
            &                    ( r_Press(k_r-1)*rm_MixRatio(k_r-1,k_m)/r_Temp(k_r-1) +    &
            &                    r_Press(k_r  )*rm_MixRatio(k_r  ,k_m)/r_Temp(k_r  ) )*   &
            &                    zm_MixRatio(k_r,k_m)
        end if
      end do !k_m in k_r
    end do ! k_r

    ! test
    !do k_r = 1, kmax
    !  print *, zw_OptDepA(k_r,:)
    !end do
    ! test

    do k_r = 1, kmax
      do k_m = 1, mmax
        if(zm_MixRatio(k_r,k_m) .ne. 0.0_8) then
          zw_OptDepA(k_r,:) = zw_OptDepA(k_r,:) + &
            &                  ( &
            &                    (rmw_CAbsorpK(k_r-1,k_m,:)                            )* &
            &                    m_MolWt(k_m)*r_Press(k_r-1)*rm_MixRatio(k_r-1,k_m)/r_Temp(k_r-1) + &
            &                    (rmw_CAbsorpK(k_r  ,k_m,:)                            )* &
            &                    m_MolWt(k_m)*r_Press(k_r  )*rm_MixRatio(k_r  ,k_m)/r_Temp(k_r  ) )/&
            &                    ( r_Press(k_r-1)*rm_MixRatio(k_r-1,k_m)/r_Temp(k_r-1) +    &
            &                      r_Press(k_r  )*rm_MixRatio(k_r  ,k_m)/r_Temp(k_r  ) )*   &
            &                    zm_MixRatio(k_r,k_m)
        end if
      end do !k_m in k_r
    end do ! k_r

    ! test
    !do k_r = 1, kmax
    !  print *, zw_OptDepA(k_r,:)
    !end do
    ! test

    do k_r = 1, kmax
      do k_w = 0, wmax
        if( (zw_OptDepA(k_r,k_w)+zw_OptDepS(k_r,k_w)) .ne. 0.0_8 ) then  ! 2016/10/18
          zw_SingleScatA(k_r,k_w) = zw_OptDepS(k_r,k_w)/(zw_OptDepA(k_r,k_w)+zw_OptDepS(k_r,k_w))
        else
          zw_SingleScatA(k_r,k_w) = 0.0_8
        end if
      end do ! k_w in k_r
    end do ! k_r  

    rw_OptDepTOA(kmax,:) = 0.0_8
    do k_r = kmax-1, 0, -1
      rw_OptDepTOA(k_r,:) = rw_OptDepTOA(k_r+1,:) + &
        &              ( zw_OptDepA(k_r+1,:) + zw_OptDepS(k_r+1,:))*&
        &              ( r_Press(k_r) - r_Press(k_r+1))/        &
        &              ( Grav*z_MolAve(k_r+1) )
    end do

    ! test
    !do k_r = 1, kmax
    !  print *, rw_OptDepTOA(k_r,:)
    !  print *, z_MolAve(k_r)
    !  print *, (r_Press(k_r-1) - r_Press(k_r))/(Grav*z_MolAve(k_r))
    !end do
    ! test

  end subroutine KtoTauSSA

end module convert
