!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Program init
!
!   * Developers: Morikawa Yasuhiro
!   * Version: $Id: init.f90,v 1.2 2005/11/29 14:45:06 yukiko Exp $
!   * Tag Name: $Name: dcpam3-20051207 $
!   * Change History: 
!
!== Overview
!
!Generate Various Initial Data
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!
!
                                                                 !=end
program init
                                                                 !=begin
  !== Dependency
  !----- 構造体参照モジュール -----
  use type_mod, only: STRING, INTKIND, REKIND, DBKIND
  use nmlfile_mod, only: nmlfile_init, nmlfile_open, nmlfile_close
  use constants_mod, only: constants_init, pi, R0

  !----- 格子点取得モジュール -----
  use grid_3d_mod,         only: grid_3d_init, im, jm, km, grid_3d_end
  use grid_wavenumber_mod, only: grid_wavenumber_init, nm, grid_wavenumber_end

  !----- 座標データ生成モジュール -----
  use axis_type_mod, only: AXISINFO
  use axis_x_mod, only: axis_x_init, axis_x_weight, axis_x_spectral, axis_x_end
  use axis_y_mod, only: axis_y_init, axis_y_weight, axis_y_spectral, axis_y_end
  use axis_z_mod, only: axis_z_init, axis_z_sigmahalf_manual, axis_z_end

  !----- データI/Oモジュール -----
  use io_gt4_out_mod, only : io_gt4_out_init   , io_gt4_out_SetDims, &
       &                     io_gt4_out_SetVars, io_gt4_out_Put    , &
       &                     io_gt4_out_end

  !----- 時刻管理モジュール -----
  use time_mod, only: time_init, tvar, ttype, tname, tunit, time_end

  !----- SPMODEL モジュール  -----
  use spml_mod,  only: spml_init, xya_wa, wa_Div_xya_xya &
       &              ,wa_LaplaInv_wa, wa_xya, xya_GradLat_wa, xya_GradLon_wa

  !----- デバッグ・汎用ツール -----
  use dc_trace,  only: SetDebug, DbgMessage, BeginSub, EndSub, DataDump
  use dc_string, only: toChar, StriEq, LChar, StrHead
  use dc_message,only: MessageNotify
                                                                 !=end
  implicit none
                                                                 !=begin
  !== NAMELIST
  !
  !初期値の種類の設定を行なう。
  !現在 condition に与えて有効なのは以下の値である。
  !
  !  * (({ rigid body rotation })) 
  !    * 剛体回転流を与える。風速の最大値は VelLonMax_rbr に与える。
  !      * これを選択した場合、強制的に VorDiv_Priority は .false. 
  !        に設定される。
  !
  !  * (({ convex of surface pressure }))
  !    * 地表面気圧の「山」を与える。位置とサイズは、度数で与える場合は
  !      Lon_Center_Dig, Lat_Center_Dig, LonLatRadius_Dig を用い、
  !      ラジアンで与える場合は Lon_Center_Rad, Lat_Center_Rad,
  !      LonLatRadius_Rad を用いる (ただし、Rad_Priority を .true.
  !      にする必要がある)。 最大値は PsMax で与える。
  !
  !  * (({ convex of temperature }))
  !    * 温度の「山」を与える。位置とサイズは、度数で与える場合は
  !      Lon_Center_Dig, Lat_Center_Dig, LonLatRadius_Dig を用い、
  !      ラジアンで与える場合は Lon_Center_Rad, Lat_Center_Rad,
  !      LonLatRadius_Rad を用いる (ただし、Rad_Priority を .true.
  !      にする必要がある)。 最大値は TempMax で与える。
  !
  !  * (({ agcm5.3 default value }))
  !    * 等温静止大気. 温度場に擾乱を与える. 既定値は以下. 
  !      QVapAve   = 1.0d-10 
  !      TempAve   = 250.0d0 
  !      PsAve     = 1.0d5 
  !      TPRTRB     = 0.1 
  !
  !  * その他
  !    * 等温無風の初期値を与える。
  !
  !VelLonAve, VelLonAve, VelLatAve, VorAve, DivAve, TempAve, QVapAve, PsAve
  !には、それぞれ平均値を与える。
  !
  !デフォルトでは風速から渦度発散を生成するが、
  !VorDiv_Priority を .true. にする事で、渦度発散から風速を生成する。
  !
  character(STRING) :: condition = ''      ! 初期値の種類

  real(DBKIND)      :: VelLonAve = 0.0d0   ! 速度経度成分平均値
  real(DBKIND)      :: VelLatAve = 0.0d0   ! 速度緯度成分平均値
  real(DBKIND)      :: VorAve    = 0.0d0   ! 渦度平均値
  real(DBKIND)      :: DivAve    = 0.0d0   ! 発散平均値
  real(DBKIND)      :: TempAve   = 273.0d0 ! 温度平均値
  real(DBKIND)      :: QVapAve   = 0.0d0   ! 比湿平均値
  real(DBKIND)      :: PsAve     = 1.0d5   ! 地表面圧力平均値

  logical           :: VorDiv_Priority = .false. ! 渦度発散から風速を生成

  !  !for 'rigid body rotation'
  real(DBKIND)      :: VelLonMax_rbr = 1.0d2    ! Maximum of 'VelLon'

  !  !for 'convex of surface pressure' or 'convex of temperature'
  real(DBKIND)      :: LonLat_Radius_Deg  = 20.0   ! 半径 (度数)
  real(DBKIND)      :: LonLat_Radius_Rad  = 0.349  ! 半径 (ラジアン)

  real(DBKIND)      :: Lat_Center_Deg = 45.0  ! 緯度の中心位置 (度数)
  real(DBKIND)      :: Lon_Center_Deg = 100.0 ! 経度の中心位置 (度数)
  real(DBKIND)      :: Lat_Center_Rad = 0.785 ! 緯度の中心位置 (ラジアン)
  real(DBKIND)      :: Lon_Center_Rad = 1.745 ! 経度の中心位置 (ラジアン)

  logical           :: Rad_Priority   = .false. ! ラジアン表記を優先

  !  !for 'convex of surface pressure'
  real(DBKIND)      :: PsMax   = -200.0d2       ! Maximum of 'Ps'

  !  !for 'convex of temperature'
  real(DBKIND)      :: TempMax = 10.0d0         ! Maximum of 'Temp'

  !  !for 'agcm5.3 default value'
  real(DBKIND)      :: TPRTRB     = 0.1         ! 温度擾乱最大値

  namelist /init_nml/ &
       & condition           , & ! 初期値の種類

       & VelLonAve           , & ! 速度経度成分平均値
       & VelLatAve           , & ! 速度緯度成分平均値
       & VorAve              , & ! 渦度平均値        
       & DivAve              , & ! 発散平均値        
       & TempAve             , & ! 温度平均値        
       & QVapAve             , & ! 比湿平均値        
       & PsAve               , & ! 地表面圧力平均値  

       & VorDiv_Priority     , & ! 渦度発散から風速を生成

       & VelLonMax_rbr       , & ! 速度経度成分 (剛体回転流用)

       & LonLat_Radius_Deg   , & ! 半径 (度数)
       & LonLat_Radius_Rad   , & ! 半径 (ラジアン)

       & Lat_Center_Deg      , & ! 緯度の中心位置 (度数)
       & Lon_Center_Deg      , & ! 経度の中心位置 (度数)
       & Lat_Center_Rad      , & ! 緯度の中心位置 (ラジアン)
       & Lon_Center_Rad      , & ! 経度の中心位置 (ラジアン)

       & Rad_Priority        , & ! ラジアン表記を優先

       & PsMax               , & ! Maximum of 'Ps'

       & TempMax                 ! Maximum of 'Temp'
                                                                 !=end

  !For Generate AXES
  type(AXISINFO)::          &
       & x_Lon            , & ! 経度座標
       & x_Lon_Weight     , & ! 経度重み座標
       & y_Lat            , & ! 緯度座標
       & y_Lat_Weight     , & ! 緯度重み座標
       & z_Sigma          , & ! σレベル(整数)座標
       & r_Sigma              ! σレベル(半整数)座標

  !For Generate Initial Data
  real(DBKIND), pointer     ::  &
       & xyz_VelLon(:,:,:)    , & ! 格子点データ(速度経度成分)
       & xyz_VelLat(:,:,:)    , & ! 格子点データ(速度緯度成分)
       & xyz_Vor(:,:,:)       , & ! 格子点データ(渦度)
       & xyz_Div(:,:,:)       , & ! 格子点データ(発散)
       & xyz_Temp(:,:,:)      , & ! 格子点データ(温度)
       & xyz_QVap(:,:,:)      , & ! 格子点データ(比湿)

       & xyz_Ps(:,:,:)               ! 格子点データ(地表面気圧)

  !Axes data (for Use as 3-dimensional data)
  real(DBKIND), allocatable :: &
       & xyz_Lon(:,:,:)      , & ! 経度座標 (ラジアン表記)
       & xyz_Lat(:,:,:)      , & ! 緯度座標 (ラジアン表記)
       & xyz_Sigma(:,:,:)    , & ! σレベル(整数)座標  
       & xyr_Sigma(:,:,:)        ! σレベル(半整数)座標
  real(DBKIND) :: RadDegFact     ! ラジアンを度数に変換する係数

  !for Generate Velocity from Vorticity and Divergence
  real(DBKIND), allocatable :: &
       & wz_Psi(:,:)         , & ! スペクトル(流線関数)
       & wz_Chi(:,:)             ! スペクトル(ポテンシャル)

  !for nmlfile_mod
  integer(INTKIND)   :: nmlstat, nmlunit
  logical            :: nmlreadable

  !----- 作業用内部変数 -----
  integer(INTKIND)   :: i,j,k

  character(STRING),  parameter:: version = &
       & '$Id: init.f90,v 1.2 2005/11/29 14:45:06 yukiko Exp $'
  character(STRING),  parameter:: tagname = '$Name: dcpam3-20051207 $'
  character(STRING),  parameter:: subname = "init"

  !-------------------------------------------------------------------
  !   Set Debug Mode
  !-------------------------------------------------------------------
