!= 速度から渦度発散を出力するプログラム
!
!= A program that generates vorticity and divergence from velocity
!
! Authors:: Yasuhiro MORIKAWA
! Version:: $Id: convert_uv2vordiv.f90,v 1.1 2008-04-21 10:28:49 morikawa Exp $
! Tag Name:: $Name: dcpam4-20080609-1 $
! Copyright:: Copyright (C) GFD Dennou Club, 2007. All rights reserved.
! License:: See COPYRIGHT[link:../../../COPYRIGHT]
!
! Note that Japanese and English are described in parallel.
!
! 入力ファイルの "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-20080609-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= --output=' )
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