m_fstr.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
! If new header is supported, change according to following method.    !
!  1) Increase FSTR_CTRL_HEADER_NUMBER                                 !
!  2) Add new header name to fstr_ctrl_header_names                    !
!  3) Describe new function to get parameters from control file        !
!     in fstr_ctrl.f90                                                 !
!  4) Describe new subroutine to set values of the parameter           !
!     in fstr_setup.f90                                                !
!  5) If initial values are necessary, set the value                   !
!     in subroutine fstr_setup_init in fstr_setup.f90                  !
module m_fstr
  use hecmw
  use m_common_struct
  use m_step
  use m_out
  use m_timepoint
  use m_elemact
  use mmechgauss
  use mcontactdef
  use m_fstr_contact_smoothing
 
  implicit none
 
  public
 
  integer(kind=kint), parameter :: kyes = 1
  integer(kind=kint), parameter :: kno  = 0
  integer(kind=kint), parameter :: kon  = 1
  integer(kind=kint), parameter :: koff = 0
 
  integer(kind=kint), parameter :: kstprecheck = 0
  integer(kind=kint), parameter :: kststatic   = 1
  integer(kind=kint), parameter :: ksteigen    = 2
  integer(kind=kint), parameter :: kstheat     = 3
  integer(kind=kint), parameter :: kstdynamic  = 4
  !integer(kind=kint), parameter :: kstNLSTATIC  =   5
  integer(kind=kint), parameter :: kststaticeigen = 6
  integer(kind=kint), parameter :: kstnzprof      = 7
 
  integer(kind=kint), parameter :: ksmcg       = 1
  integer(kind=kint), parameter :: ksmbicgstab = 2
  integer(kind=kint), parameter :: ksmgmres    = 3
  integer(kind=kint), parameter :: ksmgpbicg   = 4
  integer(kind=kint), parameter :: ksmgmresr   = 5
  integer(kind=kint), parameter :: ksmgmresren = 6
  integer(kind=kint), parameter :: ksmdirect   = 101
 
  integer(kind=kint), parameter :: knsmnewton      = 1
  integer(kind=kint), parameter :: knsmquasinewton = 2
 
  integer(kind=kint), parameter :: kcaslagrange = 1
  integer(kind=kint), parameter :: kcaalagrange = 2
 
  integer(kind=kint), parameter :: kbcffstr     =   0  ! BC described in fstr control file (default)
  integer(kind=kint), parameter :: kbcfnastran  =   1  ! nastran file
 
  integer(kind=kint), parameter :: kbcinitial   =  1
  integer(kind=kint), parameter :: kbctransit   =  2
 
  integer(kind=kint), parameter :: restart_outlast = 1
  integer(kind=kint), parameter :: restart_outall  = 2
 
  integer(kind=kint), parameter :: kel341fi     =  1
  integer(kind=kint), parameter :: kel341sesns  =  2
 
  integer(kind=kint), parameter :: kel361fi     =  1
  integer(kind=kint), parameter :: kel361bbar   =  2
  integer(kind=kint), parameter :: kel361ic     =  3
  integer(kind=kint), parameter :: kel361fbar   =  4
 
  integer(kind=kint), parameter :: kfloadtype_node = 1
  integer(kind=kint), parameter :: kfloadtype_surf = 2
 
  integer(kind=kint), parameter :: kfloadcase_re = 1
  integer(kind=kint), parameter :: kfloadcase_im = 2
 
  integer(kind=kint), parameter :: kitrcontinue = 0
  integer(kind=kint), parameter :: kitrconverged = 1
  integer(kind=kint), parameter :: kitrdiverged = 2
  integer(kind=kint), parameter :: kitrfloatingerror = 3
 
 
  integer(kind = kint) :: myrank
  integer(kind = kint) :: nprocs
 
  logical :: paracontactflag = .false.
 
  character(len=HECMW_FILENAME_LEN) :: cntfilname
  character(len=HECMW_FILENAME_LEN) :: restartfilname
 
  integer(kind=kint), parameter :: ilog = 16 ! log
  integer(kind=kint), parameter :: ista = 17 ! status
  integer(kind=kint), parameter :: iutb = 18 ! utable
  integer(kind=kint), parameter :: imsg = 51 ! message (myrank == 0 only)
  integer(kind=kint), parameter :: idbg = 52 ! debug
  integer(kind=kint), parameter :: ifvs = 53 ! visual.ini file
  integer(kind=kint), parameter :: ineu = 54 ! neutral file (heat)
  integer(kind=kint), parameter :: iresout = 100 ! ~110, keeping for result output file
 
  integer(kind=kint) :: sviarray(100)
  real(kind=kreal)   :: svrarray(100)
 
  integer(kind=kint), pointer :: iecho
  integer(kind=kint), pointer :: iresult
  integer(kind=kint), pointer :: ivisual
  integer(kind=kint), pointer :: ineutral  ! flag for femap neutral file
  integer(kind=kint), pointer :: irres     ! flag for restart, read
  integer(kind=kint), pointer :: iwres     ! flag for restart, write
  integer(kind=kint), pointer :: nrres     ! position of restart read
  integer(kind=kint), pointer :: nprint    ! interval of write
 
  integer(kind=kint), parameter :: kopss_solution = 1
  integer(kind=kint), parameter :: kopss_material = 2
  integer(kind=kint)            :: opsstype = kopss_solution ! output stress/strain type
 
 
  real(kind=kreal), pointer :: ref_temp
 
  real(kind=kreal)   :: dt    ! /=fstr_param%dtime
  real(kind=kreal)   :: etime ! /=fstr_param%etime
  integer(kind=kint) :: itmax
  real(kind=kreal)   :: eps   ! /=fstr_param%eps
 
  type tinitialcondition
     character(len=HECMW_FILENAME_LEN)          :: cond_name
     integer                    :: node_elem
     integer                    :: grpid
     integer, pointer           :: intval(:)     => null()      
     real(kind=kreal), pointer  :: realval(:)    => null()      
  end type
  type( tinitialcondition ), pointer, save :: g_initialcnd(:) => null()
 
  type fstr_param
    integer(kind=kint) :: solution_type
    integer(kind=kint) :: solver_method
    integer(kind=kint) :: nlsolver_method
    logical            :: nlgeom
 
    integer(kind=kint) :: analysis_n
    real(kind=kreal), pointer  :: dtime(:) 
    real(kind=kreal), pointer  :: etime(:) 
    real(kind=kreal), pointer  :: dtmin(:) 
    real(kind=kreal), pointer  :: delmax(:)
    integer(kind=kint), pointer:: itmax(:)
    real(kind=kreal), pointer  :: eps(:)   
    real(kind=kreal)           :: ref_temp 
    integer(kind=kint)         :: timepoint_id
 
    integer(kind=kint) :: fg_echo
    integer(kind=kint) :: fg_result
    integer(kind=kint) :: fg_visual
 
    integer(kind=kint) :: fg_neutral
    integer(kind=kint) :: fg_irres
    integer(kind=kint) :: fg_iwres
    integer(kind=kint) :: nrres
    integer(kind=kint) :: nprint
 
    integer(kind=kint) :: n_node
    integer(kind=kint) :: nn_internal
    integer(kind=kint), pointer :: global_local_id(:,:)
 
    integer( kind=kint ) :: fg_couple
    integer( kind=kint ) :: fg_couple_type
    integer( kind=kint ) :: fg_couple_first
    integer( kind=kint ) :: fg_couple_window
 
    integer( kind=kint ) :: restart_out_type
    integer( kind=kint ) :: restart_version
 
    integer( kind=kint ) :: contact_algo
    integer( kind=kint ) :: augiter
    type(tcontactparam), pointer :: contactparam(:)
    type(tcontactinterference), pointer :: contact_if(:)
 
    type(tparamautoinc), pointer :: ainc(:)
    type(time_points), pointer :: timepoints(:)
  end type fstr_param
 
  type( fstr_param ),target :: fstrpr
 
  type fstr_solid_physic_val
    real(kind=kreal), pointer :: stress(:) => null()   
    real(kind=kreal), pointer :: strain(:) => null()   
    real(kind=kreal), pointer :: mises(:) => null()    
 
    real(kind=kreal), pointer :: pstress(:) => null()   
    real(kind=kreal), pointer :: pstrain(:) => null()   
    real(kind=kreal), pointer :: pstress_vect(:,:) => null()  
    real(kind=kreal), pointer :: pstrain_vect(:,:) => null()  
 
    real(kind=kreal), pointer :: estress(:) => null()  
    real(kind=kreal), pointer :: estrain(:) => null()  
    real(kind=kreal), pointer :: emises(:) => null()   
    real(kind=kreal), pointer :: eplstrain(:) => null()   
 
    real(kind=kreal), pointer :: epstress(:) => null()  
    real(kind=kreal), pointer :: epstrain(:) => null()  
    real(kind=kreal), pointer :: epstress_vect(:,:) => null()  
    real(kind=kreal), pointer :: epstrain_vect(:,:) => null()  
    real(kind=kreal), pointer :: enqm(:) => null()     
 
 
    type(fstr_solid_physic_val), pointer :: layer(:) => null()   
    type(fstr_solid_physic_val), pointer :: plus => null()   
    type(fstr_solid_physic_val), pointer :: minus => null()  
  end type fstr_solid_physic_val
 
  type fstr_solid
    integer(kind=kint) :: file_type  ! kbcfFSTR or kbcfNASTRAN
    integer(kind=kint) :: statictype ! 1:Total, 2:Updated, 3:Infinitesimal
    integer(kind=kint) :: nstep_tot
 
    type(step_info), pointer       :: step_ctrl(:)  =>null()   
    type(t_output_ctrl),  pointer  :: output_ctrl(:)=>null()   
 
    integer(kind=kint) :: boundary_ngrp_tot
    integer(kind=kint), pointer :: boundary_ngrp_grpid  (:)  =>null()
    integer(kind=kint), pointer :: boundary_ngrp_id     (:)  =>null()
    integer(kind=kint), pointer :: boundary_ngrp_type   (:)  =>null()
    integer(kind=kint), pointer :: boundary_ngrp_amp    (:)  =>null()
    real(kind=kreal), pointer   :: boundary_ngrp_val(:)  =>null()
    integer(kind=kint), pointer :: boundary_ngrp_istot  (:)  =>null()
    integer(kind=kint) :: boundary_ngrp_rot
    integer(kind=kint), pointer :: boundary_ngrp_rotid     (:) =>null()
    integer(kind=kint), pointer :: boundary_ngrp_centerid  (:) =>null()
 
    integer(kind=kint) :: velocity_type
    integer(kind=kint) :: velocity_ngrp_tot
    integer(kind=kint), pointer :: velocity_ngrp_grpid  (:)  =>null()
    integer(kind=kint), pointer :: velocity_ngrp_id     (:)  =>null()
    integer(kind=kint), pointer :: velocity_ngrp_type   (:)  =>null()
    integer(kind=kint), pointer :: velocity_ngrp_amp    (:)  =>null()
    real(kind=kreal), pointer   :: velocity_ngrp_val(:)  =>null()
 
    integer(kind=kint) :: acceleration_type
    integer(kind=kint) :: acceleration_ngrp_tot
    integer(kind=kint), pointer :: acceleration_ngrp_grpid  (:)  =>null()
    integer(kind=kint), pointer :: acceleration_ngrp_id     (:)  =>null()
    integer(kind=kint), pointer :: acceleration_ngrp_type   (:)  =>null()
    integer(kind=kint), pointer :: acceleration_ngrp_amp    (:)  =>null()
    real(kind=kreal), pointer   :: acceleration_ngrp_val(:)  =>null()
 
    integer(kind=kint) :: cload_ngrp_tot
    integer(kind=kint), pointer :: cload_ngrp_grpid     (:) =>null()
    integer(kind=kint), pointer :: cload_ngrp_id        (:)
    integer(kind=kint), pointer :: cload_ngrp_dof       (:)
    integer(kind=kint), pointer :: cload_ngrp_amp       (:)
    real(kind=kreal), pointer   :: cload_ngrp_val(:)
    integer(kind=kint) :: cload_ngrp_rot
    integer(kind=kint), pointer :: cload_ngrp_rotid     (:) =>null()
    integer(kind=kint), pointer :: cload_ngrp_centerid  (:) =>null()
 
    integer(kind=kint) :: dload_ngrp_tot
    integer(kind=kint) :: dload_follow
    integer(kind=kint), pointer :: dload_ngrp_grpid     (:) =>null()
    integer(kind=kint), pointer :: dload_ngrp_id        (:)
    integer(kind=kint), pointer :: dload_ngrp_lid       (:)
    integer(kind=kint), pointer :: dload_ngrp_amp       (:)
    real(kind=kreal), pointer   :: dload_ngrp_params(:,:)
 
    integer(kind=kint) :: temp_ngrp_tot
    integer(kind=kint) :: temp_irres
    integer(kind=kint) :: temp_tstep
    integer(kind=kint) :: temp_interval
    integer(kind=kint) :: temp_rtype      ! type of reading result; 1: step-based; 2: time-based
    real(kind=kreal)   :: temp_factor
    integer(kind=kint), pointer :: temp_ngrp_grpid     (:) =>null()
    integer(kind=kint), pointer :: temp_ngrp_id        (:)
    real(kind=kreal), pointer   :: temp_ngrp_val(:)
 
    integer(kind=kint) :: spring_ngrp_tot
    integer(kind=kint), pointer :: spring_ngrp_grpid    (:) =>null()
    integer(kind=kint), pointer :: spring_ngrp_id       (:)
    integer(kind=kint), pointer :: spring_ngrp_dof      (:)
    integer(kind=kint), pointer :: spring_ngrp_amp      (:)
    real(kind=kreal), pointer   :: spring_ngrp_val(:)
 
    type(telemact) :: elemact
 
    integer( kind=kint ) :: couple_ngrp_tot
    integer( kind=kint ),pointer :: couple_ngrp_id(:)
 
    integer(kind=kint) :: maxn_gauss
 
    real(kind=kreal), pointer :: stress(:)    
    real(kind=kreal), pointer :: strain(:)    
    real(kind=kreal), pointer :: mises(:)     
 
    real(kind=kreal), pointer :: pstress(:)   
    real(kind=kreal), pointer :: pstrain(:)   
    real(kind=kreal), pointer :: pstress_vect(:,:)   
    real(kind=kreal), pointer :: pstrain_vect(:,:)   
 
    real(kind=kreal), pointer :: estress(:)   
    real(kind=kreal), pointer :: estrain(:)   
    real(kind=kreal), pointer :: emises(:)    
    real(kind=kreal), pointer :: eplstrain(:)    
 
    real(kind=kreal), pointer :: epstress(:)   
    real(kind=kreal), pointer :: epstrain(:)   
    real(kind=kreal), pointer :: epstress_vect(:,:)   
    real(kind=kreal), pointer :: epstrain_vect(:,:)   
 
    real(kind=kreal), pointer :: tnstrain(:)   
    real(kind=kreal), pointer :: testrain(:)   
 
    real(kind=kreal), pointer :: yield_ratio(:)    
 
    real(kind=kreal), pointer :: enqm(:)      
    real(kind=kreal), pointer :: reaction(:)    
 
    real(kind=kreal), pointer :: cont_nforce(:)  
    real(kind=kreal), pointer :: cont_fric(:)    
    real(kind=kreal), pointer :: cont_relvel(:)  
    real(kind=kreal), pointer :: cont_state(:)   
    integer(kind=kint), pointer :: cont_sgrp_id(:)
    real(kind=kreal), pointer :: cont_area(:)    
    real(kind=kreal), pointer :: cont_ntrac(:)   
    real(kind=kreal), pointer :: cont_ftrac(:)   
    real(kind=kreal), pointer :: embed_nforce(:)  
 
    type(fstr_solid_physic_val), pointer :: solid=>null()     
    type(fstr_solid_physic_val), pointer :: shell=>null()     
    type(fstr_solid_physic_val), pointer :: beam =>null()     
 
    integer(kind=kint) :: restart_nout
    integer(kind=kint) :: restart_nin
    type(step_info)    :: step_ctrl_restart
 
    integer(kind=kint) :: max_lyr
    integer(kind=kint) :: is_33shell
    integer(kind=kint) :: is_33beam
    integer(kind=kint) :: is_heat
    integer(kind=kint) :: max_ncon_stf
    integer(kind=kint) :: max_ncon
    integer(kind=kint), pointer :: is_rot(:) => null()
    integer(kind=kint) :: elemopt361
    logical            :: is_smoothing_active
    real(kind=kreal)   :: factor(2)      
    integer(kind=kint) :: nrstat_i(10)
    real(kind=kreal)   :: nrstat_r(10)     
    integer(kind=kint) :: autoinc_stat
    integer(kind=kint) :: cutback_stat
 
    real(kind=kreal), pointer :: gl(:)           
    real(kind=kreal), pointer :: gl0(:)           
    real(kind=kreal), pointer :: eforce(:)           
    real(kind=kreal), pointer :: qforce(:)           
    real(kind=kreal), pointer :: qforce_bak(:)           
    real(kind=kreal), pointer :: unode(:)      => null()   
    real(kind=kreal), pointer :: unode_bak(:)  => null()   
    real(kind=kreal), pointer :: dunode(:)     => null()   
    real(kind=kreal), pointer :: ddunode(:)    => null()   
    real(kind=kreal), pointer :: temperature(:)=> null()   
    real(kind=kreal), pointer :: temp_bak(:)   => null()
    real(kind=kreal), pointer :: last_temp(:)  => null()
 
    type( telement ), pointer :: elements(:)   =>null()  
    type( tmaterial ),pointer :: materials(:)  =>null()  
    integer                   :: n_contacts
    type( tcontact ), pointer :: contacts(:)   =>null()  
    integer                   :: n_embeds
    type( tcontact ), pointer :: embeds(:)   =>null()  
    integer                   :: n_fix_mpc
    real(kind=kreal), pointer :: mpc_const(:)  =>null()  
    type(tsection), pointer   :: sections(:)   =>null()  
 
    ! for cutback
    ! ####################### Notice #######################
    ! # If you add new variables to store analysis status, #
    ! # - backup variables with postfix "_bkup" here       #
    ! # - backup process to module m_fstr_Cutback          #
    ! # must be added if necessary.                        #
    ! ######################################################
    real(kind=kreal), pointer :: unode_bkup(:)     => null() 
    real(kind=kreal), pointer :: qforce_bkup(:)    => null() 
    real(kind=kreal), pointer :: last_temp_bkup(:) => null()
    type( telement ), pointer :: elements_bkup(:)  =>null()  
    type( tcontact ), pointer :: contacts_bkup(:)  =>null()  
    type( tcontact ), pointer :: embeds_bkup(:)  =>null()  
  end type fstr_solid
 
  type fstr_heat
    integer(kind=kint) :: is_steady
    real(kind=kreal)   :: beta
    logical :: is_iter_max_limit
 
    integer(kind=kint) :: steptot
    integer(kind=kint) :: restart_nout
    real(kind=kreal), pointer :: step_dltime(:), step_eetime(:)
    real(kind=kreal), pointer :: step_delmin(:), step_delmax(:)
    integer(kind=kint) :: timepoint_id
 
    integer(kind=kint) :: materialtot
    integer(kind=kint), pointer :: rhotab(:), cptab(:), condtab(:)
    real(kind=kreal), pointer :: rho(:,:), rhotemp(:,:)
    real(kind=kreal), pointer :: cp(:,:),  cptemp(:,:)
    real(kind=kreal), pointer :: cond(:,:),condtemp(:,:)
 
    real(kind=kreal), pointer :: rhofunca(:,:),  rhofuncb(:,:)
    real(kind=kreal), pointer :: cpfunca(:,:),  cpfuncb(:,:)
    real(kind=kreal), pointer :: condfunca(:,:),condfuncb(:,:)
 
    integer(kind=kint) :: amplitudetot
    integer(kind=kint), pointer :: ampltab(:)
    real(kind=kreal), pointer :: ampl(:,:), ampltime(:,:)
    real(kind=kreal), pointer :: amplfunca(:,:),  amplfuncb(:,:)
 
    real(kind=kreal), pointer :: temp0(:)
    real(kind=kreal), pointer :: tempc(:)
    real(kind=kreal), pointer :: temp(:)
 
    integer(kind=kint) :: t_fix_tot
    integer(kind=kint), pointer :: t_fix_node(:)
    integer(kind=kint), pointer :: t_fix_ampl(:)
    real(kind=kreal), pointer :: t_fix_val(:)
 
    integer(kind=kint) :: q_nod_tot
    integer(kind=kint), pointer :: q_nod_node(:)
    integer(kind=kint), pointer :: q_nod_ampl(:)
    real(kind=kreal), pointer :: q_nod_val(:)
 
    integer(kind=kint) :: q_vol_tot
    integer(kind=kint), pointer :: q_vol_elem(:)
    integer(kind=kint), pointer :: q_vol_ampl(:)
    real(kind=kreal), pointer :: q_vol_val(:)
 
    integer(kind=kint) :: q_suf_tot
    integer(kind=kint), pointer :: q_suf_elem(:)
    integer(kind=kint), pointer :: q_suf_ampl(:)
    integer(kind=kint), pointer :: q_suf_surf(:)
    real(kind=kreal), pointer :: q_suf_val(:)
 
    integer(kind=kint) :: r_suf_tot
    integer(kind=kint), pointer :: r_suf_elem(:)
    integer(kind=kint), pointer :: r_suf_ampl(:,:)
    integer(kind=kint), pointer :: r_suf_surf(:)
    real(kind=kreal), pointer :: r_suf_val(:,:)
 
    integer(kind=kint) :: h_suf_tot
    integer(kind=kint), pointer :: h_suf_elem(:)
    integer(kind=kint), pointer :: h_suf_ampl(:,:)
    integer(kind=kint), pointer :: h_suf_surf(:)
    real(kind=kreal), pointer :: h_suf_val(:,:)
 
    integer(kind=kint)  :: wl_tot
    type(tweldline), pointer :: weldline(:) => null()
 
    type(telemact) :: elemact
  end type fstr_heat
 
  type fstr_dynamic
    integer(kind=kint) :: idx_eqa       ! implicit or explicit
    integer(kind=kint) :: idx_resp      ! time history or steady-state harmonic response analysis
 
    integer(kind=kint) :: n_step        ! total step number of analysis
    real(kind=kreal)   :: t_start       ! start time of analysis
    real(kind=kreal)   :: t_curr        ! current time of analysis
    real(kind=kreal)   :: t_end         ! end time of analysis
    real(kind=kreal)   :: t_delta       ! time increment
    integer(kind=kint) :: restart_nout  ! output interval of restart file
    ! (if  .gt.0) restart file write
    ! (if  .lt.0) restart file read and write
    integer(kind=kint) :: restart_nin  !input number of restart file
 
    real(kind=kreal)   :: gamma         ! Newmark-beta parameter gamma
    real(kind=kreal)   :: beta          ! Newmark-beta parameter beta
 
    integer(kind=kint) :: idx_mas       ! mass matrix type
 
    integer(kind=kint) :: idx_dmp      ! damping type
    real(kind=kreal)   :: ray_m         ! Rayleigh damping parameter Rm
    real(kind=kreal)   :: ray_k         ! Rayleigh damping parameter Rk
 
    logical           :: varinitialize  ! initialization flag
 
    integer(kind=kint) :: nout           ! output interval of result
    integer(kind=kint) :: ngrp_monit     ! node of monitoring result
    integer(kind=kint) :: nout_monit     ! output interval of result monitoring
    integer(kind=kint) :: i_step         ! step number
    integer(kind=kint) :: iout_list(6)   ! 0:not output  1:output
    ! iout_list(1): displacement
    ! iout_list(2): velocity
    ! iout_list(3): acceleration
    ! iout_list(4): reaction force
    ! iout_list(5): strain
    ! iout_list(6): stress
 
    real(kind=kreal), pointer :: disp(:,:)
    real(kind=kreal), pointer :: vel(:,:)
    real(kind=kreal), pointer :: acc(:,:)
 
    real(kind=kreal) :: kineticenergy
    real(kind=kreal) :: strainenergy
    real(kind=kreal) :: totalenergy
 
    real(kind=kreal), pointer :: vec1(:)
    real(kind=kreal), pointer :: vec2(:)
    real(kind=kreal), pointer :: vec3(:)
 
    integer(kind=kint) :: dynamic_iw4        =   204
    integer(kind=kint) :: dynamic_iw5        =   205
    integer(kind=kint) :: dynamic_iw6        =   206
    integer(kind=kint) :: dynamic_iw7        =   207
    integer(kind=kint) :: dynamic_iw8        =   208
    integer(kind=kint) :: dynamic_iw9        =   209
    integer(kind=kint) :: dynamic_iw10       =   210
  end type fstr_dynamic
 
  type fstr_freqanalysis
    integer(kind=kint)                :: fload_ngrp_tot
    integer(kind=kint), pointer       :: fload_ngrp_grpid(:) => null()
    integer(kind=kint), pointer       :: fload_ngrp_id(:)    => null()
    integer(kind=kint), pointer       :: fload_ngrp_type(:)  => null()
    integer(kind=kint), pointer       :: fload_ngrp_dof(:)   => null()
    real(kind=kreal), pointer         :: fload_ngrp_valre(:) => null()
    real(kind=kreal), pointer         :: fload_ngrp_valim(:) => null()
    character(len=HECMW_FILENAME_LEN) :: eigenlog_filename
    integer(kind=kint)                :: start_mode
    integer(kind=kint)                :: end_mode
  end type fstr_freqanalysis
 
  type fstr_freqanalysis_data
    integer(kind=kint)        :: nummode
    integer(kind=kint)        :: numnodedof
    real(kind=kreal), pointer :: eigomega(:)    => null()
    real(kind=kreal), pointer :: eigvector(:,:) => null()
    real(kind=kreal)           :: rayalpha, raybeta
  end type fstr_freqanalysis_data
 
  type fstr_eigen
    integer(kind=kint)  :: nget      ! Solved eigen value number (default:5)
    integer(kind=kint)  :: maxiter   ! Max. Lcz iterations (default:60)
    integer(kind=kint)  :: iter      ! Max. Lcz iterations (default:60)
    real   (kind=kreal) :: sigma     ! 0.0
    real   (kind=kreal) :: tolerance ! Lcz tolerance (default:1.0e-8)
    real   (kind=kreal) :: totalmass
    real   (kind=kreal), pointer :: eigval(:)
    real   (kind=kreal), pointer :: eigvec(:,:)
    real   (kind=kreal), pointer :: filter(:)
    real   (kind=kreal), pointer :: mass(:)
    real   (kind=kreal), pointer :: effmass(:)
    real   (kind=kreal), pointer :: partfactor(:)
    logical :: is_free = .false.
  end type fstr_eigen
 
  type fstr_couple
    integer( kind=kint )   :: dof              ! == 3
    integer( kind=kint )   :: ndof             ! total dof (coupled_node_n*dof)
    integer( kind=kint )   :: coupled_node_n
    integer, pointer       :: coupled_node(:)  ! local node id sent from revocap
    real( kind=8 ),pointer :: trac(:)          ! input  (x,y,z,x,y,z ... )
    real( kind=8 ),pointer :: disp(:)          ! output (x,y,z,x,y,z ... )
    real( kind=8 ),pointer :: velo(:)          ! output (x,y,z,x,y,z ... )
    real( kind=8 ),pointer :: accel(:)         ! output (x,y,z,x,y,z ... )
    integer( kind=kint ),pointer :: index(:)   ! size:total node num.
  end type fstr_couple
 
  type tweldline
    integer(kind=kint) :: egrpid
    real( kind=kreal ) :: i
    real( kind=kreal ) :: u
    real( kind=kreal ) :: coe
    real( kind=kreal ) :: v
    integer(kind=kint) :: xyz
    real(kind=kreal)   :: n1, n2
    real(kind=kreal)   :: distol
    real(kind=kreal)   :: tstart
  end type tweldline
 
  type tsection
    !integer              :: mat_ID
    !integer              :: iset
    !integer              :: orien_ID
    !real(kind=kreal)     :: thickness
    integer              :: elemopt341
    !integer              :: elemopt342
    !integer              :: elemopt351
    !integer              :: elemopt352
    integer              :: elemopt361
    !integer              :: elemopt362
  end type tsection
 
