dc_units.f90

Go to the documentation of this file.

! -*- mode: f90; coding: utf-8 -*-
!-----------------------------------------------------------------------
! Copyright (c) 2000-2026 Gtool Development Group. All rights reserved.
!-----------------------------------------------------------------------
!>
!> @author    Youhei SASAKI, Eizi TOYODA, Yasuhiro MORIKAWA
!> @copyright Copyright (C) GFD Dennou Club, 2000-2026. All rights reserved. <br/>
!>            License is BSD-2-Clause. see [COPYRIGHT](@ref COPYRIGHT) in detail
!> @en
!> @brief Unit system processing module
!> @details
!> This module provides procedures for unit system processing.
!> It includes a derived type UNITS for representing physical units and
!> operators for unit arithmetic (multiplication, division, addition).
!>
!> @section units_operators Operators and Procedures
!>
!> | Operator/Procedure | Description                                      |
!> |--------------------|--------------------------------------------------|
!> | assignment(=)      | Build UNITS from string or convert to string     |
!> | operator(*)        | Multiply two UNITS                               |
!> | operator(/)        | Divide two UNITS                                 |
!> | operator(+)        | Add two UNITS (with compatibility check)         |
!> | clear              | Initialize UNITS variable                        |
!> | deallocate         | Deallocate UNITS variable                        |
!> | add_okay           | Check if two UNITS are compatible for addition   |
!>
!> @enden
!>
!> @ja
!> @brief 単位系処理用モジュール
!> @details
!> このモジュールは単位系処理のためのプログラム群を提供します.
!> 物理単位を表現する派生型 UNITS と、単位の算術演算
!> (乗算、除算、加算) のための演算子を含みます.
!>
!> @section units_operators_ja 演算子・手続一覧
!>
!> | 演算子/手続 | 説明                                           |
!> |-------------|------------------------------------------------|
!> | assignment(=) | 文字列から UNITS を構築または文字列へ変換    |
!> | operator(*)   | 2つの UNITS を乗算                           |
!> | operator(/)   | 2つの UNITS を除算                           |
!> | operator(+)   | 2つの UNITS を加算 (互換性チェック付き)      |
!> | clear         | UNITS 変数の初期化                           |
!> | deallocate    | UNITS 変数の解放                             |
!> | add_okay      | 2つの UNITS が加算互換かチェック             |
!>
!> @endja
!>
 
module dc_units
  !
 
  use dc_types, only: dp, token, string
  implicit none
 
  private
  public:: units
  public:: clear, deallocate, add_okay
  public:: assignment(=), operator(*), operator(/), operator(+)
  private:: units_simplify
 
  !> @en
  !> @brief Derived type for representing physical units
  !> @details
  !> A physical unit is represented as a product of named unit elements
  !> with powers (exponents), plus an optional scaling factor and offset.
  !> For example, "m/s2" would have name=["m","s"], power=[1,-2], factor=1.
  !> @enden
  !>
  !> @ja
  !> @brief 物理単位を表現する派生型
  !> @details
  !> 物理単位は、名前付き単位要素のべき乗（指数）の積と、
  !> オプションのスケール係数およびオフセットとして表現されます.
  !> 例えば "m/s2" は name=["m","s"], power=[1,-2], factor=1 となります.
  !> @endja
  type units
    real(dp):: factor             !< @en Scale factor @enden @ja スケール係数 @endja
    integer:: nelems              !< @en Number of unit elements @enden @ja 単位要素の数 @endja
    character(TOKEN), pointer:: name(:)  !< @en Names of unit elements @enden @ja 単位要素の名前 @endja
    character(TOKEN):: offset     !< @en Offset string (e.g., "@273.15") @enden @ja オフセット文字列 @endja
    real(dp), pointer:: power(:)  !< @en Powers (exponents) of unit elements @enden @ja 単位要素のべき乗 @endja
  end type units
 
  interface clear
    module procedure dcunitsclear
  end interface
 
  interface deallocate
    module procedure dcunitsdeallocate
  end interface
 
  interface assignment(=)
    module procedure dcunitsbuild, dcunitstostring
  end interface
 
  interface operator(*)
    module procedure dcunitsmul
  end interface
 
  interface operator(/)
    module procedure dcunitsdiv
  end interface
 
  interface operator(+)
    module procedure dcunitsadd
  end interface
 