!  call SetDebug

  call BeginSub(subname)

  !----------------------------------------------------------------
  !   Version identifier
  !----------------------------------------------------------------
  call DbgMessage('%c :: %c', c1=trim(version), c2=trim(tagname))

  !-------------------------------------------------------------------
  !   NAMELIST file Setting
  !-------------------------------------------------------------------
  call nmlfile_init('init.nml')

  !----------------------------------------------------------------
  !   Read init_nml
  !----------------------------------------------------------------
  call nmlfile_open(nmlunit, nmlreadable)
  if (nmlreadable) then
     read(nmlunit, nml=init_nml, iostat=nmlstat)
     call DbgMessage('Stat of NAMELIST init_nml Input is <%d>', &
          &           i=(/nmlstat/))
     write(0, nml=init_nml)
  else
     call DbgMessage('Not Read NAMELIST init_nml')
     call MessageNotify('W', subname, &
          & 'Can not Read NAMELIST init_nml. Force Use Default Value.')
  end if
  call nmlfile_close

  !-------------------------------------------------------------------
  !   Initialize Dependent Modules
  !-------------------------------------------------------------------
  call constants_init        ! 定数の設定
  call grid_3d_init          ! 空間格子点数の取得
  call grid_wavenumber_init  ! 波数格子点数の取得
  call axis_x_init           ! 座標軸データの生成
  call axis_y_init           ! 座標軸データの生成
  call axis_z_init           ! 座標軸データの生成
  call io_gt4_out_init       ! データ出力の初期設定
  call constants_init        ! 物理定数の取得
  call spml_init             ! SPMODEL の初期設定

  !-------------------------------------------------------------------
  !  軸データ生成
  !-------------------------------------------------------------------
