!= 速度から渦度発散を出力するプログラム
!
!= A program that generates vorticity and divergence from velocity
!
! Authors::   Yasuhiro MORIKAWA
! Version::   $Id: convert_uv2vordiv.f90,v 1.2 2008-06-14 04:07:42 morikawa Exp $
! Tag Name::  $Name: dcpam4-20080624-1 $
! Copyright:: Copyright (C) GFD Dennou Club, 2007. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!
! <b>Note that Japanese and English are described in parallel.</b>
!
! 入力ファイルの "U" と "V"  という変数の情報を参照し, 
! "Vor" と "Div" に変換して出力ファイルへと出力します. 
! 入力ファイルには軸データ "lon", "lon_weight", "lat", "lat_weight", 
! "sig", "sigm" が格納されている
! 事を想定します. 
!
! Variables "U" and "V" in an input file are converted into 
! "Vor" and "Div" . The velocities are output to an output file. 
! It is expected that axes data "lon", "lon_weight", "lat", "lat_weight", 
! "sig", "sigm" are stored in input file.
!

program convert_uv2vordiv
  use constants, only: CONST, Create, Get
  use dc_types, only: DP, STRING
  use dc_message, only: MessageNotify
  use dyn_spectral_as83, only: DYNSPAS83, DynSpAsCreate, DynSpAsClose, &
    & DynSpAsPutLine, DynSpAsInitialized, UV2VorDiv
  use dc_string, only: StoA
  use dc_args, only: ARGS, Open, HelpMsg, Option, Debug, Help, Strict, Close
  use gt4_history, only: HistoryGet, HistoryGetPointer
  use gt4_history, only: GT_HISTORY, &
    & HistoryCreate, HistoryAddVariable, HistoryPut, HistoryClose, &
    & HistoryAddAttr
  implicit none

  !---------------------------------------------------------
  !  実験の表題, モデルの名称, 所属機関名
  !  Title of a experiment, name of model, sub-organ
  !---------------------------------------------------------
  character(*), parameter:: title = &
    & 'convert_uv2vordiv $Name: dcpam4-20080624-1 $ :: ' // &
    & 'Test program of "dyn_spectral_as83" module'
  character(*), parameter:: source = &
    & 'dcpam4 ' // &
    & '(See http://www.gfd-dennou.org/library/dcpam)'
  character(*), parameter:: institution = &
    & 'GFD Dennou Club (See http://www.gfd-dennou.org)'

  !-------------------------------------------------------------------
  !  格子点数・最大全波数
  !  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

  !-----------------------------------------------------------------
  !  物理定数等の設定
  !  Configure physical constants etc.
  !-----------------------------------------------------------------
  type(CONST):: const_earth
  real(DP):: PI         ! $ \pi $ .    円周率.         Circular constant
  real(DP):: RPlanet    ! $ a $ .      惑星半径.       Radius of planet
  real(DP):: Omega      ! $ \Omega $ . 回転角速度.     Angular velocity
