!= Module BasicEnvInit
!
! Authors::   SUGIYAMA Koichiro, ODAKA Masatsugu
! Version::   $Id: initialdata_yamasaki1983.f90,v 1.1 2011-04-26 04:08:10 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview 
!
!デフォルトの基本場を設定するための変数参照型モジュール
!   * BasicEnvFile_init: 基本場の値を netCDF ファイルから取得
!   * BasicEnvCalc_Init: 基本場の情報を Namelist から取得して値を計算
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!
module initialdata_yamasaki1983
  !
  !デフォルトの基本場を設定するためのサブルーチン. 
  !基本場を計算し, BasicSet モジュールの値を初期化する. 
  !
  !コンパイルの順序の問題から, 基本場の値(hogeBasicZ な変数)を
  !計算する部分をBasicSet モジュールから切り離している. 
  !ECCM 始め, BasicSet 自体に依存するが hogeBasicZ は use しない
  !外部サブルーチンを利用するためである. 
  !

  !モジュール読み込み
  use dc_types,   only: STRING, DP
  use dc_message, only: MessageNotify
  use gridset,  only: kmin,          &!配列の Z 方向の下限
    &                 kmax,          &!配列の Z 方向の上限
    &                 nz,            &!配列の Z 方向の数
    &                 na,            &!凝縮成分の数
    &                 z_Z            !スカラー格子点での高度

  use basicset, only: BasicSetArray_Init, &!
    &                 PressBasis,    &!温位の基準圧力
    &                 GasRDry,       &!乾燥成分の定圧比熱
    &                 CpDry,         &!乾燥成分の定圧比熱
    &                 TempSfc,       &!地表面温度
    &                 PressSfc,      &!地表面圧力
    &                 GasRDry         !気体定数
  use chemcalc, only: SvapPress       !飽和蒸気圧

    
  !暗黙の型宣言禁止
  implicit none

  !デフォルトは private
  private

  real(DP), save  :: DryHeight = 0.0d0

  !初期化だけ公開
  public basicenv_yamasaki1983_set

contains