!!$  call axis_x_weight(x_Lon_Weight)    ! 経度座標重みデータ取得
  call axis_x_spectral(x_Lon)         ! 経度座標データ取得
!!$  call axis_y_weight(y_Lat_Weight)    ! 緯度座標重みデータ取得
  call axis_y_spectral(y_Lat)         ! 緯度座標データ取得
  call axis_z_sigmahalf_manual &
       & ( z_Sigma           , &      ! 整数σレベル座標データ取得
       &   r_Sigma  )                 ! 半整数σレベル座標データ取得

  !-------------------------------------------------------------------
  !   出力用の軸データ設定
  !-------------------------------------------------------------------
  call io_gt4_out_SetDims(x_Lon)        ! 経度座標重みデータ取得
!!$  call io_gt4_out_SetDims(x_Lon_Weight) ! 経度座標データ取得
  call io_gt4_out_SetDims(y_Lat)        ! 緯度座標重みデータ取得
!!$  call io_gt4_out_SetDims(y_Lat_Weight) ! 緯度座標データ取得
  call io_gt4_out_SetDims(z_Sigma)      ! 整数σレベル座標データ取得
  call io_gt4_out_SetDims(r_Sigma)  ! 半整数σレベル座標データ取得

  !-------------------------------------------------------------------
  !   出力用の変数データ設定
  !-------------------------------------------------------------------
  call io_gt4_out_SetVars('VelLon')
  call io_gt4_out_SetVars('VelLat')
  call io_gt4_out_SetVars('Vor')
  call io_gt4_out_SetVars('Div')
  call io_gt4_out_SetVars('Temp')
  call io_gt4_out_SetVars('QVap')
  call io_gt4_out_SetVars('Ps')

  !-------------------------------------------------------------------
  !   Allocate Variable
  !-------------------------------------------------------------------
  allocate( xyz_VelLon(im,jm,km) )
  allocate( xyz_VelLat(im,jm,km) )
  allocate( xyz_Vor(im,jm,km)  )
  allocate( xyz_Div(im,jm,km)  )
  allocate( xyz_Temp(im,jm,km) )
  allocate( xyz_QVap(im,jm,km) )
  allocate( xyz_Ps(im,jm,km)   )

  allocate( xyz_Lon(im,jm,km) )
  allocate( xyz_Lat(im,jm,km) )
  allocate( xyz_Sigma(im,jm,km) )
  allocate( xyr_Sigma(im,jm,km+1) )

  allocate( wz_Psi((nm+1)*(nm+1), km) )
  allocate( wz_Chi((nm+1)*(nm+1), km) )

  !-------------------------------------------------------------------
  !   Set Average Value
  !-------------------------------------------------------------------
  xyz_VelLon = VelLonAve
  xyz_VelLat = VelLatAve
  xyz_Vor    = VorAve
  xyz_Div    = DivAve
  xyz_Temp   = TempAve
  xyz_QVap   = QVapAve
  xyz_Ps     = PsAve

  !-------------------------------------------------------------------
  !  Create 3-dimentional axis data
  !-------------------------------------------------------------------
  if (  StrHead( 'radians', trim(LChar(x_Lon%axisinfo%units)) ) .or.&
       & StrHead( 'rad.', trim(LChar(x_Lon%axisinfo%units)) ) ) then
     RadDegFact = 1.
  else
     RadDegFact = 180. / pi
  end if

  do i = 1, im
     xyz_Lon(i,:,:) = x_Lon%a_Dim(i) / RadDegFact
  end do

  if (  StrHead( 'radians', trim(LChar(y_Lat%axisinfo%units)) ) .or.&
       & StrHead( 'rad.', trim(LChar(y_Lat%axisinfo%units)) ) ) then
     RadDegFact = 1.
  else
     RadDegFact = 180. / pi
  end if
  do j = 1, jm
     xyz_Lat(:,j,:) = y_Lat%a_Dim(j) / RadDegFact
  end do

  do k = 1, km
     xyz_Sigma(:,:,k) = z_Sigma%a_Dim(k)
  end do

  do k = 1, km + 1
     xyr_Sigma(:,:,k) = r_Sigma%a_Dim(k)
  end do