!!$  real(DP):: Grav       ! $ g $ .      重力加速度.     Gravitational acceleration
  real(DP):: Cp         ! $ C_p $ .    大気定圧比熱.   Specific heat of air at constant pressure
  real(DP):: RAir       ! $ R $ .      大気気体定数.   Gas constant of air
  real(DP):: EpsV       ! $ \epsilon_v $ .  水蒸気分子量比. Molecular weight ratio of water vapor
  integer:: VisOrder    ! 超粘性の次数.  Order of hyper-viscosity
  real(DP):: EFoldTime  ! 最大波数に対する e-folding time. E-folding time for maximum wavenumber
  real(DP):: DelTime

  !---------------------------------------------------------
  !  軸データ
  !  Axes data
  !---------------------------------------------------------
  real(DP), pointer:: x_Lon (:)
                              ! 経度. Longitude
  real(DP), pointer:: x_Lon_Weight(:)
                   ! 経度座標重み. 
                   ! Weight of longitude
  real(DP), pointer:: y_Lat (:)
                              ! 緯度. Latitude
  real(DP), pointer:: y_Lat_Weight(:)
                   ! 緯度座標重み. 
                   ! Weight of latitude
  real(DP), pointer:: z_Sigma (:)
                              ! $ \sigma $ レベル (整数). 
                              ! Full $ \sigma $ level
  real(DP), pointer:: r_Sigma (:)
                              ! $ \sigma $ レベル (半整数). 
                              ! Half $ \sigma $ level

  !---------------------------------------------------------
  !  物理量
  !  Physical values
  !---------------------------------------------------------
  real(DP), allocatable:: xyz_Vor (:,:,:)
                              ! $ \zeta (t-\Delta t) $ . 渦度. Vorticity
  real(DP), allocatable:: xyz_Div (:,:,:)
                              ! $ D (t-\Delta t) $ .     発散. Divergence
  real(DP), allocatable:: xyz_U (:,:,:)
                              ! $ u $ . 東西風速. Zonal wind
  real(DP), allocatable:: xyz_V (:,:,:)
                              ! $ v $ . 南北風速. Meridional wind

  !---------------------------------------------------------
  !  入出力データ
  !  Input/Output data
  !---------------------------------------------------------
  type(GT_HISTORY):: gthist
  character(STRING):: input_data
  character(STRING):: output_data

  !---------------------------------------------------------
  !  作業変数
  !  Work variables
  !---------------------------------------------------------
  type(ARGS):: arg            ! コマンドライン引数. 
                              ! Command line arguments
!!$  logical:: OPT_namelist      ! -N, --namelist オプションの有無. 
!!$                              ! Existence of '-N', '--namelist' option
!!$  character(STRING):: VAL_namelist
!!$                              ! -N, --namelist オプションの値. 
!!$                              ! Value of '-N', '--namelist' option
  logical:: OPT_input         ! -I, --input オプションの有無. 
                              ! Existence of '-I', '--input' option
  character(STRING):: VAL_input
                              ! -I, --input オプションの値. 
                              ! Value of '-I', '--input' option
  logical:: OPT_output        ! -O, --output オプションの有無. 
                              ! Existence of '-O', '--output' option
  character(STRING):: VAL_output
                              ! -O, --output オプションの値. 
                              ! Value of '-O', '--output' option

  type(DYNSPAS83):: dyn_sp_as00
  logical:: err

  character(*), parameter:: subname = 'uv2vordiv'
continue

  !---------------------------------------------------------
  !  コマンドライン引数の処理
  !  Command line arguments handling
  !---------------------------------------------------------
  call Open( arg )
  call HelpMsg( arg, 'Title', title )
  call HelpMsg( arg, 'Usage', &
    & './convert_uv2vordiv [Options] --input=<file> --output=<file>' )
  call HelpMsg( arg, 'Source', source )
  call HelpMsg( arg, 'Institution', institution )
