fstr_frequency_analysis.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
module fstr_frequency_analysis
 
  use m_fstr
  use m_fstr_creatematrix_and_dampingforce
  use m_fstr_addbc
  use m_fstr_spring
  use m_fstr_eig_setmass
  use m_fstr_nodalkinematics
  use m_fstr_update
  use m_fstr_residual, only: fstr_update_reaction_spc
  use m_fstr_nodalstress
  use m_make_result
  use m_hecmw2fstr_mesh_conv
  use hecmw_setup_util, only: node_global_to_local
  use hecmw_ebc_defer
  use m_solve_lineq !for the quasi-static displacement of !FACCELERATION
 
  implicit none
 
contains
 
  subroutine fstr_solve_frequency_analysis(hecMESH, hecMAT, fstrSOLID, fstrEIG,  &
      fstrDYNAMIC, fstrPARAM, fstrFREQ, hecLagMAT)
    !C
    !C-- global variable
    !C
    type(hecmwst_local_mesh)             :: hecMESH
    type(hecmwst_matrix)                 :: hecMAT
    type(fstr_eigen)                     :: fstrEIG
    type(fstr_solid)                     :: fstrSOLID
    type(fstr_param)                     :: fstrPARAM
    type(fstr_dynamic)                   :: fstrDYNAMIC
    type(fstr_freqanalysis)              :: fstrFREQ
    type(hecmwst_matrix_lagrange)        :: hecLagMAT
 
    !C
    !C-- local variable
    !C
    integer(kind=kint)            :: numnode, numelm, startmode, endmode, nummode, ndof, im, vistype
    integer(kind=kint)            :: numfreq, idnode, numdisp, imonit, nmonit, monitnode(1)
    integer(kind=kint)            :: freqiout(3)
    integer(kind=kint)            :: ierr
    integer(kind=kint), parameter :: ilogin = 9056
    real(kind=kreal), allocatable :: eigenvalue(:), loadvecre(:), loadvecim(:)
    real(kind=kreal), allocatable :: bjre(:), bjim(:)
    real(kind=kreal), allocatable :: dispre(:), dispim(:), velre(:), velim(:), accre(:), accim(:)
    real(kind=kreal), allocatable :: usre(:), usim(:)           !quasi-static displacement of !FACCELERATION
    real(kind=kreal), allocatable :: totloadre(:), totloadim(:) !!FLOAD plus the load of !FACCELERATION
    real(kind=kreal)              :: freq, omega, val, dx, dy, dz, f_start, f_end
    real(kind=kreal)              :: t_start, t_end, time, dxi, dyi, dzi
    complex(kind=kreal)           :: eiwt  !exp(i*omega*time)
    logical                       :: isacc !.true. if !FACCELERATION is given
    type(fstr_freqanalysis_data)  :: freqData
 
    numnode   = hecmesh%n_node
    numelm    = hecmesh%n_elem
    ndof      = hecmesh%n_dof
    startmode = fstrfreq%start_mode
    endmode   = fstrfreq%end_mode
    nummode   = endmode - startmode +1
 
    freqdata%numMode    = nummode
    freqdata%numNodeDOF = numnode*ndof
    allocate(freqdata%eigOmega(nummode))
    allocate(freqdata%eigVector(numnode*ndof, nummode))
 
    call setupfreqparam(fstrdynamic, f_start, f_end, numfreq, freqdata%rayAlpha, freqdata%rayBeta, idnode, vistype, freqiout)
    call setupoutputitem(fstrsolid, freqiout)
 
    monitnode(1) = idnode
    imonit       = 0
    if( node_global_to_local(hecmesh, monitnode, 1) == 1 ) then
      if( monitnode(1) <= hecmesh%nn_internal ) imonit = monitnode(1)  !monitor is reported by the rank owning it
    end if
    nmonit = 0
    if( imonit > 0 ) nmonit = 1
    call hecmw_allreduce_i1(hecmesh, nmonit, hecmw_sum)
    if( nmonit == 0 ) then
      if( myrank == 0 ) then
        write(*,   *) "Error: monitor node not found:", idnode
        write(ilog,*) "Error: monitor node not found:", idnode
      end if
      call hecmw_abort( hecmw_comm_get_comm() )
    end if
 
    if( myrank == 0 ) then
      write(*,*) "Rayleigh alpha:", freqdata%rayAlpha
      write(*,*) "Rayleigh beta:", freqdata%rayBeta
      write(ilog,*) "Rayleigh alpha:", freqdata%rayAlpha
      write(ilog,*) "Rayleigh beta:", freqdata%rayBeta
    end if
 
 
    allocate(eigenvalue(nummode))
    allocate(loadvecre(numnode*ndof))
    allocate(loadvecim(numnode*ndof))
    allocate(bjre(nummode))
    allocate(bjim(nummode))
    allocate(dispre(numnode*ndof))
    allocate(dispim(numnode*ndof))
    allocate(velre(numnode*ndof))
    allocate(velim(numnode*ndof))
    allocate(accre(numnode*ndof))
    allocate(accim(numnode*ndof))
    allocate(usre(numnode*ndof))
    allocate(usim(numnode*ndof))
    allocate(totloadre(numnode*ndof))
    allocate(totloadim(numnode*ndof))
 
    loadvecre(:) = 0.0d0
    loadvecim(:) = 0.0d0
    usre(:)      = 0.0d0
    usim(:)      = 0.0d0
 
    if( myrank == 0 ) then
      write(*,*) "--frequency analysis--"
      write(*,   *) "read from=", trim(fstrfreq%eigenlog_filename)
      write(ilog,*) "read from=", trim(fstrfreq%eigenlog_filename)
      write(*,   *) "start mode=", startmode
      write(ilog,*) "start mode=", startmode
      write(*,   *) "end mode=", endmode
      write(ilog,*) "end mode=", endmode
    end if
    open(unit=ilogin, file=trim(fstrfreq%eigenlog_filename), status="OLD", action="READ", iostat=ierr)
    if( ierr /= 0 ) then
      if( myrank == 0 ) write(*,*) "Error: cannot open eigenlog file: ", trim(fstrfreq%eigenlog_filename)
      call hecmw_abort( hecmw_comm_get_comm() )
    endif
    call read_eigen_values(ilogin, startmode, endmode, eigenvalue, freqdata%eigOmega)
    !call read_eigen_vector(ilogin, startmode, endmode, ndof, numnode, freqData%eigVector)
    close(ilogin)
 
    call read_eigen_vector_res(hecmesh, startmode, endmode, ndof, numnode, freqdata%eigVector)
 
    call extract_surf2node(hecmesh, fstrfreq, loadvecre, loadvecim)
    call assemble_nodeload(hecmesh, fstrfreq, ndof, loadvecre, loadvecim)
 
    !acceleration input (relative motion method)
    !  the displacement is u = -us/omega^2 + (relative response by the modes), where us is
    !  the quasi-static displacement with the acceleration amplitude given as the displacement
    !  of the excited DOFs.
    !  without !FACCELERATION, the analysis follows the same path as before.
    isacc = fstrfreq%FACC_ngrp_tot > 0
    if( isacc ) then
      !the modes must be fixed at the excited DOFs
      call checkaccelerationdof(hecmesh, fstrfreq, freqdata)
      if( myrank == 0 ) write(*,*) "calc quasi-static displacement of acceleration input"
      !us for the real part (LOAD CASE=1) and the imaginary part (LOAD CASE=2) of the amplitude
      call calcquasistaticdisp(hecmesh, hecmat, fstrsolid, fstrparam, fstrfreq, heclagmat, kfloadcase_re, usre)
      call calcquasistaticdisp(hecmesh, hecmat, fstrsolid, fstrparam, fstrfreq, heclagmat, kfloadcase_im, usim)
    end if
 
    if( myrank == 0 ) write(*,*) "calc mass matrix"
    call calcmassmatrix(fstrparam, hecmesh, hecmat, fstrsolid, fstreig, heclagmat)
    if( myrank == 0 ) write(*,*) "scale eigenvector"
    call scaleeigenvector(hecmesh, fstreig, ndof*numnode, nummode, freqdata%eigVector)
 
    if( myrank == 0 ) then
      write(*,   *) "start frequency:", f_start
      write(ilog,*) "start frequency:", f_start
      write(*,   *) "end frequency:", f_end
      write(ilog,*) "end frequency:", f_end
      write(*   ,*) "number of the sampling points", numfreq
      write(ilog,*) "number of the sampling points", numfreq
      write(*   ,*) "monitor nodeid=", idnode
      write(ilog,*) "monitor nodeid=", idnode
    end if
 
    do im=1, numfreq
      freq = (f_end-f_start)/dble(numfreq)*dble(im) + f_start
      omega = 2.0d0 * 3.14159265358979d0 * freq
 
      !load of the relative response: !FLOAD plus the load of !FACCELERATION, which depends on omega
      totloadre(:) = loadvecre(:)
      totloadim(:) = loadvecim(:)
      if( isacc ) call addaccelerationload(fstreig, freqdata, omega, usre, usim, totloadre, totloadim)
      call calcfreqcoeff(hecmesh, freqdata, totloadre, totloadim, omega, bjre, bjim)
      call calcdispvector(freqdata, bjre, bjim, dispre, dispim)
      call calcvelvector(freqdata, omega, bjre, bjim, velre, velim)
      call calcaccvector(freqdata, omega, bjre, bjim, accre, accim)
      !absolute response: add the quasi-static part
      !  disp: -us/omega^2, vel: -i*us/omega, acc: us
      if( isacc ) then
        call addquasistaticvector(dcmplx(-1.0d0/omega**2, 0.0d0), usre, usim, dispre, dispim)
        call addquasistaticvector(dcmplx(0.0d0, -1.0d0/omega), usre, usim, velre, velim)
        call addquasistaticvector(dcmplx(1.0d0, 0.0d0), usre, usim, accre, accim)
      end if
 
      if( imonit > 0 ) then
        dx  = sqrt(dispre(3*(imonit-1)+1)**2 + dispim(3*(imonit-1)+1)**2)
        dy  = sqrt(dispre(3*(imonit-1)+2)**2 + dispim(3*(imonit-1)+2)**2)
        dz  = sqrt(dispre(3*(imonit-1)+3)**2 + dispim(3*(imonit-1)+3)**2)
        val = sqrt(dx**2 + dy**2 + dz**2)
        write(*,    *) freq, "[Hz] : ", val
        write(ilog, *) freq, "[Hz] : ", val
      end if
 
      if( myrank == 0 ) then
        if(fstrparam%fg_result==1) then
          write(*,   *) freq, "[Hz] : ", im, ".res"
          write(ilog,*) freq, "[Hz] : ", im, ".res"
        end if
        if(fstrparam%fg_visual==1 .and. vistype==1) then
          write(*,   *) freq, "[Hz] : ", im, ".vis"
          write(ilog,*) freq, "[Hz] : ", im, ".vis"
        end if
      end if
      call output_result(hecmesh, hecmat, fstrsolid, fstrdynamic, fstrparam, &
        'fstrRES', 'frequency', freq, im, &
        fstrparam%fg_result==1, fstrparam%fg_visual==1 .and. vistype==1, &
        dispre, dispim, velre, velim, accre, accim)
    end do
 
    call setupdynaparam(fstrdynamic, t_start, t_end, freq, numdisp)
    if( myrank == 0 ) then
      write(*,   *) "start time:", t_start
      write(ilog,*) "start time:", t_start
      write(*,   *) "end time:", t_end
      write(ilog,*) "end time:", t_end
      write(*,   *) "frequency:", freq
      write(ilog,*) "frequency:", freq
      write(*,   *) "node id:", idnode
      write(ilog,*) "node id:", idnode
      write(*,   *) "num disp:", numdisp
      write(ilog,*) "num disp:", numdisp
    end if
 
    !the quasi-static part of the displacement and the velocity is divided by omega
    if( isacc .and. freq <= 0.0d0 ) then
      if( myrank == 0 ) then
        write(*,   *) "Error: !FACCELERATION requires a positive frequency for the time history output"
        write(ilog,*) "Error: !FACCELERATION requires a positive frequency for the time history output"
      end if
      call hecmw_abort( hecmw_comm_get_comm() )
    end if
 
    omega = 2.0d0 * 3.14159265358979d0 * freq
    !load of the relative response at the frequency of the time history output
    totloadre(:) = loadvecre(:)
    totloadim(:) = loadvecim(:)
    if( isacc ) call addaccelerationload(fstreig, freqdata, omega, usre, usim, totloadre, totloadim)
    call calcfreqcoeff(hecmesh, freqdata, totloadre, totloadim, omega, bjre, bjim)
 
    do im=1, numdisp
      time = (t_end-t_start)/dble(numdisp)*dble(im-1) + t_start
      call calcdispvectortime(freqdata, time, omega, bjre, bjim, dispre, dispim)
      call calcvelvectortime(freqdata, time, omega, bjre, bjim, velre, velim)
      call calcaccvectortime(freqdata, time, omega, bjre, bjim, accre, accim)
      !absolute response: add the quasi-static part multiplied by exp(i*omega*time)
      if( isacc ) then
        eiwt = exp(dcmplx(0.0d0, omega*time))
        call addquasistaticvector(eiwt*dcmplx(-1.0d0/omega**2, 0.0d0), usre, usim, dispre, dispim)
        call addquasistaticvector(eiwt*dcmplx(0.0d0, -1.0d0/omega), usre, usim, velre, velim)
        call addquasistaticvector(eiwt, usre, usim, accre, accim)
      end if
      if( imonit > 0 ) then
        dx  = dispre(3*(imonit-1)+1)
        dy  = dispre(3*(imonit-1)+2)
        dz  = dispre(3*(imonit-1)+3)
        dxi = dispim(3*(imonit-1)+1)
        dyi = dispim(3*(imonit-1)+2)
        dzi = dispim(3*(imonit-1)+3)
      end if
 
      if( myrank == 0 ) then
        if(fstrparam%fg_result==1) then
          write(*,   *) "time=", time, " : ", im, ".res"
          write(ilog,*) "time=", time, " : ", im, ".res"
        end if
        if(fstrparam%fg_visual==1 .and. vistype==2) then
          write(*,   *) "time=", time, " : ", im, ".vis"
          write(ilog,*) "time=", time, " : ", im, ".vis"
        end if
      end if
      call output_result(hecmesh, hecmat, fstrsolid, fstrdynamic, fstrparam, &
        'fstrDYNA', 'TOTALTIME', time, im, &
        fstrparam%fg_result==1, fstrparam%fg_visual==1 .and. vistype==2, &
        dispre, dispim, velre, velim, accre, accim)
    end do
 
    deallocate(freqdata%eigOmega)
    deallocate(freqdata%eigVector)
    deallocate(eigenvalue)
    deallocate(loadvecre)
    deallocate(loadvecim)
    deallocate(bjre)
    deallocate(bjim)
    deallocate(dispre)
    deallocate(dispim)
    deallocate(velre)
    deallocate(velim)
    deallocate(accre)
    deallocate(accim)
    deallocate(usre)
    deallocate(usim)
    deallocate(totloadre)
    deallocate(totloadim)
 
  end subroutine
 
  subroutine read_eigen_values(logfile, startmode, endmode, eigenvalue, anglfreq)
    !---- args
    integer(kind=kint), intent(in)  :: logfile
    integer(kind=kint), intent(in)  :: startmode
    integer(kind=kint), intent(in)  :: endmode
    real(kind=kreal), intent(inout) :: eigenvalue(:) !intend(endmode-startmode+1)
    real(kind=kreal), intent(inout) :: anglfreq(:)
    !---- vals
    integer(kind=kint)           :: im, endflag, id
    character(len=HECMW_MSG_LEN) :: line
    real(kind=kreal)             :: freq
    !---- body
 
    rewind(logfile)
    endflag = 0
    !Find eigenvalue header.
    do
      read(logfile, '(A80)', err=119) line
      if(trim(adjustl(line)) == "NO.  EIGENVALUE  FREQUENCY   (HZ)        X           Y           Z           X") then
        endflag = 1
        exit
      end if
    end do
    read(logfile, '(A80)') line
    !read eigenvalue
    do im=1, startmode-1
      read(logfile, '(A80)') line
    end do
    do im=1, (endmode-startmode+1)
      read(logfile, '(i5,3e12.4,a)', err=119) id, eigenvalue(im), anglfreq(im), freq, line
    end do
    return
 
    !error handling
    119 write(*,*) "Error to find eigenvalue information from logfile"
    write(ilog,*) "Error to find eigenvalue information from logfile"
    stop
  end subroutine
 
  subroutine read_eigen_vector(logfile, startmode, endmode, numdof, numnode, eigenvector)
    !---- args
    integer(kind=kint), intent(in)  :: logfile
    integer(kind=kint), intent(in)  :: startmode
    integer(kind=kint), intent(in)  :: endmode
    integer(kind=kint), intent(in)  :: numdof
    integer(kind=kint), intent(in)  :: numnode
    real(kind=kreal), intent(inout) :: eigenvector(:, :) !intend (numdof*NN,nmode)
    !---- vals
    integer(kind=kint)           :: im, in, gblid, j, idx
    real(kind=kreal)             :: vec(6)
    character(len=HECMW_MSG_LEN) :: line
    !---- body
 
    rewind(logfile)
    !Find first eigenvector header
    do im=1, startmode-1
      do
        read(logfile, '(a80)', err=119, end=119) line
        if(line(1:9) == " Mode No.") then
          exit ! goto next mode
        end if
      end do
    end do
 
    !read eigenvector
    do im=1, (endmode-startmode+1)
      !find header
      do
        read(logfile, '(a80)', err=119, end=119) line
        if(line(1:9) == " Mode No.") then
          exit !find eigenmode
        end if
      end do
 
      read(logfile, '(a80)', err=119) line
      read(logfile, '(a80)', err=119) line
      read(logfile, '(a80)', err=119) line
      read(logfile, '(a80)', err=119) line
      !read eigenvector component
      do in=1, numnode
        select case(numdof)
          case(2)
            read(logfile, '(i10,2e12.4)', err=119) gblid, (vec(j), j=1,2)
 
          case(3)
            !read(logfile, '(i10,3e12.4)', ERR=119) gblid, (vec(j), j=1,3)
            read(logfile, '(i10,3e16.8)', err=119) gblid, (vec(j), j=1,3)
 
          case(6)
            read(logfile, '(i10,6e12.4)', err=119) gblid, (vec(j), j=1,6)
          case default
            !error
            goto 119
        end select
 
        do j=1, numdof
          idx = (in-1)*numdof + j
          eigenvector(idx,im) = vec(j)
        end do
      end do
    end do
    return
 
    !error handling
    119 write(*,*) "Error to find eigenvector from logfile"
    write(ilog,*) "Error to find eigenvector from logfile"
    stop
 
  end subroutine
 
  subroutine read_eigen_vector_res(hecMESH, startmode, endmode, numdof, numnode, eigenvector)
    !---- args
    type(hecmwst_local_mesh), intent(in) :: hecMESH
    integer(kind=kint), intent(in)       :: startmode
    integer(kind=kint), intent(in)       :: endmode
    integer(kind=kint), intent(in)       :: numdof
    integer(kind=kint), intent(in)       :: numnode
    real(kind=kreal), intent(inout)      :: eigenvector(:, :) !intend (numdof*NN,nmode)
    !---- vals
    integer(kind=kint), parameter :: compidx = 1 !Component index of displacement
    integer(kind=kint)            :: imode, idx, ind, a, b, nallcomp, j
    type(hecmwst_result_data)     :: eigenres
    character(len=HECMW_NAME_LEN) :: name
    !---- body
 
    name = 'result-in'
    do imode=startmode, endmode
      call nullify_result_data(eigenres)
      call hecmw_result_read_by_name(hecmesh, name, imode, eigenres)
 
      nallcomp = 0
      do ind=1,eigenres%nn_component
        nallcomp = nallcomp + eigenres%nn_dof(ind)
      end do
 
      idx = imode - startmode + 1
      do ind=1, numnode
        do j=1, numdof
          a = (ind-1)*nallcomp + j  !src vector index
          b = (ind-1)*numdof   + j
          eigenvector(b,idx) = eigenres%node_val_item(a)
        end do
      end do
      call free_result_data(eigenres)
    end do
 
  contains
 
    subroutine free_result_data(res)
      !---- args
      type(hecmwst_result_data) :: res
      !---- vals
      !---- body
      if(associated(res%nn_dof))        deallocate(res%nn_dof)
      if(associated(res%ne_dof))        deallocate(res%ne_dof)
      if(associated(res%node_label))    deallocate(res%node_label)
      if(associated(res%elem_label))    deallocate(res%elem_label)
      if(associated(res%node_val_item)) deallocate(res%node_val_item)
      if(associated(res%elem_val_item))  deallocate(res%elem_val_item)
    end subroutine
 
    subroutine nullify_result_data(res)
      !---- args
      type(hecmwst_result_data) :: res
      !---- vals
      !---- body
      nullify(res%nn_dof)
      nullify(res%ne_dof)
      nullify(res%node_label)
      nullify(res%elem_label)
      nullify(res%node_val_item)
      nullify(res%elem_val_item)
    end subroutine
 
 
  end subroutine
 
  !
  ! The displacement, the velocity and the acceleration of the frequency response
  ! analysis are selected by iout_list of !DYNAMIC, and the other items follow
  ! !OUTPUT_RES as in the static and the transient dynamic analysis
  !
  subroutine setupoutputitem(fstrSOLID, iout)
    !---- args
    type(fstr_solid), intent(inout) :: fstrSOLID
    integer(kind=kint), intent(in)  :: iout(3)
    !---- vals
    integer(kind=kint) :: i
    !---- body
 
    do i=3, 4
      fstrsolid%output_ctrl(i)%outinfo%on(1)  = (iout(1) == 1)
      fstrsolid%output_ctrl(i)%outinfo%on(15) = (iout(2) == 1)
      fstrsolid%output_ctrl(i)%outinfo%on(16) = (iout(3) == 1)
    end do
  end subroutine
 
  !
  ! Set one phase of the response to fstrSOLID and fstrDYNAMIC, and recover the
  ! strain and the stress of it with the routines of the static analysis.
  ! The incremental form keeps the incompatible mode of the 361 element consistent
  ! with the displacement, as the material is linear elastic here.
  !
  subroutine updatesolidresult(hecMESH, hecMAT, fstrSOLID, fstrDYNAMIC, fstrPARAM, istp, disp, vel, acc)
    !---- args
    type(hecmwst_local_mesh)       :: hecMESH
    type(hecmwst_matrix)           :: hecMAT
    type(fstr_solid)               :: fstrSOLID
    type(fstr_dynamic)             :: fstrDYNAMIC
    type(fstr_param)               :: fstrPARAM
    integer(kind=kint), intent(in) :: istp
    real(kind=kreal), intent(in)   :: disp(:) !intend (numnodeDOF)
    real(kind=kreal), intent(in)   :: vel(:) !intend (numnodeDOF)
    real(kind=kreal), intent(in)   :: acc(:) !intend (numnodeDOF)
    !---- vals
    integer(kind=kint) :: ndof, ntotal, idx
    !---- body
 
    ndof   = hecmesh%n_dof
    ntotal = ndof*hecmesh%n_node
 
    idx = 1
    if( fstrdynamic%idx_eqa == 1 .and. istp > 0 ) idx = 2
    fstrdynamic%DISP(1:ntotal,idx) = disp(1:ntotal)
    fstrdynamic%VEL (1:ntotal,idx) = vel(1:ntotal)
    fstrdynamic%ACC (1:ntotal,idx) = acc(1:ntotal)
 
    hecmat%X(1:ntotal) = disp(1:ntotal) - fstrsolid%unode(1:ntotal) - fstrsolid%dunode(1:ntotal)
    call fstr_apply_solution_increment( hecmesh, fstrsolid, ndof, hecmat%X )
    call fstr_updatenewton( hecmesh, hecmat, fstrsolid, fstrparam, 0.0d0, 0.0d0, 1 )
    call fstr_update_reaction_spc( 1, hecmesh, fstrsolid )
 
    if( ndof == 2 ) then
      call fstr_nodalstress2d( hecmesh, fstrsolid )
    else if( ndof == 6 ) then
      call fstr_nodalstress6d( hecmesh, fstrsolid )
    else
      call fstr_nodalstress3d( hecmesh, fstrsolid )
    end if
  end subroutine
 
  !
  ! The response is complex, so the real part and the imaginary part are written as
  ! two sets of the same items, distinguished by the label suffix. Both are recovered
  ! once and shared by the result file and the visualization data.
  !
  subroutine output_result(hecMESH, hecMAT, fstrSOLID, fstrDYNAMIC, fstrPARAM, nameid, glabel, gvalue, &
      istp, wres, wvis, dispRe, dispIm, velRe, velIm, accRe, accIm)
    !---- args
    type(hecmwst_local_mesh)       :: hecMESH
    type(hecmwst_matrix)           :: hecMAT
    type(fstr_solid)               :: fstrSOLID
    type(fstr_dynamic)             :: fstrDYNAMIC
    type(fstr_param)               :: fstrPARAM
    character(len=*), intent(in)   :: nameid
    character(len=*), intent(in)   :: glabel
    real(kind=kreal), intent(in)   :: gvalue
    integer(kind=kint), intent(in) :: istp
    logical, intent(in)            :: wres, wvis
    real(kind=kreal), intent(in)   :: dispre(:), dispim(:) !intend (numnodeDOF)
    real(kind=kreal), intent(in)   :: velre(:), velim(:) !intend (numnodeDOF)
    real(kind=kreal), intent(in)   :: accre(:), accim(:) !intend (numnodeDOF)
    !---- vals
    type(hecmwst_result_data)       :: visre, visim
    character(len=HECMW_HEADER_LEN) :: header
    character(len=HECMW_MSG_LEN)    :: comment
    character(len=HECMW_NAME_LEN)   :: label, resid
    real(kind=kreal)                :: gval(1)
    !---- body
 
    if( .not. (wres .or. wvis) ) return
 
    if( wres ) then
      header='*fstrresult'
      comment='frequency_result'
      call hecmw_result_init(hecmesh, istp, header, comment)
 
      label   = glabel
      gval(1) = gvalue
      call hecmw_result_add(hecmw_result_dtype_global, 1, label, gval)
    end if
 
    call updatesolidresult(hecmesh, hecmat, fstrsolid, fstrdynamic, fstrparam, istp, dispre, velre, accre)
    if( wres ) call fstr_write_result_add(hecmesh, fstrsolid, istp, '_real', fstrdynamic)
    if( wvis ) call fstr_make_result(hecmesh, fstrsolid, visre, istp, gvalue, fstrdynamic, '_real')
 
    call updatesolidresult(hecmesh, hecmat, fstrsolid, fstrdynamic, fstrparam, istp, dispim, velim, accim)
    if( wres ) call fstr_write_result_add(hecmesh, fstrsolid, istp, '_imag', fstrdynamic)
    if( wvis ) call fstr_make_result(hecmesh, fstrsolid, visim, istp, gvalue, fstrdynamic, '_imag')
 
    if( wres ) then
      resid = nameid
      call hecmw_result_write_by_name(resid)
      call hecmw_result_finalize()
    end if
 
    if( wvis ) then
      call fstr_merge_result(hecmesh, visre, visim)
      call hecmw_result_free(visim)
      call fstr2hecmw_mesh_conv(hecmesh)
      call hecmw_visualize_init
      call hecmw_visualize( hecmesh, visre, istp )
      call hecmw2fstr_mesh_conv(hecmesh)
      call hecmw_result_free(visre)
    end if
    return
  end subroutine
 
  subroutine extract_surf2node(hecMESH, freqData, loadvecRe, loadvecIm)
    !---- args
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    type(fstr_freqanalysis), intent(in)  :: freqData
    real(kind=kreal), intent(inout)      :: loadvecre(:) !intend(numnode*ndof)
    real(kind=kreal), intent(inout)      :: loadvecim(:) !intend(numnode*ndof)
    !---- vals
    integer(kind=kint), parameter   :: MAXNODE = 100
    integer(kind=kint)              :: sgrpID, is, ie, ic, nsurf, ic_type
    integer(kind=kint)              :: nn, dof_index, ndof
    integer(kind=kint)              :: i, j, k, l, isn, nsize
    integer(kind=kint)              :: iwk(60), nodLOCAL(20)
    real(kind=kreal)                :: vect(60), xx(20), yy(20), zz(20), forcere(3), forceim(3)
    !---- body
 
    ndof = 3
    do i=1,freqdata%FLOAD_ngrp_tot
      if(freqdata%FLOAD_ngrp_TYPE(i) == kfloadtype_surf) then  !FLOAD type=surface
        sgrpid     = freqdata%FLOAD_ngrp_ID(i)
        dof_index  = freqdata%FLOAD_ngrp_DOF(i)
        forcere(:) = 0.0d0
        forceim(:) = 0.0d0
        forcere(dof_index) = freqdata%FLOAD_ngrp_valre(i)
        forceim(dof_index) = freqdata%FLOAD_ngrp_valim(i)
 
        is = hecmesh%surf_group%grp_index(sgrpid-1) + 1
        ie = hecmesh%surf_group%grp_index(sgrpid)
        do j=is, ie
          ic      = hecmesh%surf_group%grp_item(2*j-1)
          nsurf   = hecmesh%surf_group%grp_item(2*j)
          ic_type = hecmesh%elem_type(ic)
          nn  = hecmw_get_max_node(ic_type)
          isn = hecmesh%elem_node_index(ic-1)
          do k=1, nn
            nodlocal(k) = hecmesh%elem_node_item(isn+k)
            xx(k) = hecmesh%node(3*nodlocal(k)-2)
            yy(k) = hecmesh%node(3*nodlocal(k)-1)
            zz(k) = hecmesh%node(3*nodlocal(k)  )
            do l=1, ndof
              iwk(ndof*(k-1)+l) = ndof*(nodlocal(k)-1)+l
            end do
          end do
 
          call dl_c3_freq(ic_type, nn, xx, yy, zz, nsurf, forcere, vect, nsize)
          do k=1,nsize
            loadvecre(iwk(k)) = loadvecre(iwk(k)) + vect(k)
          end do
 
          call dl_c3_freq(ic_type, nn, xx, yy, zz, nsurf, forceim, vect, nsize)
          do k=1,nsize
            loadvecim(iwk(k)) = loadvecim(iwk(k)) + vect(k)
          end do
        end do
      end if
    end do
 
    return
  end subroutine
 
  subroutine dl_c3_freq(ETYPE, NN, XX, YY, ZZ, LTYPE, force, VECT, nsize)
    !---- args
    integer(kind=kint), intent(in)    :: ETYPE !--solid element type
    integer(kind=kint), intent(in)    :: NN !--node num
    integer(kind=kint), intent(in)    :: LTYPE !--solid element face
    real(kind=kreal), intent(in)      :: xx(:) !--node x pos
    real(kind=kreal), intent(in)      :: yy(:) !--node y pos
    real(kind=kreal), intent(in)      :: zz(:) !--node z pos
    real(kind=kreal), intent(in)      :: force(3) !--node surfforce
    real(kind=kreal), intent(inout)   :: vect(:)
    integer(kind=kint), intent(inout) :: nsize
    !---- vals
    integer(kind=kint), parameter :: NDOF = 3
    real(kind=kreal)              :: wg, area
    integer(kind=kint)            :: NOD(NN)
    real(kind=kreal)              :: elecoord(3, nn), localcoord(3), normal(3)
    real(kind=kreal)              :: h(nn)
    integer(kind=kint)            :: I, IG2, NSUR, SURTYPE
    !---- body
 
    call getsubface( etype, ltype, surtype, nod )
    nsur = getnumberofnodes( surtype )
 
    do i=1,nsur
      elecoord(1,i)=xx(nod(i))
      elecoord(2,i)=yy(nod(i))
      elecoord(3,i)=zz(nod(i))
    end do
    nsize         = nn*ndof
    vect(1:nsize) = 0.0d0
    do ig2=1,numofquadpoints( surtype )
      call getquadpoint( surtype, ig2, localcoord(1:2) )
      call getshapefunc( surtype, localcoord(1:2), h(1:nsur) )
 
      wg=getweight( surtype, ig2 )
      normal=surfacenormal( surtype, nsur, localcoord(1:2), elecoord(:,1:nsur) )
      area=dsqrt( dot_product(normal, normal) )
      do i=1,nsur
        vect(3*nod(i)-2)=vect(3*nod(i)-2)+area*wg*h(i)*force(1)
        vect(3*nod(i)-1)=vect(3*nod(i)-1)+area*wg*h(i)*force(2)
        vect(3*nod(i)  )=vect(3*nod(i)  )+area*wg*h(i)*force(3)
      end do
    end do
  end subroutine
 
  subroutine assemble_nodeload(hecMESH, freqData, numdof, loadvecRe, loadvecIm)
    !---- args
    type(hecmwst_local_mesh), intent(in) :: hecMESH
    type(fstr_freqanalysis), intent(in)  :: freqData
    integer(kind=kint), intent(in)       :: numdof
    real(kind=kreal), intent(inout)      :: loadvecre(:)
    real(kind=kreal), intent(inout)      :: loadvecim(:)
    !---- vals
    integer(kind=kint) :: i, ig, is, ie, in, nodeid, dof_index
 
    !---- body
 
    do i=1, freqdata%FLOAD_ngrp_tot
      if(freqdata%FLOAD_ngrp_TYPE(i) == kfloadtype_node) then
        ig = freqdata%FLOAD_ngrp_ID(i)
        is = hecmesh%node_group%grp_index(ig-1) + 1
        ie = hecmesh%node_group%grp_index(ig)
        do in=is, ie
          nodeid    = hecmesh%node_group%grp_item(in)
          dof_index = freqdata%FLOAD_ngrp_DOF(i)
          loadvecre((nodeid-1)*numdof + dof_index) = loadvecre((nodeid-1)*numdof + dof_index) + freqdata%FLOAD_ngrp_valre(i)
          loadvecim((nodeid-1)*numdof + dof_index) = loadvecim((nodeid-1)*numdof + dof_index) + freqdata%FLOAD_ngrp_valim(i)
        end do
      end if
    end do
    return
  end subroutine
 
  subroutine calcmassmatrix(fstrPARAM, hecMESH, hecMAT, fstrSOLID, fstrEIG, hecLagMAT)
    !---- args
    type(fstr_param), intent(in)         :: fstrPARAM
    type(hecmwst_local_mesh), intent(in) :: hecMESH
    type(hecmwst_matrix), intent(inout)  :: hecMAT
    type(fstr_solid), intent(inout)      :: fstrSOLID
    type(fstr_eigen), intent(inout)      :: fstrEIG
    type(hecmwst_matrix_lagrange), intent(inout) :: hecLagMAT
    !---- vals
    type(hecmwst_ebc) :: hecEBC
    !---- body
 
 
    fstrsolid%dunode = 0.d0
    call fstr_creatematrix_and_dampingforce( hecmesh, hecmat, fstrsolid, 0.d0, 0.d0 )
    call hecmw_ebc_init(hecmat, hecebc)
    call fstr_addbc(1, hecmesh, hecmat, fstrsolid, fstrparam, heclagmat, 2, hecebc)
    call hecmw_ebc_apply(hecmesh, hecmat, hecebc)
    call hecmw_ebc_finalize(hecebc)
 
    call setmass(fstrsolid, hecmesh, hecmat, fstreig)
 
  end subroutine
 
  !
  ! The response to the acceleration input is the quasi-static displacement us, which
  ! satisfies K us = 0 with the input amplitude on the excited DOFs and zero on the other
  ! constrained DOFs, plus the relative response expanded by the modes fixed at the
  ! excited DOFs. us is computed for the part of the amplitude selected by loadcase.
  !
  subroutine calcquasistaticdisp(hecMESH, hecMAT, fstrSOLID, fstrPARAM, fstrFREQ, hecLagMAT, loadcase, us)
    !---- args
    type(hecmwst_local_mesh)             :: hecMESH
    type(hecmwst_matrix)                 :: hecMAT
    type(fstr_solid)                     :: fstrSOLID
    type(fstr_param)                     :: fstrPARAM
    type(fstr_freqanalysis), intent(in)  :: fstrFREQ
    type(hecmwst_matrix_lagrange)        :: hecLagMAT
    integer(kind=kint), intent(in)       :: loadcase
    real(kind=kreal), intent(inout)      :: us(:) !intend (numnodeDOF)
    !---- vals
    type(hecmwst_local_mesh), pointer :: hecmeshmpc
    type(hecmwst_matrix), pointer     :: hecMATmpc
    type(hecmwst_ebc)                 :: hecEBC
    integer(kind=kint)                :: ndof, ntotal, i, ig, ik, in, idof, idofS, idofE
    real(kind=kreal)                  :: val(fstrfreq%FACC_ngrp_tot) !amplitude of each entry
    !---- body
 
    ndof   = hecmesh%n_dof
    ntotal = ndof*hecmesh%n_node
    us(:)  = 0.0d0
 
    !amplitude of each entry for the part selected by loadcase
    if( loadcase == kfloadcase_re ) then
      val(:) = fstrfreq%FACC_ngrp_valre(1:fstrfreq%FACC_ngrp_tot)
    else
      val(:) = fstrfreq%FACC_ngrp_valim(1:fstrfreq%FACC_ngrp_tot)
    end if
    !us = 0 without solving if this part of the amplitude is zero
    if( maxval(abs(val)) <= 0.0d0 ) return
 
    !the following steps are the same as those of the linear static analysis (fstr_Newton)
    call hecmw_mpc_mat_init(hecmesh, hecmat, hecmeshmpc, hecmatmpc)
 
    !stiffness matrix at the initial state
    fstrsolid%dunode = 0.d0
    call fstr_creatematrix_and_dampingforce( hecmesh, hecmat, fstrsolid, 0.d0, 0.d0 )
    !!SPRING: the same stiffness as the eigen analysis (fstr_solve_eigen)
    call fstr_addspring(1, hecmesh, hecmat, fstrsolid)
    !no external force
    call hecmw_mat_clear_b(hecmat)
    !!BOUNDARY: iter=2 gives the constrained DOFs with zero displacement
    call hecmw_ebc_init(hecmat, hecebc)
    call fstr_addbc(1, hecmesh, hecmat, fstrsolid, fstrparam, heclagmat, 2, hecebc)
    !!FACCELERATION: the amplitude is given as the displacement of the excited DOFs.
    !it overwrites the value of !BOUNDARY on the same DOF.
    do i=1, fstrfreq%FACC_ngrp_tot
      ig    = fstrfreq%FACC_ngrp_ID(i)
      !decode 10*dof_start + dof_end
      idofs = fstrfreq%FACC_ngrp_type(i)/10
      idofe = fstrfreq%FACC_ngrp_type(i) - idofs*10
      !all nodes of the group, including the external nodes
      do ik=hecmesh%node_group%grp_index(ig-1)+1, hecmesh%node_group%grp_index(ig)
        in = hecmesh%node_group%grp_item(ik)
        do idof=idofs, idofe
          call hecmw_ebc_set(hecebc, in, idof, val(i))
        end do
      end do
    end do
    !apply the MPC and the prescribed displacements to the matrix and the right hand side
    call hecmw_mpc_mat_ass(hecmesh, hecmat, hecmeshmpc, hecmatmpc)
    call hecmw_mpc_trans_rhs(hecmesh, hecmat, hecmatmpc)
    call hecmw_ebc_apply(hecmeshmpc, hecmatmpc, hecebc)
 
    !solve K us = 0
    hecmatmpc%Iarray(98) = 1 !need symbolic setup of the solver
    hecmatmpc%Iarray(97) = 1 !need numerical setup of the solver
    hecmatmpc%X = 0.0d0
    call solve_lineq(hecmeshmpc, hecmatmpc)
    call hecmw_mpc_tback_sol(hecmesh, hecmat, hecmatmpc)
    !values of the external nodes
    call hecmw_update_r(hecmesh, hecmat%X, hecmat%NP, ndof)
    us(1:ntotal) = hecmat%X(1:ntotal)
 
    call hecmw_ebc_finalize(hecebc)
    call hecmw_mpc_mat_finalize(hecmesh, hecmat, hecmeshmpc, hecmatmpc)
  end subroutine
 
  !
  ! The modes must vanish on the excited DOFs, i.e. the eigen analysis must constrain them.
  ! Otherwise the relative response would move the excited DOFs and their acceleration
  ! would differ from the input. Stop with an error in that case.
  !
  subroutine checkaccelerationdof(hecMESH, fstrFREQ, freqData)
    !---- args
    type(hecmwst_local_mesh), intent(in)     :: hecmesh
    type(fstr_freqanalysis), intent(in)      :: fstrFREQ
    type(fstr_freqanalysis_data), intent(in) :: freqData
    !---- vals
    real(kind=kreal), parameter :: tol = 1.0d-8 !relative to the largest component of the mode
    integer(kind=kint)          :: ndof, i, ig, ik, in, idof, idofS, idofE, imode, nerr
    real(kind=kreal)            :: vmax
    !---- body
 
    ndof = hecmesh%n_dof
    nerr = 0
    do imode=1, freqdata%numMode
      !largest component of the mode over all ranks, used as the scale of the tolerance
      vmax = maxval(abs(freqdata%eigVector(1:hecmesh%nn_internal*ndof, imode)))
      call hecmw_allreduce_r1(hecmesh, vmax, hecmw_max)
      do i=1, fstrfreq%FACC_ngrp_tot
        ig    = fstrfreq%FACC_ngrp_ID(i)
        !decode 10*dof_start + dof_end (DOFs beyond ndof are not checked)
        idofs = fstrfreq%FACC_ngrp_type(i)/10
        idofe = min(fstrfreq%FACC_ngrp_type(i) - idofs*10, ndof)
        do ik=hecmesh%node_group%grp_index(ig-1)+1, hecmesh%node_group%grp_index(ig)
          in = hecmesh%node_group%grp_item(ik)
          !each node is checked by the rank that owns it
          if( in > hecmesh%nn_internal ) cycle
          do idof=idofs, idofe
            !count the excited DOFs where the mode does not vanish
            if( abs(freqdata%eigVector(ndof*(in-1)+idof, imode)) > tol*vmax ) nerr = nerr + 1
          end do
        end do
      end do
    end do
    !stop on all ranks if any rank has found an error
    call hecmw_allreduce_i1(hecmesh, nerr, hecmw_sum)
    if( nerr > 0 ) then
      if( myrank == 0 ) then
        write(*,   *) "Error: the DOFs of !FACCELERATION must be constrained in the eigen analysis"
        write(ilog,*) "Error: the DOFs of !FACCELERATION must be constrained in the eigen analysis"
      end if
      call hecmw_abort( hecmw_comm_get_comm() )
    end if
  end subroutine
 
  subroutine scaleeigenvector(hecMESH, fstrEIG, ntotaldof, nmode, eigenvector)
    !---- args
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    type(fstr_eigen), intent(in)    :: fstrEIG
    integer(kind=kint), intent(in)  :: ntotaldof
    integer(kind=kint), intent(in)  :: nmode
    real(kind=kreal), intent(inout) :: eigenvector(:, :)
    !---- vals
    integer(kind=kint) :: imode, idof, nintdof
    real(kind=kreal)   :: mas
    !---- body
 
    nintdof = hecmesh%nn_internal*hecmesh%n_dof
 
    do imode=1,nmode
      mas = 0.0d0
      do idof=1,nintdof
        mas = mas + fstreig%mass(idof)*eigenvector(idof,imode)**2
      end do
      call hecmw_allreduce_r1(hecmesh, mas, hecmw_sum)
      do idof=1,ntotaldof
        eigenvector(idof,imode) =  eigenvector(idof,imode) / sqrt(mas)
      end do
    end do
  end subroutine
 
  subroutine checkorthvector(fstrEIG, eigenvector, imode, jmode, prod)
    !---- args
    type(fstr_eigen), intent(in)    :: fstreig
    real(kind=kreal), intent(in)    :: eigenvector(:, :)
    integer(kind=kint), intent(in)  :: imode
    integer(kind=kint), intent(in)  :: jmode
    real(kind=kreal), intent(inout) :: prod
    !---- vals
    integer(kind=kint) :: idof, s
    !---- body
    s = size(eigenvector(:,1))
    prod = 0.0d0
 
    do idof=1, s
      prod = prod + eigenvector(idof,imode)*fstreig%mass(idof)*eigenvector(idof,jmode)
    end do
    return
  end subroutine
 
  subroutine writeoutvector(im, vector)
    !---- args
    integer(kind=kint), intent(in) :: im
    real(kind=kreal), intent(in)   :: vector(:)
    !---- vals
    integer(kind=kint) :: i, s
    !---- body
    s = size(vector)
    do i=1, s
      if(i == 1) then
        write(*,'("eigenvec",i2.2,":[",e12.5,", ")') im, vector(i)
      else if(i /= s) then
        write(*,'(e12.5,", ")') vector(i)
      else
        write(*,'(e12.5,"];")') vector(i)
      end if
    end do
    write(*,*)
    return
  end subroutine
 
  subroutine calcdotproduct(hecMESH, a, b, c)
    !---- args
    type(hecmwst_local_mesh), intent(in) :: hecMESH
    real(kind=kreal), intent(in)    :: a(:)
    real(kind=kreal), intent(in)    :: b(:)
    real(kind=kreal), intent(inout) :: c
    !---- vals
    integer(kind=kint) :: nintdof
    !---- body
    nintdof = hecmesh%nn_internal*hecmesh%n_dof  !external nodes are owned by another rank
    c = dot_product(a(1:nintdof), b(1:nintdof))
    call hecmw_allreduce_r1(hecmesh, c, hecmw_sum)
    return
  end subroutine
 
  subroutine calcfreqcoeff(hecMESH, freqData, loadRe, loadIm, inpOmega, bjRe, bjIm)
    !---- args
    type(hecmwst_local_mesh), intent(in)     :: hecMESH
    type(fstr_freqanalysis_data), intent(in) :: freqData
    real(kind=kreal), intent(in)             :: loadre(:) !intend (numNodeDOF)
    real(kind=kreal), intent(in)             :: loadim(:) !intend (numNodeDOF)
    real(kind=kreal), intent(in)             :: inpomega
    real(kind=kreal), intent(inout)          :: bjre(:) !intend (numMode)
    real(kind=kreal), intent(inout)          :: bjim(:) !intend (numMode)
    !---- vals
    integer(kind=kint) :: imode
    real(kind=kreal)   :: ujfr, ujfi, a, b, alp, beta
    !---- body
 
    alp  = freqdata%rayAlpha
    beta = freqdata%rayBeta
 
    do imode=1, freqdata%numMode
      call calcdotproduct(hecmesh, freqdata%eigVector(:,imode), loadre, ujfr)
      call calcdotproduct(hecmesh, freqdata%eigVector(:,imode), loadim, ujfi)
 
      a = ujfr*(freqdata%eigOmega(imode)**2 - inpomega**2) + ujfi*(alp + beta*freqdata%eigOmega(imode)**2)*inpomega
      b = (freqdata%eigOmega(imode)**2 - inpomega**2)**2 + ((alp + beta*freqdata%eigOmega(imode)**2)*inpomega)**2
      bjre(imode) = a / b
 
      a = ujfi*(freqdata%eigOmega(imode)**2 -inpomega**2) - ujfr*(alp + beta*freqdata%eigOmega(imode)**2)*inpomega
      b = (freqdata%eigOmega(imode)**2 - inpomega**2)**2 + ((alp + beta*freqdata%eigOmega(imode)**2)*inpomega)**2
      bjim(imode) = a / b
    end do
    return
  end subroutine
 
  !
  ! The inertia and the mass proportional damping of the quasi-static motion load the
  ! relative response as -M us (1 - i alpha/omega). The stiffness proportional part
  ! vanishes on the free DOFs since K us = 0 there.
  !
  subroutine addaccelerationload(fstrEIG, freqData, omega, usRe, usIm, loadRe, loadIm)
    !---- args
    type(fstr_eigen), intent(in)             :: fstreig
    type(fstr_freqanalysis_data), intent(in) :: freqData
    real(kind=kreal), intent(in)             :: omega
    real(kind=kreal), intent(in)             :: usre(:) !intend (numNodeDOF)
    real(kind=kreal), intent(in)             :: usim(:) !intend (numNodeDOF)
    real(kind=kreal), intent(inout)          :: loadre(:) !intend (numNodeDOF)
    real(kind=kreal), intent(inout)          :: loadim(:) !intend (numNodeDOF)
    !---- vals
    integer(kind=kint) :: n
    real(kind=kreal)   :: c
    !---- body
 
    n = freqdata%numNodeDOF
    c = freqdata%rayAlpha / omega
    !fstrEIG%mass is the lumped mass (diagonal) computed by calcMassMatrix
    !real and imaginary parts of -M (usRe + i usIm) (1 - i c)
    loadre(1:n) = loadre(1:n) - fstreig%mass(1:n)*(usre(1:n) + c*usim(1:n))
    loadim(1:n) = loadim(1:n) - fstreig%mass(1:n)*(usim(1:n) - c*usre(1:n))
  end subroutine
 
  !
  ! Add coef*us (complex) to the vector vec (complex), i.e. the quasi-static part of
  ! the displacement, the velocity or the acceleration.
  !
  subroutine addquasistaticvector(coef, usRe, usIm, vecRe, vecIm)
    !---- args
    complex(kind=kreal), intent(in) :: coef
    real(kind=kreal), intent(in)    :: usre(:) !intend (numNodeDOF)
    real(kind=kreal), intent(in)    :: usim(:) !intend (numNodeDOF)
    real(kind=kreal), intent(inout) :: vecre(:) !intend (numNodeDOF)
    real(kind=kreal), intent(inout) :: vecim(:) !intend (numNodeDOF)
    !---- body
 
    !(coefRe + i coefIm) (usRe + i usIm)
    vecre(:) = vecre(:) + dble(coef)*usre(:) - aimag(coef)*usim(:)
    vecim(:) = vecim(:) + dble(coef)*usim(:) + aimag(coef)*usre(:)
  end subroutine
 
  subroutine calcdispvector(freqData, bjRe, bjIm, dispRe, dispIm)
    !---- args
    type(fstr_freqanalysis_data), intent(in) :: freqdata
    real(kind=kreal), intent(in)             :: bjre(:) !intend (numMode)
    real(kind=kreal), intent(in)             :: bjim(:) !intend (numMode)
    real(kind=kreal), intent(inout)          :: dispre(:) !intend (numNodeDOF)
    real(kind=kreal), intent(inout)          :: dispim(:) !intend (numNodeDOF)
    !---- vals
    integer(kind=kint) :: imode
    !---- body
 
    dispre(:) = 0.0d0
    dispim(:) = 0.0d0
 
    do imode=1, freqdata%numMode
      dispre(:) = dispre(:) + bjre(imode)*freqdata%eigVector(:,imode)
      dispim(:) = dispim(:) + bjim(imode)*freqdata%eigVector(:,imode)
    end do
    return
  end subroutine
 
  subroutine calcvelvector(freqData, omega, bjRe, bjIm, velRe, velIm)
    !---- args
    type(fstr_freqanalysis_data), intent(in) :: freqData
    real(kind=kreal), intent(in)             :: omega
    real(kind=kreal), intent(in)             :: bjre(:) !intend (numMode)
    real(kind=kreal), intent(in)             :: bjim(:) !intend (numMode)
    real(kind=kreal), intent(inout)          :: velre(:) !intend (numNodeDOF)
    real(kind=kreal), intent(inout)          :: velim(:) !intend (numNodeDOF)
    !---- vals
    integer(kind=kint) :: imode
    !---- body
 
    velre(:) = 0.0d0
    velim(:) = 0.0d0
 
    do imode=1, freqdata%numMode
      velre(:) = velre(:) - omega * bjim(imode) * freqdata%eigVector(:,imode)
      velim(:) = velim(:) + omega * bjre(imode) * freqdata%eigVector(:,imode)
    end do
  end subroutine
 
  subroutine calcaccvector(freqData, omega, bjRe, bjIm, accRe, accIm)
    !---- args
    type(fstr_freqanalysis_data), intent(in) :: freqData
    real(kind=kreal), intent(in)             :: omega
    real(kind=kreal), intent(in)             :: bjre(:) !intend (numMode)
    real(kind=kreal), intent(in)             :: bjim(:) !intend (numMode)
    real(kind=kreal), intent(inout)          :: accre(:) !intend (numNodeDOF)
    real(kind=kreal), intent(inout)          :: accim(:) !intend (numNodeDOF)
    !---- vals
    integer(kind=kint) :: imode
    !---- body
 
    accre(:) = 0.0d0
    accim(:) = 0.0d0
 
    do imode=1, freqdata%numMode
      accre(:) = accre(:) - omega**2 * bjre(imode) * freqdata%eigVector(:,imode)
      accim(:) = accim(:) - omega**2 * bjim(imode) * freqdata%eigVector(:,imode)
    end do
 
  end subroutine
 
  subroutine setupfreqparam(fstrDYNAMIC, f_start, f_end, numfreq, raym, rayk, idnode, vistype, ioutl)
    !---- args
    type(fstr_dynamic), intent(in)    :: fstrDYNAMIC
    real(kind=kreal), intent(inout)   :: f_start
    real(kind=kreal), intent(inout)   :: f_end
    integer(kind=kint), intent(inout) :: numfreq
    real(kind=kreal), intent(inout)   :: raym
    real(kind=kreal), intent(inout)   :: rayk
    integer(kind=kint), intent(inout) :: idnode
    integer(kind=kint), intent(inout) :: vistype
    integer(kind=kint), intent(inout) :: ioutl(3)
    !---- vals
 
    !---- body
    f_start    = fstrdynamic%t_start
    f_end      = fstrdynamic%t_end
    numfreq    = fstrdynamic%n_step
    raym       = fstrdynamic%ray_m
    rayk       = fstrdynamic%ray_k
    idnode     = fstrdynamic%nout_monit
    vistype    = fstrdynamic%ngrp_monit
    ioutl(1:3) = fstrdynamic%iout_list(1:3)
    return
  end subroutine
 
  subroutine calcdispvectortime(freqData, time, omega, bjRe, bjIm, dispRe, dispIm)
    !---- args
    type(fstr_freqanalysis_data), intent(in) :: freqData
    real(kind=kreal), intent(in)             :: time
    real(kind=kreal), intent(in)             :: omega
    real(kind=kreal), intent(in)             :: bjre(:) !intend (numMode)
    real(kind=kreal), intent(in)             :: bjim(:) !intend (numMode)
    real(kind=kreal), intent(inout)          :: dispre(:) !intend (numNodeDOF)
    real(kind=kreal), intent(inout)          :: dispim(:) !intend (numNodeDOF)
    !---- vals
    integer(kind=kint)  :: imode, idf, s
    complex(kind=kreal) :: a, b, c
    !---- body
 
    dispre(:) = 0.0d0
    dispim(:) = 0.0d0
    a = exp(dcmplx(0.0d0, omega*time))
 
    do imode=1, freqdata%numMode
      s = size(freqdata%eigvector(:,imode))
      b = dcmplx(bjre(imode), bjim(imode)) * a
      do idf=1, s
        c = b*dcmplx(freqdata%eigVector(idf,imode), 0.0d0)
        dispre(idf) = dispre(idf) + dble(c)
        dispim(idf) = dispim(idf) + imag(c)
      end do
    end do
    return
  end subroutine
 
  subroutine calcvelvectortime(freqData, time, omega, bjRe, bjIm, velRe, velIm)
    !---- args
    type(fstr_freqanalysis_data), intent(in) :: freqData
    real(kind=kreal), intent(in)             :: time
    real(kind=kreal), intent(in)             :: omega
    real(kind=kreal), intent(in)             :: bjre(:) !intend (numMode)
    real(kind=kreal), intent(in)             :: bjim(:) !intend (numMode)
    real(kind=kreal), intent(inout)          :: velre(:) !intend (numNodeDOF)
    real(kind=kreal), intent(inout)          :: velim(:) !intend (numNodeDOF)
    !---- vals
    integer(kind=kint)  :: imode, idf, s
    complex(kind=kreal) :: a, b, c
    !---- body
 
    velre(:) = 0.0d0
    velim(:) = 0.0d0
    a = dcmplx(0.0d0, 1.0d0)*dcmplx(omega, 0.0d0)*exp(dcmplx(0.0d0, omega*time))
 
    do imode=1, freqdata%numMode
      s = size(freqdata%eigvector(:,imode))
      b = dcmplx(bjre(imode), bjim(imode)) * a
      do idf=1, s
        c = b*dcmplx(freqdata%eigVector(idf,imode), 0.0d0)
        velre(idf) = velre(idf) + dble(c)
        velim(idf) = velim(idf) + imag(c)
      end do
    end do
    return
  end subroutine
 
  subroutine calcaccvectortime(freqData, time, omega, bjRe, bjIm, accRe, accIm)
    !---- args
    type(fstr_freqanalysis_data), intent(in) :: freqData
    real(kind=kreal), intent(in)             :: time
    real(kind=kreal), intent(in)             :: omega
    real(kind=kreal), intent(in)             :: bjre(:) !intend (numMode)
    real(kind=kreal), intent(in)             :: bjim(:) !intend (numMode)
    real(kind=kreal), intent(inout)          :: accre(:) !intend (numNodeDOF)
    real(kind=kreal), intent(inout)          :: accim(:) !intend (numNodeDOF)
    !---- vals
    integer(kind=kint)  :: imode, idf, s
    complex(kind=kreal) :: a, b, c
    !---- body
 
    accre(:) = 0.0d0
    accim(:) = 0.0d0
    a = dcmplx(-1.0d0, 0.0d0)*dcmplx(omega**2, 0.0d0)*exp(dcmplx(0.0d0, omega*time))
 
    do imode=1, freqdata%numMode
      s = size(freqdata%eigvector(:,imode))
      b = dcmplx(bjre(imode), bjim(imode)) * a
      do idf=1, s
        c = b*dcmplx(freqdata%eigVector(idf,imode), 0.0d0)
        accre(idf) = accre(idf) + dble(c)
        accim(idf) = accim(idf) + imag(c)
      end do
    end do
    return
  end subroutine
 
  subroutine setupdynaparam(fstrDYNAMIC, t_start, t_end, dynafreq, numdisp)
    !---- args
    type(fstr_dynamic), intent(in)    :: fstrDYNAMIC
    real(kind=kreal), intent(inout)   :: t_start
    real(kind=kreal), intent(inout)   :: t_end
    real(kind=kreal), intent(inout)   :: dynafreq
    integer(kind=kint), intent(inout) :: numdisp
    !---- vals
    !---- body
    t_start  = fstrdynamic%gamma
    t_end    = fstrdynamic%beta
    dynafreq = fstrdynamic%t_delta
    numdisp  = fstrdynamic%nout
    return
  end subroutine
 
end module