contains
 
  subroutine fstr_nullify_fstr_param( P )
    implicit none
    type( fstr_param ) :: P
 
    nullify( p%dtime )
    nullify( p%etime )
    nullify( p%dtmin )
    nullify( p%delmax )
    nullify( p%itmax )
    nullify( p%eps )
    nullify( p%global_local_ID)
    nullify( p%timepoints )
  end subroutine fstr_nullify_fstr_param
 
  subroutine fstr_nullify_fstr_solid( S )
    implicit none
    type( fstr_solid ) :: S
 
    nullify( s%BOUNDARY_ngrp_ID )
    nullify( s%BOUNDARY_ngrp_type )
    nullify( s%BOUNDARY_ngrp_amp )
    nullify( s%BOUNDARY_ngrp_val)
    nullify( s%BOUNDARY_ngrp_rotID )
    nullify( s%BOUNDARY_ngrp_centerID )
    nullify( s%CLOAD_ngrp_ID )
    nullify( s%CLOAD_ngrp_DOF )
    nullify( s%CLOAD_ngrp_amp )
    nullify( s%CLOAD_ngrp_rotID )
    nullify( s%CLOAD_ngrp_centerID )
    nullify( s%CLOAD_ngrp_val )
    nullify( s%DLOAD_ngrp_ID )
    nullify( s%DLOAD_ngrp_LID )
    nullify( s%DLOAD_ngrp_amp )
    nullify( s%DLOAD_ngrp_params )
    nullify( s%TEMP_ngrp_ID )
    nullify( s%TEMP_ngrp_val )
    nullify( s%SPRING_ngrp_ID )
    nullify( s%SPRING_ngrp_DOF )
    nullify( s%SPRING_ngrp_amp )
    nullify( s%SPRING_ngrp_val )
    nullify( s%STRESS )
    nullify( s%STRAIN )
    nullify( s%MISES )
    nullify( s%PSTRESS )
    nullify( s%PSTRAIN )
    nullify( s%PSTRESS_VECT )
    nullify( s%PSTRAIN_VECT )
    nullify( s%REACTION )
    nullify( s%ESTRESS )
    nullify( s%ESTRAIN )
    nullify( s%EMISES )
    nullify( s%EPLSTRAIN )
    nullify( s%EPSTRESS )
    nullify( s%EPSTRAIN )
    nullify( s%EPSTRESS_VECT )
    nullify( s%EPSTRAIN_VECT )
    nullify( s%ENQM )
    nullify( s%GL          )
    nullify( s%GL0         )
    nullify( s%QFORCE      )
    nullify( s%VELOCITY_ngrp_ID )
    nullify( s%VELOCITY_ngrp_type )
    nullify( s%VELOCITY_ngrp_amp )
    nullify( s%VELOCITY_ngrp_val )
    nullify( s%ACCELERATION_ngrp_ID )
    nullify( s%ACCELERATION_ngrp_type )
    nullify( s%ACCELERATION_ngrp_amp )
    nullify( s%ACCELERATION_ngrp_val )
    nullify( s%COUPLE_ngrp_ID )
  end subroutine fstr_nullify_fstr_solid
 
  subroutine fstr_nullify_fstr_heat( H )
    implicit none
    type( fstr_heat ) :: H
 
    nullify( h%STEP_DLTIME )
    nullify( h%STEP_EETIME )
    nullify( h%STEP_DELMIN )
    nullify( h%STEP_DELMAX )
    nullify( h%RHO )
    nullify( h%RHOtemp  )
    nullify( h%CP )
    nullify( h%CPtemp  )
    nullify( h%COND )
    nullify( h%CONDtemp  )
    nullify( h%RHOtab )
    nullify( h%CPtab )
    nullify( h%CONDtab )
    nullify( h%RHOfuncA )
    nullify( h%RHOfuncB )
    nullify( h%CPfuncA  )
    nullify( h%CPfuncB )
    nullify( h%CONDfuncA  )
    nullify( h%CONDfuncB )
    nullify( h%AMPL )
    nullify( h%AMPLtime  )
    nullify( h%AMPLtab )
    nullify( h%AMPLfuncA )
    nullify( h%AMPLfuncB )
    nullify( h%TEMP0 )
    nullify( h%TEMPC )
    nullify( h%TEMP  )
    nullify( h%T_FIX_node )
    nullify( h%T_FIX_ampl )
    nullify( h%T_FIX_val )
    nullify( h%Q_NOD_node )
    nullify( h%Q_NOD_ampl )
    nullify( h%Q_NOD_val )
    nullify( h%Q_VOL_elem )
    nullify( h%Q_VOL_ampl )
    nullify( h%Q_VOL_val )
    nullify( h%Q_SUF_elem )
    nullify( h%Q_SUF_ampl )
    nullify( h%Q_SUF_surf )
    nullify( h%Q_SUF_val )
    nullify( h%R_SUF_elem )
    nullify( h%R_SUF_ampl )
    nullify( h%R_SUF_surf )
    nullify( h%R_SUF_val )
    nullify( h%H_SUF_elem )
    nullify( h%H_SUF_ampl )
    nullify( h%H_SUF_surf )
    nullify( h%H_SUF_val )
  end subroutine fstr_nullify_fstr_heat
 
  subroutine fstr_nullify_fstr_dynamic( DY )
    implicit none
    type( fstr_dynamic ) :: DY
 
    nullify( dy%DISP )
    nullify( dy%VEL  )
    nullify( dy%ACC  )
    nullify( dy%VEC1 )
    nullify( dy%VEC2 )
    nullify( dy%VEC3 )
  end subroutine fstr_nullify_fstr_dynamic
 
  subroutine fstr_nullify_fstr_freqanalysis( f )
    implicit none
    type( fstr_freqanalysis ), intent(inout) :: f
 
    f%FLOAD_ngrp_tot = 0
    nullify( f%FLOAD_ngrp_GRPID )
    nullify( f%FLOAD_ngrp_ID )
    nullify( f%FLOAD_ngrp_TYPE )
    nullify( f%FLOAD_ngrp_DOF )
    nullify( f%FLOAD_ngrp_valre )
    nullify( f%FLOAD_ngrp_valim )
  end subroutine fstr_nullify_fstr_freqanalysis
 
  subroutine fstr_nullify_fstr_eigen( E )
    implicit none
    type( fstr_eigen ) :: E
 
    nullify( e%mass )
  end subroutine fstr_nullify_fstr_eigen
 
  subroutine fstr_nullify_fstr_couple( C )
    implicit none
    type( fstr_couple ) :: C
 
    nullify( c%coupled_node )
    nullify( c%trac )
    nullify( c%velo )
    nullify( c%index )
  end subroutine fstr_nullify_fstr_couple
 
  subroutine fstr_mat_init( hecMAT )
    implicit none
    type(hecmwst_matrix) :: hecMAT
 
    hecmat%Iarray(1) =  100    ! = nier
    hecmat%Iarray(2) =    1    ! = method
    hecmat%Iarray(3) =    1    ! = precond
    hecmat%Iarray(4) =    0    ! = nset
    hecmat%Iarray(5) =    1    ! = iterpremax
    hecmat%Iarray(6) =   10    ! = nrest
    hecmat%Iarray(7) =    0    ! = scaling
    hecmat%Iarray(21)=  kno    ! = iterlog
    hecmat%Iarray(22)=  kno    ! = timelog
    hecmat%Iarray(31)=    0    ! = dumptype
    hecmat%Iarray(32)=    0    ! = dumpexit
    hecmat%Iarray(33)=    0    ! = usejad
    hecmat%Iarray(34)=   10    ! = ncolor_in
    hecmat%Iarray(13)=    0    ! = mpc_method
    hecmat%Iarray(14)=    0    ! = estcond
    hecmat%Iarray(35)=    3    ! = maxrecycle_precond
    hecmat%Iarray(41)=    0    ! = solver_opt1
    hecmat%Iarray(42)=    0    ! = solver_opt2
    hecmat%Iarray(43)=    0    ! = solver_opt3
    hecmat%Iarray(44)=    0    ! = solver_opt4
    hecmat%Iarray(45)=    0    ! = solver_opt5
    hecmat%Iarray(46)=    0    ! = solver_opt6
 
    hecmat%Rarray(1) =  1.0e-8 ! = resid
    hecmat%Rarray(2) =  1.0    ! = sigma_diag
    hecmat%Rarray(3) =  0.0    ! = sigma
    hecmat%Rarray(4) =  0.1    ! = thresh
    hecmat%Rarray(5) =  0.1    ! = filter
    hecmat%Rarray(11)=  1.0e+4 ! = penalty
 
    hecmat%Iarray(96) =   0    ! nrecycle_precond
    hecmat%Iarray(97) = kyes   ! flag_numfact
    hecmat%Iarray(98) = kyes   ! flag_symbfact
    hecmat%Iarray(99) = kyes   ! indirect method
  end subroutine fstr_mat_init
 
  subroutine hecmat_init( hecMAT )
    implicit none
    type( hecmwst_matrix ) :: hecMAT
    integer ::  ndof, nn, ierror
    ndof = hecmat%NDOF
    nn = ndof*ndof
    allocate (hecmat%AL(nn*hecmat%NPL)        ,stat=ierror )
    if( ierror /= 0 ) then
      write(*,*) "##ERROR : not enough memory"
      write(idbg,*) 'stop due to allocation error'
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm() )
    end if
    allocate (hecmat%AU(nn*hecmat%NPU)        ,stat=ierror )
    if( ierror /= 0 ) then
      write(*,*) "##ERROR : not enough memory"
      write(idbg,*) 'stop due to allocation error'
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm() )
    end if
    allocate (hecmat%B(ndof*hecmat%NP)          ,stat=ierror )
    if( ierror /= 0 ) then
      write(*,*) "##ERROR : not enough memory"
      write(idbg,*) 'stop due to allocation error'
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm() )
    end if
    hecmat%B(:)=0.d0
    allocate (hecmat%D(nn*hecmat%NP)          ,stat=ierror )
    if( ierror /= 0 ) then
      write(*,*) "##ERROR : not enough memory"
      write(idbg,*) 'stop due to allocation error'
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm() )
    end if
    allocate (hecmat%X(ndof*hecmat%NP)          ,stat=ierror )
    if( ierror /= 0 ) then
      write(*,*) "##ERROR : not enough memory"
      write(idbg,*) 'stop due to allocation error'
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm() )
    end if