!!$  call DataDump('x_Lon1', x_Lon%a_Dim, strlen=70)
!!$  call DataDump('x_Lon2', x_Lon%a_Dim * RadDegFact , strlen=70)
!!$  call DataDump('x_Lon3', xyz_Lon / RadDegFact , strlen=70)
!!$  call DataDump('xyz_Lon', xyz_Lon, strlen=70)
!!$  call DataDump('xyz_Lat', xyz_Lat, strlen=70)
!!$  call DataDump('xyz_Sigma', xyz_Sigma, strlen=70)

  !-------------------------------------------------------------------
  !   Set Various Conditions
  !-------------------------------------------------------------------
  if (  StriEq( 'rigid body rotation', trim(condition) )  ) then
     ! 剛体回転流
     xyz_VelLon = xyz_VelLon + VelLonMax_rbr * cos( xyz_Lat )
     VorDiv_Priority = .false.

     call MessageNotify('M', subname, &
          & 'Selected <%c>.', c1=trim(LChar(condition))  )

!!$     call DataDump('xyz_Lat', xyz_Lat, strlen=80)
!!$     call DataDump('y_Lat', y_Lat%a_Dim, strlen=80)
!!$     call DataDump('xyz_VelLon', xyz_VelLon, strlen=80)

  elseif (  StriEq( 'convex of surface pressure', &
       &             trim(condition) )  ) then

     ! 度数をラジアンに変換するため
     RadDegFact = 180. / pi
     if (.not. Rad_Priority) then
        LonLat_Radius_Rad = LonLat_Radius_Deg / RadDegFact
        Lat_Center_Rad    = Lat_Center_Deg    / RadDegFact
        Lon_Center_Rad    = Lon_Center_Deg    / RadDegFact
     end if

     xyz_Ps = xyz_Ps &
          & + PsMax &
          &   * ( 1.0d0 &
          &        - sqrt(  &
          &               min(  LonLat_Radius_Rad**2 ,              &
          &                     ( xyz_Lon - Lon_Center_Rad )**2     &
          &                       + ( xyz_Lat - Lat_Center_Rad )**2 &
          &                  )     &
          &              ) / LonLat_Radius_Rad  &
          &      )

     call MessageNotify('M', subname, &
          & 'Selected <%c>.', c1=trim(LChar(condition))  )

  elseif (  StriEq( 'convex of temperature', &
       &             trim(condition) )  ) then

     ! 度数をラジアンに変換するため
     RadDegFact = 180. / pi
     if (.not. Rad_Priority) then
        LonLat_Radius_Rad = LonLat_Radius_Deg / RadDegFact
        Lat_Center_Rad    = Lat_Center_Deg    / RadDegFact
        Lon_Center_Rad    = Lon_Center_Deg    / RadDegFact
     end if

     xyz_Temp = xyz_Temp &
          & + TempMax &
          &   * ( 1.0D0 &
          &        - sqrt(  &
          &               min(  LonLat_Radius_Rad**2 ,              &
          &                     ( xyz_Lon - Lon_Center_Rad )**2     &
          &                       + ( xyz_Lat - Lat_Center_Rad )**2 &
          &                  )     &
          &              ) / LonLat_Radius_Rad  &
          &      )