!!------------------------------------------------------------------------------!!!
  subroutine initialdata_yamasaki1983_set( z_Temp, z_Press, za_MolFr )
    !
    !== 概要
    !  * Yamasaki, 1983 の温度と相対湿度の観測値で基本場を作成する
    !    * 温度の基本場
    !      * 観測データをNetCDFファイル化したものから値を読み込む
    !        * 読み込んだ値を線形補間して温度の基本場を作成する
    !    * 湿度の基本場
    !      * 相対湿度は subroutine HUM で作成済み
    !        * ここでは相対湿度をモル比に変換して湿度の基本場を作成する 
    !    * 気圧の基本場
    !      * 湿潤成分と乾燥成分の分子量差を考慮した静水圧平衡から計算する
    !

    implicit none
    
    real(DP), intent(out):: z_Press(kmin:kmax)!圧力
    real(DP), intent(out):: z_Temp(kmin:kmax) !温度
    real(DP), intent(out):: za_MolFr(kmin:kmax, 1:na) !モル比
    real(DP)             :: Ob_Alt(1:89)            !観測データに対応する高度
    real(DP)             :: Ob_Hum(1:88)            !相対湿度の観測値(1次元)
    real(DP)             :: Ob_Temp(1:89)           !NetCDF から読み込んだ気温の観測値(1次元)
    real(DP)             :: z_Hum(kmin:kmax)        !相対湿度の基本場の値(%)
    integer              :: k, s, i, m, b

    integer :: mmin = 1
    integer :: mmax = 89
    integer :: ID = 6    ! 水の物質番号


    ! 観測データ. 
    ! なお, Yamsaki(1983)にある元データはデータ間隔がばらばらであったので
    ! 250m間隔に直したものを使用している (by 今関)
    !
    do m = mmin, mmax
      Ob_Alt(m) = 25.0d1*(m-1)
    end do

    Ob_Temp = (/ &
      & 300.00, 298.00, 296.00, 294.22, 293.15, 291.08, 289.80, 288.55, 287.30, 286.50, &
      & 284.78, 283.49, 282.16, 280.81, 279.45, 278.07, 276.70, 275.33, 273.95, 272.58, &
      & 271.20, 269.83, 268.45, 267.06, 265.63, 264.19, 262.70, 261.20, 259.68, 258.11, &
      & 256.55, 254.86, 253.18, 251.41, 249.60, 247.76, 245.83, 243.89, 241.90, 240.30, &
      & 237.90, 235.90, 233.90, 231.90, 229.90, 227.90, 225.90, 223.91, 221.97, 220.03, &
      & 218.09, 216.16, 214.24, 212.42, 210.64, 208.86, 207.37, 205.96, 204.55, 203.35, &
      & 202.48, 201.60, 200.73, 200.36, 200.12, 199.89, 199.65, 199.91, 200.30, 200.70, &
      & 201.09, 201.55, 202.30, 203.05, 203.80, 204.55, 205.30, 206.05, 206.80, 207.55, &
      & 208.30, 209.05, 209.69, 210.31, 210.93, 211.54, 212.16, 212.77, 212.99 /)
    
    Ob_Hum = (/ &
      & 83.00, 88.85, 89.70, 91.70, 93.70, 94.35, 95.00, 94.35, 93.70, 91.85, 90.00, &
      & 87.50, 85.00, 82.90, 80.80, 78.80, 76.70, 74.35, 72.00, 71.00, 70.00, 68.85, &
      & 67.70, 66.85, 66.00, 65.60, 65.20, 64.75, 64.30, 63.90, 63.50, 63.10, 62.70, &
      & 62.25, 61.80, 61.40, 61.00, 61.00, 61.00, 61.35, 61.70, 61.85, 62.00, 62.15, &
      & 62.30, 62.65, 63.00, 63.00, 63.00, 63.35, 63.70, 64.50, 65.30, 66.60, 67.90, &
      & 70.10, 72.30, 73.30, 74.30, 74.60, 74.90, 68.50, 62.10, 54.05, 46.00, 41.45, &
      & 36.90, 33.25, 29.60, 24.65, 19.70, 15.75, 11.80, 9.00, 6.20, 5.25, 4.30, 3.35, &
      & 2.40, 3.35, 4.30, 2.95, 1.60, 1.05, 0.50, 0.45, 0.40, 0.35 /)


    ! 線形補完
    !
    ! Ob_Hum を線形補間して z_Hum を作成する
    ! Ob_Alt(m) < z_Z(k) < Ob_Alt(m+1),
    ! Ob_Alt = z_Z(k),
    ! Ob_Alt < z_Z(k)
    ! の3つで場合分け
    !
    do k = 1, nz
      do m = mmin, mmax
        ! z_Z が Ob_altitude のある区間に挟まれるとき,
        ! その区間の Obaltitude(m), Obaltitude(m+1) を結ぶ
        ! 直線で Ob_TempZ を線形補間する
        if (Ob_Alt(m) /=  z_Z(k) .AND. z_Z(k) > Ob_Alt(m) .AND. z_Z(k) < Ob_Alt(m+1)) then

          z_Temp(k) = Ob_Temp(m) &
            & + ( ( Ob_Temp(m+1) - Ob_Temp(m) ) / ( Ob_Alt(m+1) - Ob_Alt(m) ) ) &
            &   * ( z_Z(k) - Ob_Alt(m) )

          z_Hum(k) = Ob_Hum(m) &
            & + ( ( Ob_Hum(m+1) - Ob_Hum(m) ) / ( Ob_Alt(m+1) - Ob_Alt(m) ) ) &
            &   * (z_Z(k) - Ob_Alt(m))
          
        ! z_Z(k) = Ob_Alt(m) の時は, Ob_TempZ(m) = z_Temp(k) である
        else if (Ob_Alt(m) == z_Z(k)) then
          z_Temp(k)  = Ob_Temp(m)
          z_Hum(k) = Ob_Hum(m)

        ! z_Z(k) > Ob_altitude では観測データが無いので等温大気にする          
        else if (Ob_Alt(m).lt.z_Z(k)) then
          z_Temp(k)  = z_Temp(k-1)
          z_Hum(k) = z_Hum(k-1)
        end if
      end do
    end do

    ! 初期化
    z_Press = 1.0d-60

    ! 地表面気圧と温度から大気最下層の気圧を求める
    ! 静水圧の式 dP/dz = - \rho * g を使用
    z_Press(RegZMin+1)  = &
      & PressSfc *((TempSfc / z_Temp(RegZMin+1)) ** (- CpDry /  GasRDry))

    ! 大気最下層のモル比を計算する
    za_MolFr(1,1) = &
      & SvapPress( ID, z_Temp(1)) * ( z_Hum(1) / 100.0d0 ) / z_Press(1)

    ! 大気最下層の値を元に, 気圧と湿度の基本場を求める
    ! 気圧は湿潤成分と乾燥成分の分子量差を考慮した静水圧平衡から計算する
    ! 相対湿度はモル比に変換する
    ! 以下のif文のDryHeightは基本場に乾燥した層を入れる場合に必要な値
    ! NAMELISTで鉛直領域以上の値を指定しておけば乾燥した層は入らない

    ! 湿度の基本場の計算
    do k = 1, nz
      za_MolFr(k,1) = SvapPress( ID, z_Temp(k)) * (z_Hum(k) / 100 ) / z_Press(k)
      
      z_Press(k+1) =                                                  &
        &  z_Press(k) * ((z_Temp(k) / z_Temp(k+1)) ** (- CpDry / GasRDry)) 
    end do

  end subroutine INITIALDATA_YAMASAKI1983_SET
  
end module Initialdata_yamasaki1983


