!= 初期値生成プログラムのサンプル
!
!= Initial data generation sample program
!
! Authors::   Yasuhiro MORIKAWA
! Version::   $Id: init_sample.f90,v 1.1.1.1 2007/06/21 03:23:58 morikawa Exp $
! Tag Name::  $Name: dcpam4-20070730-2 $
! Copyright:: Copyright (C) GFD Dennou Club, 2007. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!

program init_sample
  !
  ! 初期値生成のための実行プログラムのサンプルです.
  !
  ! This is sample executable program for initial data generation.
  !

  !---------------------------------------------------------
  !  初期値生成
  !  Generate initial data
  !---------------------------------------------------------
  use initial_data, only: INIDAT, Create, Get, Close, PutLine

  !---------------------------------------------------------
  !  座標変数設定
  !  Configure coordinate variables
  !---------------------------------------------------------

  !-------------------------------------
  !  力学過程
  !  Dynamical core
  use dyn_spectral_as83, only: DYNSPAS83, Create, Close, GetAxes

  !-------------------------------------
  !  物理定数
  !  Physical constants
  use constants, only: CONST, Create

  !---------------------------------------------------------
  !  データ I/O
  !  Data I/O
  !---------------------------------------------------------
  use gt4_history, only: GT_HISTORY, &
    & HistoryCreate, HistoryAddVariable, HistoryPut, HistoryClose, &
    & HistoryAddAttr

  !---------------------------------------------------------
  !  汎用ユーティリティ
  !  Common utilities
  !---------------------------------------------------------
  use dc_types, only: DP, STRING, TOKEN
  use dc_args, only: ARGS, Open, Debug, Help, Strict, Close
  use dc_trace, only: DbgMessage, BeginSub, EndSub
  use dc_message,only: MessageNotify
  use dc_string, only: StoA
  use dc_clock, only: CLOCK, Create, Close, Start, Stop, Result, &
    & Predict, operator(+)
  implicit none


  !-------------------------------------------------------------------
  !  格子点数・最大全波数
  !  Grid points and maximum truncated wavenumber
  !-------------------------------------------------------------------
  integer:: nmax ! 最大全波数. 
                 ! Maximum truncated wavenumber
  integer:: imax ! 経度格子点数. 
                 ! Number of grid points in longitude
  integer:: jmax ! 緯度格子点数. 
                 ! Number of grid points in latitude
  integer:: kmax ! 鉛直層数. 
                 ! Number of vertical level

  character(TOKEN):: init_nc = 'init.nc'
                              ! 初期値データ netCDF ファイル
                              ! NetCDF file for initial data
  character(TOKEN):: init_time_unit = 'minutes'
  real:: init_timeB = -20.0
                              ! 初期値データ ( $ t-\Delta t $ ) の時刻
                              ! Time of initial data ( $ t-\Delta t $ )
  real:: init_timeN = 0.0
                              ! 初期値データ ( $ t $ ) の時刻
                              ! Time of initial data ( $ t $ )

  !---------------------------------------------------------
  !  物理定数
  !  Physical constants
  !---------------------------------------------------------
  real(DP):: PI               ! $ \pi $ .    円周率.         Circular constant


  !---------------------------------------------------------
  !  配列の定義
  !  Declaration of array
  !---------------------------------------------------------

  !-------------------------------------
  !  座標変数
  !  Coordinate variables
  real(DP), allocatable:: x_Lon (:) ! 経度. Longitude
  real(DP), allocatable:: y_Lat (:) ! 緯度. Latitude
  real(DP), allocatable:: z_Sigma (:)
                              ! $ \sigma $ レベル (整数). 
                              ! Full $ \sigma $ level
  real(DP), allocatable:: r_Sigma (:)
                              ! $ \sigma $ レベル (半整数). 
                              ! Half $ \sigma $ level
  real(DP), allocatable:: z_DelSigma (:)
                              ! $ \Delta \sigma $ (整数). 
                              ! $ \Delta \sigma $ (Full)

  !-------------------------------------
  !  予報変数
  !  Prediction variables
  real(DP), allocatable:: xyz_VorB (:,:,:)
                              ! $ \zeta (t-\Delta t) $ . 渦度. Vorticity
  real(DP), allocatable:: xyz_DivB (:,:,:)
                              ! $ D (t-\Delta t) $ .     発散. Divergence
  real(DP), allocatable:: xyz_TempB (:,:,:)
                              ! $ T (t-\Delta t) $ .     温度. Temperature
  real(DP), allocatable:: xyz_QVapB (:,:,:)
                              ! $ q (t-\Delta t) $ .     比湿. Specific humidity
  real(DP), allocatable:: xy_PsB (:,:)
                              ! $ P_s (t-\Delta t) $ .   地表面気圧. Surface pressure

  real(DP), allocatable:: xyz_VorN (:,:,:)
                              ! $ \zeta (t) $ . 渦度. Vorticity
  real(DP), allocatable:: xyz_DivN (:,:,:)
                              ! $ D (t) $ .     発散. Divergence
  real(DP), allocatable:: xyz_TempN (:,:,:)
                              ! $ T (t) $ .     温度. Temperature
  real(DP), allocatable:: xyz_QVapN (:,:,:)
                              ! $ q (t) $ .     比湿. Specific humidity
  real(DP), allocatable:: xy_PsN (:,:)
                              ! $ P_s (t) $ .   地表面気圧. Surface pressure

  !---------------------------------------------------------
  !  作業変数
  !  Work variables
  !---------------------------------------------------------
  type(ARGS):: arg            ! コマンドライン引数. 
                              ! Command line arguments
  type(CONST):: const_earth   ! 物理定数. Physical constants.
  type(DYNSPAS83):: dyn       ! 力学過程 (座標変数設定用)
                              ! Dynamical core (For configuration of coordinate variables)
  type(INIDAT):: ini_dat      ! 初期値データ生成
                              ! Generation of initial data
  type(GT_HISTORY):: gthist_init
                              ! 初期値データ出力. 
                              ! Output of initial data
  type(CLOCK):: clk_setup, clk_histput
                              ! CPU 時間モニター. 
                              ! CPU time monitor

  !-------------------------------------
  !  データ出力用メタデータ
  !  Meta data for output data
  character(*), parameter:: title = 'DCPAM initial data'
  character(*), parameter:: source = 'dcpam4 (see http://www.gfd-dennou.org/library/dcpam)'
  character(*), parameter:: institution = 'GFD Dennou Club'

  character(*), parameter:: version = &
    & '$Name: dcpam4-20070730-2 $' // &
    & '$Id: init_sample.f90,v 1.1.1.1 2007/06/21 03:23:58 morikawa Exp $'
  character(STRING),  parameter:: subname = "init"

