phys/complex_number_module.F

   1 !------------------------------------------------------------------------!
   2 !  The Community Multiscale Air Quality (CMAQ) system software is in     !
   3 !  continuous development by various groups and is based on information  !
   4 !  from these groups: Federal Government employees, contractors working  !
   5 !  within a United States Government contract, and non-Federal sources   !
   6 !  including research institutions.  These groups give the Government    !
   7 !  permission to use, prepare derivative works of, and distribute copies !
   8 !  of their work in the CMAQ system to the public and to permit others   !
   9 !  to do so.  The United States Environmental Protection Agency          !
  10 !  therefore grants similar permission to use the CMAQ system software,  !
  11 !  but users are requested to provide copies of derivative works or      !
  12 !  products designed to operate in the CMAQ system to the United States  !
  13 !  Government without restrictions as to use by others.  Software        !
  14 !  that is used with the CMAQ system but distributed under the GNU       !
  15 !  General Public License or the GNU Lesser General Public License is    !
  16 !  subject to their copyright restrictions.                              !
  17 !------------------------------------------------------------------------!
  18
  19 ! complex number general function
  20
  21 ! Revision History:
  22 !   2012/07/31 David Wong Original version
  23 !   2012/10/22 David Wong Added treatment to avoid division by 0
  24 !   2013/11/20 David Wong modified the way to compute
  25 !              max(abs(c_div_cc%real_part), min_val) in subroutine
  26 !              c_div_cc to satisfy absoft compiler requirement
  27 !   2015/12/16 David Wong renamed argument list for c_add_cr, c_add_rc,
  28 !              c_sub_cr and c_sub_rc to avoid gfortran not able to
  29 !              distinguish those routines in the interface block issue
  30 !   2016/02/23 David Wong extracted the entire module and put it in a
  31 !              file alone.
  32
  33         module complex_number_module
  34
  35         implicit none
  36
  37 !       integer, parameter :: loc_real_precision = selected_real_kind(p=16, r=60)
  38         integer, parameter :: loc_real_precision = 8
  39
  40         real (kind=loc_real_precision), parameter, private :: min_val = 1.0e-30_loc_real_precision
  41
  42         type complex_number
  43           real(kind=loc_real_precision) :: real_part, imag_part
  44         end type complex_number
  45
  46         interface c_add
  47           module procedure c_add_cc,    &      ! z1 + z2
  48                            c_add_cr,    &      ! z1 + num, where num is a real number
  49                            c_add_rc            ! num + z1, where num is a real number
  50         end interface
  51
  52         interface c_sub
  53           module procedure c_sub_cc,    &      ! z1 - z2
  54                            c_sub_cr,    &      ! z1 - num, where num is a real number
  55                            c_sub_rc            ! num - z1, where num is a real number
  56         end interface
  57
  58         interface c_mul
  59           module procedure c_mul_cc,    &      ! z1 * z2
  60                            c_mul_rc            ! num * z1, where num is a real number
  61         end interface
  62
  63         interface c_div
  64           module procedure c_div_cc,    &      ! z1 / z2
  65                            c_div_rc            ! num / z1, where num is a real number
  66         end interface
  67
  68         contains
  69
  70 ! --------------------------------------------------------------------------
  71         type (complex_number) function c_set (x, y)
  72
  73 ! initialize a complex number
  74
  75           real(kind=loc_real_precision), intent(in) :: x, y
  76
  77           character (len = 80) :: str
  78
  79           write (str, *) x
  80           read(str, *) c_set%real_part
  81           write (str, *) y
  82           read(str, *) c_set%imag_part
  83
  84         end function c_set
  85
  86 ! --------------------------------------------------------------------------
  87         type (complex_number) function c_add_cc (z1, z2)
  88
  89           type (complex_number), intent(in) :: z1, z2
  90
  91           c_add_cc%real_part = z1%real_part + z2%real_part
  92           c_add_cc%imag_part = z1%imag_part + z2%imag_part
  93
  94         end function c_add_cc
  95
  96 ! --------------------------------------------------------------------------
  97         type (complex_number) function c_add_cr (z3, num1)
  98
  99           type (complex_number),         intent(in) :: z3
 100           real(kind=loc_real_precision), intent(in) :: num1
 101
 102           c_add_cr%real_part = z3%real_part + num1
 103           c_add_cr%imag_part = z3%imag_part
 104
 105         end function c_add_cr
 106
 107 ! --------------------------------------------------------------------------
 108         type (complex_number) function c_add_rc (num2, z4)
 109
 110           type (complex_number),         intent(in) :: z4
 111           real(kind=loc_real_precision), intent(in) :: num2
 112
 113           c_add_rc%real_part = z4%real_part + num2
 114           c_add_rc%imag_part = z4%imag_part
 115
 116         end function c_add_rc
 117
 118 ! --------------------------------------------------------------------------
 119         type (complex_number) function c_sub_cc (z1, z2)
 120
 121           type (complex_number), intent(in) :: z1, z2
 122
 123           c_sub_cc%real_part = z1%real_part - z2%real_part
 124           c_sub_cc%imag_part = z1%imag_part - z2%imag_part
 125
 126         end function c_sub_cc
 127
 128 ! --------------------------------------------------------------------------
 129         type (complex_number) function c_sub_cr (z3, num1)
 130
 131           type (complex_number),         intent(in) :: z3
 132           real(kind=loc_real_precision), intent(in) :: num1
 133
 134           c_sub_cr%real_part = z3%real_part - num1
 135           c_sub_cr%imag_part = z3%imag_part
 136
 137         end function c_sub_cr
 138
 139 ! --------------------------------------------------------------------------
 140         type (complex_number) function c_sub_rc (num2, z4)
 141
 142           type (complex_number),         intent(in) :: z4
 143           real(kind=loc_real_precision), intent(in) :: num2
 144
 145           c_sub_rc%real_part = num2 - z4%real_part
 146           c_sub_rc%imag_part = - z4%imag_part
 147
 148         end function c_sub_rc
 149
 150 ! --------------------------------------------------------------------------
 151         type (complex_number) function c_mul_cc (z1, z2)
 152
 153           type (complex_number), intent(in) :: z1, z2
 154
 155           c_mul_cc%real_part =   z1%real_part * z2%real_part    &
 156                                - z1%imag_part * z2%imag_part
 157           c_mul_cc%imag_part =   z1%real_part * z2%imag_part    &
 158                                + z1%imag_part * z2%real_part
 159
 160         end function c_mul_cc
 161
 162 ! --------------------------------------------------------------------------
 163         type (complex_number) function c_mul_rc (x, z1)
 164
 165           type (complex_number),         intent(in) :: z1
 166           real(kind=loc_real_precision), intent(in) :: x
 167
 168           c_mul_rc%real_part = z1%real_part * x
 169           c_mul_rc%imag_part = z1%imag_part * x
 170
 171         end function c_mul_rc
 172
 173 ! --------------------------------------------------------------------------
 174         type (complex_number) function c_div_cc (z1, z2)
 175
 176           type (complex_number), intent(in) :: z1, z2
 177
 178           real(kind=loc_real_precision) :: denom
 179           real(kind=loc_real_precision) :: temp(2)
 180
 181           denom = 1.0 / (  z2%real_part * z2%real_part &
 182                          + z2%imag_part * z2%imag_part)
 183
 184           c_div_cc%real_part = (  z1%real_part * z2%real_part          &
 185                                 + z1%imag_part * z2%imag_part) * denom
 186
 187           temp(1) = abs(c_div_cc%real_part)
 188           temp(2) = min_val
 189
 190           c_div_cc%real_part = sign(maxval(temp), c_div_cc%real_part)
 191
 192           c_div_cc%imag_part = (  z1%imag_part * z2%real_part          &
 193                                 - z1%real_part * z2%imag_part) * denom
 194
 195           temp(1) = abs(c_div_cc%imag_part)
 196
 197           c_div_cc%imag_part = sign(maxval(temp), c_div_cc%imag_part)
 198
 199         end function c_div_cc
 200
 201 ! --------------------------------------------------------------------------
 202         type (complex_number) function c_div_rc (num, z1)
 203
 204 ! compute 1 / z1
 205
 206           real(kind=loc_real_precision), intent(in) :: num
 207           type (complex_number),         intent(in) :: z1
 208
 209           real(kind=loc_real_precision) :: denom, temp
 210
 211           temp = z1%real_part * z1%real_part + z1%imag_part * z1%imag_part
 212           temp = sign(max(abs(temp), min_val), temp)
 213
 214           denom = num / temp
 215           c_div_rc%real_part = z1%real_part * denom
 216           c_div_rc%imag_part = -1.0 * z1%imag_part * denom
 217
 218         end function c_div_rc
 219
 220 ! --------------------------------------------------------------------------
 221         type (complex_number) function c_sin (z1)
 222
 223 ! compute sin of a complex number
 224
 225           type (complex_number), intent(in) :: z1
 226
 227           c_sin%real_part = sin(z1%real_part) * cosh(z1%imag_part)
 228           c_sin%imag_part = cos(z1%real_part) * sinh(z1%imag_part)
 229
 230         end function c_sin
 231
 232 ! --------------------------------------------------------------------------
 233         type (complex_number) function c_cos (z1)
 234
 235           type (complex_number), intent(in) :: z1
 236
 237           c_cos%real_part = cos(z1%real_part) * cosh(z1%imag_part)
 238           c_cos%imag_part = -1.0 * sin(z1%real_part) * sinh(z1%imag_part)
 239
 240         end function c_cos
 241
 242 ! --------------------------------------------------------------------------
 243         real(kind=loc_real_precision) function c_abs (z1)
 244
 245 ! computer absolute value of a complex number
 246
 247           type (complex_number), intent(in) :: z1
 248
 249           c_abs = sqrt(z1%real_part**2 + z1%imag_part**2)
 250
 251         end function c_abs
 252
 253         end module complex_number_module