!= ニュートン冷却モジュール
!
! Authors::   SUGIYAMA Ko-ichiro
! Version::   $Id: damping.f90,v 1.14 2014/11/07 06:46:45 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!

module NewtonCoolGeneral
  !
  ! ニュートン冷却用モジュール (一般的)
  ! 指定された高度領域を指定された時定数でダンピングする. 
  !

  !モジュール読み込み
  use dc_types, only : DP

  !暗黙の型宣言禁止
  implicit none

  !private 属性を指定
  private 
  
  !関数には public 属性を指定
  public NewtonCoolGeneral_Init
  public NewtonCoolGeneral_forcing
  
  !変数定義
  real(DP), save, private, allocatable :: xyz_Nwtn(:,:,:) !ダンピング係数

contains 
  
!!!------------------------------------------------------------------------!!!
  subroutine NewtonCoolGeneral_Init
    !
    ! ニュートン冷却 初期化モジュール
    ! 
    use dc_types,   only: DP
    use dc_iounit,  only: FileOpen
    use dc_message, only: MessageNotify
    use gtool_historyauto, only: HistoryAutoAddVariable
    use gridset, only: imin,       &! x 方向の配列の下限
      &                imax,       &! x 方向の配列の上限
      &                jmin,       &! y 方向の配列の下限
      &                jmax,       &! y 方向の配列の上限
      &                kmin,       &! z 方向の配列の下限
      &                kmax         ! z 方向の配列の上限
    use axesset, only: z_Z          !Z 座標軸(スカラー格子点)
    use namelist_util,  &
      &          only: namelist_filename
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    integer    :: unit ! 装置番号
    integer    :: k
    real(DP)   :: EFTime = 1.0d10   ! e-folding time
    real(DP)   :: AltMin = 0.0d0    ! 高度領域 (下端)  
    real(DP)   :: AltMax = 1.0d10   ! 高度領域 (上端)

    ! NAMELIST から取得
    !
    NAMELIST /newton_cool_nml/ EFTime, AltMin, AltMax

    call FileOpen(unit, file=namelist_filename, mode='r')
    read(unit, NML=newton_cool_nml)
    close(unit)

    ! 配列の用意
    !
    allocate( xyz_Nwtn(imin:imax,jmin:jmax,kmin:kmax) )
    xyz_Nwtn = 0.0d0

    ! 指定された高度を一様に冷却
    !
    do k = kmin, kmax
      if ( AltMin <= z_Z(k) .AND. z_Z(k) <= AltMax ) then 
        xyz_Nwtn(:,:,k) = 1.0d0 / EFTime
      end if
    end do

    !-----------------------------------------------------------------    
    ! 値の確認
    !
    call MessageNotify( "M", "Damping_init", "EFTime = %f", d=(/EFTime/) )
    call MessageNotify( "M", "Damping_init", "AltMin = %f", d=(/AltMin/) )
    call MessageNotify( "M", "Damping_init", "AltMax = %f", d=(/AltMax/) )  

    !-----------------------------------------------------------------    
    ! 出力
    !
    call HistoryAutoAddVariable(                                 &
      & varname='DPTempDtNwtn',                                  &
      & dims=(/'x','y','z','t'/),                                &
      & longname='Newton cooling term of potential temperature', &
      & units='K.s-1',                                           &
      & xtype='double')

  end subroutine NewtonCoolGeneral_Init


  subroutine NewtonCoolGeneral_forcing(   &
    & xyz_PTempBl,                        & !(in) 
    & xyz_DPTempDt )                        !(inout)

    use dc_types, only : DP
    use gtool_historyauto, only: HistoryAutoPut
    use timeset, only: TimeN        ! 現在の時刻
    use gridset, only: imin,       &! x 方向の配列の下限
      &                imax,       &! x 方向の配列の上限
      &                jmin,       &! y 方向の配列の下限
      &                jmax,       &! y 方向の配列の上限
      &                kmin,       &! z 方向の配列の下限
      &                kmax,       &! z 方向の配列の上限
      &                nx,         &! x 方向の物理領域の上限
      &                ny,         &! y 方向の物理領域の上限
      &                nz           ! z 方向の物理領域の上限

    !暗黙の型宣言禁止
    implicit none

    real(DP), intent(in)    :: xyz_PTempBl(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(inout) :: xyz_DPTempDt(imin:imax, jmin:jmax, kmin:kmax)
    real(DP)                :: xyz_Damp(imin:imax, jmin:jmax, kmin:kmax)

    !--------------------------------------------------------
    ! 擾乱場の値をゼロに戻すような強制を与える 
    !
    xyz_Damp = - xyz_Nwtn * xyz_PTempBl 
    xyz_DPTempDt = xyz_DPTempDt + xyz_Damp
    
    call HistoryAutoPut(TimeN, 'DPTempDtNwtn', xyz_Damp(1:nx,1:ny,1:nz))
    
  end subroutine NewtonCoolGeneral_forcing
  
end module NewtonCoolGeneral
