Path: | main/convert_vordiv2uv.f90 |
Last Update: | Sat Jun 14 13:07:42 +0900 2008 |
Authors: | Yasuhiro MORIKAWA |
Version: | $Id: convert_vordiv2uv.f90,v 1.3 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 |
Note that Japanese and English are described in parallel.
入力ファイルの "Vor" と "Div" という変数の情報を参照し, "U" と "V" に変換して出力ファイルへと出力します. 入力ファイルには軸データ "lon", "lon_weight", "lat", "lat_weight", "sig", "sigm" が格納されている 事を想定します.
Variables "Vor" and "Div" in an input file are converted into "U" and "V". 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.
Main Program : |
program convert_vordiv2uv 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, VorDiv2UV 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_vordiv2uv $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 = 'vordiv2uv' continue !--------------------------------------------------------- ! コマンドライン引数の処理 ! Command line arguments handling !--------------------------------------------------------- call Open( arg ) call HelpMsg( arg, 'Title', title ) call HelpMsg( arg, 'Usage', './convert_vordiv2uv [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, PI = PI, RPlanet = RPlanet, Omega = Omega, Cp = Cp, RAir = RAir, EpsV = EpsV, VisOrder = VisOrder, EFoldTime = EFoldTime ) ! (out) !--------------------------------------------------------- ! 軸データ取得 ! Get axes data !--------------------------------------------------------- call HistoryGetPointer ( file = input_data, varname = 'lon', array = x_Lon ) ! (out) call HistoryGetPointer ( file = input_data, varname = 'lon_weight', array = x_Lon_Weight ) ! (out) call HistoryGetPointer ( file = input_data, varname = 'lat', array = y_Lat ) ! (out) call HistoryGetPointer ( file = input_data, varname = 'lat_weight', array = y_Lat_Weight ) ! (out) call HistoryGetPointer ( file = input_data, varname = 'sig', array = z_Sigma ) ! (out) call HistoryGetPointer ( file = input_data, varname = 'sigm', 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, nmax = nmax, imax = imax, jmax = jmax, kmax = kmax, PI = PI, RPlanet = RPlanet, Omega = Omega, Cp = Cp, RAir = RAir, EpsV = EpsV, VisOrder = VisOrder, EFoldTime = EFoldTime, DelTIme = DelTime ) ! (in) call DynSpAsPutLine( dyn_sp_as00 ) ! (in) !--------------------------------------------------------- ! 渦度発散データ読み込み ! Load vorticity and divergence data !--------------------------------------------------------- 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 HistoryGet ( file = input_data, varname = 'Vor', array = xyz_Vor ) ! (out) call HistoryGet ( file = input_data, varname = 'Div', array = xyz_Div ) ! (out) !--------------------------------------------------------- ! 速度データへ変換 ! Convert into velocities !--------------------------------------------------------- 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 VorDiv2UV( dyn_sp_as = dyn_sp_as00, xyz_Vor = xyz_Vor, xyz_Div = xyz_Div, xyz_U = xyz_U, xyz_V = xyz_V ) ! (out) !--------------------------------------------------------- ! データ出力 ! Output data !--------------------------------------------------------- call HistoryCreate( history = gthist, file = output_data, title = 'U and V data converted from "' // trim(input_data) // '"', source = source, institution = institution, dims = StoA('lon', 'lat', 'sig', 'sigm'), dimsizes = (/imax, jmax, kmax, kmax + 1/), longnames = StoA('longitude', 'latitude', 'sigma at layer midpoints', 'sigma at layer end-points (half level))'), units = StoA('degree_east', 'degree_north', '1', '1') ) ! (in) call HistoryPut( history = gthist, varname = 'lon', array = x_Lon / PI * 180.0_DP ) ! (in) call HistoryPut( history = gthist, varname = 'lat', array = y_Lat / PI * 180.0_DP ) ! (in) call HistoryPut( history = gthist, varname = 'sig', array = z_Sigma ) ! (in) call HistoryPut( history = gthist, varname = 'sigm', array = r_Sigma ) ! (in) call HistoryAddAttr( history = gthist, varname = 'lon', attrname = 'standard_name', value = 'longitude' ) ! (in) call HistoryAddAttr( history = gthist, varname = 'lat', attrname = 'standard_name', value = 'latitude' ) ! (in) call HistoryAddAttr( history = gthist, varname = 'sig', attrname = 'standard_name', value = 'atmosphere_sigma_coordinate' ) ! (in) call HistoryAddAttr( history = gthist, varname = 'sigm', attrname = 'standard_name', value = 'atmosphere_sigma_coordinate' ) ! (in) call HistoryAddAttr( history = gthist, varname = 'sig', attrname = 'positive', value = 'down' ) ! (in) call HistoryAddAttr( history = gthist, varname = 'sigm', attrname = 'positive', value = 'down' ) ! (in) call HistoryAddVariable( history = gthist, varname = 'lon_weight', dims = StoA('lon'), longname = 'weight for integration in longitude', units = 'radian', xtype = 'double' ) ! (in) call HistoryAddAttr( history = gthist, varname = 'lon', attrname = 'gt_calc_weight', value = 'lon_weight' ) ! (in) call HistoryPut( history = gthist, varname = 'lon_weight', array = x_Lon_Weight ) ! (in) call HistoryAddVariable( history = gthist, varname = 'lat_weight', dims = StoA('lat'), longname = 'weight for integration in latitude', units = 'radian', xtype = 'double' ) ! (in) call HistoryAddAttr( history = gthist, varname = 'lat', attrname = 'gt_calc_weight', value = 'lat_weight' ) ! (in) call HistoryPut( history = gthist, varname = 'lat_weight', array = y_Lat_Weight ) ! (in) call HistoryAddVariable( history = gthist, varname = 'U', dims = StoA('lon', 'lat', 'sig'), longname = 'eastward wind', units = 'm/s', xtype = 'double' ) ! (in) call HistoryAddAttr( history = gthist, varname = 'U', attrname = 'standard_name', value = 'eastward_wind' ) ! (in) call HistoryAddVariable( history = gthist, varname = 'V', dims = StoA('lon', 'lat', 'sig'), longname = 'northward wind', units = 'm/s', xtype = 'double' ) ! (in) call HistoryAddAttr( history = gthist, varname = 'V', attrname = 'standard_name', value = 'northward_wind' ) ! (in) call HistoryPut( history = gthist, varname = 'U', array = xyz_U ) ! (in) call HistoryPut( history = gthist, varname = 'V', array = xyz_V ) ! (in) call HistoryClose( history = gthist ) ! (inout) !--------------------------------------------------------- ! 終了処理 ! Termination !--------------------------------------------------------- call DynSpAsClose( dyn_sp_as00 ) ! (inout) call DynSpAsPutLine( dyn_sp_as00 ) end program convert_vordiv2uv