contains
 
  !> @en
  !> @brief Simplify unit representation by combining like terms
  !> @details
  !> Internal subroutine that combines unit elements with the same name
  !> by adding their powers.
  !> @param[in,out] u UNITS variable to simplify
  !> @param[in] name Array of unit element names
  !> @param[in] power Array of unit element powers
  !> @enden
  !>
  !> @ja
  !> @brief 同類項を結合して単位表現を簡約化
  !> @details
  !> 同じ名前の単位要素のべき乗を加算して結合する内部サブルーチン.
  !> @param[in,out] u 簡約化する UNITS 変数
  !> @param[in] name 単位要素の名前の配列
  !> @param[in] power 単位要素のべき乗の配列
  !> @endja
  subroutine units_simplify(u, name, power)
    type(UNITS), intent(inout):: u
    character(*), intent(in):: name(u%nelems)
    real(DP), intent(in):: power(u%nelems)
    integer:: i, n, j, onazi
    integer:: table(u%nelems)
 
    if (u%nelems < 1) return
    table(:) = 0
    n = 0
    do, i = 1, u%nelems
      if (name(i) == '') cycle
      onazi = 0
      do, j = 1, i - 1
        if (name(j) == name(i)) then
          onazi = j
        endif
      enddo
      if (onazi > 0) then
        table(i) = table(onazi)
      else
        n = n + 1
        table(i) = n
      endif
    enddo
    allocate(u%name(n), u%power(n))
    u%power = 0.0_dp
    do, i = 1, u%nelems
      if (table(i) == 0) cycle
      u%name(table(i)) = name(i)
      u%power(table(i)) = u%power(table(i)) + power(i)
    enddo
    u%nelems = n
  end subroutine units_simplify
 
  !> @en
  !> @brief Multiply two UNITS
  !> @param[in] u1 First UNITS variable
  !> @param[in] u2 Second UNITS variable
  !> @return Product of two UNITS
  !> @enden
  !>
  !> @ja
  !> @brief 2つの UNITS を乗算
  !> @param[in] u1 1番目の UNITS 変数
  !> @param[in] u2 2番目の UNITS 変数
  !> @return 2つの UNITS の積
  !> @endja
  type(units) function dcunitsmul(u1, u2) result(result)
    type(units), intent(in):: u1, u2
    integer:: n
    character(TOKEN), allocatable:: name(:)
    real(dp), allocatable:: power(:)
    result%factor = u1%factor * u2%factor
    result%nelems = u1%nelems + u2%nelems
    result%offset = ""
    n = result%nelems
    if (n == 0) then
      nullify(result%name, result%power)
      return
    endif
    allocate(name(n), power(n))
    name = (/u1%name, u2%name/)
    power = (/u1%power, u2%power/)
    call units_simplify(result, name, power)
    deallocate(name, power)
  end function dcunitsmul
 
  !> @en
  !> @brief Divide two UNITS
  !> @param[in] u1 Dividend UNITS variable
  !> @param[in] u2 Divisor UNITS variable
  !> @return Quotient of two UNITS
  !> @enden
  !>
  !> @ja
  !> @brief 2つの UNITS を除算
  !> @param[in] u1 被除数の UNITS 変数
  !> @param[in] u2 除数の UNITS 変数
  !> @return 2つの UNITS の商
  !> @endja
  type(units) function dcunitsdiv(u1, u2) result(result)
    type(units), intent(in):: u1, u2
    integer:: n, n1
    character(TOKEN), allocatable:: name(:)
    real(dp), allocatable:: power(:)
    if (abs(u2%factor) < tiny(u2%factor)) then
      result%factor = sign(u1%factor, 1.0_dp) * &
        & sign(u2%factor, 1.0_dp) * &
        & huge(1.0_dp)
    else
      result%factor = u1%factor / u2%factor
    endif
    result%nelems = u1%nelems + u2%nelems
    result%offset = ""
    n = result%nelems
    if (n == 0) then
      nullify(result%name, result%power)
      return
    endif
    allocate(name(n), power(n))
    n1 = u1%nelems
    if (n1 >= 1) then
      name(1:n1) = u1%name(1:n1)
      power(1:n1) = u1%power(1:n1)
    endif
    n1 = n1 + 1
    if (n >= n1) then
      name(n1:n) = u2%name(1:u2%nelems)
      power(n1:n) = -u2%power(1:u2%nelems)
    endif
    call units_simplify(result, name, power)
    deallocate(name, power)
  end function dcunitsdiv
 
  !> @en
  !> @brief Add two UNITS (with dimension compatibility check)
  !> @details
  !> Addition is only valid when both units have the same dimensions.
  !> If dimensions don't match, result%offset is set to "MISMATCH".
  !> @param[in] u1 First UNITS variable
  !> @param[in] u2 Second UNITS variable
  !> @return Sum of two UNITS (or error if incompatible)
  !> @enden
  !>
  !> @ja
  !> @brief 2つの UNITS を加算 (次元互換性チェック付き)
  !> @details
  !> 加算は両単位が同じ次元を持つ場合のみ有効です.
  !> 次元が一致しない場合、result%offset に "MISMATCH" が設定されます.
  !> @param[in] u1 1番目の UNITS 変数
  !> @param[in] u2 2番目の UNITS 変数
  !> @return 2つの UNITS の和 (非互換の場合はエラー)
  !> @endja
  type(units) function dcunitsadd(u1, u2) result(result)
    type(units), intent(in):: u1, u2
    type(units):: x
    result%offset = u1%offset
    result%nelems = u1%nelems
    result%factor = u1%factor + u2%factor
    x = u1 / u2
    if (x%nelems == 0) then
      nullify(result%name, result%power)   
      return
    endif
    if (all(abs(x%power(1:result%nelems)) < tiny(0.0_dp))) then
      allocate(result%name(result%nelems), result%power(result%nelems))
      result%name = u1%name
      result%power = u1%power
      return
    endif
    result%factor = 0.0
    result%nelems = -1
    result%offset = "MISMATCH"
    nullify(result%name, result%power)   
  end function dcunitsadd
 
  !> @en
  !> @brief Check if two UNITS are compatible for addition
  !> @details
  !> Two units are compatible if they have the same dimensions
  !> (i.e., the same unit elements with the same powers).
  !> @param[in] u1 First UNITS variable
  !> @param[in] u2 Second UNITS variable
  !> @return .true. if compatible, .false. otherwise
  !> @enden
  !>
  !> @ja
  !> @brief 2つの UNITS が加算互換かチェック
  !> @details
  !> 2つの単位は同じ次元を持つ場合（すなわち、同じ単位要素と
  !> 同じべき乗を持つ場合）に互換となります.
  !> @param[in] u1 1番目の UNITS 変数
  !> @param[in] u2 2番目の UNITS 変数
  !> @return 互換なら .true., そうでなければ .false.
  !> @endja
  logical function add_okay(u1, u2) result(result)
    type(units), intent(in):: u1, u2
    type(units):: x
    character(STRING):: debug
    call clear(x)
    x = u1 / u2
    debug = u1
    debug = u2
    debug = x
    if (x%nelems == 0) then
      result = .true.
    else if (all(abs(x%power(1:x%nelems)) < tiny(0.0_dp))) then
      result = .true.
    else
      result = .false.
    endif
    call deallocate(x)
  end function add_okay
 
  !> @en
  !> @brief Initialize UNITS variable
  !> @details
  !> Nullifies pointers and sets factor to 1.0 and nelems to 0.
  !> @param[in,out] u UNITS variable to initialize
  !> @enden
  !>
  !> @ja
  !> @brief UNITS 変数を初期化
  !> @details
  !> ポインタを nullify し、factor を 1.0、nelems を 0 に設定します.
  !> @param[in,out] u 初期化する UNITS 変数
  !> @endja
  subroutine dcunitsclear(u)
    type(units), intent(inout):: u
    nullify(u%name)
    nullify(u%power)
    u%factor = 1.0_dp
    u%offset = ""
    u%nelems = 0
  end subroutine dcunitsclear
 
  !> @en
  !> @brief Deallocate UNITS variable
  !> @details
  !> Deallocates allocated arrays and resets to initial state.
  !> @param[in,out] u UNITS variable to deallocate
  !> @enden
  !>
  !> @ja
  !> @brief UNITS 変数を解放
  !> @details
  !> 割り当てられた配列を解放し、初期状態にリセットします.
  !> @param[in,out] u 解放する UNITS 変数
  !> @endja
  subroutine dcunitsdeallocate(u)
    type(units), intent(inout):: u
    if (associated(u%name)) deallocate(u%name)
    if (associated(u%power)) deallocate(u%power)
    u%factor = 1.0_dp
    u%offset = ""
    u%nelems = 0
  end subroutine dcunitsdeallocate
 
  !> @en
  !> @brief Convert UNITS to string representation
  !> @details
  !> Generates a string representation of the UNITS variable.
  !> For example, UNITS with name=["m","s"], power=[1,-2], factor=1
  !> produces "m s-2".
  !> @param[out] string Output string representation
  !> @param[in] u UNITS variable to convert
  !> @enden
  !>
  !> @ja
  !> @brief UNITS を文字列表現に変換
  !> @details
  !> UNITS 型変数から文字列表現を生成します.
  !> 例えば、name=["m","s"], power=[1,-2], factor=1 の UNITS は
  !> "m s-2" を生成します.
  !> @param[out] string 出力する文字列表現
  !> @param[in] u 変換する UNITS 変数
  !> @endja
  subroutine dcunitstostring(string, u)
    character(*), intent(out):: string
    type(units), intent(in):: u
    integer:: i, ip, npower
    character(TOKEN):: buffer
    character:: mul = '.'
    real(DP), parameter:: allowed = epsilon(1.0_dp) * 16.0
 
    if (u%nelems < 0) then
      string = 'error from ' // u%offset
      return 
    endif
 
    write(buffer, "(1pg20.12)") u%factor
    string = buffer
    if (u%nelems < 1) return
 
    if (abs(u%factor - 1.0) < allowed) then
      string = ""
    else if (abs(u%factor + 1.0) < allowed) then
      string = "-"
    endif
 
    ip = len_trim(string) + 1
    do, i = 1, u%nelems
      npower = nint(u%power(i))
      if (abs(1.0 - u%power(i)) < allowed) then
        buffer = " "
      else if (abs(npower - u%power(i)) < allowed) then
        write(buffer, "(i10)") npower
        buffer = adjustl(buffer)
      else
        write(buffer, "(1pg10.3)") u%power(i)
        buffer = adjustl(buffer)
      endif
      if (buffer == '0') cycle
      string = trim(string) // mul // trim(u%name(i)) // trim(buffer)
    enddo
    if (ip <= len(string)) string(ip:ip) = ' '
    if (string(1:1) == " ") string = adjustl(string)
    if (u%offset /= "") then
      string = trim(string) // '@' // trim(u%offset)
    endif
  end subroutine dcunitstostring
 
  !> @en
  !> @brief Build UNITS from string (constructor)
  !> @details
  !> Parses a unit string and builds a UNITS variable.
  !> Supports multiplication (.), division (/), exponents (**),
  !> parentheses, and offset notation (@).
  !>
  !> Examples:
  !> - "m/s" -> meters per second
  !> - "kg.m.s-2" -> Newton (kilogram meter per second squared)
  !> - "K@273.15" -> Kelvin with offset
  !>
  !> @param[out] u UNITS variable to build
  !> @param[in] cunits Unit string to parse
  !> @enden
  !>
  !> @ja
  !> @brief 文字列から UNITS を構築 (コンストラクタ)
  !> @details
  !> 単位文字列をパースして UNITS 変数を構築します.
  !> 乗算 (.)、除算 (/)、べき乗 (**)、括弧、オフセット表記 (@)
  !> をサポートします.
  !>
  !> 例:
  !> - "m/s" -> メートル毎秒
  !> - "kg.m.s-2" -> ニュートン (キログラム・メートル毎秒毎秒)
  !> - "K@273.15" -> オフセット付きケルビン
  !>
  !> @param[out] u 構築する UNITS 変数
  !> @param[in] cunits パースする単位文字列
  !> @endja
  subroutine dcunitsbuild(u, cunits)
    use dcunits_com
    type(units), intent(out):: u
    character(*), intent(in):: cunits
 
    ! 構築中の情報、乗算対象の列として保持する。
    ! これは shift オペレータ付き単位を乗算できないことを示す。
    type elem_units
      character(TOKEN):: name
      real(DP):: power, factor
    end type elem_units
    type(elem_units), target:: ustack(100)
    integer:: ui = 1
 
    ! 構文単位が占める乗算対象の stack における最初の添字
    type paren_t
      real(DP):: factor
      integer:: factor_exp
      logical:: factor_inv
      integer:: power_exp
      integer:: paren_exp
    end type paren_t
    type(paren_t):: pstack(50)
    integer:: pi = 1
 
    ! パーサの状態遷移
    integer, parameter:: Y_INIT = 1, y_number = 2, y_name = 3, &
      & y_nx = 4, y_ni = 5, y_mul = 6, y_shift = 7
    integer:: yparse_status = y_init
 
    ! 字句
    integer:: ltype
    integer:: ivalue(5)
    real(DP):: dvalue
    character(TOKEN):: cvalue
    ! その他
    integer:: i
        
    ! --- 実行部 ---
    ! 初期化
    if (associated(u%name)) deallocate(u%name)
    if (associated(u%power)) deallocate(u%power)
    u%nelems = 0
    u%offset = ""
    u%factor = 1.0_dp
    if (cunits == "") return
    call dcunitssetline(cunits)
    call ustack_clear
    call pstack_clear
    yparse_status = y_init
 
    do
      call dcunitsgettoken(ltype, ivalue, dvalue, cvalue)
      select case(ltype)
      case (s_integer)
        select case(yparse_status)
        case (y_init, y_mul)
          pstack(pi)%factor = pstack(pi)%factor * ivalue(1)
          yparse_status = y_number
        case (y_name, y_nx)
          i = pstack(pi)%power_exp
          ustack(i:ui)%power = ustack(i:ui)%power * ivalue(1) 
          call power_next
          yparse_status = y_ni
        case (y_shift)
          u%offset = cvalue
        case default
          call error
        end select
      case (s_real)
        select case(yparse_status)
        case (y_init, y_mul)
          pstack(pi)%factor = pstack(pi)%factor * dvalue
          yparse_status = y_number
        case (y_name, y_nx)
          i = pstack(pi)%power_exp
          ustack(i:ui)%power = ustack(i:ui)%power * dvalue
          call power_next
          yparse_status = y_ni
        case (y_shift)
          u%offset = cvalue
        case default
          call error
        end select
      case (s_text)
        select case(yparse_status)
        case (y_init, y_number, y_mul)
          ustack(ui)%name = cvalue
          yparse_status = y_name
        case (y_name, y_ni)
          call ustack_grow
          call power_next
          ustack(ui)%name = cvalue
          yparse_status = y_name
        case (y_shift)
          u%offset = cvalue
        case default
          call error
        end select
      case (s_exponent)
        select case(yparse_status)
        case (y_name)
          yparse_status = y_nx
        case default
          call error
        end select
      case (s_multiply)
        select case(yparse_status)
        case (y_number, y_name)
          call factor_next
          yparse_status = y_mul
        case default
          call error
        end select
      case (s_divide)
        select case(yparse_status)
        case (y_number, y_name)
          call factor_next
          pstack(pi)%factor_inv = .true.
          yparse_status = y_mul
        case default
          call error
        end select
      case (s_shift)
        if (yparse_status == y_nx) call cancel_exp
        call units_finalize
        yparse_status = y_shift
      case (s_openpar)
        if (yparse_status == y_nx) call cancel_exp
        call pstack_push
      case (s_closepar)
        call units_finalize
        call pstack_pop
      case (s_eof)
        exit
      case default
        call error
      end select
    enddo
    if (yparse_status == y_nx) call cancel_exp
    call units_finalize
 
    u%nelems = ui
    u%factor = product(ustack(1:ui)%factor)
    call units_simplify(u, ustack(1:ui)%name, ustack(1:ui)%power)
 
  contains
 
    subroutine cancel_exp
      print *, "DCUnitsBuild: syntax error, operator(**) ignored"
    end subroutine cancel_exp
 
    subroutine error
      print *, "DCUnitsBuild: unexpected token <", &
        & trim(cvalue), "> ignored"
    end subroutine error
 
    subroutine power_next
      ! power_exp の終了処理
      call ustack_grow
      pstack(pi)%power_exp = ui
    end subroutine power_next
 
    subroutine factor_next
      ! factor_exp の終了処理
      real(DP):: factor
      i = pstack(pi)%factor_exp
      factor = product(ustack(i:ui)%factor) * pstack(pi)%factor
      if (pstack(pi)%factor_inv) then
        ustack(i:ui)%power = -ustack(i:ui)%power
        factor = 1.0_dp / factor
      endif
      ustack(i)%factor = factor
      ustack(i+1:ui)%factor = 1.0_dp
      call power_next
      pstack(pi)%factor = 1.0_dp
      pstack(pi)%factor_exp = ui
    end subroutine factor_next
 
    subroutine units_finalize
      call factor_next
    end subroutine units_finalize
 
    subroutine ustack_clear
      ui = 0
      call ustack_grow
    end subroutine ustack_clear
 
    subroutine ustack_grow
      if (ui >= size(ustack)) stop 'DCUnitsBuild: too many elements'
      ui = ui + 1
      ustack(ui)%name = ""
      ustack(ui)%factor = 1.0_dp
      ustack(ui)%power = 1.0_dp
    end subroutine ustack_grow
 
    subroutine pstack_clear
      pi = 0
      call pstack_push
    end subroutine pstack_clear
 
    subroutine pstack_push
      if (pi >= size(pstack)) stop 'DCUnitsBuild: too many parens'
      pi = pi + 1
      call ustack_grow
      pstack(pi)%factor_exp = ui
      pstack(pi)%factor = 1.0_dp
      pstack(pi)%factor_inv = .false.
      pstack(pi)%power_exp = ui
      pstack(pi)%paren_exp = ui
    end subroutine pstack_push
 
    subroutine pstack_pop
      ! factor_exp の終了処理
      call power_next
      pi = pi - 1
    end subroutine pstack_pop
 
  end subroutine dcunitsbuild
 
  !> @namespace dc_units
end module dc_units