!= 初期値データファイル生成主プログラム
!
!= dcpam-primitive main program for generation of initial data file
!
! Authors::   Yoshiyuki O. Takahashi, Yasuhiro MORIKAWA
!             modified by Shin-ichi Takehiro for primitive model.
! Version::   $Id: init_data_primitive_PST2004.f90,v 1.1.1.1 2010-08-17 05:24:50 takepiro Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

program init_data_primitive
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 初期値データファイルを生成します. 
  !
  ! Initial data file is created. 
  !

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

!!$  ! 初期値データ (リスタートデータ) 提供
!!$  ! Prepare initial data (restart data)
!!$  !
!!$  use initial_data, only: SetInitData

  ! リスタートデータ入出力
  ! Restart data input/output
  !
  use restart_file_io, only: RestartFileOutPut

  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 倍精度実数型. Double precision. 
    &                 STRING, &  ! 文字列.       Strings. 
    &                 TOKEN      ! キーワード.   Keywords. 

  ! 宣言文 ; Declaration statements
  !
  implicit none

  ! 予報変数 (ステップ $ t-\Delta t $ , $ t $ , $ t+\Delta t $ )
  ! Prediction variables  (Step $ t-\Delta t $ , $ t $ , $ t+\Delta t $ )
  !
  real(DP), allocatable:: xyz_U (:,:,:)
                              ! $ u $ .   東西風速. Eastward wind
  real(DP), allocatable:: xyz_V (:,:,:)
                              ! $ v $ .   南北風速. Northward wind
  real(DP), allocatable:: xyz_Temp (:,:,:)
                              ! $ T $ .   温度. Temperature
  real(DP), allocatable:: xy_Ps (:,:)
                              ! $ p_s $ . 地表面気圧. Surface pressure

  ! 作業変数
  ! Work variables
  !

  ! 実行文 ; Executable statement
  !

  ! 主プログラムの初期化 (内部サブルーチン)
  ! Initialization for the main program (Internal subroutine)
  !
  call MainInit

  ! 初期値データの作成
  ! Generate initial data
  !
  call SetInitData( &
    & xyz_U, xyz_V, xyz_Temp, xy_Ps )  ! (out)

  ! 初期値データの出力
  ! Output initial data
  !
  call RestartFileOutput( &
    & xyz_U, xyz_V, xyz_Temp, xy_Ps )  ! (in)

  ! 主プログラムの終了処理 (内部サブルーチン)
  ! Termination for the main program (Internal subroutine)
  !
  call MainTerminate

