!= ニュートン冷却モジュール
!
! 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 NewtonCoolCrisp1989
  !
  ! ニュートン冷却用モジュール (金星用)
  !
  ! Crisp (1989) Fig.12 が大元であるが, 実際には Takagi and Matsuda
  ! (2005) の Fig.2 から値を読み取った. 簡単のために高度と緩和時間は比
  ! 例しているものとした. 高度 120 km までの緩和時間の鉛直プロファイル
  ! が必要だったため.
  !
  !  高度 (km)      緩和時間 (day)
  !  120              0.1
  !   95              1.0
  !   65             10.0
  !   35            100.0
  !

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

  !暗黙の型宣言禁止
  implicit none

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

contains 
  
!!!------------------------------------------------------------------------!!!
  subroutine NewtonCoolCrisp1989_Init
    !
    ! ニュートン冷却 初期値モジュール
    ! 
    use mpi_wrapper, only: myrank
    use dc_types,   only: DP
    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

    !変数定義
    real(DP), parameter   :: EFTimeLog0 =  2.0d0  ! dlog10( 100.0 ) @ 35 km
    real(DP), parameter   :: EFTimeLog1 = -1.0d0  ! dlog10( 0.1 )   @ 125 Km
    real(DP), parameter   :: AltMin =  35.0d3
    real(DP), parameter   :: AltMax = 125.0d3
    real(DP), parameter   :: DayTime = 86400.0d0
    real(DP)              :: EFTimeLog
    real(DP)              :: DEFTimeLogDZ
    integer               :: k

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

    ! ダンピング係数を決める
    !
    DEFTimeLogDZ = ( EFTimeLog1 - EFTimeLog0 ) / ( AltMax - AltMin )

    do k = kmin, kmax
      EFTimeLog = EFTimeLog0 + DEFTimeLogDZ * ( z_Z(k) - AltMin )
      xyz_Nwtn(:,:,k) = 1.d0 / ( 10.0d0 ** EFTimeLog ) / DayTime

      if (myrank == 0) then 
         write(*,*) k, z_Z(k), EFTimeLog, 10.0d0 ** EFTimeLog, xyz_Nwtn(1,1,k)
      end if
    end do

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


  subroutine NewtonCoolCrisp1989_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 NewtonCoolCrisp1989_forcing
  
end module NewtonCoolCrisp1989