!!!     xyz_VelLon = xyz_VelLon &
!!!          & + TempMax &
!!!          &   * ( 1.0D0 &
!!!          &        - sqrt(  &
!!!          &               min(  LonLat_Radius_Rad**2 ,              &
!!!          &                     ( xyz_Lon - Lon_Center_Rad )**2     &
!!!          &                       + ( xyz_Lat - Lat_Center_Rad )**2 &
!!!          &                  )     &
!!!          &              ) / LonLat_Radius_Rad  &
!!!          &      )
!!!
!!!     xyz_VelLat = xyz_VelLat &
!!!          & + TempMax &
!!!          &   * ( 1.0D0 &
!!!          &        - sqrt(  &
!!!          &               min(  LonLat_Radius_Rad**2 ,              &
!!!          &                     ( xyz_Lon - Lon_Center_Rad )**2     &
!!!          &                       + ( xyz_Lat - Lat_Center_Rad )**2 &
!!!          &                  )     &
!!!          &              ) / LonLat_Radius_Rad  &
!!!          &      )
!!!
!!!     xyz_QVap = xyz_QVap &
!!!          & + TempMax &
!!!          &   * ( 1.0D0 &
!!!          &        - sqrt(  &
!!!          &               min(  LonLat_Radius_Rad**2 ,              &
!!!          &                     ( xyz_Lon - Lon_Center_Rad )**2     &
!!!          &                       + ( xyz_Lat - Lat_Center_Rad )**2 &
!!!          &                  )     &
!!!          &              ) / LonLat_Radius_Rad  &
!!!          &      )

     call MessageNotify('M', subname, &
          & 'Selected <%c>.', c1=trim(LChar(condition))  )

  elseif (  StriEq( 'agcm5.3 default value', &
       &             trim(condition) )  ) then

     do k = 1, km
        do i = 1, im
           do j = 1, jm
              xyz_Temp(i,j,k) = xyz_Temp(i,j,k) &
                   &  + SIN ( REAL( i*j*(km-k) ) / REAL( im*jm*km )*10.0 ) &
                   &   * TPRTRB
           end do
        end do
     end do

     call MessageNotify('M', subname, &
          & 'Selected <%c>.', c1=trim(LChar(condition))  )
     
  else
     call MessageNotify('M', subname, &
          & 'Selected Isothermal and Nowind (Default).'  )
  end if


  !-------------------------------------------------------------------
  !   Generate Vorticity and Divergence from Velocity
  !-------------------------------------------------------------------
  if ( .not. VorDiv_Priority ) then
     xyz_Vor = &
          & xya_wa(                                            &
          &   wa_Div_xya_xya( xyz_VelLat , - xyz_VelLon ) / R0 &
          & )

     xyz_Div = &
          & xya_wa(                                            &
          &   wa_Div_xya_xya( xyz_VelLon ,   xyz_VelLat ) / R0 &
          & )
  end if


  !-------------------------------------------------------------------
  !   Generate Velocity from Vorticity and Divergence
  !-------------------------------------------------------------------
  if ( VorDiv_Priority ) then
     wz_Psi = wa_LaplaInv_wa(  wa_xya( xyz_Vor )  ) * R0**2
     wz_Chi = wa_LaplaInv_wa(  wa_xya( xyz_Div )  ) * R0**2

     xyz_VelLon = (  xya_GradLon_wa( wz_Chi ) &
          &                - xya_GradLat_wa( wz_Psi )  ) / R0

     xyz_VelLat = (  xya_GradLon_wa( wz_Psi ) &
          &                + xya_GradLat_wa( wz_Chi )  ) / R0
  end if


  !-------------------------------------------------------------------
  !   データの出力 (t-Δt)
  !-------------------------------------------------------------------
  call io_gt4_out_Put('VelLon' , real(xyz_VelLon(:,:,:), REKIND) )
  call io_gt4_out_Put('VelLat' , real(xyz_VelLat(:,:,:), REKIND) )
  call io_gt4_out_Put('Vor'  , real(xyz_Vor(:,:,:) , REKIND) )
  call io_gt4_out_Put('Div'  , real(xyz_Div(:,:,:) , REKIND) )
  call io_gt4_out_Put('Temp' , real(xyz_Temp(:,:,:), REKIND) )
  call io_gt4_out_Put('QVap' , real(xyz_QVap(:,:,:), REKIND) )
  call io_gt4_out_Put('Ps'   , real(xyz_Ps(:,:,1)  , REKIND) )


  !-------------------------------------------------------------------
  !   データの出力 (t)
  !-------------------------------------------------------------------
  call io_gt4_out_Put('VelLon' , real(xyz_VelLon(:,:,:), REKIND) )
  call io_gt4_out_Put('VelLat' , real(xyz_VelLat(:,:,:), REKIND) )
  call io_gt4_out_Put('Vor'  , real(xyz_Vor(:,:,:) , REKIND) )
  call io_gt4_out_Put('Div'  , real(xyz_Div(:,:,:) , REKIND) )
  call io_gt4_out_Put('Temp' , real(xyz_Temp(:,:,:), REKIND) )
  call io_gt4_out_Put('QVap' , real(xyz_QVap(:,:,:), REKIND) )
  call io_gt4_out_Put('Ps'   , real(xyz_Ps(:,:,1)  , REKIND) )

  !-------------------------------------------------------------------
  !   終了処理
  !-------------------------------------------------------------------
  call grid_3d_end          ! 空間格子点数の取得
  call grid_wavenumber_end  ! 波数格子点数の取得
  call axis_x_end           ! 座標軸データの生成
  call axis_y_end           ! 座標軸データの生成
  call axis_z_end           ! 座標軸データの生成
  call io_gt4_out_end       ! データ出力の初期設定

  call EndSub(subname)
end program init