!!$  call Option( arg, StoA('-N', '--namelist'), &
!!$    & OPT_namelist, VAL_namelist, help = "NAMELIST filename" )
  call Option( arg, StoA('-I', '--input'), &
    & OPT_input, VAL_input, help = "input data file" )
  call Option( arg, StoA('-O', '--output'), &
    & OPT_output, VAL_output, help = "output data file" )
  call Debug( arg ) ; call Help( arg ) ; call Strict( arg, severe = .true. )
  call Close( arg )

  !---------------------------------------------------------
  !  入出力データ
  !  Input/Output data
  !---------------------------------------------------------
  input_data = VAL_input
  output_data = VAL_output

  if ( trim(input_data) == '' ) then
    call MessageNotify('E', subname, 'Specify input data like "--input=VorDiv.nc"')
  end if

  if ( trim(output_data) == '' ) then
    call MessageNotify('E', subname, 'Specify output data like "--output=UV.nc"')
  end if

  !---------------------------------------------------------
  !  物理定数の準備
  !  Prepare physical constants
  !---------------------------------------------------------
  call Create( const_earth ) ! (inout)

  DelTime = 480.0_DP

  call Get( constant = const_earth, &              ! (inout)
    & PI = PI, RPlanet = RPlanet, Omega = Omega, & ! (out)
    & Cp = Cp, RAir = RAir, EpsV = EpsV, &         ! (out)
    & VisOrder = VisOrder, EFoldTime = EFoldTime ) ! (out)

  !---------------------------------------------------------
  !  軸データ取得
  !  Get axes data
  !---------------------------------------------------------
  call HistoryGetPointer ( &
    & file = input_data, varname = 'lon', & ! (in)
    & array = x_Lon )                       ! (out)

  call HistoryGetPointer ( &
    & file = input_data, varname = 'lon_weight', & ! (in)
    & array = x_Lon_Weight )                       ! (out)

  call HistoryGetPointer ( &
    & file = input_data, varname = 'lat', & ! (in)
    & array = y_Lat )                       ! (out)

  call HistoryGetPointer ( &
    & file = input_data, varname = 'lat_weight', & ! (in)
    & array = y_Lat_Weight )                       ! (out)

  call HistoryGetPointer ( &
    & file = input_data, varname = 'sig', & ! (in)
    & array = z_Sigma )                     ! (out)

  call HistoryGetPointer ( &
    & file = input_data, varname = 'sigm', & ! (in)
    & array = r_Sigma )                      ! (out)

  !-------------------------------------------------------------------
  !  格子点数・最大全波数の設定
  !  Configure grid points and maximum truncated wavenumber
  !-------------------------------------------------------------------
  imax = size ( x_Lon )
  jmax = size ( y_Lat )
  kmax = size ( z_Sigma )

  !-------------------------------------
  !  いい加減な最大全波数自動設定
  !  Irresponsible auto-setting of maximum truncated wavenumber
  if ( imax < 32 .and. jmax < 16 ) then
    call MessageNotify('E', subname, '"imax=%d" and "jmax=%d" is too small', &
      & i = (/imax, jmax/) )
  elseif ( imax < 64 .and. jmax < 32 ) then
    nmax = 10
  elseif ( imax < 128 .and. jmax < 64 ) then
    nmax = 21
  elseif ( imax < 192 .and. jmax < 96 ) then
    nmax = 42
  else
    nmax = 63
  end if

  call MessageNotify('M', subname, &
    & 'Resolution is nmax=%d imax=%d, jmax=%d, kmax=%d', &
    & i = (/nmax, imax, jmax, kmax/) )

  !---------------------------------------------------------
  !  初期設定
  !  Initialization
  !---------------------------------------------------------
  call DynSpAsCreate( dyn_sp_as = dyn_sp_as00, &    ! (inout)
    & nmax = nmax, &                                ! (in)
    & imax = imax, jmax = jmax, kmax = kmax, &      ! (in)
    & PI = PI, RPlanet = RPlanet, &                 ! (in)
    & Omega = Omega, Cp = Cp, RAir = RAir, &        ! (in)
    & EpsV = EpsV, &                                ! (in)
    & VisOrder = VisOrder, EFoldTime = EFoldTime, & ! (in)
    & DelTIme = DelTime )                           ! (in)
  call DynSpAsPutLine( dyn_sp_as00 ) ! (in)


  !---------------------------------------------------------
  !  渦度発散データ読み込み
  !  Load vorticity and divergence data
  !---------------------------------------------------------
  allocate( xyz_U(0:imax-1, 0:jmax-1, 0:kmax-1) )
  allocate( xyz_V(0:imax-1, 0:jmax-1, 0:kmax-1) )

  call HistoryGet ( &
    & file = input_data, varname = 'U', & ! (in)
    & array = xyz_U )                     ! (out)

  call HistoryGet ( &
    & file = input_data, varname = 'V', & ! (in)
    & array = xyz_V )                     ! (out)

  !---------------------------------------------------------
  !  速度データへ変換
  !  Convert into velocities
  !---------------------------------------------------------
  allocate( xyz_Vor(0:imax-1, 0:jmax-1, 0:kmax-1) )
  allocate( xyz_Div(0:imax-1, 0:jmax-1, 0:kmax-1) )

  call UV2VorDiv( dyn_sp_as = dyn_sp_as00, &  ! (inout)
    & xyz_U   = xyz_U,   xyz_V   = xyz_V, &   ! (in)
    & xyz_Vor = xyz_Vor, xyz_Div = xyz_Div )  ! (out)

  !---------------------------------------------------------
  !  データ出力
  !  Output data
  !---------------------------------------------------------
  call HistoryCreate( &
    & history = gthist,   &                                 ! (inout)
    & file = output_data, &                                 ! (in)
    & title = 'Vor and Div data converted from "' // &
    &   trim(input_data) // '"', &                          ! (in)
    & source = source, institution = institution, &         ! (in)
    & dims = StoA('lon', 'lat', 'sig', 'sigm'), &           ! (in)
    & dimsizes = (/imax, jmax, kmax, kmax + 1/), &          ! (in)
    & longnames = &
    &   StoA('longitude', 'latitude', &
    &        'sigma at layer midpoints', &
    &        'sigma at layer end-points (half level))'), &  ! (in)
    & units = StoA('degree_east', 'degree_north', &
    &              '1', '1') )                              ! (in)

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

  call HistoryAddAttr( &
    & history = gthist, &                            ! (inout)
    & varname = 'lon', attrname = 'standard_name', & ! (in)
    & value = 'longitude' )                          ! (in)
  call HistoryAddAttr( &
    & history = gthist, &                            ! (inout)
    & varname = 'lat', attrname = 'standard_name', & ! (in)
    & value = 'latitude' )                           ! (in)

  call HistoryAddAttr( &
    & history = gthist, &                              ! (inout)
    & varname = 'sig', attrname = 'standard_name', &   ! (in)
    & value = 'atmosphere_sigma_coordinate' )          ! (in)
  call HistoryAddAttr( &
    & history = gthist, &                              ! (inout)
    & varname = 'sigm', attrname = 'standard_name', &  ! (in)
    & value = 'atmosphere_sigma_coordinate' )          ! (in)
  call HistoryAddAttr( &
    & history = gthist, &                         ! (inout)
    & varname = 'sig', attrname = 'positive', &   ! (in)
    & value = 'down' )                            ! (in)
  call HistoryAddAttr( &
    & history = gthist, &                         ! (inout)
    & varname = 'sigm', attrname = 'positive', &  ! (in)
    & value = 'down' )                            ! (in)

  call HistoryAddVariable( &
    & history = gthist, &                                 ! (inout)
    & varname = 'lon_weight', &                           ! (in)
    & dims = StoA('lon'), &                               ! (in)
    & longname = 'weight for integration in longitude', & ! (in)
    & units = 'radian', xtype = 'double' )                ! (in)
  call HistoryAddAttr( &
    & history = gthist, &                             ! (inout)
    & varname = 'lon', attrname = 'gt_calc_weight', & ! (in)
    & value = 'lon_weight' )                          ! (in)
  call HistoryPut( &
    & history = gthist, &                            ! (inout)
    & varname = 'lon_weight', array = x_Lon_Weight ) ! (in)

  call HistoryAddVariable( &
    & history = gthist, &                                ! (inout)
    & varname = 'lat_weight', &                          ! (in)
    & dims = StoA('lat'), &                              ! (in)
    & longname = 'weight for integration in latitude', & ! (in)
    & units = 'radian', xtype = 'double' )               ! (in)
  call HistoryAddAttr( &
    & history = gthist, &                             ! (inout)
    & varname = 'lat', attrname = 'gt_calc_weight', & ! (in)
    & value = 'lat_weight' )                          ! (in)
  call HistoryPut( &
    & history = gthist, &                            ! (inout)
    & varname = 'lat_weight', array = y_Lat_Weight ) ! (in)

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

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

  call HistoryPut( &
    & history = gthist, &                ! (inout)
    & varname = 'Vor', array = xyz_Vor ) ! (in)
  call HistoryPut( &
    & history = gthist, &                ! (inout)
    & varname = 'Div', array = xyz_Div ) ! (in)

  call HistoryClose( history = gthist ) ! (inout)

  !---------------------------------------------------------
  !  終了処理
  !  Termination
  !---------------------------------------------------------
  call DynSpAsClose( dyn_sp_as00 ) ! (inout)
  call DynSpAsPutLine( dyn_sp_as00 )

end program convert_uv2vordiv