contains

  !-------------------------------------------------------------------

  subroutine MainInit
    !
    ! 主プログラムの初期化手続き. 
    !
    ! Initialization procedure for the main program. 
    !

    ! MPI
    !
    use mpi_wrapper, only : MPIWrapperInit

    use dc_message, only: MessageNotify

    ! コマンドライン引数処理
    ! Command line option parser
    !
    use option_parser, only: OptParseInit

    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: NmlutilInit

    ! 時刻管理
    ! Time control
    !
    use timeset, only: TimesetInit
                              ! ステップ $ t $ の時刻. Time of step $ t $. 

    ! 出力ファイルの基本情報管理
    ! Management basic information for output files
    ! 
    use fileset, only: FilesetInit

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: GridsetInit, &
      &                imax, & ! 経度格子点数. 
                               ! Number of grid points in longitude
      &                jmax, & ! 緯度格子点数. 
                               ! Number of grid points in latitude
      &                kmax    ! 鉛直層数. 
                               ! Number of vertical level

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: ConstantsInit

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: AxessetInit

    ! リスタートデータ入出力
    ! Restart data input/output
    !
    use restart_file_io, only: RestartFileOpen, RestartFileGet

    ! 宣言文 ; Declaration statements
    !
    implicit none
    character(*), parameter:: prog_name = 'init_data'
                              ! 主プログラム名. 
                              ! Main program name
    character(*), parameter:: version = &
      & '$Name:  $' // &
      & '$Id: init_data_primitive_PST2004.f90,v 1.1.1.1 2010-08-17 05:24:50 takepiro Exp $'
                              ! 主プログラムのバージョン
                              ! Main program version
    character(STRING):: brief
                              ! 実行ファイルの簡潔な説明
                              ! Brief account of executable file

    ! 実行文 ; Executable statement
    !

    ! Initialize MPI
    !
    call MPIWrapperInit


    brief = 'Initial data generation'
    call MessageNotify( 'M', prog_name, 'Run: %c', c1 = trim(brief) )
    call MessageNotify( 'M', prog_name, '-- version = %c', c1 = trim(version) )

    ! コマンドライン引数処理
    ! Command line option parser
    !
    call OptParseInit(prog_name, brief)

    ! NAMELIST ファイル名入力
    ! Input NAMELIST file name
    !
    call NmlutilInit

    ! 時刻管理
    ! Time control
    !
    call TimesetInit

    ! 出力ファイルの基本情報管理
    ! Management basic information for output files
    ! 
    call FilesetInit

    ! 格子点設定
    ! Grid points settings
    !
    call GridsetInit

    ! 物理定数設定
    ! Physical constants settings
    !
    call ConstantsInit

    ! 座標データ設定
    ! Axes data settings
    !
    call AxessetInit

    ! 予報変数の割付
    ! Allocation of prediction variables
    !
    allocate( xyz_U    (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_V    (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_Temp (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xy_Ps    (0:imax-1, 1:jmax) )

    ! 初期値データ出力
    ! Initial data output
    !
    call RestartFileOpen( flag_init_data = .true. ) ! (in) optional

  end subroutine MainInit

  !-------------------------------------------------------------------

  subroutine MainTerminate
    !
    ! 主プログラムの終了処理手続き. 
    !
    ! Termination procedure for the main program. 
    !

    ! MPI
    !
    use mpi_wrapper, only : MPIWrapperFinalize

    ! 時刻管理
    ! Time control
    !
    use timeset, only: TimesetClose

    ! リスタートデータ入出力
    ! Restart data input/output
    !
    use restart_file_io, only: RestartFileClose

    ! 宣言文 ; Declaration statements
    !
    implicit none

    ! 実行文 ; Executable statement
    !

    ! リスタートデータファイルクローズ
    ! Close restart data input
    !
    call RestartFileClose

    ! 時刻管理終了処理
    ! Termination of time control
    !
    call TimesetClose

    ! 予報変数の割付解除
    ! Deallocation of prediction variables
    !
    deallocate( xyz_U    )
    deallocate( xyz_V    )
    deallocate( xyz_Temp )
    deallocate( xy_Ps    )

    ! Finalize MPI
    !
    call MPIWrapperFinalize


  end subroutine MainTerminate

  !-------------------------------------------------------------------

  subroutine SetInitData( &
    & xyz_U, xyz_V, xyz_Temp, xy_Ps )  ! (out)
    !
    ! 初期データの設定
    !
    ! Setting initial values
    !

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: &
      &                imax, & ! 経度格子点数. 
                               ! Number of grid points in longitude
      &                jmax, & ! 緯度格子点数. 
                               ! Number of grid points in latitude
      &                kmax    ! 鉛直層数. 
                               ! Number of vertical level

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: PI, Omega, RPlanet, GasRDry

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: x_Lon, y_Lat, z_Sigma, y_Lat_Weight

    use wa_module, only : a_AvrLat_ya

    real(DP), intent(OUT) :: xyz_U (0:,:,:)
                                    ! $ u $ .   東西風速. Eastward wind
    real(DP), intent(OUT) :: xyz_V (0:,:,:)
                                    ! $ v $ .   南北風速. Northward wind
    real(DP), intent(OUT) :: xyz_Temp (0:,:,:)
                                    ! $ T $ .   温度. Temperature
    real(DP), intent(OUT) :: xy_Ps (0:,:)     
                                    ! $ p_s $ . 地表面気圧. Surface pressure

    integer, parameter :: nmax=100
    integer   :: i, j, k, n
    real(DP)  :: lat, dlat

    real(DP)  :: z_Z(1:kmax)
    real(DP)  :: z_F(1:kmax)
    real(DP)  :: z_dFdZ(1:kmax)
    real(DP)  :: z_Temp(1:kmax)

    real(DP)  :: yz_Temp(1:jmax,1:kmax)

    real(DP), parameter :: u0  = 50.0d0
    real(DP), parameter :: H   = 7.34d3
    real(DP), parameter :: z0  = 22.0d3
    real(DP), parameter :: z1  = 30.0d3
    real(DP), parameter :: dz0 = 5.0d3
    real(DP), parameter :: p0  = 1.0d5

    real(DP), parameter :: That    = 1.0d0
    real(DP), parameter :: lambda0 = 0.0d0
    real(DP), parameter :: phi0    = PI/4
    real(DP), parameter :: alpha   = 1.0d0/3.0d0
    real(DP), parameter :: beta    = 1.0d0/6.0d0
    
    xyz_V    = 0.0D0
    xy_Ps    = p0

    z_Z    = - H*log(z_Sigma)
    z_F    = 0.5D0 * ( 1 - tanh((z_Z-z0)/dz0)**3 ) * sin(PI*z_Z/z1)
    !
    !  zonal velocity profile
    ! 
    do k=1,kmax
       do j=1,jmax
          do i=0,imax-1
             if ( y_Lat(j) > 0.0d0 ) then
                xyz_U(i,j,k) = u0 * sin(PI*sin(y_Lat(j))**2)**3 * z_F(k)
             else
                xyz_U(i,j,k) = 0.0d0
             endif
          end do
       end do
    end do

    !
    !  temprature profile
    ! 
    z_dFdZ =-1.5d0 * tanh((z_Z-z0)/dz0)**2 / cosh((z_Z-z0)/dz0)**2 / dz0 &
         &   * sin(PI*z_Z/z1) &
         & + 0.5D0 * (1 - tanh((z_Z-z0)/dz0)**3) &
         &   * cos(PI*z_Z/z1) * (PI/z1)

    do k=1,kmax
       yz_Temp(1,k) = 0.0D0
       do j=1,jmax-1
          yz_Temp(j+1,k) = yz_Temp(j,k)
          dlat =(y_Lat(j+1)-y_Lat(j))/(nmax+1)
          do n=0,nmax
             lat = y_Lat(j) + dlat * (n+0.5)/(nmax+1)
             if ( lat > 0.0d0 ) then
                yz_Temp(j+1,k) = yz_Temp(j+1,k) &
                     -H/GasRDry                 &
                     * (   RPlanet * 2*Omega*sin(lat) &
                         + 2 * u0*sin(PI*sin(lat)**2)**3*z_F(k) * tan(lat) ) &
                     * u0*sin(PI*sin(lat)**2)**3*z_dFdz(k) * dlat
             else
                yz_Temp(j+1,k) = yz_Temp(j+1,k)
             endif
          enddo
       enddo
    enddo

    z_Temp = a_AvrLat_ya(yz_Temp)

    do k=1,kmax
       do j=1,jmax
          do i=0,imax-1
             xyz_Temp(i,j,k) = yz_Temp(j,k)-z_Temp(k) + TempUSStandard(z_Z(k))
          enddo
       end do
    end do

    !
    ! adding temperature disturbance
    !
    do k=1,kmax
       do j=1,jmax
          do i=0,imax-1
             if ( x_Lon(i) <= PI ) then
                xyz_Temp(i,j,k) = xyz_Temp(i,j,k) &
                     + That * cosh((x_Lon(i)-lambda0)/alpha)**(-2) &
                            * cosh((y_Lat(j)-phi0)/beta)**(-2) 
             else
                xyz_Temp(i,j,k) = xyz_Temp(i,j,k) &
                     + That * cosh((x_Lon(i)-lambda0-2*PI)/alpha)**(-2) &
                            * cosh((y_Lat(j)-phi0)/beta)**(-2) 
             endif
          enddo
       end do
    end do

  end subroutine SetInitData

  function TempUSStandard(z)
    !
    ! アメリカ標準大気(COESA 1976)温度の計算
    !
    ! Compute U.S. Standard Atmosphere temperature
    ! 
    use dc_message, only: MessageNotify

    real(DP), intent(IN) :: z               ! $ z $ . 高度  Height
    real(DP)             :: TempUSStandard  
                         ! 米標準大気温度 U.S. Standard Atmosphere temperature

    if ( z >= 0  .AND. z < 11.0d3 ) then
       TempUSStandard = 288.15 - 6.5d-3*z 
    else if ( z >= 11.0d3 .AND. z < 20.0d3 ) then
       TempUSStandard = 288.15 - 6.5d-3*11.0d3
    else if ( z >= 20.0d3 .AND. z < 32.0d3 ) then
       TempUSStandard = 288.15 - 6.5d-3*11.0d3 + 1.0d-3*(z-20.0d3)
    else if ( z >= 32.0d3 .AND. z < 47.0d3 ) then
       TempUSStandard = 288.15 - 6.5d-3*11.0d3 + 1.0d-3*(32.0d3-20.0d3) &
            & + 2.8d-3 * (z-32.0d3)
    else if ( z >= 47.0d3 .AND. z < 51.0d3 ) then
       TempUSStandard = 288.15 - 6.5d-3*11.0d3 + 1.0d-3*(32.0d3-20.0d3) &
            & + 2.8d-3 * (47.0d3-32.0d3)
    else if ( z >= 51.0d3 .AND. z < 71.0d3 ) then
       TempUSStandard = 288.15 - 6.5d-3*11.0d3 + 1.0d-3*(32.0d3-20.0d3) &
            & + 2.8d-3 * (47.0d3-32.0d3) - 2.8d3 * (z-51.0d3)
    else if ( z >= 71.0d3 .AND. z < 80.0d3 ) then
       TempUSStandard = 288.15 - 6.5d-3*11.0d3 + 1.0d-3*(32.0d3-20.0d3) &
            & + 2.8d-3 * (47.0d3-32.0d3) - 2.8d3 * (71.0d3-51.0d3) &
            & - 2.0d-3 * (z-71.0d3)
    else if ( z >= 80.0d3 ) then
       TempUSStandard = 288.15 - 6.5d-3*11.0d3 + 1.0d-3*(32.0d3-20.0d3) &
            & + 2.8d-3 * (47.0d3-32.0d3) - 2.8d3 * (71.0d3-51.0d3) &
            & - 2.0d-3 * (80.0d3-71.0d3)
    else
       call MessageNotify('E','TempUSStandard','z must be larger than zero')
    endif
  end function TempUSStandard

end program init_data_primitive
