!= デフォルトの基本場を設定するための変数参照型モジュール
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: basicset.f90,v 1.18 2014/07/08 01:05:32 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2007. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!

module basicset_3d
  !
  != デフォルトの基本場を設定するための変数参照型モジュール
  !
  ! hogeBZ となる配列は, 外部から値を入力することで初期化を行うことに注意. 
  !

  !モジュール読み込み
  use dc_types,   only:  DP
  use gridset,    only:  imin,       &! 配列の X 方向の下限
    &                    imax,       &! 配列の X 方向の上限
    &                    jmin,       &! 配列の Z 方向の下限
    &                    jmax,       &! 配列の Z 方向の上限
    &                    kmin,       &! 配列の Z 方向の下限
    &                    kmax,       &! 配列の Z 方向の上限
    &                    ncmax        ! 化学種の数  
  use composition, only: GasNum,     &!気体の数
    &                    IdxG,       &!気体の配列添え字
    &                    MolWtWet     !凝縮成分の分子量
  use axesset,  only:    dz

  !暗黙の型宣言禁止
  implicit none

  !save 属性
  private

  !Public Interface
  real(DP), allocatable, save, public :: xyz_DensBasicZ(:,:,:)    !密度
  real(DP), allocatable, save, public :: pyz_DensBasicZ(:,:,:)    !密度
  real(DP), allocatable, save, public :: xqz_DensBasicZ(:,:,:)    !密度
  real(DP), allocatable, save, public :: xyr_DensBasicZ(:,:,:)    !密度
  real(DP), allocatable, save, public :: xyz_PressBasicZ(:,:,:)   !無次元圧力
  real(DP), allocatable, save, public :: xyz_ExnerBasicZ(:,:,:)   !無次元圧力
  real(DP), allocatable, save, public :: xyz_TempBasicZ(:,:,:)    !温度
  real(DP), allocatable, save, public :: xyz_PotTempBasicZ(:,:,:)   !温位
  real(DP), allocatable, save, public :: xyr_PotTempBasicZ(:,:,:)   !温位
  real(DP), allocatable, save, public :: xyz_VPotTempBasicZ(:,:,:)  !仮温位
  real(DP), allocatable, save, public :: pyz_VPotTempBasicZ(:,:,:)  !仮温位
  real(DP), allocatable, save, public :: xqz_VPotTempBasicZ(:,:,:)  !仮温位
  real(DP), allocatable, save, public :: xyr_VPotTempBasicZ(:,:,:)  !仮温位
  real(DP), allocatable, save, public :: xyz_VelSoundBasicZ(:,:,:)!音速
  real(DP), allocatable, save, public :: xyz_VelSW(:,:,:)   !音速
  real(DP), allocatable, save, public :: xyza_MixRtBasicZ(:,:,:,:) !凝縮成分混合比
  real(DP), allocatable, save, public :: xyz_EffMolWtBasicZ(:,:,:)!分子量効果
  real(DP), allocatable, save, public :: xyz_MixRtBasicZPerMolWt(:,:,:)
                              !基本場の混合比 / 分子量 の総和
  real(DP), allocatable, save, public :: xyr_MixRtBasicZPerMolWt(:,:,:)
                              !基本場の混合比 / 分子量 の総和
  real(DP), allocatable, save, public :: xyz_MixRtBasicZ(:,:,:)
                              !基本場の混合比の総和
  real(DP), allocatable, save, public :: xyr_MixRtBasicZ(:,:,:)
                              !基本場の混合比の総和
  real(DP), allocatable, save, public :: xyr_DPTempBZDz(:,:,:)   !温位の鉛直微分
  real(DP), allocatable, save, public :: xyrf_DQMixBZDz(:,:,:,:) !混合比の鉛直微分

  public basicset_init

contains