continue
  !-------------------------------------------------------------------
  !  コマンドライン引数の取得
  !  Get command line arguments
  !-------------------------------------------------------------------
  call Open(arg)
  call Debug(arg); call Help(arg) ; call Strict(arg)
  call Close(arg)

  call BeginSub(subname, version=version)

  !-------------------------------------------------------------------
  !  CPU 時間モニターの初期設定
  !  Configure the settings for CPU time monitor
  !-------------------------------------------------------------------
  call Create( &
    &         clk_setup, &    ! (out)
    &         'Setup')        ! (in)
  call Create( &
    &         clk_histput, &  ! (out)
    &         'HistoryPut')   ! (in)

  call Start(clk_setup) ! (inout)

  !-------------------------------------------------------------------
  !  格子点数の設定
  !  Configure the grid points and max wave number
  !-------------------------------------------------------------------
  imax = 64
  jmax = 32
  kmax = 12
  nmax = 21

  !-------------------------------------------------------------------
  !  初期値データ出力の設定
  !  Configure the settings for initial data generation
  !-------------------------------------------------------------------
  call Create( ini_dat, &       ! (inout)
    &  imax, jmax, kmax )       ! (in)

  !-------------------------------------------------------------------
  !  物理定数の設定
  !  Configure the physical constants
  !-------------------------------------------------------------------
  call Create(const_earth)
  PI = const_earth % PI

  !-------------------------------------------------------------------
  !  力学過程の設定 (緯度経度データ, 鉛直レベル設定のため)
  !  Configure the settings for dynamical core
  !  (for configuration of the data of latitude and longitude 
  !  and vertical level)
  !-------------------------------------------------------------------
  call Create( dyn_sp_as = dyn, &                                ! (inout)
    &  nmax = nmax, imax = imax, jmax = jmax, kmax = kmax, &     ! (in)
    &  PI = const_earth % PI, RPlanet = const_earth % RPlanet, & ! (in)
    &  Grav = const_earth % Grav, Omega = const_earth % Omega, & ! (in)
    &  Cp = const_earth % Cp, RAir = const_earth % RAir, &       ! (in)
    &  EpsVT = const_earth % EpsVT, &                            ! (in)
    &  VisOrder = const_earth % VisOrder, &                      ! (in)
    &  EFoldTime = const_earth % EFoldTime, &                    ! (in)
    &  DelTime = real(init_timeN - init_timeB, DP) )             ! (in)

  !-------------------------------------------------------------------
  !  緯度経度データ, 鉛直レベルの設定
  !  (リスタートファイル, ヒストリファイル出力用)
  !  Configure the data of latitude and longitude and vertical level
  !  for output of restart file and history files
  !-------------------------------------------------------------------
  allocate( x_Lon(0:imax-1) )
  allocate( y_Lat(0:jmax-1) )
  allocate( z_Sigma(0:kmax-1) )
  allocate( r_Sigma(0:kmax) )
  allocate( z_DelSigma(0:kmax-1) )

  call GetAxes( dyn_sp_as = dyn, &                 ! (inout)
    & x_Lon = x_Lon,     y_Lat = y_Lat, &          ! (out)
    & z_Sigma = z_Sigma, r_Sigma = r_Sigma, &      ! (out)
    & z_DelSigma = z_DelSigma )                    ! (out)

  !-------------------------------------------------------------------
  !  予報変数の割付
  !  Allocations of prediction variables
  !-------------------------------------------------------------------
  allocate( xyz_VorB(0:imax-1, 0:jmax-1, 0:kmax-1) )
  allocate( xyz_DivB(0:imax-1, 0:jmax-1, 0:kmax-1) )
  allocate( xyz_TempB(0:imax-1, 0:jmax-1, 0:kmax-1) )
  allocate( xyz_QVapB(0:imax-1, 0:jmax-1, 0:kmax-1) )
  allocate( xy_PsB(0:imax-1, 0:jmax-1) )

  allocate( xyz_VorN(0:imax-1, 0:jmax-1, 0:kmax-1) )
  allocate( xyz_DivN(0:imax-1, 0:jmax-1, 0:kmax-1) )
  allocate( xyz_TempN(0:imax-1, 0:jmax-1, 0:kmax-1) )
  allocate( xyz_QVapN(0:imax-1, 0:jmax-1, 0:kmax-1) )
  allocate( xy_PsN(0:imax-1, 0:jmax-1) )

  !-------------------------------------------------------------------
  !  初期値データの取得
  !  Get initial data
  !-------------------------------------------------------------------
  call Get( ini_dat, &                                    ! (inout)
    & xyz_VorB, xyz_DivB, xyz_TempB, xyz_QVapB, xy_PsB, & ! (out)
    & xyz_VorN, xyz_DivN, xyz_TempN, xyz_QVapN, xy_PsN )  ! (out)

  call Stop(clk_setup)   ! (inout)

  !-------------------------------------------------------------------
  !  初期値ファイルへのデータ出力設定
  !  Configure the settings for initial data output
  !-------------------------------------------------------------------
  call Start(clk_histput)   ! (inout)

  call HistoryCreate( &
    & file = init_nc, title = title, &                    ! (in)
    & source = source, institution = institution, &       ! (in)
    & dims = StoA('lon', 'lat', 'sig', 'sigm', 'time'), & ! (in)
    & dimsizes = (/imax, jmax, kmax, kmax + 1, 0/), &     ! (in)
    & longnames = StoA('longitude', 'latitude', &
    &                  'sigma at layer midpoints', &
    &                  'sigma at layer end-points (half level)', &
    &                  'time'), &                        ! (in)
    & units = StoA('degree_east', 'degree_north', &
    &              '1', '1', init_time_unit), &          ! (in)
    & origin = init_timeB, &                             ! (in)
    & interval = init_timeN - init_timeB, &              ! (in)
    & history = gthist_init)                             ! (out)

  call HistoryPut( &
    & varname = 'lon', array = x_Lon / PI * 180.0_DP, & ! (in)
    & history = gthist_init)            ! (inout)
  call HistoryPut( &
    & varname = 'lat', array = y_Lat / PI * 180.0_DP, & ! (in)
    & history = gthist_init)            ! (inout)
  call HistoryPut( &
    & varname = 'sig', array = z_Sigma, & ! (in)
    & history = gthist_init)              ! (inout)
  call HistoryPut( &
    & varname = 'sigm', array = r_Sigma, & ! (in)
    & history = gthist_init)               ! (inout)

  call HistoryAddAttr( &
    & varname = 'lon', attrname = 'standard_name', & ! (in)
    & value = 'longitude', &                         ! (in)
    & history = gthist_init)                         ! (inout)
  call HistoryAddAttr( &
    & varname = 'lat', attrname = 'standard_name', & ! (in)
    & value = 'latitude', &                          ! (in)
    & history = gthist_init)                         ! (inout)
  call HistoryAddAttr( &
    & varname = 'sig', attrname = 'standard_name', &   ! (in)
    & value = 'atmosphere_sigma_coordinate', &         ! (in)
    & history = gthist_init)                           ! (inout)
  call HistoryAddAttr( &
    & varname = 'sigm', attrname = 'standard_name', &  ! (in)
    & value = 'atmosphere_sigma_coordinate', &         ! (in)
    & history = gthist_init)                           ! (inout)
  call HistoryAddAttr( &
    & varname = 'time', attrname = 'standard_name', & ! (in)
    & value = 'time', &                               ! (in)
    & history = gthist_init)                          ! (inout)
  call HistoryAddAttr( &
    & varname = 'sig', attrname = 'positive', &   ! (in)
    & value = 'down', &                           ! (in)
    & history = gthist_init)                      ! (inout)
  call HistoryAddAttr( &
    & varname = 'sigm', attrname = 'positive', &  ! (in)
    & value = 'down', &                           ! (in)
    & history = gthist_init)                      ! (inout)

  call HistoryAddVariable( &
    & varname = 'Vor', &                            ! (in)
    & dims = StoA('lon', 'lat', 'sig', 'time'), &   ! (in)
    & longname = 'vorticity', &                     ! (in)
    & units = 's-1', xtype = 'double',          &   ! (in)
    & history = gthist_init)                        ! (inout)
  call HistoryAddAttr(&
    & varname = 'Vor', attrname = 'standard_name', & ! (in)
    & value = 'atmosphere_relative_vorticity', &     ! (in)
    & history = gthist_init)                         ! (inout)

  call HistoryAddVariable( &
    & varname = 'Div', &                            ! (in)
    & dims = StoA('lon', 'lat', 'sig', 'time'), &   ! (in)
    & longname = 'divergence', &                    ! (in)
    & units = 's-1', xtype = 'double',          &   ! (in)
    & history = gthist_init)                        ! (inout)
  call HistoryAddAttr(&
    & varname = 'Div', attrname = 'standard_name', & ! (in)
    & value = 'divergence_of_wind', &                ! (in)
    & history = gthist_init)                         ! (inout)

  call HistoryAddVariable( &
    & varname = 'Temp', &                           ! (in)
    & dims = StoA('lon', 'lat', 'sig', 'time'), &   ! (in)
    & longname = 'temperature', &                   ! (in)
    & units = 'K', xtype = 'double',          &     ! (in)
    & history = gthist_init)                        ! (inout)
  call HistoryAddAttr(&
    & varname = 'Temp', attrname = 'standard_name', & ! (in)
    & value = 'air_temperature', &                    ! (in)
    & history = gthist_init)                          ! (inout)

  call HistoryAddVariable( &
    & varname = 'QVap', &                           ! (in)
    & dims = StoA('lon', 'lat', 'sig', 'time'), &   ! (in)
    & longname = 'specific humidity', &             ! (in)
    & units = '1', xtype = 'double',          &     ! (in)
    & history = gthist_init)                        ! (inout)
  call HistoryAddAttr(&
    & varname = 'QVap', attrname = 'standard_name', & ! (in)
    & value = 'specific_humidity', &                  ! (in)
    & history = gthist_init)                          ! (inout)

  call HistoryAddVariable( &
    & varname = 'Ps', &                       ! (in)
    & dims = StoA('lon', 'lat', 'time'), &    ! (in)
    & longname = 'surface pressure', &        ! (in)
    & units = 'Pa', xtype = 'double', &       ! (in)
    & history = gthist_init)                  ! (inout)
  call HistoryAddAttr(&
    & varname = 'Ps', attrname = 'standard_name', & ! (in)
    & value = 'surface_air_pressure', &             ! (in)
    & history = gthist_init)                        ! (inout)

  !----------------------------------------------------------------
  !  初期値ファイルへのデータ出力
  !  Initial data output
  !----------------------------------------------------------------
  call HistoryPut( &
    & varname = 'Vor', array = xyz_VorB, &  ! (in)
    & history = gthist_init)                ! (inout)
  call HistoryPut( &
    & varname = 'Div', array = xyz_DivB, &  ! (in)
    & history = gthist_init)                ! (inout)
  call HistoryPut( &
    & varname = 'Temp', array = xyz_TempB, &  ! (in)
    & history = gthist_init)                  ! (inout)
  call HistoryPut( &
    & varname = 'QVap', array = xyz_QVapB, &  ! (in)
    & history = gthist_init)                  ! (inout)
  call HistoryPut( &
    & varname = 'Ps', array = xy_PsB, &  ! (in)
    & history = gthist_init)             ! (inout)

  call HistoryPut( &
    & varname = 'Vor', array = xyz_VorN, &  ! (in)
    & history = gthist_init)                ! (inout)
  call HistoryPut( &
    & varname = 'Div', array = xyz_DivN, &  ! (in)
    & history = gthist_init)                ! (inout)
  call HistoryPut( &
    & varname = 'Temp', array = xyz_TempN, &  ! (in)
    & history = gthist_init)                  ! (inout)
  call HistoryPut( &
    & varname = 'QVap', array = xyz_QVapN, &  ! (in)
    & history = gthist_init)                  ! (inout)
  call HistoryPut( &
    & varname = 'Ps', array = xy_PsN, &  ! (in)
    & history = gthist_init)             ! (inout)

  !----------------------------------------------------------------
  !  初期値ファイルへのデータ出力の終了処理
  !  Terminate initial data output
  !----------------------------------------------------------------
  call HistoryClose( history = gthist_init )  ! (inout)

  call Stop(clk_histput)   ! (inout)

  !----------------------------------------------------------------
  !  ファイル出力に関してメッセージを表示
  !  Print message of file output
  !----------------------------------------------------------------
  call MessageNotify( 'M', subname, &
    & 'Initial data file "%c" is generated. (time=%*r)', &
    & r=(/init_timeB, init_timeN/), n=(/2/), c1=trim(init_nc) )

  !----------------------------------------------------------------
  !  CPU 時間の総計を表示
  !  Print total CPU time
  !----------------------------------------------------------------
  call Result( &
    & clks = (/clk_setup, clk_histput/), &       ! (in)
    & total_auto = .true.)                       ! (in)

  call EndSub(subname)
end program init_sample