#ifdef _OPENACC
    allocate (hecmat%A(nn*(hecmat%NPA))        ,stat=ierror )
    if( ierror /= 0 ) then
      write(*,*) "##ERROR : not enough memory"
      write(idbg,*) 'stop due to allocation error'
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm() )
    endif
#endif
    hecmat%D  = 0.0d0
    hecmat%AL = 0.0d0
    hecmat%AU = 0.0d0
    hecmat%B  = 0.0d0
    hecmat%X  = 0.0d0
#ifdef _OPENACC
    hecmat%A  = 0.0d0
#endif
  end subroutine hecmat_init
 
  subroutine hecmat_finalize( hecMAT )
    implicit none
    type( hecmwst_matrix ) :: hecMAT
    integer ::  ndof, nn, ierror
    ndof = hecmat%NDOF
    nn = ndof*ndof
    if( associated(hecmat%AL) ) then
      deallocate(hecmat%AL                  ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(hecmat%AU) ) then
      deallocate(hecmat%AU                  ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(hecmat%B) ) then
      deallocate(hecmat%B                   ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(hecmat%D) ) then
      deallocate(hecmat%D                   ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(hecmat%X) ) then
      deallocate(hecmat%X                   ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(hecmat%A) ) then
      deallocate(hecmat%A                   ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
  end subroutine hecmat_finalize
 
  subroutine fstr_param_init( fstrPARAM, hecMESH )
    implicit none
    type(fstr_param) :: fstrPARAM
    type(hecmwst_local_mesh) :: hecMESH
    integer(kind=kint) :: i
    external fstr_sort_index
 
    fstrparam%solution_type = kststatic
    fstrparam%nlgeom        = .false.
    fstrparam%solver_method = ksmcg
    fstrparam%nlsolver_method = knsmnewton
 
    !!STATIC !HEAT
    fstrparam%analysis_n = 0
    fstrparam%ref_temp   = 0
 
    ! output control
    fstrparam%fg_echo       = koff
    fstrparam%fg_result     = koff
    fstrparam%fg_visual     = koff
 
    ! for heat ...
    fstrparam%fg_neutral    = koff
    fstrparam%fg_irres      = kno
    fstrparam%fg_iwres      = kno
    fstrparam%nrres         = 1
    fstrparam%nprint        = 1
 
    ! for couple
    fstrparam%fg_couple      = 0
    fstrparam%fg_couple_type = 0
    fstrparam%fg_couple_first= 0
    fstrparam%fg_couple_window= 0
 
    ! for restart control
    fstrparam%restart_version = 5
 
    ! for contact analysis
    fstrparam%contact_algo = kcaslagrange  ! default: Standard Lagrange
    fstrparam%augiter = 2                  ! default augmentation iteration for ALagrange
 
    ! index table for global node ID sorting
 
    fstrparam%n_node = hecmesh%n_node;
    fstrparam%nn_internal = hecmesh%nn_internal;
    allocate( fstrparam%global_local_ID(2,hecmesh%nn_internal))
    do i = 1, hecmesh%nn_internal
      fstrparam%global_local_ID(1,i) = hecmesh%global_node_ID(i)
      fstrparam%global_local_ID(2,i) = i
    end do
    call fstr_sort_index( fstrparam%global_local_ID, hecmesh%nn_internal )
  end subroutine fstr_param_init
 
  logical function fstr_isboundaryactive( fstrSOLID, nbc, cstep )
    type(fstr_solid)    :: fstrsolid
    integer, intent(in) :: nbc
    integer, intent(in) :: cstep
    fstr_isboundaryactive = .true.
    if( .not. associated(fstrsolid%step_ctrl) ) return
    if( cstep>fstrsolid%nstep_tot ) return
    fstr_isboundaryactive = isboundaryactive( nbc, fstrsolid%step_ctrl(cstep) )
  end function
 
  logical function fstr_isloadactive( fstrSOLID, nbc, cstep )
    type(fstr_solid)    :: fstrsolid
    integer, intent(in) :: nbc
    integer, intent(in) :: cstep
    fstr_isloadactive = .true.
    if( cstep > 0 ) then
      if( .not. associated(fstrsolid%step_ctrl) ) return
      if( cstep>fstrsolid%nstep_tot ) return
      fstr_isloadactive = isloadactive( nbc, fstrsolid%step_ctrl(cstep) )
    else
      fstr_isloadactive = isloadactive( nbc, fstrsolid%step_ctrl_restart )
    endif
  end function
 
  logical function fstr_iscontactactive( fstrSOLID, nbc, cstep )
    type(fstr_solid)     :: fstrsolid
    integer, intent(in) :: nbc
    integer, intent(in) :: cstep
    fstr_iscontactactive = .true.
    if( .not. associated(fstrsolid%step_ctrl) ) return
    if( cstep>fstrsolid%nstep_tot ) return
    fstr_iscontactactive = iscontactactive( nbc, fstrsolid%step_ctrl(cstep) )
  end function
 
  logical function fstr_isembedactive( fstrSOLID, nbc, cstep )
    type(fstr_solid)     :: fstrsolid
    integer, intent(in) :: nbc
    integer, intent(in) :: cstep
    fstr_isembedactive = .true.
    if( .not. associated(fstrsolid%step_ctrl) ) return
    if( cstep>fstrsolid%nstep_tot ) return
    fstr_isembedactive = iscontactactive( nbc, fstrsolid%step_ctrl(cstep) )
  end function
 
  logical function fstr_iselemactivationactive( fstrSOLID, nbc, cstep )
    type(fstr_solid)    :: fstrsolid
    integer, intent(in) :: nbc
    integer, intent(in) :: cstep
    fstr_iselemactivationactive = .true.
    if( .not. associated(fstrsolid%step_ctrl) ) return
    if( cstep>fstrsolid%nstep_tot ) return
    fstr_iselemactivationactive = iselemactivationactive( nbc, fstrsolid%step_ctrl(cstep) )
  end function
 
  subroutine get_coordsys( cdsys_ID, hecMESH, fstrSOLID, coords )
    integer, intent(in)             :: cdsys_ID
    type(hecmwst_local_mesh)       :: hecMESH
    type(fstr_solid), intent(inout) :: fstrSOLID
    real(kind=kreal), intent(out)   :: coords(3,3)
    integer :: ik
 
    coords = 0.d0
    if( cdsys_id>0 ) then
      if( iscoordneeds(g_localcoordsys(cdsys_id)) ) then
        coords=g_localcoordsys(cdsys_id)%CoordSys
      else
        ik=g_localcoordsys(cdsys_id)%node_ID(1)
        coords(1,:)= hecmesh%node(3*ik-2:3*ik)+fstrsolid%unode(3*ik-2:3*ik)  &
          + fstrsolid%dunode(3*ik-2:3*ik)
        ik=g_localcoordsys(cdsys_id)%node_ID(2)
        coords(2,:)= hecmesh%node(3*ik-2:3*ik)+fstrsolid%unode(3*ik-2:3*ik)  &
          + fstrsolid%dunode(3*ik-2:3*ik)
        ik=g_localcoordsys(cdsys_id)%node_ID(3)
        if(ik>0) coords(3,:)= hecmesh%node(3*ik-2:3*ik)+fstrsolid%unode(3*ik-2:3*ik)  &
          + fstrsolid%dunode(3*ik-2:3*ik)
      endif
    endif
  end subroutine get_coordsys
 
  subroutine fstr_set_current_config_to_mesh(hecMESH,fstrSOLID,coord)
    implicit none
    type(hecmwst_local_mesh), intent(inout) :: hecMESH
    type (fstr_solid), intent(in) :: fstrSOLID
    real(kind=kreal), intent(out) :: coord(:)
    integer(kind=kint) :: i
    if(hecmesh%n_dof == 4) return
    do i = 1, hecmesh%nn_internal*min(hecmesh%n_dof,3)
      coord(i) = hecmesh%node(i)
      hecmesh%node(i) = coord(i)+fstrsolid%unode(i)+fstrsolid%dunode(i)
    enddo
  end subroutine fstr_set_current_config_to_mesh
 
  subroutine fstr_recover_initial_config_to_mesh(hecMESH,fstrSOLID,coord)
    implicit none
    type(hecmwst_local_mesh), intent(inout) :: hecMESH
    type (fstr_solid), intent(in) :: fstrSOLID
    real(kind=kreal), intent(in) :: coord(:)
    integer(kind=kint) :: i
    if(hecmesh%n_dof == 4) return
    do i = 1, hecmesh%nn_internal*min(hecmesh%n_dof,3)
      hecmesh%node(i) = coord(i)
    enddo
  end subroutine fstr_recover_initial_config_to_mesh
 
  subroutine fstr_solid_phys_zero(phys)
    implicit none
    type(fstr_solid_physic_val), pointer :: phys
    phys%STRAIN  = 0.0d0
    phys%STRESS  = 0.0d0
    phys%MISES   = 0.0d0
    phys%ESTRAIN = 0.0d0
    phys%ESTRESS = 0.0d0
    phys%EMISES  = 0.0d0
    phys%EPLSTRAIN  = 0.0d0
    phys%ENQM    = 0.0d0
  end subroutine fstr_solid_phys_zero
 
  subroutine fstr_solid_phys_clear(fstrSOLID)
    implicit none
    type (fstr_solid)         :: fstrSOLID
    integer(kind=kint) :: i
 
    if (associated(fstrsolid%SOLID)) then
      call fstr_solid_phys_zero(fstrsolid%SOLID)
    end if
    if (associated(fstrsolid%SHELL)) then
      call fstr_solid_phys_zero(fstrsolid%SHELL)
      do i=1,fstrsolid%max_lyr
        call fstr_solid_phys_zero(fstrsolid%SHELL%LAYER(i)%PLUS)
        call fstr_solid_phys_zero(fstrsolid%SHELL%LAYER(i)%MINUS)
      end do
    end if
    if (associated(fstrsolid%BEAM)) then
      call fstr_solid_phys_zero(fstrsolid%BEAM)
    end if
  end subroutine fstr_solid_phys_clear
 
  subroutine table_amp(hecMESH, fstrSOLID, cstep, jj_n_amp, time, f_t)
    type ( hecmwST_local_mesh ), intent(in) :: hecMESH
    type ( fstr_solid         ), intent(in) :: fstrSOLID
    integer(kind=kint), intent(in)          :: cstep
    integer(kind=kint), intent(in)          :: jj_n_amp
    real(kind=kreal), intent(in)            :: time       
    real(kind=kreal), intent(out)           :: f_t        
 
    integer(kind=kint) :: i
    integer(kind=kint) :: jj1, jj2
    integer(kind=kint) :: s1, s2, flag
    real(kind=kreal) :: t_1, t_2, t_t, f_1, f_2, tincre
 
    s1 = 0; s2 = 0
 
    if( jj_n_amp <= 0 ) then  ! Amplitude not defined
      if(myrank == 0) then
        write(imsg,*) 'internal error: amplitude table not defined'
      endif
      call hecmw_abort( hecmw_comm_get_comm() )
    endif
 
    tincre = fstrsolid%step_ctrl( cstep )%initdt
    jj1 = hecmesh%amp%amp_index(jj_n_amp - 1)
    jj2 = hecmesh%amp%amp_index(jj_n_amp)
 
    jj1 = jj1 + 2
    t_t = time-fstrsolid%step_ctrl(cstep)%starttime
 
    !      if(jj2 .eq. 0) then
    !         f_t = 1.0
    if(t_t .gt. hecmesh%amp%amp_table(jj2)) then
      f_t = hecmesh%amp%amp_val(jj2)
    else if(t_t .le. hecmesh%amp%amp_table(jj2)) then
      flag=0
      do i = jj1, jj2
        if(t_t .le. hecmesh%amp%amp_table(i)) then
          s2 = i
          s1 = i - 1
          flag = 1
        endif
        if( flag == 1 ) exit
      end do
 
      t_2 = hecmesh%amp%amp_table(s2)
      t_1 = hecmesh%amp%amp_table(s1)
      f_2 = hecmesh%amp%amp_val(s2)
      f_1 = hecmesh%amp%amp_val(s1)
      if( t_2-t_1 .lt. 1.0e-20) then
        if(myrank == 0) then
          write(imsg,*) 'stop due to t_2-t_1 <= 0'
        endif
        call hecmw_abort( hecmw_comm_get_comm())
      endif
      f_t = ((t_2*f_1 - t_1*f_2) + (f_2 - f_1)*t_t) / (t_2 - t_1)
    endif
 
  end subroutine table_amp
 
end module m_fstr