!!!-----------------------------------------------------------------!!!
  subroutine basicset_init(                                     &
    &   xyz_Press, xyz_Exner, xyz_Temp, xyz_PotTemp, xyz_Dens,    & 
    &   xyz_VelSound, xyza_MixRt, xyz_EffMolWt &
    & )
    !
    != 基本場の値を外部から取得する. 
    !
    ! dry の場合は混合比や乾燥成分と凝結成分の存在比を陽に使わないので,
    ! これらの変数は optional にしている. 
    !
    
    !入力変数
    real(DP), intent(in) :: xyz_Press(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in) :: xyz_Exner(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in) :: xyz_Temp(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in) :: xyz_PotTemp(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in) :: xyz_Dens(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in) :: xyz_VelSound(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in), optional :: xyza_MixRt(imin:imax,jmin:jmax,kmin:kmax,ncmax)
    real(DP), intent(in), optional :: xyz_EffMolWt(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyza_MixRtBasicZPerMolWt(imin:imax,jmin:jmax,kmin:kmax,GasNum)
    integer              :: s, k

    ! 配列の初期化
    call basicset_array_init
    
    !値の入力
    xyz_PressBasicZ    = xyz_Press
    xyz_ExnerBasicZ    = xyz_Exner
    xyz_TempBasicZ     = xyz_Temp
    xyz_PotTempBasicZ    = xyz_PotTemp
    xyz_DensBasicZ     = xyz_Dens
    xyz_VelSoundBasicZ = xyz_VelSound
    xyz_VelSW      = xyz_VelSound  ! 配列名を短くするために. 

    if (present(xyz_EffMolWt)) xyz_EffMolWtBasicZ = xyz_EffMolWt

    if (present(xyza_MixRt)) then 
      xyza_MixRtBasicZ = xyza_MixRt
      xyza_MixRtBasicZPerMolWt = 0.0d0

      do s = 1, GasNum
        xyza_MixRtBasicZPerMolWt(:,:,:,s) = &
          & xyza_MixRtBasicZ(:,:,:,IdxG(s)) / MolWtWet(IdxG(s))
      end do

      xyz_MixRtBasicZPerMolWt = sum(xyza_MixRtBasicZPerMolWt, 4) 
      xyz_MixRtBasicZ         = sum(xyza_MixRtBasicZ, 4) 

      do k = kmin, kmax-1
        xyr_MixRtBasicZPerMolWt(:,:,k) =        &
          & (                               &
          &  + xyz_MixRtBasicZPerMolWt(:,:,k)   &
          &  + xyz_MixRtBasicZPerMolWt(:,:,k+1) &
          & ) * 5.0d-1

        xyr_MixRtBasicZ(:,:,k) =           &
          & (                              &
          &  + xyz_MixRtBasicZ(:,:,k)          &
          &  + xyz_MixRtBasicZ(:,:,k+1)        &
          & ) * 5.0d-1
      end do

!      xyr_MixRtBasicZPerMolWt = xyr_xyz( sum(xyza_MixRtBasicZPerMolWt, 4) )
!      xyr_MixRtBasicZ         = xyr_xyz( sum(xyza_MixRtBasicZ, 4) )     
    end if

    xyz_VPotTempBasicZ = xyz_PotTempBasicZ / xyz_EffMolWtBasicZ
 
    do k = kmin, kmax-1    
      xyr_PotTempBasicZ(:,:,k) = &
        & ( xyz_PotTempBasicZ(:,:,k) + xyz_PotTempBasicZ(:,:,k+1) ) * 5.0d-1
      
      xyr_VPotTempBasicZ(:,:,k) = &
        & ( xyz_VPotTempBasicZ(:,:,k) + xyz_VPotTempBasicZ(:,:,k+1) ) * 5.0d-1
      
      xyr_DensBasicZ(:,:,k) = &
        & ( xyz_DensBasicZ(:,:,k) + xyz_DensBasicZ(:,:,k+1) ) * 5.0d-1
    end do
   
!    xyr_PotTempBasicZ = xyr_xyz( xyz_PotTempBasicZ )
!    xyr_VPotTempBasicZ = xyr_xyz( xyz_VPotTempBasicZ )
!    xyr_DensBasicZ   = xyr_xyz( xyz_Dens )

    ! 水平一様なので, 平均操作は必要ない
    !
    pyz_DensBasicZ   = xyz_Dens
    xqz_DensBasicZ   = xyz_Dens
    pyz_VPotTempBasicZ = xyz_VPotTempBasicZ 
    xqz_VPotTempBasicZ = xyz_VPotTempBasicZ 

    ! 鉛直微分
    !
    xyr_DPTempBZDz = xyr_dz_xyz( xyz_PotTempBasicZ ) 
    xyrf_DQMixBZDz = xyrf_dz_xyzf( xyza_MixRtBasicZ )
  
  end subroutine basicset_init


  subroutine basicset_array_init
    !
    ! *BasicZ な配列の初期化
    !

    allocate( & 
      & xyz_DensBasicZ(imin:imax,jmin:jmax,kmin:kmax), &
      & pyz_DensBasicZ(imin:imax,jmin:jmax,kmin:kmax), &
      & xqz_DensBasicZ(imin:imax,jmin:jmax,kmin:kmax), &
      & xyr_DensBasicZ(imin:imax,jmin:jmax,kmin:kmax), &
      & xyz_PressBasicZ(imin:imax,jmin:jmax,kmin:kmax), &
      & xyz_ExnerBasicZ(imin:imax,jmin:jmax,kmin:kmax), &
      & xyz_TempBasicZ(imin:imax,jmin:jmax,kmin:kmax), &
      & xyz_PotTempBasicZ (imin:imax,jmin:jmax,kmin:kmax), &
      & xyr_PotTempBasicZ (imin:imax,jmin:jmax,kmin:kmax), &
      & xyz_VPotTempBasicZ (imin:imax,jmin:jmax,kmin:kmax), &
      & pyz_VPotTempBasicZ (imin:imax,jmin:jmax,kmin:kmax), &
      & xqz_VPotTempBasicZ (imin:imax,jmin:jmax,kmin:kmax), &
      & xyr_VPotTempBasicZ (imin:imax,jmin:jmax,kmin:kmax), &
      & xyz_VelSoundBasicZ(imin:imax,jmin:jmax,kmin:kmax), &
      & xyz_VelSW(imin:imax,jmin:jmax,kmin:kmax), &
      & xyza_MixRtBasicZ(imin:imax,jmin:jmax,kmin:kmax,ncmax), &
      & xyz_EffMolWtBasicZ(imin:imax,jmin:jmax,kmin:kmax)  &
      & )
    allocate( &
      & xyr_MixRtBasicZPerMolWt(imin:imax,jmin:jmax,kmin:kmax),           &
      & xyz_MixRtBasicZPerMolWt(imin:imax,jmin:jmax,kmin:kmax),           &
      & xyr_MixRtBasicZ(imin:imax,jmin:jmax,kmin:kmax),                   &
      & xyz_MixRtBasicZ(imin:imax,jmin:jmax,kmin:kmax)                    &
      & )
    allocate( xyr_DPTempBZDz(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyrf_DQMixBZDz(imin:imax,jmin:jmax,kmin:kmax,ncmax) )

    ! 値の確定
    xyz_DensBasicZ     = 0.0d0
    xyz_PressBasicZ    = 0.0d0   
    xyz_ExnerBasicZ    = 0.0d0   
    xyz_TempBasicZ     = 0.0d0   
    xyz_PotTempBasicZ    = 0.0d0
    xyz_VPotTempBasicZ   = 0.0d0
    xyz_VelSoundBasicZ = 0.0d0
    xyza_MixRtBasicZ    = 0.0d0
    xyz_EffMolWtBasicZ = 1.0d0 ! dry の場合は必ず 1.0

    xyr_MixRtBasicZPerMolWt = 0.0d0
    xyz_MixRtBasicZPerMolWt = 0.0d0
    xyr_MixRtBasicZ = 0.0d0
    xyz_MixRtBasicZ = 0.0d0

    xyr_DPTempBZDz = 0.0d0
    xyrf_DQMixBZDz = 0.0d0
    
  end subroutine basicset_array_init


  function xyr_dz_xyz(xyz_Var)
    ! 半整数格子点上の 1 階微分を計算する
    
    real(DP),intent(in) :: xyz_Var(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP)            :: xyr_dz_xyz(imin:imax,jmin:jmax,kmin:kmax)
    integer             :: kz
    
    ! 初期化
    ! * 0 割りを防ぐためにはマシンイプシロン値を用いるべき
    !
    xyr_dz_xyz = 0.0d0
    
    ! 1 階微分の計算
    ! * r_dz_z を用いて計算
    !
    do kz = kmin+1, kmax-2
       xyr_dz_xyz(:,:,kz) = &
            &   ((xyz_Var(:,:,kz+1) - xyz_Var(:,:,kz))*9.0d0/(8.0d0*dz)) &
            & - ((xyz_Var(:,:,kz+2) - xyz_Var(:,:,kz-1))/(24.0d0*dz))
    end do
    
  end function xyr_dz_xyz


  function xyrf_dz_xyzf(xyzf_Var)
    ! 半整数格子点上の 1 階微分を計算する
    
    real(DP),intent(in) :: xyzf_Var(imin:imax,jmin:jmax,kmin:kmax,ncmax) 
    real(DP)            :: xyrf_dz_xyzf(imin:imax,jmin:jmax,kmin:kmax,ncmax)
    integer             :: kz
    
    ! 初期化
    ! * 0 割りを防ぐためにはマシンイプシロン値を用いるべき
    !
    xyrf_dz_xyzf = 0.0d0
    
    ! 1 階微分の計算
    ! * r_dz_z を用いて計算
    !
    do kz = kmin+1, kmax-2
       xyrf_dz_xyzf(:,:,kz,:) = &
            &   ((xyzf_Var(:,:,kz+1,:) - xyzf_Var(:,:,kz,:))*9.0d0/(8.0d0*dz)) &
            & - ((xyzf_Var(:,:,kz+2,:) - xyzf_Var(:,:,kz-1,:))/(24.0d0*dz))
    end do
    
  end function xyrf_dz_xyzf
  
end module basicset_3d
