fstr_setup.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 m_fstr_setup
  use m_fstr
  use hecmw_setup_util, only &
    : node_grp_name_to_id_ex &
    , elem_grp_name_to_id_ex &
    , surf_grp_name_to_id_ex &
    , hecmw_streqr
  use fstr_setup_util
  use fstr_ctrl_common
  use fstr_ctrl_static
  use fstr_ctrl_heat
  use fstr_ctrl_eigen
  use fstr_ctrl_dynamic
  use fstr_ctrl_material
  use mcontactparam
  use m_static_get_prop
  use m_out
  use m_step
  use m_utilities
  use fstr_ctrl_util_f
  implicit none
 
  type fstr_param_pack
    type(hecmwst_local_mesh), pointer :: mesh
    type(fstr_param), pointer         :: param
    type(fstr_solid), pointer         :: solid
    type(fstr_heat), pointer          :: heat
    type(fstr_eigen), pointer           :: eigen
    type(fstr_dynamic), pointer       :: dyn
    type(fstr_couple), pointer        :: cpl
    type(fstr_freqanalysis), pointer  :: freq
  end type fstr_param_pack
 
contains
 
  !=============================================================================!
  !=============================================================================!
  subroutine fstr_setup( cntl_filename, hecMESH, fstrPARAM,  &
      fstrSOLID, fstrEIG, fstrHEAT, fstrDYNAMIC, fstrCPL, fstrFREQ )
    use mmaterial
    character(len=HECMW_FILENAME_LEN) :: cntl_filename, input_filename
    type(hecmwst_local_mesh),target :: hecMESH
    type(fstr_param),target   :: fstrPARAM
    type(fstr_solid),target   :: fstrSOLID
    type(fstr_eigen),target     :: fstrEIG
    type(fstr_heat),target    :: fstrHEAT
    type(fstr_dynamic),target :: fstrDYNAMIC
    type(fstr_couple),target  :: fstrCPL
    type(fstr_freqanalysis), target :: fstrFREQ
 
    integer(kind=kint) :: ctrl, ctrl_list(20), ictrl
    type(fstr_param_pack) :: P
 
    integer, parameter :: MAXOUTFILE = 10
    double precision, parameter :: dpi = 3.14159265358979323846d0
 
    integer(kind=kint) :: version, result, visual, femap, n_totlyr
    integer(kind=kint) :: rcode, n, i, cid, nout, ierror, cparam_id
    character(len=HECMW_NAME_LEN) :: header_name
    real(kind=kreal) :: ee, pp, rho, alpha, thick, alpha_over_mu
    real(kind=kreal) :: beam_radius,                          &
      beam_angle1, beam_angle2, beam_angle3,&
      beam_angle4, beam_angle5, beam_angle6
    logical          :: isOK
    type(t_output_ctrl) :: outctrl
    type(tshellmat),pointer :: shmat(:)
    character(len=HECMW_FILENAME_LEN) :: logfileNAME, mName, mName2, mName3
 
    ! counters
    integer(kind=kint) :: c_solution, c_solver, c_nlsolver, c_step, c_echo, c_amplitude
    integer(kind=kint) :: c_static, c_boundary, c_cload, c_dload, c_temperature, c_reftemp, c_spring, c_elemact
    integer(kind=kint) :: c_heat, c_fixtemp, c_cflux, c_dflux, c_sflux, c_film, c_sfilm, c_radiate, c_sradiate
    integer(kind=kint) :: c_eigen, c_contact, c_contactparam, c_embed, c_contact_if
    integer(kind=kint) :: c_dynamic, c_velocity, c_acceleration
    integer(kind=kint) :: c_fload, c_facc, c_eigenread
    integer(kind=kint) :: c_couple, c_material
    integer(kind=kint) :: c_mpc, c_weldline, c_initial
    integer(kind=kint) :: c_istep, c_localcoord, c_section
    integer(kind=kint) :: c_elemopt, c_aincparam, c_timepoints, c_cnvparam
    integer(kind=kint) :: c_output, islog
    integer(kind=kint) :: k
    integer(kind=kint) :: ss_nslav, ss_nmas_vis, ss_nmas_glob, ss_bad   ! MORTAR=YES visibility guard scratch
    integer(kind=kint) :: cache = 1
 
    write( logfilename, '(i5,''.log'')' ) myrank
 
    ! packaging
    p%MESH   => hecmesh
    p%PARAM  => fstrparam
    p%SOLID  => fstrsolid
    p%EIGEN  => fstreig
    p%HEAT   => fstrheat
    p%DYN    => fstrdynamic
    p%CPL    => fstrcpl
    p%FREQ   => fstrfreq
 
    fstrparam%contact_algo = kcaalagrange
 
    c_solution = 0; c_solver   = 0; c_nlsolver = 0; c_step   = 0; c_output = 0; c_echo = 0; c_amplitude = 0
    c_static   = 0; c_boundary = 0; c_cload  = 0; c_dload = 0; c_temperature = 0; c_reftemp = 0; c_spring = 0;
    c_elemact = 0;
    c_heat     = 0; c_fixtemp  = 0; c_cflux  = 0; c_dflux = 0; c_sflux = 0
    c_film     = 0; c_sfilm    = 0; c_radiate= 0; c_sradiate = 0
    c_eigen    = 0; c_contact  = 0; c_contactparam = 0; c_embed  = 0; c_contact_if = 0
    c_dynamic  = 0; c_velocity = 0; c_acceleration = 0
    c_couple   = 0; c_material = 0; c_section =0
    c_mpc      = 0; c_weldline = 0; c_initial = 0
    c_istep    = 0; c_localcoord = 0
    c_fload    = 0; c_facc = 0; c_eigenread = 0
    c_elemopt  = 0;
    c_aincparam= 0; c_timepoints = 0; c_cnvparam = 0
 
    ctrl_list = 0
    ictrl = 1
    ctrl  = fstr_ctrl_open( cntl_filename )
    if( ctrl < 0 ) then
      write(*,*) '### Error: Cannot open FSTR control file : ', cntl_filename
      write(ilog,*) '### Error: Cannot open FSTR control file : ', cntl_filename
      stop hecmw_exit_input
    end if
 
    version =0
    do
      rcode = fstr_ctrl_get_c_h_name( ctrl, header_name, hecmw_name_len )
      if(     header_name == '!VERSION' ) then
        rcode = fstr_ctrl_get_data_ex( ctrl, 1, 'i ', version )
      else if(     header_name == '!SOLUTION' ) then
        c_solution = c_solution + 1
        call fstr_setup_solution( ctrl, c_solution, p )
      else if(     header_name == '!NONLINEAR_SOLVER' ) then
        c_nlsolver = c_nlsolver + 1
        call fstr_setup_nonlinear_solver( ctrl, c_nlsolver, p )
      else if( header_name == '!SOLVER' ) then
        c_solver = c_solver + 1
        call fstr_setup_solver( ctrl, c_solver, p )
      else if( header_name == '!ISTEP' ) then
        c_istep = c_istep + 1
      else if( header_name == '!STEP' ) then
        if( version==0 ) then
          c_step = c_step + 1
          call fstr_setup_step( ctrl, c_step, p )
        else
          c_istep = c_istep + 1
        endif
      else if( header_name == '!WRITE' ) then
        call fstr_ctrl_get_output( ctrl, outctrl, islog, result, visual, femap )
        if( visual==1 ) p%PARAM%fg_visual= 1
        if( result==1 ) p%PARAM%fg_result = 1
        c_output = c_output+1
      else if( header_name == '!ECHO' ) then
        c_echo = c_echo + 1
        call fstr_setup_echo( ctrl, c_echo, p )
      else if( header_name == '!RESTART' ) then
        call fstr_setup_restart( ctrl, nout, p%PARAM%restart_version )
        fstrsolid%restart_nout= nout
        fstrdynamic%restart_nout= nout
        fstrheat%restart_nout= nout
      else if( header_name == '!ORIENTATION' ) then
        c_localcoord = c_localcoord + 1
      else if( header_name == '!AUTOINC_PARAM' ) then
        c_aincparam = c_aincparam + 1
      else if( header_name == '!CONVERG_PARAM' ) then
        c_cnvparam = c_cnvparam + 1
      else if( header_name == '!TIME_POINTS' ) then
        c_timepoints = c_timepoints + 1
      else if( header_name == '!OUTPUT_SSTYPE' ) then
        call fstr_setup_output_sstype( ctrl, p )
      else if( header_name == '!INITIAL_CONDITION' ) then
        c_initial = c_initial + 1
      else if( header_name == '!AMPLITUDE' ) then
        c_amplitude = c_amplitude + 1
        call fstr_setup_amplitude( ctrl, p )
      else if( header_name == '!ELEMENT_ACTIVATION' ) then
        c_elemact = c_elemact + 1
        call fstr_setup_element_activation( ctrl, c_elemact, p )
 
        !--------------- for static -------------------------
 
      else if( header_name == '!STATIC' ) then
        c_static = c_static + 1
        call fstr_setup_static( ctrl, c_static, p )
      else if( header_name == '!BOUNDARY' ) then
        c_boundary = c_boundary + 1
        call fstr_setup_boundary( ctrl, c_boundary, p )
      else if( header_name == '!CLOAD' ) then
        c_cload = c_cload + 1
        call fstr_setup_cload( ctrl, c_cload, p )
        n = fstr_ctrl_get_data_line_n( ctrl )
      else if( header_name == '!DLOAD' ) then
        c_dload = c_dload + 1
        call fstr_setup_dload( ctrl, c_dload, p )
      else if( header_name == '!CONTACT_ALGO' ) then
        call fstr_setup_contactalgo( ctrl, p )
      else if( header_name == '!CONTACT' ) then
        n = fstr_ctrl_get_data_line_n( ctrl )
        c_contact = c_contact + n
      else if( header_name == '!EMBED' ) then
        n = fstr_ctrl_get_data_line_n( ctrl )
        c_embed = c_embed + n
      else if( header_name == '!CONTACT_PARAM' ) then
        c_contactparam = c_contactparam + 1
      else if( header_name == '!CONTACT_INTERFERENCE' ) then
        n = fstr_ctrl_get_data_line_n( ctrl )
        c_contact_if = c_contact_if + n
      else if( header_name == '!MATERIAL' ) then
        c_material = c_material + 1
      else if( header_name == '!TEMPERATURE' ) then
        c_temperature = c_temperature + 1
        call fstr_setup_temperature( ctrl, c_temperature, p )
      else if( header_name == '!SPRING' ) then
        c_spring = c_spring + 1
        call fstr_setup_spring( ctrl, c_spring, p )
      else if( header_name == '!REFTEMP' ) then
        c_reftemp = c_reftemp + 1
        call fstr_setup_reftemp( ctrl, c_reftemp, p )
 
        !--------------- for heat -------------------------
 
      else if( header_name == '!HEAT' ) then
        c_heat = c_heat + 1
      else if( header_name == '!FIXTEMP' ) then
        c_fixtemp = c_fixtemp + 1
        call fstr_setup_fixtemp( ctrl, c_fixtemp, p )
      else if( header_name == '!CFLUX' ) then
        c_cflux = c_cflux + 1
        call fstr_setup_cflux( ctrl, c_cflux, p )
      else if( header_name == '!DFLUX' ) then
        c_dflux = c_dflux + 1
        call fstr_setup_dflux( ctrl, c_dflux, p )
      else if( header_name == '!SFLUX' ) then
        c_sflux = c_sflux + 1
        call fstr_setup_sflux( ctrl, c_sflux, p )
      else if( header_name == '!FILM' ) then
        c_film = c_film + 1
        call fstr_setup_film( ctrl, c_film, p )
      else if( header_name == '!SFILM' ) then
        c_sfilm = c_sfilm + 1
        call fstr_setup_sfilm( ctrl, c_sfilm, p )
      else if( header_name == '!RADIATE' ) then
        c_radiate = c_radiate + 1
        call fstr_setup_radiate( ctrl, c_radiate, p )
      else if( header_name == '!SRADIATE' ) then
        c_sradiate = c_sradiate + 1
        call fstr_setup_sradiate( ctrl, c_sradiate, p )
      else if( header_name == '!WELD_LINE' ) then
        c_weldline = c_weldline + 1
 
        !--------------- for eigen -------------------------
 
      else if( header_name == '!EIGEN' ) then
        c_eigen = c_eigen + 1
        call fstr_setup_eigen( ctrl, c_eigen, p )
 
        !--------------- for dynamic -------------------------
 
      else if( header_name == '!DYNAMIC' ) then
        c_dynamic = c_dynamic + 1
        call fstr_setup_dynamic( ctrl, c_eigen, p )
      else if( header_name == '!VELOCITY' ) then
        c_velocity = c_velocity + 1
        call fstr_setup_velocity( ctrl, c_eigen, p )
      else if( header_name == '!ACCELERATION' ) then
        c_acceleration = c_acceleration + 1
        call fstr_setup_acceleration( ctrl, c_eigen, p )
      else if( header_name == '!FLOAD' ) then
        c_fload = c_fload + 1
        call fstr_setup_fload( ctrl , c_fload, p )
      else if( header_name == '!FACCELERATION' ) then
        ! acceleration input of the frequency response analysis
        c_facc = c_facc + 1
        call fstr_setup_facceleration( ctrl, c_facc, p )
      else if( header_name == '!EIGENREAD' ) then
        c_eigenread = c_eigenread + 1
        call fstr_setup_eigenread( ctrl, c_eigenread, p )
 
        !--------------- for couple -------------------------
 
      else if( header_name == '!COUPLE' ) then
        c_couple = c_couple + 1
        call fstr_setup_couple( ctrl, c_couple, p )
 
        !--------------- for mpc -------------------------
 
      else if( header_name == '!MPC' ) then
        c_mpc = c_mpc + 1
        call fstr_setup_mpc( ctrl, c_mpc, p )
 
        !--------------------- for input -------------------------
 
      else if( header_name == '!INCLUDE' ) then
        ctrl_list(ictrl) = ctrl
        input_filename   = ""
        ierror = fstr_ctrl_get_param_ex( ctrl, 'INPUT ', '# ', 0, 'F', input_filename )
        ctrl   = fstr_ctrl_open( input_filename )
        if( ctrl < 0 ) then
          write(*,*) '### Error: Cannot open FSTR control file : ', input_filename
          write(ilog,*) '### Error: Cannot open FSTR control file : ', input_filename
          stop hecmw_exit_input
        end if
        ictrl = ictrl + 1
        cycle
 
        !--------------------- END -------------------------
 
      else if( header_name == '!END' ) then
        exit
      end if
 
      ! next
      if( fstr_ctrl_seek_next_header(ctrl) == 0 )then
        if( ictrl == 1 )then
          exit
        else
          ierror= fstr_ctrl_close( ctrl )
          ictrl = ictrl - 1
          ctrl  = ctrl_list(ictrl)
          if( fstr_ctrl_seek_next_header(ctrl) == 0 ) exit
        endif
      endif
    end do
 
    ! -----
    fstrsolid%n_contacts = c_contact
    if( c_contact>0 ) then
      allocate( fstrsolid%contacts( c_contact ) )
      ! convert SURF_SURF contact to NODE_SURF contact
      call fstr_convert_contact_type( p%MESH )
    endif
    call fstr_convert_rotation_dof( p )
    fstrsolid%n_embeds = c_embed
    if( c_embed>0 )  allocate( fstrsolid%embeds( c_embed ) )
    if( c_weldline>0 ) allocate( fstrheat%weldline( c_weldline ) )
    if( c_initial>0 ) allocate( g_initialcnd( c_initial ) )
    if( c_istep>0 ) then
      allocate( fstrsolid%step_ctrl( c_istep ) )
      do i=1, c_istep
        call init_stepinfo( fstrsolid%step_ctrl(i) )
        if( p%PARAM%solution_type==kstdynamic ) then
          fstrsolid%step_ctrl(i)%num_substep = fstrdynamic%n_step
          fstrsolid%step_ctrl(i)%initdt      = fstrdynamic%t_delta
          fstrsolid%step_ctrl(i)%elapsetime  = dble(fstrdynamic%n_step) * fstrdynamic%t_delta
          fstrsolid%step_ctrl(i)%mindt       = fstrdynamic%t_delta
          fstrsolid%step_ctrl(i)%maxdt       = fstrdynamic%t_delta
        endif
      end do
    endif
    if( c_localcoord>0 ) allocate( g_localcoordsys(c_localcoord) )
    allocate( fstrparam%ainc(0:c_aincparam) )
    do i=0,c_aincparam
      call init_aincparam( fstrparam%ainc(i) )
    end do
    allocate( fstrparam%cnvparam(0:c_cnvparam) )
    do i=0,c_cnvparam
      call init_convergparam( fstrparam%cnvparam(i) )
    end do
    if( c_timepoints>0 ) allocate( fstrparam%timepoints(c_timepoints) )
    allocate( fstrparam%contactparam(0:c_contactparam) )
    do i=0,c_contactparam
      call init_contactparam( fstrparam%contactparam(i) )
    end do
    if( c_contact_if>0 )then 
      allocate( fstrparam%contact_if( c_contact_if ) )
      do i=1,c_contact_if
        call init_contact_if( fstrparam%contact_if(i) )
      end do
    end if
 
    p%SOLID%is_33shell = 0
    p%SOLID%is_33beam  = 0
 
    do i=1,hecmesh%n_elem_type
      n =  hecmesh%elem_type_item(i)
      if (n == 781 .or. n == 761)then
        p%SOLID%is_33shell = 1
      elseif (n == 641)then
        p%SOLID%is_33beam  = 1
      endif
    enddo
 
    n = c_material
    if( hecmesh%material%n_mat>n ) n= hecmesh%material%n_mat
    if( n==0 ) stop "material property not defined!"
    allocate( fstrsolid%materials( n ) )
    do i = 1, n
      call initmaterial(fstrsolid%materials(i))
    enddo
 
    ! for section control
    allocate( fstrsolid%sections(hecmesh%section%n_sect) )
    do i=1,hecmesh%section%n_sect
      ! set default 361 element formulation
      if( p%PARAM%solution_type==kststatic .or. p%PARAM%solution_type==kstdynamic ) then
        if( p%PARAM%nlgeom ) then
          fstrsolid%sections(i)%elemopt361 = kel361fbar
        else
          fstrsolid%sections(i)%elemopt361 = kel361ic
        end if
      else if( p%PARAM%solution_type==ksteigen ) then
        fstrsolid%sections(i)%elemopt361 = kel361ic
      else if( p%PARAM%solution_type==kststaticeigen ) then
        fstrsolid%sections(i)%elemopt361 = kel361fbar
      else
        fstrsolid%sections(i)%elemopt361 = kel361fi
      end if
      fstrsolid%sections(i)%elemopt341 = kel341fi
      fstrsolid%sections(i)%elemopt611 = kel611euler
      ! INTERFACE sections are skipped in the loop below and keep this thickness
      fstrsolid%sections(i)%thickness = 1.0d0
    enddo
 
    if( hecmesh%section%n_sect >0 ) then
      do i=1,hecmesh%section%n_sect
        if( hecmesh%section%sect_type(i) == 4 ) cycle
        cid = hecmesh%section%sect_mat_ID_item(i)
        if( cid>n ) stop "Error in material property definition!"
        if( fstrparam%nlgeom .or. fstrparam%solution_type==kststaticeigen ) &
          fstrsolid%materials(cid)%nlgeom_flag = 1
        nullify(shmat)
        call fstr_get_prop(hecmesh,shmat,i,ee,pp,rho,alpha,thick,&
          n_totlyr,alpha_over_mu, &
          beam_radius,beam_angle1,beam_angle2,beam_angle3, &
          beam_angle4,beam_angle5,beam_angle6)
        fstrsolid%sections(i)%thickness = thick
        fstrsolid%materials(cid)%name = hecmesh%material%mat_name(cid)
        fstrsolid%materials(cid)%variables(m_youngs)=ee
        fstrsolid%materials(cid)%variables(m_poisson)=pp
        fstrsolid%materials(cid)%variables(m_density)=rho
        fstrsolid%materials(cid)%variables(m_expansion)=alpha
        fstrsolid%materials(cid)%variables(m_alpha_over_mu)= alpha_over_mu
        fstrsolid%materials(cid)%variables(m_beam_radius)=beam_radius
        fstrsolid%materials(cid)%variables(m_beam_angle1)=beam_angle1
        fstrsolid%materials(cid)%variables(m_beam_angle2)=beam_angle2
        fstrsolid%materials(cid)%variables(m_beam_angle3)=beam_angle3
        fstrsolid%materials(cid)%variables(m_beam_angle4)=beam_angle4
        fstrsolid%materials(cid)%variables(m_beam_angle5)=beam_angle5
        fstrsolid%materials(cid)%variables(m_beam_angle6)=beam_angle6
        fstrsolid%materials(cid)%mtype = elastic
        if( hecmesh%section%sect_type(i) == 2 ) then
          fstrsolid%materials(cid)%totallyr = n_totlyr
          fstrsolid%materials(cid)%shell_var => shmat
        endif
      enddo
    endif
 
    allocate( fstrsolid%output_ctrl( 4 ) )
    call fstr_init_outctrl(fstrsolid%output_ctrl(1))
    fstrsolid%output_ctrl( 1 )%filename = trim(logfilename)
    fstrsolid%output_ctrl( 1 )%filenum = ilog
    call fstr_init_outctrl(fstrsolid%output_ctrl(2))
    call fstr_init_outctrl(fstrsolid%output_ctrl(3))
    call fstr_init_outctrl(fstrsolid%output_ctrl(4))
 
    ! -----
    rcode = fstr_ctrl_rewind( ctrl )
 
    c_istep    = 0
    c_heat     = 0
    c_material = 0
    c_output = 0
    c_contact  = 0
    c_contactparam  = 0
    c_contact_if  = 0
    c_embed = 0
    c_initial = 0
    c_localcoord = 0
    c_section = 0
    fstrheat%WL_tot = 0
    c_elemopt = 0
    c_aincparam = 0
    c_timepoints = 0
    c_cnvparam = 0
    fstrsolid%elemopt361 = 0
    fstrsolid%AutoINC_stat = 0
    fstrsolid%CutBack_stat = 0
    fstrsolid%NRstat_i(:) = 0
    fstrsolid%NRstat_r(:) = 0.d0
    ictrl = 1
    do
      rcode = fstr_ctrl_get_c_h_name( ctrl, header_name, hecmw_name_len )
 
      if( header_name == '!ORIENTATION' ) then
        c_localcoord = c_localcoord + 1
        if( fstr_setup_orientation( ctrl, hecmesh, c_localcoord, g_localcoordsys(c_localcoord) )/=0 ) then
          write(*,*) '### Error: Fail in read in ORIENTATION definition : ', c_localcoord
          write(ilog,*) '### Error: Fail in read in ORIENTATION definition : ', c_localcoord
          stop hecmw_exit_input
        endif
 
        ! ----- CONTACT condition setting
      elseif( header_name == '!CONTACT' ) then
        n = fstr_ctrl_get_data_line_n( ctrl )
        if( .not. fstr_ctrl_get_contact( ctrl, n, fstrsolid%contacts(c_contact+1:c_contact+n)   &
            ,ee, pp, rho, alpha, p%PARAM%contact_algo, mname, k ) ) then
          write(*,*) '### Error: Fail in read in contact condition : ', c_contact
          write(ilog,*) '### Error: Fail in read in contact condition : ', c_contact
          stop hecmw_exit_input
        endif
        cparam_id = 0
        do i=1,size(fstrparam%contactparam)-1
          if( hecmw_streqr( fstrparam%contactparam(i)%name, mname ) ) then
            cparam_id = i; exit
          endif
        enddo
        ! initialize contact condition (cdotp and mut are now obsolete - penalty is in tContact)
        if( rho>0.d0 ) cgn = rho
        if( alpha>0.d0 ) cgt = alpha
        do i=1,n
          fstrsolid%contacts(c_contact+i)%smoothing = k
          if( .not. fstr_contact_check( fstrsolid%contacts(c_contact+i), p%MESH ) ) then
            write(*,*) '### Error: Inconsistence in contact and surface definition : ' , i+c_contact
            write(ilog,*) '### Error: Inconsistence in contact and surface definition : ', i+c_contact
            stop hecmw_exit_model
          else
            isok = fstr_contact_init( fstrsolid%contacts(c_contact+i), p%MESH, fstrparam%contactparam(cparam_id))
            ! !CONTACT_ALGO, FRICTION_CONE=FOLLOW gives up the symmetry of the ALagrange
            ! friction terms; see getContactStiffness_Alag
            fstrsolid%contacts(c_contact+i)%symmetric = .not. p%PARAM%fric_cone_follow
            fstrsolid%contacts(c_contact+i)%eps_fric_band = p%PARAM%eps_fric_band
          endif
          if( fstrsolid%contacts(c_contact+i)%method == contacts2s ) then
            ! Only the static path is wired for MORTAR=YES. The implicit dynamic solver calls
            ! the shared scan / matrix connectivity / stiffness routines, so the pair would be
            ! dispatched to the mortar code, but it never calls fstr_begin_lambda_txn /
            ! fstr_commit_lambda_txn and the multipliers would not advance over the substeps.
            ! The explicit solver applies the NODE-SURF scan to every pair regardless of the
            ! method. Both would be silent, so this must be an error.
            ! ELEMCHECK and NZPROF run no analysis at all (fstr_input_precheck returns before
            ! the solver dispatch and never touches the contact pairs), so they are let through.
            if( p%PARAM%solution_type /= kststatic .and. &
              & p%PARAM%solution_type /= kstprecheck .and. &
              & p%PARAM%solution_type /= kstnzprof ) then
              write(*,*)    '### Error: MORTAR=YES is supported only in static analysis : ', i+c_contact
              write(ilog,*) '### Error: MORTAR=YES is supported only in static analysis : ', i+c_contact
              write(*,*)    '           Specify !SOLUTION, TYPE=STATIC or TYPE=NLSTATIC.'
              write(ilog,*) '           Specify !SOLUTION, TYPE=STATIC or TYPE=NLSTATIC.'
              stop hecmw_exit_model
            endif
            ! Restart is not supported. The restart file keeps the node-level contact state
            ! only (contacts(i)%states over contacts(i)%slave), so the segment states in
            ! slave_surf and the converged multipliers are not saved. A resumed run would
            ! start with every segment free and without the warm start, and the initial scan
            ! is skipped on restart, so the difference would be silent.
            if( fstrsolid%restart_nout < 0 ) then
              write(*,*)    '### Error: MORTAR=YES does not support restart : ', i+c_contact
              write(ilog,*) '### Error: MORTAR=YES does not support restart : ', i+c_contact
              write(*,*)    '           The contact state of mortar pairs is not saved in the restart file.'
              write(ilog,*) '           The contact state of mortar pairs is not saved in the restart file.'
              stop hecmw_exit_model
            endif
            if( p%PARAM%contact_algo == kcaslagrange ) then
              write(*,*)    '### Error: MORTAR=YES is not supported with !CONTACT_ALGO TYPE=SLAGRANGE : ', i+c_contact
              write(ilog,*) '### Error: MORTAR=YES is not supported with !CONTACT_ALGO TYPE=SLAGRANGE : ', i+c_contact
              stop hecmw_exit_model
            endif
            if( fstrsolid%contacts(c_contact+i)%surf_id1_sgrp <= 0 ) then
              write(*,*)    '### Error: MORTAR=YES requires !CONTACT_PAIR TYPE=SURF-SURF : ', i+c_contact
              write(ilog,*) '### Error: MORTAR=YES requires !CONTACT_PAIR TYPE=SURF-SURF : ', i+c_contact
              stop hecmw_exit_model
            endif
            ! Only finite sliding is implemented. With any other INTERACTION the per-integration
            ! point tracking is skipped, and with it the activation of the slave segment, so the
            ! assembly sees every segment as free and no contact force is produced at all
            ! (a silent no-op). INTERACTION defaults to SSLID, so this must be an error.
            if( fstrsolid%contacts(c_contact+i)%algtype /= contactfslid ) then
              write(*,*)    '### Error: MORTAR=YES is supported only with INTERACTION=FSLID : ', i+c_contact
              write(ilog,*) '### Error: MORTAR=YES is supported only with INTERACTION=FSLID : ', i+c_contact
              write(*,*)    '           Specify INTERACTION=FSLID on !CONTACT (the default is SSLID).'
              write(ilog,*) '           Specify INTERACTION=FSLID on !CONTACT (the default is SSLID).'
              stop hecmw_exit_model
            endif
            if( fstrsolid%contacts(c_contact+i)%smoothing /= kcsnone ) then
              write(*,*)    '### Error: MORTAR=YES is not supported with SMOOTHING= : ', i+c_contact
              write(ilog,*) '### Error: MORTAR=YES is not supported with SMOOTHING= : ', i+c_contact
              stop hecmw_exit_model
            endif
            ! Visibility guard: the mortar integral of a slave segment needs every master face
            ! it can project onto, so a rank owning slave faces must see the entire master
            ! surface (owned + ghost), which !PARTITION, CONTACT_OWNER=SLAVE provides. Verified
            ! by comparing the visible master count with the global owned-only count.
            ! Serial (np=1) skips the block, so the numerical path is unchanged.
            if( hecmw_comm_get_size() > 1 ) then
              ss_nslav     = size( fstrsolid%contacts(c_contact+i)%slave_surf )
              ss_nmas_vis  = size( fstrsolid%contacts(c_contact+i)%master )
              ss_nmas_glob = fstrsolid%contacts(c_contact+i)%n_master_owned
              call hecmw_allreduce_i1( p%MESH, ss_nmas_glob, hecmw_sum )
              ! a slave-owning rank must see every master face; flag incomplete ghost
              if( ss_nslav > 0 .and. ss_nmas_vis /= ss_nmas_glob ) then
                ss_bad = 1
              else
                ss_bad = 0
              endif
              call hecmw_allreduce_i1( p%MESH, ss_bad, hecmw_max )
              if( ss_bad /= 0 ) then
                write(*,*)    '### Error: MORTAR=YES contact pair has incomplete master ghost : ', i+c_contact
                write(ilog,*) '### Error: MORTAR=YES contact pair has incomplete master ghost : ', i+c_contact
                write(*,*)    '           A rank owning slave faces cannot see the full master surface.'
                write(ilog,*) '           A rank owning slave faces cannot see the full master surface.'
                write(*,*)    '           Re-partition the mesh with !PARTITION, CONTACT_OWNER=SLAVE,'
                write(ilog,*) '           Re-partition the mesh with !PARTITION, CONTACT_OWNER=SLAVE,'
                write(*,*)    '           or run serial (np=1).'
                write(ilog,*) '           or run serial (np=1).'
                stop hecmw_exit_model
              endif
            endif
          endif
        enddo
        c_contact = c_contact+n
 
        ! ----- EMBED condition setting
      elseif( header_name == '!EMBED' ) then
        n = fstr_ctrl_get_data_line_n( ctrl )
        if( .not. fstr_ctrl_get_embed( ctrl, n, fstrsolid%embeds(c_embed+1:c_embed+n), mname, k ) ) then
          write(*,*) '### Error: Fail in read in embed condition : ', c_embed
          write(ilog,*) '### Error: Fail in read in embed condition : ', c_embed
          stop hecmw_exit_input
        endif
        cparam_id = 0
        do i=1,size(fstrparam%contactparam)-1
          if( hecmw_streqr( fstrparam%contactparam(i)%name, mname ) ) then
            cparam_id = i; exit
          endif
        enddo
        do i=1,n
          fstrsolid%embeds(c_embed+i)%smoothing = k
          if( .not. fstr_contact_check( fstrsolid%embeds(c_embed+i), p%MESH ) ) then
            write(*,*) '### Error: Inconsistence in contact and surface definition : ' , i+c_embed
            write(ilog,*) '### Error: Inconsistence in contact and surface definition : ', i+c_embed
            stop hecmw_exit_model
          else
            isok = fstr_embed_init( fstrsolid%embeds(c_embed+i), p%MESH, fstrparam%contactparam(cparam_id))
          endif
        enddo
        c_embed = c_embed+n
 
      else if( header_name == '!ISTEP'  ) then
        c_istep = c_istep+1
        if( .not. fstr_ctrl_get_istep( ctrl, hecmesh, fstrsolid%step_ctrl(c_istep), mname, mname2, mname3 ) ) then
          write(*,*) '### Error: Fail in read in step definition : ' , c_istep
          write(ilog,*) '### Error: Fail in read in step definition : ', c_istep
          stop hecmw_exit_input
        endif
        call fstr_setup_step_convergence( fstrsolid%step_ctrl(c_istep), mname3, fstrparam%cnvparam, c_istep )
        if( associated(fstrparam%timepoints) ) then
          do i=1,size(fstrparam%timepoints)
            if( hecmw_streqr( fstrparam%timepoints(i)%name, mname ) ) then
              fstrsolid%step_ctrl(c_istep)%timepoint_id = i; exit
            endif
          enddo
        endif
        if( associated(fstrparam%ainc) ) then
          do i=1,size(fstrparam%ainc)
            if( hecmw_streqr( fstrparam%ainc(i)%name, mname2 ) ) then
              fstrsolid%step_ctrl(c_istep)%AincParam_id = i; exit
            endif
          enddo
        endif
      else if( header_name == '!STEP' .and. version>=1 ) then
        c_istep = c_istep+1
        if( .not. fstr_ctrl_get_istep( ctrl, hecmesh, fstrsolid%step_ctrl(c_istep), mname, mname2, mname3 ) ) then
          write(*,*) '### Error: Fail in read in step definition : ' , c_istep
          write(ilog,*) '### Error: Fail in read in step definition : ', c_istep
          stop hecmw_exit_input
        endif
        call fstr_setup_step_convergence( fstrsolid%step_ctrl(c_istep), mname3, fstrparam%cnvparam, c_istep )
        ! For DYNAMIC fixed-increment: keep the !DYNAMIC time increment while preserving !STEP duration.
        ! fstr_ctrl_get_ISTEP unconditionally sets initdt=1/num_substep which is wrong for DYNAMIC.
        ! Only override initdt/mindt/maxdt; keep elapsetime and num_substep as-is.
        if( p%PARAM%solution_type==kstdynamic .and. &
          & fstrsolid%step_ctrl(c_istep)%inc_type == stepfixedinc ) then
          fstrsolid%step_ctrl(c_istep)%initdt      = fstrdynamic%t_delta
          fstrsolid%step_ctrl(c_istep)%mindt       = fstrdynamic%t_delta
          fstrsolid%step_ctrl(c_istep)%maxdt       = fstrdynamic%t_delta
        endif
        if( associated(fstrparam%timepoints) ) then
          do i=1,size(fstrparam%timepoints)
            if( hecmw_streqr( fstrparam%timepoints(i)%name, mname ) ) then
              fstrsolid%step_ctrl(c_istep)%timepoint_id = i; exit
            endif
          enddo
        endif
        if( associated(fstrparam%ainc) ) then
          do i=1,size(fstrparam%ainc)-1
            if( hecmw_streqr( fstrparam%ainc(i)%name, mname2 ) ) then
              fstrsolid%step_ctrl(c_istep)%AincParam_id = i; exit
            endif
          enddo
        endif
 
      else if( header_name == '!HEAT'  ) then
        c_heat = c_heat + 1
        call fstr_setup_heat( ctrl, c_heat, p )
 
      else if( header_name == '!WELD_LINE'  ) then
        fstrheat%WL_tot = fstrheat%WL_tot+1
        if( fstr_ctrl_get_weldline( ctrl, hecmesh, hecmw_name_len, fstrheat%weldline(fstrheat%WL_tot) )/=0 ) then
          write(*,*) '### Error: Fail in read in Weld Line definition : ' , fstrheat%WL_tot
          write(ilog,*) '### Error: Fail in read in Weld Line definition : ', fstrheat%WL_tot
          stop hecmw_exit_input
        endif
 
      else if( header_name == '!INITIAL_CONDITION' .or. header_name == '!INITIAL CONDITION' ) then
        c_initial = c_initial+1
        if( fstr_setup_initial( ctrl, g_initialcnd(c_initial), p%MESH )/=0 ) then
           write(*,*) '### Error: Fail in read in INITIAL CONDITION definition : ' ,c_initial
           write(ilog,*) '### Error: Fail in read in INITIAL CONDITION definition : ', c_initial
           stop hecmw_exit_input
        endif
 
      else if( header_name == '!SECTION'  ) then
        c_section = c_section+1
        if( fstr_ctrl_get_section( ctrl, hecmesh, fstrsolid%sections )/=0 ) then
          write(*,*) '### Error: Fail in read in SECTION definition : ' , c_section
          write(ilog,*) '### Error: Fail in read in SECTION definition : ', c_section
          stop hecmw_exit_input
        endif
 
      else if( header_name == '!ELEMOPT'  ) then
        c_elemopt = c_elemopt+1
        if( fstr_ctrl_get_elemopt( ctrl, fstrsolid%elemopt361 )/=0 ) then
          write(*,*) '### Error: Fail in read in ELEMOPT definition : ' , c_elemopt
          write(ilog,*) '### Error: Fail in read in ELEMOPT definition : ', c_elemopt
          stop hecmw_exit_input
        endif
 
        !== following material properties ==
      else if( header_name == '!MATERIAL' ) then
        c_material = c_material+1
        if( fstr_ctrl_get_material( ctrl, mname )/=0 ) then
          write(*,*) '### Error: Fail in read in material definition : ' , c_material
          write(ilog,*) '### Error: Fail in read in material definition : ', c_material
          stop hecmw_exit_input
        endif
        cid = 0
        if(cache < hecmesh%material%n_mat) then
          if(hecmw_streqr( hecmesh%material%mat_name(cache), mname ))then
            cid = cache
            cache = cache + 1
          endif
        endif
        if(cid == 0)then
          do i=1,hecmesh%material%n_mat
            if( hecmw_streqr( hecmesh%material%mat_name(i), mname ) ) then
              cid = i
              cache = i + 1
              exit
            endif
          enddo
        endif
        if(cid == 0)then
          write(*,*) '### Error: Fail in read in material definition : ' , c_material
          write(ilog,*) '### Error: Fail in read in material definition : ', c_material
          stop hecmw_exit_input
        endif
        fstrsolid%materials(cid)%name = hecmesh%material%mat_name(cid)
        if(c_material>hecmesh%material%n_mat) call initmaterial( fstrsolid%materials(cid) )
 
      else if( header_name == '!ELASTIC' ) then
        if( c_material >0 ) then
          if( fstr_ctrl_get_elasticity( ctrl,                                        &
              fstrsolid%materials(cid)%mtype,       &
              fstrsolid%materials(cid)%nlgeom_flag, &
              fstrsolid%materials(cid)%variables,   &
              fstrsolid%materials(cid)%dict)/=0 ) then
            write(*,*) '### Error: Fail in read in elasticity definition : ' , cid
            write(ilog,*) '### Error: Fail in read in elasticity definition : ', cid
            stop hecmw_exit_input
          endif
        endif
      else if( header_name == '!PLASTIC' ) then
        if( cid >0 ) then
          if( fstr_ctrl_get_plasticity( ctrl,                                        &
              fstrsolid%materials(cid)%mtype,       &
              fstrsolid%materials(cid)%nlgeom_flag, &
              fstrsolid%materials(cid)%variables,   &
              fstrsolid%materials(cid)%table,   &
              fstrsolid%materials(cid)%dict)/=0 ) then
            write(*,*) '### Error: Fail in read in plasticity definition : ' , cid
            write(ilog,*) '### Error: Fail in read in plasticity definition : ', cid
            stop hecmw_exit_input
          endif
        endif
      else if( header_name == '!HYPERELASTIC' ) then
        if( cid >0 ) then
          if( fstr_ctrl_get_hyperelastic( ctrl,                                        &
              fstrsolid%materials(cid)%mtype,       &
              fstrsolid%materials(cid)%nlgeom_flag, &
              fstrsolid%materials(cid)%variables )/=0 ) then
            write(*,*) '### Error: Fail in read in elasticity definition : ' , cid
            write(ilog,*) '### Error: Fail in read in elasticity definition : ', cid
            stop hecmw_exit_input
          endif
        endif
      else if( header_name == '!VISCOELASTIC' ) then
        if( cid >0 ) then
          if( fstr_ctrl_get_viscoelasticity( ctrl,                                        &
              fstrsolid%materials(cid)%mtype,       &
              fstrsolid%materials(cid)%nlgeom_flag, &
              fstrsolid%materials(cid)%dict)/=0 ) then
            write(*,*) '### Error: Fail in read in plasticity definition : ' , cid
            write(ilog,*) '### Error: Fail in read in plasticity definition : ', cid
            stop hecmw_exit_input
          endif
        endif
      else if( header_name == '!TRS' ) then
        if( cid >0 ) then
          if( fstrsolid%materials(cid)%mtype/=viscoelastic ) then
            write(*,*) '### WARNING: TRS can only be defined for viscoelastic material! It is ignored! '
            write(ilog,*) '### WARNING: TRS can only be defined for viscoelastic material! It is ignored! '
          else
            if( fstr_ctrl_get_trs( ctrl, fstrsolid%materials(cid)%mtype, fstrsolid%materials(cid)%variables)/=0 ) then
              write(*,*) '### Error: Fail in read in TRS definition : ' , cid
              write(ilog,*) '### Error: Fail in read in TRS definition : ', cid
              stop hecmw_exit_input
            endif
          endif
        endif
      else if( header_name == '!CREEP' ) then
        if( cid >0 ) then
          if( fstr_ctrl_get_viscoplasticity( ctrl,                                        &
              fstrsolid%materials(cid)%mtype,       &
              fstrsolid%materials(cid)%nlgeom_flag, &
              fstrsolid%materials(cid)%dict)/=0 ) then
            write(*,*) '### Error: Fail in read in plasticity definition : ' , cid
            write(ilog,*) '### Error: Fail in read in plasticity definition : ', cid
            stop hecmw_exit_input
          endif
        endif
      else if( header_name == '!DENSITY' ) then
        if( cid >0 ) then
          if( fstr_ctrl_get_density( ctrl, fstrsolid%materials(cid)%variables )/=0 ) then
            write(*,*) '### Error: Fail in read in density definition : ' , cid
            write(ilog,*) '### Error: Fail in read in density definition : ', cid
            stop hecmw_exit_input
          endif
        endif
      else if( header_name == '!EXPANSION_COEF' .or. header_name == '!EXPANSION_COEFF' .or. &
          header_name == '!EXPANSION') then
        if( cid >0 ) then
          if( fstr_ctrl_get_expansion_coeff( ctrl, fstrsolid%materials(cid)%variables, &
              fstrsolid%materials(cid)%dict)/=0 )  then
            write(*,*) '### Error: Fail in read in expansion coefficient definition : ' , cid
            write(ilog,*) '### Error: Fail in read in expansion coefficient definition : ', cid
            stop hecmw_exit_input
          endif
        endif
      else if( header_name == '!DAMPING') then
        if( cid >0 ) then
          if( fstr_ctrl_get_rayleigh_damping( ctrl, fstrsolid%materials(cid)%variables, &
              fstrsolid%materials(cid)%is_elem_Rayleigh_damping)/=0 )  then
            write(*,*) '### Error: Fail in read in damping definition : ' , cid
            write(ilog,*) '### Error: Fail in read in damping definition : ', cid
            stop hecmw_exit_input
          endif
        endif
      else if( header_name == '!FLUID' ) then
        if( c_material >0 ) then
          if( fstr_ctrl_get_fluid( ctrl,                                 &
              fstrsolid%materials(cid)%mtype,       &
              fstrsolid%materials(cid)%nlgeom_flag, &
              fstrsolid%materials(cid)%variables,   &
              fstrsolid%materials(cid)%dict)/=0 ) then
            write(*,*) '### Error: Fail in read in fluid definition : ' , cid
            write(ilog,*) '### Error: Fail in read in fluid definition : ', cid
            stop hecmw_exit_input
          endif
        endif
      else if( header_name == '!SPRING_D' ) then
        if( c_material >0 ) then
          if( fstr_ctrl_get_spring_d( ctrl,                                        &
              fstrsolid%materials(cid)%mtype,       &
              fstrsolid%materials(cid)%nlgeom_flag, &
              fstrsolid%materials(cid)%variables_i,  &
              fstrsolid%materials(cid)%dict)/=0 ) then
            write(*,*) '### Error: Fail in read in spring_d definition : ' , cid
            write(ilog,*) '### Error: Fail in read in spring_d definition : ', cid
            stop hecmw_exit_input
          endif
        endif
      else if( header_name == '!SPRING_A' ) then
        if( c_material >0 ) then
          if( fstr_ctrl_get_spring_a( ctrl,                                        &
              fstrsolid%materials(cid)%mtype,       &
              fstrsolid%materials(cid)%nlgeom_flag, &
              fstrsolid%materials(cid)%variables_i,  &
              fstrsolid%materials(cid)%dict)/=0 ) then
            write(*,*) '### Error: Fail in read in spring_a definition : ' , cid
            write(ilog,*) '### Error: Fail in read in spring_a definition : ', cid
            stop hecmw_exit_input
          endif
        endif
      else if( header_name == '!DASHPOT_D' ) then
        if( c_material >0 ) then
          if( fstr_ctrl_get_dashpot_d( ctrl,                                        &
              fstrsolid%materials(cid)%mtype,       &
              fstrsolid%materials(cid)%nlgeom_flag, &
              fstrsolid%materials(cid)%variables_i,  &
              fstrsolid%materials(cid)%dict)/=0 ) then
            write(*,*) '### Error: Fail in read in spring_d definition : ' , cid
            write(ilog,*) '### Error: Fail in read in spring_d definition : ', cid
            stop hecmw_exit_input
          endif
        endif
      else if( header_name == '!DASHPOT_A' ) then
        if( c_material >0 ) then
          if( fstr_ctrl_get_dashpot_a( ctrl,                                        &
              fstrsolid%materials(cid)%mtype,       &
              fstrsolid%materials(cid)%nlgeom_flag, &
              fstrsolid%materials(cid)%variables_i,  &
              fstrsolid%materials(cid)%dict)/=0 ) then
            write(*,*) '### Error: Fail in read in spring_a definition : ' , cid
            write(ilog,*) '### Error: Fail in read in spring_a definition : ', cid
            stop hecmw_exit_input
          endif
        endif
      else if( header_name == '!USER_MATERIAL' ) then
        if( cid >0 ) then
          if( fstr_ctrl_get_usermaterial( ctrl, fstrsolid%materials(cid)%mtype,   &
              fstrsolid%materials(cid)%nlgeom_flag, fstrsolid%materials(cid)%nfstatus, &
              fstrsolid%materials(cid)%variables(101:) )/=0 ) then
            write(*,*) '### Error: Fail in read in user defined material : ' , cid
            write(ilog,*) '### Error: Fail in read in user defined material : ', cid
            stop hecmw_exit_input
          endif
        endif
 
 
        ! == Following output control ==
      else if( header_name == '!WRITE' ) then
        call fstr_ctrl_get_output( ctrl, outctrl, islog, result, visual, femap )
        if( islog == 1 ) then
          c_output=1
          outctrl%filename = trim(logfilename)
          outctrl%filenum = ilog
          call fstr_copy_outctrl(fstrsolid%output_ctrl(c_output), outctrl)
        endif
        if( femap == 1 ) then
          c_output=2
          write( outctrl%filename, '(a,i0,a)') 'utable.',myrank,'.dat'
          outctrl%filenum = iutb
          call fstr_copy_outctrl(fstrsolid%output_ctrl(c_output), outctrl)
          open( unit=outctrl%filenum, file=outctrl%filename, status='REPLACE', iostat=ierror )
          if( ierror /= 0 ) then
            write(*,*) 'Warning: cannot open output file: ', trim(outctrl%filename)
          endif
          fstrsolid%output_ctrl(c_output)%outinfo%grp_id = 1
        endif
        if( result == 1 ) then
          c_output=3
          call fstr_copy_outctrl(fstrsolid%output_ctrl(c_output), outctrl)
        endif
        if( visual == 1 ) then
          c_output=4
          call fstr_copy_outctrl(fstrsolid%output_ctrl(c_output), outctrl)
        endif
 
      else if( header_name == '!OUTPUT_RES' ) then
        c_output=3
        if( .not. fstr_ctrl_get_outitem( ctrl, hecmesh, fstrsolid%output_ctrl(c_output)%outinfo ) ) then
          write(*,*) '### Error: Fail in read in node output definition : ' , c_output
          write(ilog,*) '### Error: Fail in read in node output definition : ', c_output
          stop hecmw_exit_input
        endif
        if( fstrsolid%output_ctrl(c_output)%outinfo%grp_id_name /= 'ALL' ) then
          c_output=2
          do i=1,hecmesh%node_group%n_grp
            if( fstrsolid%output_ctrl(c_output)%outinfo%grp_id_name == hecmesh%node_group%grp_name(i) ) then
              fstrsolid%output_ctrl(c_output)%outinfo%grp_id = i; exit
            endif
          enddo
        endif
      else if( header_name == '!OUTPUT_VIS' ) then
        c_output=4
        if( .not. fstr_ctrl_get_outitem( ctrl, hecmesh, fstrsolid%output_ctrl(c_output)%outinfo ) ) then
          write(*,*) '### Error: Fail in read in element output definition : ' , c_output
          write(ilog,*) '### Error: Fail in read in element output definition : ', c_output
          stop hecmw_exit_input
        endif
        if( fstrsolid%output_ctrl(c_output)%outinfo%grp_id_name /= 'ALL' ) then
          c_output=2
          do i=1,hecmesh%node_group%n_grp
            if( fstrsolid%output_ctrl(c_output)%outinfo%grp_id_name == hecmesh%node_group%grp_name(i) ) then
              fstrsolid%output_ctrl(c_output)%outinfo%grp_id = i; exit
            endif
          enddo
        endif
      else if( header_name == '!AUTOINC_PARAM' ) then
        c_aincparam = c_aincparam + 1
        if( fstr_get_autoinc( ctrl, fstrparam%ainc(c_aincparam) ) /=0 ) then
          write(*,*) '### Error: Fail in read in AUTOINC_PARAM definition : ' , c_aincparam
          write(ilog,*) '### Error: Fail in read in AUTOINC_PARAM definition : ', c_aincparam
          stop hecmw_exit_input
        endif
      else if( header_name == '!CONVERG_PARAM' ) then
        c_cnvparam = c_cnvparam + 1
        if( fstr_ctrl_get_convergparam( ctrl, fstrparam%cnvparam(c_cnvparam) ) /=0 ) then
          write(*,*) '### Error: Fail in read in CONVERG_PARAM definition : ' , c_cnvparam
          write(ilog,*) '### Error: Fail in read in CONVERG_PARAM definition : ', c_cnvparam
          stop hecmw_exit_input
        endif
      else if( header_name == '!TIME_POINTS'  ) then
        c_timepoints = c_timepoints + 1
        if( fstr_ctrl_get_timepoints( ctrl, fstrparam%timepoints(c_timepoints) )/=0 ) then
          write(*,*) '### Error: Fail in read in TIME_POINTS definition : ' , c_timepoints
          write(ilog,*) '### Error: Fail in read in TIME_POINTS definition : ', c_timepoints
          stop hecmw_exit_input
        endif
      else if( header_name == '!CONTACT_PARAM' ) then
        c_contactparam = c_contactparam + 1
        if( fstr_ctrl_get_contactparam( ctrl, fstrparam%contactparam(c_contactparam) ) /=0 ) then
          write(*,*) '### Error: Fail in read in CONTACT_PARAM definition : ' , c_contactparam
          write(ilog,*) '### Error: Fail in read in CONTACT_PARAM definition : ', c_contactparam
          stop hecmw_exit_input
        endif
      else if( header_name == '!CONTACT_INTERFERENCE' ) then
        n = fstr_ctrl_get_data_line_n( ctrl )
        if( fstr_ctrl_get_contact_if( ctrl, n, fstrparam%contact_if(c_contact_if+1:c_contact_if+n) ) /= 0 ) then
          write(*,*) '### Error: Fail in read in CONTACT_INTERFERENCE definition : ' , c_contact_if
          write(ilog,*) '### Error: Fail in read in CONTACT_INTERFERENCE definition : ', c_contact_if
          stop hecmw_exit_input
        endif
        do i=1, n
          if( check_apply_contact_if(fstrparam%contact_if(c_contact_if+i), fstrsolid%contacts) /= 0) then
            write(*,*) '### Error:(INTERFERENCE) Inconsistence of contact_pair in CONTACTS: ' , i+c_contact_if
            write(ilog,*) '### Error:(INTERFERENCE)  Inconsistence of contact_pair in CONTACTS: ', i+c_contact_if
            stop hecmw_exit_model
          end if
        end do
        c_contact_if = c_contact_if + n
      else if( header_name == '!ULOAD' ) then
        if( fstr_ctrl_get_userload( ctrl )/=0 ) then
          write(*,*) '### Error: Fail in read in ULOAD definition : '
          write(ilog,*) '### Error: Fail in read in ULOAD definition : '
          stop hecmw_exit_input
        endif
 
      else if( header_name == '!INCLUDE' ) then
        ctrl_list(ictrl) = ctrl
        input_filename   = ""
        ierror = fstr_ctrl_get_param_ex( ctrl, 'INPUT ', '# ', 0, 'F', input_filename )
        ctrl   = fstr_ctrl_open( input_filename )
        if( ctrl < 0 ) then
          write(*,*) '### Error: Cannot open FSTR control file : ', input_filename
          write(ilog,*) '### Error: Cannot open FSTR control file : ', input_filename
          stop hecmw_exit_input
        end if
        ictrl = ictrl + 1
        cycle
 
      else if( header_name == '!END' ) then
        exit
      endif
 
      ! next
      if( fstr_ctrl_seek_next_header(ctrl) == 0 )then
        if( ictrl == 1 )then
          exit
        else
          ierror= fstr_ctrl_close( ctrl )
          ictrl = ictrl - 1
          ctrl  = ctrl_list(ictrl)
          if( fstr_ctrl_seek_next_header(ctrl) == 0 ) exit
        endif
      endif
 
    end do
 
    ! ----- material type judgement. in case of infinitive analysis, nlgeom_flag=0
    if( .not. p%PARAM%nlgeom ) then
      do i=1, c_material
        fstrsolid%materials(i)%nlgeom_flag = 0
      enddo
    endif
 
    if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 ) then
      allocate ( fstrsolid%temperature( hecmesh%n_node )      ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <FSTR_SOLID, TEMPERATURE>'
        write(idbg,*) '  rank = ', myrank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
      fstrsolid%temperature = ref_temp
      allocate ( fstrsolid%last_temp( hecmesh%n_node )      ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <FSTR_SOLID, LAST_TEMP>'
        write(idbg,*) '  rank = ', myrank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
      fstrsolid%last_temp = 0.d0
      allocate ( fstrsolid%temp_bak( hecmesh%n_node )      ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <FSTR_SOLID, TEMP_BAK>'
        write(idbg,*) '  rank = ', myrank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
      fstrsolid%temp_bak = 0.d0
    endif
 
    if( associated(fstrsolid%step_ctrl) )  then
      fstrsolid%nstep_tot = size(fstrsolid%step_ctrl)
      call setup_stepinfo_starttime( fstrsolid%step_ctrl )
    else
      if( p%PARAM%solution_type==kststatic .and. p%PARAM%nlgeom ) then
        write( *,* ) " ERROR: STEP not defined!"
        write( idbg,* ) "ERROR: STEP not defined!"
        call flush(idbg)
        call fstr_abort( hecmw_exit_model )
      endif
 
      if( myrank==0 ) write(*,*)"Step control not defined! Using default step=1"
      fstrsolid%nstep_tot = 1
      allocate( fstrsolid%step_ctrl(1) )
      call init_stepinfo( fstrsolid%step_ctrl(1) )
      if( p%PARAM%solution_type==kstdynamic ) then
        fstrsolid%step_ctrl(1)%num_substep = fstrdynamic%n_step
        fstrsolid%step_ctrl(1)%initdt      = fstrdynamic%t_delta
        fstrsolid%step_ctrl(1)%elapsetime  = dble(fstrdynamic%n_step) * fstrdynamic%t_delta
        fstrsolid%step_ctrl(1)%mindt       = fstrdynamic%t_delta
        fstrsolid%step_ctrl(1)%maxdt       = fstrdynamic%t_delta
      endif
      n =  fstrsolid%BOUNDARY_ngrp_tot
      if( n>0 ) allocate( fstrsolid%step_ctrl(1)%Boundary(n) )
      do i = 1, n
        fstrsolid%step_ctrl(1)%Boundary(i) = fstrsolid%BOUNDARY_ngrp_GRPID(i)
      enddo
      n = fstrsolid%CLOAD_ngrp_tot + fstrsolid%DLOAD_ngrp_tot + fstrsolid%TEMP_ngrp_tot + fstrsolid%SPRING_ngrp_tot
      if( n>0 ) allocate( fstrsolid%step_ctrl(1)%Load(n) )
      n = 0
      do i = 1, fstrsolid%CLOAD_ngrp_tot
        n = n + 1
        fstrsolid%step_ctrl(1)%Load(n) = fstrsolid%CLOAD_ngrp_GRPID(i)
      enddo
      do i = 1, fstrsolid%DLOAD_ngrp_tot
        n = n + 1
        fstrsolid%step_ctrl(1)%Load(n) = fstrsolid%DLOAD_ngrp_GRPID(i)
      enddo
      do i = 1, fstrsolid%TEMP_ngrp_tot
        n = n + 1
        fstrsolid%step_ctrl(1)%Load(n) = fstrsolid%TEMP_ngrp_GRPID(i)
      enddo
      do i = 1, fstrsolid%SPRING_ngrp_tot
        n = n + 1
        fstrsolid%step_ctrl(1)%Load(n) = fstrsolid%SPRING_ngrp_GRPID(i)
      enddo
      n = fstrsolid%elemact%ELEMACT_egrp_tot
      if( n>0 ) allocate( fstrsolid%step_ctrl(1)%ElemActivation(n) )
      do i = 1, n
        fstrsolid%step_ctrl(1)%ElemActivation(i) = fstrsolid%elemact%ELEMACT_egrp_GRPID(i)
      enddo
    endif
 
    call fstr_element_init( hecmesh, fstrsolid, p%PARAM%solution_type )
    call fstr_check_shell_j2_inputs(hecmesh, fstrsolid, fstrparam)
    call fstr_check_up_inputs(hecmesh, fstrsolid)
    call fstr_check_shell_gauss_output(hecmesh, fstrsolid)
    if( p%PARAM%solution_type==kststatic .or. p%PARAM%solution_type==kstdynamic .or.   &
      p%PARAM%solution_type==ksteigen  .or. p%PARAM%solution_type==kststaticeigen )  &
      call fstr_solid_alloc( hecmesh, fstrsolid )
 
    if( p%PARAM%solution_type == kstheat) then
      p%PARAM%fg_irres = fstrsolid%output_ctrl(3)%frequency
      p%PARAM%fg_iwres = fstrsolid%output_ctrl(4)%frequency
      p%HEAT%elemact = p%SOLID%elemact
    endif
 
    n_totlyr = 1
    do i=1,hecmesh%section%n_sect
      cid = hecmesh%section%sect_mat_ID_item(i)
      n  = fstrsolid%materials(cid)%totallyr
      if (n > n_totlyr)then
        n_totlyr = n
      endif
    enddo
    p%SOLID%max_lyr = n_totlyr
 
    call fstr_setup_post( ctrl, p )
    rcode = fstr_ctrl_close( ctrl )
 
  end subroutine fstr_setup
 
  subroutine fstr_check_up_inputs(hecMESH, fstrSOLID)
    type(hecmwst_local_mesh), intent(in) :: hecMESH
    type(fstr_solid), intent(in) :: fstrSOLID
    integer(kind=kint) :: i, isect, mtype
 
    do i = 1, hecmesh%n_elem
      if( hecmesh%elem_type(i) /= 361 ) cycle
      isect = hecmesh%section_ID(i)
      if( fstrsolid%sections(isect)%elemopt361 /= kel361up ) cycle
      if( .not.associated(fstrsolid%elements(i)%gausses) ) cycle
      mtype = fstrsolid%elements(i)%gausses(1)%pMaterial%mtype
      if( .not.iselastoplastic(mtype) .and. .not.iscreep(mtype) ) cycle
      write(*,*) '###ERROR### : Elastoplastic and creep materials are not supported by FORM361=UP'
      write(ilog,*) '###ERROR### : Elastoplastic and creep materials are not supported by FORM361=UP'
      call hecmw_abort(hecmw_comm_get_comm())
    enddo
  end subroutine fstr_check_up_inputs
 
  subroutine fstr_check_shell_gauss_output(hecMESH, fstrSOLID)
    use mmechgauss, only: fstr_shell_output_point_count
    type(hecmwst_local_mesh), intent(in) :: hecMESH
    type(fstr_solid), intent(in) :: fstrSOLID
    integer(kind=kint) :: i
 
    if( hecmesh%n_dof /= 6 ) return
    if( .not.any(fstrsolid%output_ctrl(3)%outinfo%on(9:10)) ) return
    do i=1,hecmesh%n_elem
      if( fstr_shell_output_point_count(fstrsolid%elements(i))>0 ) cycle
      write(*,*) '###ERROR### : ISTRAIN/ISTRESS require shell thickness-point output support; element type = ', &
        hecmesh%elem_type(i)
      call hecmw_abort(hecmw_comm_get_comm())
    enddo
  end subroutine fstr_check_shell_gauss_output
 
  subroutine fstr_check_shell_j2_inputs(hecMESH, fstrSOLID, fstrPARAM)
    use m_elastoplastic, only: checkplanestressj2material, planestressj2errormessage
    use elementinfo, only: numofshellthicknessquadpoints
    use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
    type(hecmwst_local_mesh), intent(in) :: hecMESH
    type(fstr_solid), intent(in) :: fstrSOLID
    type(fstr_param), intent(in) :: fstrPARAM
    type(tmaterial), pointer :: material
    integer(kind=kint) :: i, nn, ierr, isect, index, ilayer
    real(kind=kreal) :: elastic(2)
    logical :: missing
    character(len=160) :: message
 
    do i = 1, hecmesh%n_elem
      if( hecmesh%elem_type(i) /= 741 ) cycle
      if( .not.associated(fstrsolid%elements(i)%gausses) ) cycle
      material => fstrsolid%elements(i)%gausses(1)%pMaterial
      if( .not.iselastoplastic(material%mtype) ) cycle
      nn = hecmesh%elem_node_index(i)-hecmesh%elem_node_index(i-1)
      message = ''
      call checkplanestressj2material(material, ierr)
      if( ierr /= 0 ) then
        message = planestressj2errormessage(ierr)
      else if( .not.fstr_uses_plane_stress_j2_shell(741, nn, material) ) then
        message = 'Elastoplastic MITC4 shells require INFINITESIMAL material kinematics'
      else if( fstrparam%solution_type /= kststatic .or. .not.fstrparam%nlgeom ) then
        message = 'Elastoplastic MITC4 shells require !SOLUTION, TYPE=STATIC, NONLINEAR'
      else if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 ) then
        message = 'Temperature loading is not supported in models with elastoplastic MITC4 shells'
      endif
      if( len_trim(message) == 0 .and. associated(material%shell_var) ) then
        ! The J2 return mapping uses the common elastic table, not layer-specific constants.
        call fetch_tabledata(mc_isoelastic, material%dict, elastic, missing)
        do ilayer = 1, material%totallyr
          if( .not.ieee_is_finite(material%shell_var(ilayer)%ee) .or. &
              .not.ieee_is_finite(material%shell_var(ilayer)%pp) ) then
            message = 'Invalid elastic properties in J2 shell layer'
            exit
          endif
          if( abs(material%shell_var(ilayer)%ee-elastic(1)) > 1.0d-12*elastic(1) .or. &
              abs(material%shell_var(ilayer)%pp-elastic(2)) > 1.0d-12 ) then
            message = 'J2 shell layers must share the !ELASTIC Young modulus and Poisson ratio'
            exit
          endif
        enddo
      endif
      if( len_trim(message) == 0 ) then
        isect = hecmesh%section_ID(i)
        index = hecmesh%section%sect_I_index(isect-1)
        if( hecmesh%section%sect_I_index(isect)-index /= 1 ) then
          message = 'Elastoplastic MITC4 shells require a SHELL section with 2 thickness integration points'
        else if( hecmesh%section%sect_I_item(index+1) /= numofshellthicknessquadpoints(741) ) then
          message = 'Elastoplastic MITC4 shells require 2 thickness integration points in !SECTION'
        endif
      endif
      if( len_trim(message) == 0 ) cycle
      write(*,*) '###ERROR### : ', trim(message)
      write(ilog,*) '###ERROR### : ', trim(message)
      call hecmw_abort(hecmw_comm_get_comm())
    enddo
  end subroutine fstr_check_shell_j2_inputs
 
 
  subroutine fstr_solid_init( hecMESH, fstrSOLID )
    use m_fstr
    type(hecmwst_local_mesh),target :: hecMESH
    type(fstr_solid)                :: fstrSOLID
 
    fstrsolid%file_type  = kbcffstr
 
    fstrsolid%BOUNDARY_ngrp_tot = 0
    fstrsolid%BOUNDARY_ngrp_rot = 0
    fstrsolid%CLOAD_ngrp_tot    = 0
    fstrsolid%CLOAD_ngrp_rot    = 0
    fstrsolid%DLOAD_ngrp_tot    = 0
    fstrsolid%DLOAD_follow      = 1
    fstrsolid%TEMP_ngrp_tot     = 0
    fstrsolid%SPRING_ngrp_tot   = 0
    fstrsolid%TEMP_irres        = 0
    fstrsolid%TEMP_tstep        = 1
    fstrsolid%TEMP_interval     = 1
    fstrsolid%TEMP_rtype        = 1
    fstrsolid%TEMP_factor       = 1.d0
    fstrsolid%VELOCITY_ngrp_tot = 0
    fstrsolid%VELOCITY_ngrp_rot = 0
    fstrsolid%ACCELERATION_ngrp_tot = 0
    fstrsolid%COUPLE_ngrp_tot   = 0
 
    fstrsolid%restart_nout= 0
    fstrsolid%is_smoothing_active = .false.
    fstrsolid%has_finite_rotation_kinematics = .false.
    fstrsolid%finite_rotation_state_ready = .false.
 
  end subroutine fstr_solid_init
 
  subroutine fstr_solid_alloc( hecMESH, fstrSOLID )
    use m_fstr
    use m_fstr_finiterotationkinematics, only: fstr_has_finite_rotation_kinematics, &
      fstr_mark_finite_rotation_nodes
    type(hecmwst_local_mesh),target :: hecMESH
    type(fstr_solid)                :: fstrSOLID
 
    integer :: ndof, ntotal, ierror
 
    ndof=hecmesh%n_dof
    ntotal=ndof*hecmesh%n_node
    fstrsolid%has_finite_rotation_kinematics = &
      fstr_has_finite_rotation_kinematics( hecmesh, fstrsolid )
    fstrsolid%finite_rotation_state_ready = .false.
 
    allocate ( fstrsolid%GL( ntotal )          ,stat=ierror )
    if( ierror /= 0 ) then
      write(idbg,*) 'stop due to allocation error <FSTR_SOLID, GL>'
      write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm())
    end if
    allocate ( fstrsolid%GL0( ntotal )          ,stat=ierror )
    if( ierror /= 0 ) then
      write(idbg,*) 'stop due to allocation error <FSTR_SOLID, GL0>'
      write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm())
    end if
    allocate ( fstrsolid%EFORCE( ntotal )      ,stat=ierror )
    if( ierror /= 0 ) then
      write(idbg,*) 'stop due to allocation error <FSTR_SOLID, EFORCE>'
      write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm())
    end if
    !        allocate ( fstrSOLID%TOTAL_DISP( ntotal )  ,STAT=ierror )
    !            if( ierror /= 0 ) then
    !              write(idbg,*) 'stop due to allocation error <FSTR_SOLID, TOTAL_DISP>'
    !              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
    !              call flush(idbg)
    !              call hecmw_abort( hecmw_comm_get_comm())
    !            end if
    allocate ( fstrsolid%unode( ntotal )  ,stat=ierror )
    if( ierror /= 0 ) then
      write(idbg,*) 'stop due to allocation error <FSTR_SOLID, unode>'
      write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm())
    end if
    allocate ( fstrsolid%unode_bak( ntotal )  ,stat=ierror )
    if( ierror /= 0 ) then
      write(idbg,*) 'stop due to allocation error <FSTR_SOLID, unode>'
      write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm())
    end if
    allocate ( fstrsolid%dunode( ntotal )  ,stat=ierror )
    if( ierror /= 0 ) then
      write(idbg,*) 'stop due to allocation error <FSTR_SOLID, dunode>'
      write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm())
    end if
    allocate ( fstrsolid%ddunode( ntotal )  ,stat=ierror )
    if( ierror /= 0 ) then
      write(idbg,*) 'stop due to allocation error <FSTR_SOLID, ddunode>'
      write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm())
    end if
    if( fstrsolid%has_finite_rotation_kinematics ) then
      allocate ( fstrsolid%shell_node_mode( hecmesh%n_node )  ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <FSTR_SOLID, shell_node_mode>'
        write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
      call fstr_mark_finite_rotation_nodes( hecmesh, fstrsolid, ndof, fstrsolid%shell_node_mode )
      allocate ( fstrsolid%shell_rot_state( hecmesh%n_node )  ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <FSTR_SOLID, shell_rot_state>'
        write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
      allocate ( fstrsolid%shell_ref_triad( 9*hecmesh%n_node )  ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <FSTR_SOLID, shell_ref_triad>'
        write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
      allocate ( fstrsolid%shell_triad( 9*hecmesh%n_node )  ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <FSTR_SOLID, shell_triad>'
        write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
      allocate ( fstrsolid%shell_triad_bak( 9*hecmesh%n_node )  ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <FSTR_SOLID, shell_triad_bak>'
        write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
      allocate ( fstrsolid%shell_dtriad( 9*hecmesh%n_node )  ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <FSTR_SOLID, shell_dtriad>'
        write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
      allocate ( fstrsolid%shell_drill( hecmesh%n_node )  ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <FSTR_SOLID, shell_drill>'
        write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
      allocate ( fstrsolid%shell_drill_bak( hecmesh%n_node )  ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <FSTR_SOLID, shell_drill_bak>'
        write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
      allocate ( fstrsolid%shell_ddrill( hecmesh%n_node )  ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <FSTR_SOLID, shell_ddrill>'
        write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    end if
    allocate ( fstrsolid%QFORCE( ntotal )      ,stat=ierror )
    if( ierror /= 0 ) then
      write(idbg,*) 'stop due to allocation error <FSTR_SOLID, QFORCE>'
      write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm())
    end if
    allocate ( fstrsolid%DFORCE( ntotal )      ,stat=ierror )
    if( ierror /= 0 ) then
      write(idbg,*) 'stop due to allocation error <FSTR_SOLID, DFORCE>'
      write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm())
    end if
    allocate ( fstrsolid%QFORCE_bak( ntotal )      ,stat=ierror )
    if( ierror /= 0 ) then
      write(idbg,*) 'stop due to allocation error <FSTR_SOLID, QFORCE_bak>'
      write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm())
    end if
 
    fstrsolid%GL(:)=0.d0
    fstrsolid%GL0(:)=0.d0
    !        fstrSOLID%TOTAL_DISP(:)=0.d0
    fstrsolid%unode(:)      = 0.d0
    fstrsolid%unode_bak(:)  = 0.d0
    fstrsolid%dunode(:)     = 0.d0
    fstrsolid%ddunode(:)    = 0.d0
    if( associated(fstrsolid%shell_rot_state) ) fstrsolid%shell_rot_state(:) = 0
    if( associated(fstrsolid%shell_ref_triad) ) fstrsolid%shell_ref_triad(:) = 0.d0
    if( associated(fstrsolid%shell_triad) )     fstrsolid%shell_triad(:)     = 0.d0
    if( associated(fstrsolid%shell_triad_bak) ) fstrsolid%shell_triad_bak(:) = 0.d0
    if( associated(fstrsolid%shell_dtriad) )    fstrsolid%shell_dtriad(:)    = 0.d0
    if( associated(fstrsolid%shell_drill) )     fstrsolid%shell_drill(:)     = 0.d0
    if( associated(fstrsolid%shell_drill_bak) ) fstrsolid%shell_drill_bak(:) = 0.d0
    if( associated(fstrsolid%shell_ddrill) )    fstrsolid%shell_ddrill(:)    = 0.d0
    fstrsolid%QFORCE(:)     = 0.d0
    fstrsolid%QFORCE_bak(:) = 0.d0
    ! DFORCE is assembled only by implicit dynamic analysis, but the convergence
    ! check reads it in every analysis type.
    fstrsolid%DFORCE(:)     = 0.d0
    fstrsolid%FACTOR( 1:2 ) = 0.d0
 
    ! for MPC
    fstrsolid%n_fix_mpc = hecmesh%mpc%n_mpc
    if( fstrsolid%n_fix_mpc>0 ) then
      allocate( fstrsolid%mpc_const( fstrsolid%n_fix_mpc ) )
      fstrsolid%mpc_const(:) = hecmesh%mpc%mpc_const(:)
    endif
 
    ! initialize for linear static problems
    fstrsolid%FACTOR(2)=1.d0
    fstrsolid%FACTOR(1)=0.d0
  end subroutine fstr_solid_alloc
 
  subroutine fstr_smoothed_element_init( hecMESH, fstrSOLID )
    type(hecmwst_local_mesh),target :: hecMESH
    type(fstr_solid)                :: fstrSOLID
 
    logical, allocatable :: is_selem_list(:)
    integer :: i, isect
 
    do isect=1,hecmesh%section%n_sect
      if( fstrsolid%sections(isect)%elemopt341 == kel341sesns ) fstrsolid%is_smoothing_active = .true.
    end do
    if( .not. fstrsolid%is_smoothing_active ) return
 
    allocate(is_selem_list(hecmesh%n_elem), stat=i)
    if( i /= 0 ) then
      write(*,*) 'Allocation error: is_selem_list'
      return
    endif
    is_selem_list(:) = .false.
 
    do i=1,hecmesh%n_elem
      isect= hecmesh%section_ID(i)
      if( hecmesh%elem_type(i) /= fe_tet4n ) cycle
      if( fstrsolid%sections(isect)%elemopt341 == kel341sesns ) is_selem_list(i) = .true.
    enddo
 
    call hecmw_create_smoothing_element_connectivity(hecmesh,is_selem_list)
 
    deallocate(is_selem_list)
 
  end subroutine
 
  subroutine fstr_smoothed_element_calcmaxcon( hecMESH, fstrSOLID )
    use m_static_lib_c3d4sesns
    type(hecmwst_local_mesh),target :: hecMESH
    type(fstr_solid)                :: fstrSOLID
 
    integer :: i, nodlocal(fstrSOLID%max_ncon), iiS, nn, con_stf
 
    if( fstrsolid%max_ncon_stf > 20 ) fstrsolid%max_ncon_stf = 20   
 
    do i=1,hecmesh%elem_type_index(hecmesh%n_elem_type)
      if( hecmesh%elem_type(i) /= fe_nodesmooth_tet4n .and. hecmesh%elem_type(i) /= fe_edgesmooth_tet4n ) cycle
      iis = hecmesh%elem_node_index(i-1)
      nn = hecmesh%elem_node_index(i) - iis
      nodlocal(1:nn) = hecmesh%elem_node_item(iis+1:iis+nn)
      con_stf = return_nn_comp_c3d4_sesns(nn, nodlocal)
      if( con_stf > fstrsolid%max_ncon_stf ) fstrsolid%max_ncon_stf = con_stf
    enddo
  end subroutine
 
  subroutine fstr_element_init( hecMESH, fstrSOLID, solution_type )
    use elementinfo
    use mmechgauss
    use m_fstr_finiterotationkinematics, only: fstr_uses_finite_rotation_kinematics
    use m_fstr
    type(hecmwst_local_mesh),target :: hecMESH
    type(fstr_solid)                :: fstrSOLID
    integer(kind=kint), intent(in)  :: solution_type
 
    integer :: i, j, ng, isect, ndof, id, nn, n_elem, nthick
 
    if( hecmesh%n_elem <=0 ) then
      stop "no element defined!"
    endif
 
    fstrsolid%maxn_gauss = 0
    fstrsolid%max_ncon   = 0
 
    ! elemopt341 = kel341ES
    call fstr_smoothed_element_init( hecmesh, fstrsolid )
 
    ! number of elements
    n_elem = hecmesh%elem_type_index(hecmesh%n_elem_type)
    allocate( fstrsolid%elements(n_elem) )
    
    do i=1,n_elem
      fstrsolid%elements(i)%elemact_flag = kelact_undefined
      if( solution_type == kstheat) cycle !fstrSOLID is used only for elemact element in heat transfer analysis
 
      fstrsolid%elements(i)%etype = hecmesh%elem_type(i)
      if( hecmesh%elem_type(i)==301 ) fstrsolid%elements(i)%etype=111
      if (hecmw_is_etype_link(fstrsolid%elements(i)%etype)) cycle
      if (hecmw_is_etype_patch(fstrsolid%elements(i)%etype)) cycle
      ng = numofquadpoints( fstrsolid%elements(i)%etype )
      if( ng > fstrsolid%maxn_gauss ) fstrsolid%maxn_gauss = ng
      if(ng>0) allocate( fstrsolid%elements(i)%gausses( ng ) )
 
      isect= hecmesh%section_ID(i)
      ndof = getspacedimension( fstrsolid%elements(i)%etype )
      if (ndof == 2) then              ! why do this???
        id=hecmesh%section%sect_opt(isect)
        if( id==0 ) then
          fstrsolid%elements(i)%iset=1
        else if( id==1) then
          fstrsolid%elements(i)%iset=0
        else if( id==2) then
          fstrsolid%elements(i)%iset=2
        endif
      endif
 
      if( isect<0 .or. isect>hecmesh%section%n_sect )   &
        stop "Error in element's section definition"
      id = hecmesh%section%sect_mat_ID_item(isect)
      fstrsolid%materials(id)%cdsys_ID = hecmesh%section%sect_orien_ID(isect)
      nn = hecmesh%elem_node_index(i)-hecmesh%elem_node_index(i-1)
      do j=1,ng
        fstrsolid%elements(i)%gausses(j)%pMaterial => fstrsolid%materials(id)
        call fstr_init_gauss( fstrsolid%elements(i)%gausses( j )  )
      enddo
      nthick = 0
      if( fstr_uses_finite_rotation_kinematics( fstrsolid%elements(i)%etype, nn, &
          fstrsolid%materials(id) ) .or. &
          fstr_uses_plane_stress_j2_shell(fstrsolid%elements(i)%etype, nn, &
          fstrsolid%materials(id)) ) then
        nthick = fstr_shell_num_thickness_points( fstrsolid%elements(i)%etype )
      endif
      if( nthick > 0 ) call fstr_init_shell_layer_gausses( fstrsolid%elements(i), ng, &
        fstrsolid%materials(id)%totallyr, nthick )
 
      allocate(fstrsolid%elements(i)%equiForces(nn*ndof))
      fstrsolid%elements(i)%equiForces = 0.0d0
      if( nn > fstrsolid%max_ncon ) fstrsolid%max_ncon = nn
 
      if( hecmesh%elem_type(i)==361 ) then
        if( fstrsolid%sections(isect)%elemopt361==kel361ic ) then
          allocate( fstrsolid%elements(i)%aux(3,3) )
          fstrsolid%elements(i)%aux = 0.0d0
        endif
        fstrsolid%elements(i)%p(1) = 0.0d0
      endif
 
    enddo
 
    fstrsolid%max_ncon_stf = fstrsolid%max_ncon
    if( fstrsolid%is_smoothing_active ) call fstr_smoothed_element_calcmaxcon( hecmesh, fstrsolid )
 
    call hecmw_allreduce_i1(hecmesh,fstrsolid%maxn_gauss,hecmw_max)
  end subroutine
 
  recursive subroutine fstr_solid_physic_val_finalize( phys )
    use m_fstr
    implicit none
    type(fstr_solid_physic_val) :: phys
    integer :: i
 
    if(associated(phys%STRESS)) deallocate(phys%STRESS)
    if(associated(phys%STRAIN)) deallocate(phys%STRAIN)
    if(associated(phys%MISES)) deallocate(phys%MISES)
 
    if(associated(phys%PSTRESS)) deallocate(phys%PSTRESS)
    if(associated(phys%PSTRAIN)) deallocate(phys%PSTRAIN)
    if(associated(phys%PSTRESS_VECT)) deallocate(phys%PSTRESS_VECT)
    if(associated(phys%PSTRAIN_VECT)) deallocate(phys%PSTRAIN_VECT)
 
    if(associated(phys%ESTRESS)) deallocate(phys%ESTRESS)
    if(associated(phys%ESTRAIN)) deallocate(phys%ESTRAIN)
    if(associated(phys%EMISES)) deallocate(phys%EMISES)
    if(associated(phys%EPLSTRAIN)) deallocate(phys%EPLSTRAIN)
 
    if(associated(phys%EPSTRESS)) deallocate(phys%EPSTRESS)
    if(associated(phys%EPSTRAIN)) deallocate(phys%EPSTRAIN)
    if(associated(phys%EPSTRESS_VECT)) deallocate(phys%EPSTRESS_VECT)
    if(associated(phys%EPSTRAIN_VECT)) deallocate(phys%EPSTRAIN_VECT)
    if(associated(phys%ENQM)) deallocate(phys%ENQM)
 
    if(associated(phys%LAYER)) then
      do i=1, size(phys%LAYER)
        call fstr_solid_physic_val_finalize(phys%LAYER(i))
      end do
      deallocate(phys%LAYER)
    end if
    if(associated(phys%PLUS)) then
      call fstr_solid_physic_val_finalize(phys%PLUS)
      deallocate(phys%PLUS)
    end if
    if(associated(phys%MINUS)) then
      call fstr_solid_physic_val_finalize(phys%MINUS)
      deallocate(phys%MINUS)
    end if
  end subroutine
 
  subroutine fstr_solid_finalize( fstrSOLID )
    use mmechgauss, only: fstr_finalize_shell_layer_gausses
    implicit none
    type(fstr_solid) :: fstrSOLID
    integer :: i, j, ierror
    if( associated(fstrsolid%materials) ) then
      do j=1,size(fstrsolid%materials)
        call finalizematerial(fstrsolid%materials(j))
      enddo
      deallocate( fstrsolid%materials )
    endif
    if( .not. associated(fstrsolid%elements ) ) return
    do i=1,size(fstrsolid%elements)
      if( associated(fstrsolid%elements(i)%gausses) ) then
        do j=1,size(fstrsolid%elements(i)%gausses)
          call fstr_finalize_gauss(fstrsolid%elements(i)%gausses(j))
        enddo
        deallocate( fstrsolid%elements(i)%gausses )
      endif
      call fstr_finalize_shell_layer_gausses( fstrsolid%elements(i) )
      if(associated(fstrsolid%elements(i)%equiForces) ) then
        deallocate(fstrsolid%elements(i)%equiForces)
      endif
      if( associated(fstrsolid%elements(i)%aux) ) then
        deallocate(fstrsolid%elements(i)%aux)
      endif
    enddo
 
    deallocate( fstrsolid%elements )
    if( associated( fstrsolid%mpc_const ) ) then
      deallocate( fstrsolid%mpc_const )
    endif
    call free_stepinfo( fstrsolid%step_ctrl_restart )
    if( associated(fstrsolid%step_ctrl) ) then
      do i=1,size(fstrsolid%step_ctrl)
        call free_stepinfo( fstrsolid%step_ctrl(i) )
      enddo
      deallocate( fstrsolid%step_ctrl )
    endif
    if(associated(fstrsolid%output_ctrl) ) then
      do i=1,size(fstrsolid%output_ctrl)
        if( fstrsolid%output_ctrl(i)%filenum==iutb )  &
          close(fstrsolid%output_ctrl(i)%filenum)
      enddo
      deallocate(fstrsolid%output_ctrl)
    endif
    if( associated( fstrsolid%sections ) ) then
      deallocate( fstrsolid%sections )
    endif
 
    if( associated(fstrsolid%GL) ) then
      deallocate(fstrsolid%GL               ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, GL>'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(fstrsolid%EFORCE) ) then
      deallocate(fstrsolid%EFORCE           ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, EFORCE>'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(fstrsolid%unode) ) then
      deallocate(fstrsolid%unode       ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, unode>'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(fstrsolid%unode_bak) ) then
      deallocate(fstrsolid%unode_bak   ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, unode_bak>'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(fstrsolid%dunode) ) then
      deallocate(fstrsolid%dunode       ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, dunode>'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(fstrsolid%ddunode) ) then
      deallocate(fstrsolid%ddunode       ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, ddunode>'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(fstrsolid%shell_node_mode) ) then
      deallocate(fstrsolid%shell_node_mode       ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, shell_node_mode>'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(fstrsolid%shell_rot_state) ) then
      deallocate(fstrsolid%shell_rot_state       ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, shell_rot_state>'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(fstrsolid%shell_ref_triad) ) then
      deallocate(fstrsolid%shell_ref_triad       ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, shell_ref_triad>'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(fstrsolid%shell_triad) ) then
      deallocate(fstrsolid%shell_triad       ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, shell_triad>'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(fstrsolid%shell_triad_bak) ) then
      deallocate(fstrsolid%shell_triad_bak       ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, shell_triad_bak>'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(fstrsolid%shell_dtriad) ) then
      deallocate(fstrsolid%shell_dtriad       ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, shell_dtriad>'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(fstrsolid%shell_drill) ) then
      deallocate(fstrsolid%shell_drill       ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, shell_drill>'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(fstrsolid%shell_drill_bak) ) then
      deallocate(fstrsolid%shell_drill_bak       ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, shell_drill_bak>'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(fstrsolid%shell_ddrill) ) then
      deallocate(fstrsolid%shell_ddrill       ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, shell_ddrill>'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(fstrsolid%QFORCE) ) then
      deallocate(fstrsolid%QFORCE           ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, QFORCE>'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(fstrsolid%DFORCE) ) then
      deallocate(fstrsolid%DFORCE           ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, DFORCE>'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(fstrsolid%temperature) ) then
      deallocate(fstrsolid%temperature       ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, temperature>'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(fstrsolid%last_temp) ) then
      deallocate(fstrsolid%last_temp       ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, reftemp>'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(fstrsolid%temp_bak) ) then
      deallocate(fstrsolid%temp_bak       ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, reftemp>'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
 
    ! Allocated in in f str_setup_BOUNDARY */
    if( associated(fstrsolid%BOUNDARY_ngrp_GRPID) ) then
       deallocate(fstrsolid%BOUNDARY_ngrp_GRPID, stat=ierror)
       if( ierror /= 0 ) then
          write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, BOUNDARY_ngrp_GRPID>'
          call flush(idbg)
          call hecmw_abort( hecmw_comm_get_comm())
       end if
    endif
    if( associated(fstrsolid%BOUNDARY_ngrp_ID) ) then
       deallocate(fstrsolid%BOUNDARY_ngrp_ID, stat=ierror)
       if( ierror /= 0 ) then
          write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, BOUNDARY_ngrp_ID>'
          call flush(idbg)
          call hecmw_abort( hecmw_comm_get_comm())
       end if
    endif
    if( associated(fstrsolid%BOUNDARY_ngrp_type) ) then
       deallocate(fstrsolid%BOUNDARY_ngrp_type, stat=ierror)
       if( ierror /= 0 ) then
          write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, BOUNDARY_ngrp_type>'
          call flush(idbg)
          call hecmw_abort( hecmw_comm_get_comm())
       end if
    endif
    if( associated(fstrsolid%BOUNDARY_ngrp_val) ) then
       deallocate(fstrsolid%BOUNDARY_ngrp_val, stat=ierror)
       if( ierror /= 0 ) then
          write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, BOUNDARY_ngrp_val>'
          call flush(idbg)
          call hecmw_abort( hecmw_comm_get_comm())
       end if
    endif
    if( associated(fstrsolid%BOUNDARY_ngrp_amp) ) then
       deallocate(fstrsolid%BOUNDARY_ngrp_amp, stat=ierror)
       if( ierror /= 0 ) then
          write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, BOUNDARY_ngrp_amp>'
          call flush(idbg)
          call hecmw_abort( hecmw_comm_get_comm())
       end if
    endif
    if( associated(fstrsolid%BOUNDARY_ngrp_istot) ) then
       deallocate(fstrsolid%BOUNDARY_ngrp_istot, stat=ierror)
       if( ierror /= 0 ) then
          write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, BOUNDARY_ngrp_istot>'
          call flush(idbg)
          call hecmw_abort( hecmw_comm_get_comm())
       end if
    endif
    if( associated(fstrsolid%BOUNDARY_ngrp_rotID) ) then
       deallocate(fstrsolid%BOUNDARY_ngrp_rotID, stat=ierror)
       if( ierror /= 0 ) then
          write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, BOUNDARY_ngrp_rotID>'
          call flush(idbg)
          call hecmw_abort( hecmw_comm_get_comm())
       end if
    endif
    if( associated(fstrsolid%BOUNDARY_ngrp_centerID) ) then
       deallocate(fstrsolid%BOUNDARY_ngrp_centerID, stat=ierror)
       if( ierror /= 0 ) then
          write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, BOUNDARY_ngrp_centerID>'
          call flush(idbg)
          call hecmw_abort( hecmw_comm_get_comm())
       end if
    endif
 
    if(associated(fstrsolid%VELOCITY_ngrp_GRPID)) deallocate(fstrsolid%VELOCITY_ngrp_GRPID)
    if(associated(fstrsolid%VELOCITY_ngrp_ID)) deallocate(fstrsolid%VELOCITY_ngrp_ID)
    if(associated(fstrsolid%VELOCITY_ngrp_type)) deallocate(fstrsolid%VELOCITY_ngrp_type)
    if(associated(fstrsolid%VELOCITY_ngrp_amp)) deallocate(fstrsolid%VELOCITY_ngrp_amp)
    if(associated(fstrsolid%VELOCITY_ngrp_val)) deallocate(fstrsolid%VELOCITY_ngrp_val)
    if(associated(fstrsolid%VELOCITY_ngrp_rotID)) deallocate(fstrsolid%VELOCITY_ngrp_rotID)
    if(associated(fstrsolid%VELOCITY_ngrp_centerID)) deallocate(fstrsolid%VELOCITY_ngrp_centerID)
 
    if(associated(fstrsolid%ACCELERATION_ngrp_GRPID)) deallocate(fstrsolid%ACCELERATION_ngrp_GRPID)
    if(associated(fstrsolid%ACCELERATION_ngrp_ID)) deallocate(fstrsolid%ACCELERATION_ngrp_ID)
    if(associated(fstrsolid%ACCELERATION_ngrp_type)) deallocate(fstrsolid%ACCELERATION_ngrp_type)
    if(associated(fstrsolid%ACCELERATION_ngrp_amp)) deallocate(fstrsolid%ACCELERATION_ngrp_amp)
    if(associated(fstrsolid%ACCELERATION_ngrp_val)) deallocate(fstrsolid%ACCELERATION_ngrp_val)
 
    ! Allocated in in fstr_setup_CLOAD
    if( associated(fstrsolid%CLOAD_ngrp_GRPID) ) then
       deallocate(fstrsolid%CLOAD_ngrp_GRPID, stat=ierror)
       if( ierror /= 0 ) then
          write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, CLOAD_ngrp_GRPID>'
          call flush(idbg)
          call hecmw_abort( hecmw_comm_get_comm())
       end if
    endif
    if( associated(fstrsolid%CLOAD_ngrp_ID) ) then
       deallocate(fstrsolid%CLOAD_ngrp_ID, stat=ierror)
       if( ierror /= 0 ) then
          write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, CLOAD_ngrp_ID>'
          call flush(idbg)
          call hecmw_abort( hecmw_comm_get_comm())
       end if
    endif
    if( associated(fstrsolid%CLOAD_ngrp_DOF) ) then
       deallocate(fstrsolid%CLOAD_ngrp_DOF, stat=ierror)
       if( ierror /= 0 ) then
          write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, CLOAD_ngrp_DOF>'
          call flush(idbg)
          call hecmw_abort( hecmw_comm_get_comm())
       end if
    endif
    if( associated(fstrsolid%CLOAD_ngrp_val) ) then
       deallocate(fstrsolid%CLOAD_ngrp_val, stat=ierror)
       if( ierror /= 0 ) then
          write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, CLOAD_ngrp_val>'
          call flush(idbg)
          call hecmw_abort( hecmw_comm_get_comm())
       end if
    endif
    if( associated(fstrsolid%CLOAD_ngrp_amp) ) then
       deallocate(fstrsolid%CLOAD_ngrp_amp, stat=ierror)
       if( ierror /= 0 ) then
          write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, CLOAD_ngrp_amp>'
          call flush(idbg)
          call hecmw_abort( hecmw_comm_get_comm())
       end if
    endif
    if( associated(fstrsolid%CLOAD_ngrp_rotID) ) then
       deallocate(fstrsolid%CLOAD_ngrp_rotID, stat=ierror)
       if( ierror /= 0 ) then
          write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, CLOAD_ngrp_rotID>'
          call flush(idbg)
          call hecmw_abort( hecmw_comm_get_comm())
       end if
    endif
    if( associated(fstrsolid%CLOAD_ngrp_centerID) ) then
       deallocate(fstrsolid%CLOAD_ngrp_centerID, stat=ierror)
       if( ierror /= 0 ) then
          write(idbg,*) 'stop due to deallocation error <FSTR_SOLID, CLOAD_ngrp_centerID>'
          call flush(idbg)
          call hecmw_abort( hecmw_comm_get_comm())
       end if
    endif
 
    if(associated(fstrsolid%DLOAD_ngrp_GRPID)) deallocate(fstrsolid%DLOAD_ngrp_GRPID)
    if(associated(fstrsolid%DLOAD_ngrp_ID)) deallocate(fstrsolid%DLOAD_ngrp_ID)
    if(associated(fstrsolid%DLOAD_ngrp_LID)) deallocate(fstrsolid%DLOAD_ngrp_LID)
    if(associated(fstrsolid%DLOAD_ngrp_amp)) deallocate(fstrsolid%DLOAD_ngrp_amp)
    if(associated(fstrsolid%DLOAD_ngrp_params)) deallocate(fstrsolid%DLOAD_ngrp_params)
 
    if(associated(fstrsolid%TEMP_ngrp_GRPID)) deallocate(fstrsolid%TEMP_ngrp_GRPID)
    if(associated(fstrsolid%TEMP_ngrp_ID)) deallocate(fstrsolid%TEMP_ngrp_ID)
    if(associated(fstrsolid%TEMP_ngrp_val)) deallocate(fstrsolid%TEMP_ngrp_val)
 
    if(associated(fstrsolid%SPRING_ngrp_GRPID)) deallocate(fstrsolid%SPRING_ngrp_GRPID)
    if(associated(fstrsolid%SPRING_ngrp_ID)) deallocate(fstrsolid%SPRING_ngrp_ID)
    if(associated(fstrsolid%SPRING_ngrp_DOF)) deallocate(fstrsolid%SPRING_ngrp_DOF)
    if(associated(fstrsolid%SPRING_ngrp_amp)) deallocate(fstrsolid%SPRING_ngrp_amp)
    if(associated(fstrsolid%SPRING_ngrp_val)) deallocate(fstrsolid%SPRING_ngrp_val)
    if(associated(fstrsolid%SPRING_incremental)) deallocate(fstrsolid%SPRING_incremental)
 
    if(associated(fstrsolid%COUPLE_ngrp_ID)) deallocate(fstrsolid%COUPLE_ngrp_ID)
 
!    if(associated(fstrSOLID%STRESS)) deallocate(fstrSOLID%STRESS)
!    if(associated(fstrSOLID%STRAIN)) deallocate(fstrSOLID%STRAIN)
!    if(associated(fstrSOLID%MISES)) deallocate(fstrSOLID%MISES)
!    if(associated(fstrSOLID%PSTRESS)) deallocate(fstrSOLID%PSTRESS)
!    if(associated(fstrSOLID%PSTRAIN)) deallocate(fstrSOLID%PSTRAIN)
!    if(associated(fstrSOLID%PSTRESS_VECT)) deallocate(fstrSOLID%PSTRESS_VECT)
!    if(associated(fstrSOLID%PSTRAIN_VECT)) deallocate(fstrSOLID%PSTRAIN_VECT)
 
!    if(associated(fstrSOLID%ESTRESS)) deallocate(fstrSOLID%ESTRESS)
!    if(associated(fstrSOLID%ESTRAIN)) deallocate(fstrSOLID%ESTRAIN)
!    if(associated(fstrSOLID%EMISES)) deallocate(fstrSOLID%EMISES)
!    if(associated(fstrSOLID%EPLSTRAIN)) deallocate(fstrSOLID%EPLSTRAIN)
 
!    if(associated(fstrSOLID%EPSTRESS)) deallocate(fstrSOLID%EPLSTRAIN)
!    if(associated(fstrSOLID%EPSTRAIN)) deallocate(fstrSOLID%EPLSTRAIN)
!    if(associated(fstrSOLID%EPSTRESS_VECT)) deallocate(fstrSOLID%EPLSTRAIN)
!    if(associated(fstrSOLID%EPSTRAIN_VECT)) deallocate(fstrSOLID%EPLSTRAIN)
 
    if(associated(fstrsolid%TNSTRAIN)) deallocate(fstrsolid%TNSTRAIN)
    if(associated(fstrsolid%TESTRAIN)) deallocate(fstrsolid%TESTRAIN)
 
    if(associated(fstrsolid%YIELD_RATIO)) deallocate(fstrsolid%YIELD_RATIO)
 
!    if(associated(fstrSOLID%ENQM)) deallocate(fstrSOLID%ENQM)
    if(associated(fstrsolid%REACTION)) deallocate(fstrsolid%REACTION)
 
    if(associated(fstrsolid%CONT_NFORCE)) deallocate(fstrsolid%CONT_NFORCE)
    if(associated(fstrsolid%CONT_FRIC)) deallocate(fstrsolid%CONT_FRIC)
    if(associated(fstrsolid%CONT_RELVEL)) deallocate(fstrsolid%CONT_RELVEL)
    if(associated(fstrsolid%CONT_STATE)) deallocate(fstrsolid%CONT_STATE)
    if(associated(fstrsolid%CONT_SGRP_ID)) deallocate(fstrsolid%CONT_SGRP_ID)
    if(associated(fstrsolid%CONT_AREA)) deallocate(fstrsolid%CONT_AREA)
    if(associated(fstrsolid%CONT_NTRAC)) deallocate(fstrsolid%CONT_NTRAC)
    if(associated(fstrsolid%CONT_FTRAC)) deallocate(fstrsolid%CONT_FTRAC)
    if(associated(fstrsolid%EMBED_NFORCE)) deallocate(fstrsolid%EMBED_NFORCE)
 
    if(associated(fstrsolid%SOLID)) call fstr_solid_physic_val_finalize(fstrsolid%SOLID)
    if(associated(fstrsolid%SHELL)) call fstr_solid_physic_val_finalize(fstrsolid%SHELL)
    if(associated(fstrsolid%BEAM)) call fstr_solid_physic_val_finalize(fstrsolid%BEAM)
 
    if(associated(fstrsolid%is_rot)) deallocate(fstrsolid%is_rot)
    if(associated(fstrsolid%GL0)) deallocate(fstrsolid%GL0)
    if(associated(fstrsolid%QFORCE_bak)) deallocate(fstrsolid%QFORCE_bak)
 
    if(associated(fstrsolid%contacts)) then
!      do i=1, size(fstrSOLID%contacts)
!        call fstr_contact_finalize(fstrSOLID%contacts(i))
!      end do
      deallocate(fstrsolid%contacts)
    end if
    if(associated(fstrsolid%embeds)) then
!      do i=1, size(fstrSOLID%embeds)
!        call fstr_contact_finalize(fstrSOLID%embeds(i))
!      end do
      deallocate(fstrsolid%embeds)
    end if
    if(associated(fstrsolid%unode_bkup)) deallocate(fstrsolid%unode_bkup)
    if(associated(fstrsolid%QFORCE_bkup)) deallocate(fstrsolid%QFORCE_bkup)
    if(associated(fstrsolid%DFORCE_bkup)) deallocate(fstrsolid%DFORCE_bkup)
    if(associated(fstrsolid%last_temp_bkup)) deallocate(fstrsolid%last_temp_bkup)
    if(associated(fstrsolid%shell_triad_bkup)) deallocate(fstrsolid%shell_triad_bkup)
    if(associated(fstrsolid%shell_drill_bkup)) deallocate(fstrsolid%shell_drill_bkup)
    if(associated(fstrsolid%elements_bkup)) then
      do i=1,size(fstrsolid%elements_bkup)
        if( associated(fstrsolid%elements_bkup(i)%gausses) ) then
          do j=1,size(fstrsolid%elements_bkup(i)%gausses)
            call fstr_finalize_gauss(fstrsolid%elements_bkup(i)%gausses(j))
          enddo
          deallocate( fstrsolid%elements_bkup(i)%gausses )
        endif
        call fstr_finalize_shell_layer_gausses( fstrsolid%elements_bkup(i) )
        if(associated(fstrsolid%elements_bkup(i)%equiForces) ) then
          deallocate(fstrsolid%elements_bkup(i)%equiForces)
        endif
        if( associated(fstrsolid%elements_bkup(i)%aux) ) then
          deallocate(fstrsolid%elements_bkup(i)%aux)
        endif
      enddo
      deallocate(fstrsolid%elements_bkup)
    end if
    if(associated(fstrsolid%contacts_bkup)) then
!      do i=1, size(fstrSOLID%contacts_bkup)
!        call fstr_contact_finalize(fstrSOLID%contacts_bkup(i))
!      end do
      deallocate(fstrsolid%contacts_bkup)
    end if
    if(associated(fstrsolid%embeds_bkup)) then
!      do i=1, size(fstrSOLID%embeds_bkup)
!        call fstr_contact_finalize(fstrSOLID%embeds_bkup(i))
!      end do
      deallocate(fstrsolid%embeds_bkup)
    end if
 
  end subroutine
 
  subroutine fstr_heat_init( fstrHEAT )
    implicit none
    type(fstr_heat) :: fstrHEAT
 
    fstrheat%STEPtot     = 0
    fstrheat%MATERIALtot = 0
    fstrheat%AMPLITUDEtot= 0
    fstrheat%T_FIX_tot   = 0
    fstrheat%Q_NOD_tot   = 0
    fstrheat%Q_VOL_tot   = 0
    fstrheat%Q_SUF_tot   = 0
    fstrheat%R_SUF_tot   = 0
    fstrheat%H_SUF_tot   = 0
    fstrheat%WL_tot      = 0
    fstrheat%beta        = -1.0d0
  end subroutine fstr_heat_init
 
  subroutine fstr_eigen_init( fstrEIG )
    implicit none
    type(fstr_eigen) :: fstrEIG
 
    fstreig%nget        = 5
    fstreig%maxiter     = 60
    fstreig%iter        = 0
    !fstrEIG%sigma       = 0.0d0
    fstreig%sigma       = 0.01d0
    fstreig%tolerance   = 1.0d-6
    fstreig%totalmass   = 0.0d0
  end subroutine fstr_eigen_init
 
  subroutine fstr_dynamic_init( fstrDYNAMIC )
    use m_fstr
    type(fstr_dynamic) :: fstrDYNAMIC
    fstrdynamic%idx_eqa  = 1
    fstrdynamic%idx_resp = 1
    fstrdynamic%n_step   = 1
    fstrdynamic%t_start  = 0.0
    fstrdynamic%t_curr   = 0.0d0
    fstrdynamic%t_end    = 1.0
    fstrdynamic%t_delta  = 1.0
    fstrdynamic%gamma    = 0.5
    fstrdynamic%beta     = 0.25
    fstrdynamic%idx_mas  = kmasslumped
    fstrdynamic%idx_dmp  = 1
    fstrdynamic%ray_m    = 0.0
    fstrdynamic%ray_k    = 0.0
    fstrdynamic%restart_nout = 0
    fstrdynamic%nout         = 100
    fstrdynamic%ngrp_monit   = 0
    fstrdynamic%nout_monit   = 1
    fstrdynamic%iout_list(1) = 0
    fstrdynamic%iout_list(2) = 0
    fstrdynamic%iout_list(3) = 0
    fstrdynamic%iout_list(4) = 0
    fstrdynamic%iout_list(5) = 0
    fstrdynamic%iout_list(6) = 0
 
  end subroutine fstr_dynamic_init
 
 
  subroutine fstr_dynamic_alloc( hecMESH, fstrDYNAMIC )
    use m_fstr
    type(hecmwst_local_mesh),target :: hecMESH
    type(fstr_dynamic) :: fstrDYNAMIC
 
    integer :: ierror, ndof,nnod
 
    ndof=hecmesh%n_dof
    nnod=hecmesh%n_node
    if(fstrdynamic%idx_eqa == 11) then
      allocate( fstrdynamic%DISP(ndof*nnod,3)  ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <fstr_solve_LINEAR_DYNAMIC, DISP>'
        write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
      allocate( fstrdynamic%VEL (ndof*nnod,1)  ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <fstr_solve_LINEAR_DYNAMIC, VEL>'
        write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
      allocate( fstrdynamic%ACC (ndof*nnod,1)  ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <fstr_solve_LINEAR_DYNAMIC, ACC>'
        write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    else
      allocate( fstrdynamic%DISP(ndof*nnod,2)  ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <fstr_solve_LINEAR_DYNAMIC, DISP>'
        write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
      allocate( fstrdynamic%VEL (ndof*nnod,2)  ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <fstr_solve_LINEAR_DYNAMIC, VEL>'
        write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
      allocate( fstrdynamic%ACC (ndof*nnod,2)  ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <fstr_solve_LINEAR_DYNAMIC, ACC>'
        write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
 
 
    allocate( fstrdynamic%VEC1(ndof*nnod)    ,stat=ierror )
    if( ierror /= 0 ) then
      write(idbg,*) 'stop due to allocation error <fstr_solve_LINEAR_DYNAMIC, VEC1>'
      write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm())
    end if
    allocate( fstrdynamic%VEC2(ndof*nnod)    ,stat=ierror )
    if( ierror /= 0 ) then
      write(idbg,*) 'stop due to allocation error <fstr_solve_LINEAR_DYNAMIC, VEC2>'
      write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm())
    end if
    allocate( fstrdynamic%VEC3(ndof*nnod)    ,stat=ierror )
    if( ierror /= 0 ) then
      write(idbg,*) 'stop due to allocation error <fstr_solve_LINEAR_DYNAMIC, VEC3>'
      write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm())
    end if
 
  end subroutine fstr_dynamic_alloc
 
  subroutine fstr_dynamic_finalize( fstrDYNAMIC )
    type(fstr_dynamic) :: fstrDYNAMIC
 
    integer :: ierror
    if( associated(fstrdynamic%DISP) )        &
      deallocate( fstrdynamic%DISP    ,stat=ierror )
    if( ierror /= 0 ) then
      write(idbg,*) 'stop due to deallocation error <fstr_solve_LINEAR_DYNAMIC, DISP>'
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm())
    end if
    if( associated(fstrdynamic%VEL) )        &
      deallocate( fstrdynamic%VEL     ,stat=ierror )
    if( ierror /= 0 ) then
      write(idbg,*) 'stop due to deallocation error <fstr_solve_LINEAR_DYNAMIC, VEL>'
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm())
    end if
    if( associated(fstrdynamic%ACC) )        &
      deallocate( fstrdynamic%ACC     ,stat=ierror )
    if( ierror /= 0 ) then
      write(idbg,*) 'stop due to deallocation error <fstr_solve_LINEAR_DYNAMIC, ACC>'
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm())
    end if
    if( associated(fstrdynamic%VEC1) )        &
      deallocate( fstrdynamic%VEC1    ,stat=ierror )
    if( ierror /= 0 ) then
      write(idbg,*) 'stop due to deallocation error <fstr_solve_LINEAR_DYNAMIC, VEC1>'
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm())
    end if
    if( associated(fstrdynamic%VEC2) )        &
      deallocate( fstrdynamic%VEC2    ,stat=ierror )
    if( ierror /= 0 ) then
      write(idbg,*) 'stop due to deallocation error <fstr_solve_LINEAR_DYNAMIC, VEC2>'
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm())
    end if
    if( associated(fstrdynamic%VEC3) )        &
      deallocate( fstrdynamic%VEC3    ,stat=ierror )
    if( ierror /= 0 ) then
      write(idbg,*) 'stop due to deallocation error <fstr_solve_LINEAR_DYNAMIC, VEC3>'
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm())
    end if
 
  end subroutine
 
 
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_post_phys_alloc(phys, NDOF, n_node, n_elem)
    implicit none
    type(fstr_solid_physic_val), pointer :: phys
    integer(kind=kint) :: NDOF, n_node, n_elem, mdof, istat
    mdof = (ndof*ndof+ndof)/2;
    allocate ( phys%STRAIN  (mdof*n_node), stat=istat)
    if( istat /= 0 ) stop "Allocation error: phys%STRAIN"
    allocate ( phys%STRESS  (mdof*n_node), stat=istat)
    if( istat /= 0 ) stop "Allocation error: phys%STRESS"
    allocate ( phys%MISES   (     n_node), stat=istat)
    if( istat /= 0 ) stop "Allocation error: phys%MISES"
    allocate ( phys%ESTRAIN (mdof*n_elem), stat=istat)
    if( istat /= 0 ) stop "Allocation error: phys%ESTRAIN"
    allocate ( phys%ESTRESS (mdof*n_elem), stat=istat)
    if( istat /= 0 ) stop "Allocation error: phys%ESTRESS"
    allocate ( phys%EMISES  (     n_elem), stat=istat)
    if( istat /= 0 ) stop "Allocation error: phys%EMISES"
    allocate ( phys%EPLSTRAIN (   n_elem), stat=istat)
    if( istat /= 0 ) stop "Allocation error: phys%EPLSTRAIN"
    allocate ( phys%ENQM    (12*n_elem), stat=istat)
    if( istat /= 0 ) stop "Allocation error: phys%ENQM"
  end subroutine fstr_setup_post_phys_alloc
 
  subroutine fstr_setup_post( ctrl, P )
    implicit none
    integer(kind=kint) :: ctrl, i
    type(fstr_param_pack) :: P
    type(fstr_solid_physic_val), pointer :: phys => null()
 
    if( p%PARAM%solution_type == kststatic &
        .or. p%PARAM%solution_type == ksteigen  &
        .or. p%PARAM%solution_type == kstdynamic &
        .or. p%PARAM%solution_type == kststaticeigen ) then
      ! Memory Allocation for Result Vectors ------------
      if( p%MESH%n_dof == 6 .or. p%SOLID%is_33shell == 1 ) then
        allocate ( p%SOLID%SHELL )
        call fstr_setup_post_phys_alloc(p%SOLID%SHELL,3, p%MESH%n_node,p%MESH%n_elem)
        allocate ( p%SOLID%SHELL%LAYER(p%SOLID%max_lyr) )
        do i=1,p%SOLID%max_lyr
          allocate ( p%SOLID%SHELL%LAYER(i)%PLUS )
          allocate ( p%SOLID%SHELL%LAYER(i)%MINUS )
          call fstr_setup_post_phys_alloc(p%SOLID%SHELL%LAYER(i)%PLUS , 3, p%MESH%n_node, p%MESH%n_elem)
          call fstr_setup_post_phys_alloc(p%SOLID%SHELL%LAYER(i)%MINUS, 3, p%MESH%n_node, p%MESH%n_elem)
        enddo
        phys => p%SOLID%SHELL
      else
        allocate ( p%SOLID%SOLID )
        phys => p%SOLID%SOLID
        call fstr_setup_post_phys_alloc(phys, p%MESH%n_dof, p%MESH%n_node,  p%MESH%n_elem)
      end if
      p%SOLID%STRAIN => phys%STRAIN
      p%SOLID%STRESS => phys%STRESS
      p%SOLID%MISES  => phys%MISES
      p%SOLID%ESTRAIN => phys%ESTRAIN
      p%SOLID%ESTRESS => phys%ESTRESS
      p%SOLID%EMISES  => phys%EMISES
      p%SOLID%EPLSTRAIN  => phys%EPLSTRAIN
      p%SOLID%ENQM    => phys%ENQM
      allocate( p%SOLID%REACTION( p%MESH%n_dof*p%MESH%n_node ), stat=i )
      if( i /= 0 ) stop "Allocation error: REACTION"
      p%SOLID%REACTION = 0.d0
    end if
 
    if( p%PARAM%fg_visual == kon )then
      call fstr_setup_visualize( ctrl, p%MESH )
    end if
 
    call hecmw_barrier( p%MESH ) ! JP-7
 
    if( p%HEAT%STEPtot == 0 ) then ! No !HEAT Input
      if( p%PARAM%analysis_n == 0 ) then  ! No !STATIC Input
        call reallocate_real( p%PARAM%dtime, 1)
        call reallocate_real( p%PARAM%etime, 1)
        call reallocate_real( p%PARAM%dtmin, 1)
        call reallocate_real( p%PARAM%delmax,1)
        call reallocate_integer( p%PARAM%itmax, 1)
        call reallocate_real( p%PARAM%eps,   1)
        p%PARAM%analysis_n = 1
        p%PARAM%dtime = 0
        p%PARAM%etime = 0
        p%PARAM%dtmin = 0
        p%PARAM%delmax = 0
        p%PARAM%itmax = 20
        p%PARAM%eps = 1.0e-6
      end if
      p%HEAT%STEPtot = 1
      call reallocate_real( p%HEAT%STEP_DLTIME, 1)
      call reallocate_real( p%HEAT%STEP_EETIME, 1)
      call reallocate_real( p%HEAT%STEP_DELMIN, 1)
      call reallocate_real( p%HEAT%STEP_DELMAX, 1)
      p%HEAT%STEP_DLTIME = 0
      p%HEAT%STEP_EETIME = 0
      p%HEAT%STEP_DELMIN = 0
      p%HEAT%STEP_DELMAX = 0
    end if
  end subroutine fstr_setup_post
 
  !*****************************************************************************!
  !* GENERAL HEADERS ***********************************************************!
  !*****************************************************************************!
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_solution( ctrl, counter, P )
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
 
    rcode = fstr_ctrl_get_solution( ctrl, p%PARAM%solution_type, p%PARAM%nlgeom )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
  end subroutine fstr_setup_solution
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_nonlinear_solver( ctrl, counter, P )
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
 
    rcode = fstr_ctrl_get_nonlinear_solver( ctrl, p%PARAM%nlsolver_method )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
  end subroutine fstr_setup_nonlinear_solver
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_solver( ctrl, counter, P )
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack),target :: P
 
    integer(kind=kint) :: rcode
 
    if( counter >= 2 ) then
      write(ilog,*) '### Error : !SOLVER exists twice in FSTR control file.'
      stop hecmw_exit_input
    endif
 
    !   nier       => svIarray(1)
    !   method     => svIarray(2)
    !   precond    => svIarray(3)
    !   nset       => svIarray(4)
    !   iterpremax => svIarray(5)
    !   nrest      => svIarray(6)
    !   scaling    => svIarray(7)
    !   iterlog    => svIarray(21)
    !   timelog    => svIarray(22)
    !   steplog    => svIarray(23)
    !   loglevel   => svIarray(24)
    !   dumptype   => svIarray(31)
    !   dumpexit   => svIarray(32)
    !   usejad     => svIarray(33)
    !   ncolor_in  => svIarray(34)
    !   mpc_method => svIarray(13)
    !   estcond    => svIarray(14)
    !   contact_elim=> svIarray(15)
    !   recompute_residual => svIarray(16)
    !   method2    => svIarray(8)
    !   recyclepre => svIarray(35)
    !   matvec_impl=> svIarray(36)
    !   precond_impl=> svIarray(37)
    !   solver_opt => svIarray(41:50)
    !   nBFGS      => svIarray(60)
 
    !   resid      => svRarray(1)
    !   sigma_diag => svRarray(2)
    !   omega      => svRarray(3)
    !   thresh     => svRarray(4)
    !   filter     => svRarray(5)
    !   solver_ropt=> svRarray(41:50)
 
    rcode = fstr_ctrl_get_solver( ctrl,                      &
      sviarray(2), sviarray(3), sviarray(4), sviarray(21), sviarray(22), sviarray(23),&
      sviarray(1), sviarray(5), sviarray(6), sviarray(60), sviarray(7), &
      sviarray(31), sviarray(32), sviarray(33), sviarray(34), sviarray(13), sviarray(14), sviarray(8),&
      sviarray(35), sviarray(41:50), sviarray(15), &
      svrarray(1), svrarray(2), svrarray(3),                &
      svrarray(4), svrarray(5), svrarray(41:50), sviarray(24), &
      sviarray(36), sviarray(37), sviarray(16) )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
    if( sviarray(2) <= 100 ) then
      sviarray(99) = 1  ! indirect method
    else
      sviarray(99) = sviarray(2)-99 !2  ! direct method
    end if
 
  end subroutine fstr_setup_solver
 
  !* ----------------------------------------------------------------------------------------------- *!
  !* ----------------------------------------------------------------------------------------------- *!
 
  integer function fstr_setup_orientation( ctrl, hecMESH, cnt, coordsys )
    use hecmw_setup_util, only: node_global_to_local
    implicit none
    integer(kind=kint)         :: ctrl
    type( hecmwst_local_mesh ) :: hecmesh
    integer                    :: cnt
    type( tlocalcoordsys )     :: coordsys
 
    character(len=HECMW_NAME_LEN) :: grp_id_name(1)
 
    integer :: nid, dtype
    character(len=HECMW_NAME_LEN) :: data_fmt
    real(kind=kreal) :: fdum, xyza(3), xyzb(3), xyzc(3), ff1(3), ff2(3), ff3(3)
 
    fstr_setup_orientation = -1
 
    nid = 1
    coordsys%sys_type = 10
 
    nid = 1
    data_fmt = 'COORDINATES,NODES,LOCAL_NODES '
    if( fstr_ctrl_get_param_ex( ctrl, 'DEFINITION ', data_fmt, 0, 'P', nid )/=0 ) return
    dtype = nid-1
    coordsys%sys_type = coordsys%sys_type + dtype
 
    if( fstr_ctrl_get_param_ex( ctrl, 'NAME ',  '# ',  1, 'S', grp_id_name(1) )/= 0) return
    coordsys%sys_name = grp_id_name(1)
 
    if( dtype==0 ) then
      data_fmt = "RRRRRRrrr "
      xyzc(:) = 0.d0
      if( fstr_ctrl_get_data_ex( ctrl, 1, data_fmt, xyza(1), xyza(2),  &
        xyza(3), xyzb(1), xyzb(2), xyzb(3), xyzc(1), xyzc(2), xyzc(3) )/=0 ) return
      if( coordsys%sys_type==10 ) then
        ff1 = xyza-xyzc
        fdum = dsqrt( dot_product(ff1, ff1) )
        if( fdum==0.d0 ) return
        ff1 = ff1/fdum
        ff2 = xyzb-xyzc
        call cross_product(ff1,ff2,ff3)
        coordsys%CoordSys(1,:) = ff1
 
        fdum = dsqrt( dot_product(ff3, ff3) )
        if( fdum==0.d0 ) return
        coordsys%CoordSys(3,:) = ff3/fdum
 
        call cross_product(coordsys%CoordSys(3,:), coordsys%CoordSys(1,:), coordsys%CoordSys(2,:) )
      else
        coordsys%CoordSys(1,:) = xyza
        coordsys%CoordSys(2,:) = xyzb
      endif
 
    else
      coordsys%node_ID(3) = 0   ! global origin
      data_fmt = "IIi "
      if( fstr_ctrl_get_data_ex( ctrl, 1, data_fmt, coordsys%node_ID(1),  &
        coordsys%node_ID(2), coordsys%node_ID(3) )/=0 ) return
      ! DEFINITION=LOCAL_NODES: three local node numbers (1..number of nodes of
      ! the element) defining the coordinate system of each element. All three
      ! are mandatory here, unlike the third node of DEFINITION=NODES.
      if( dtype==2 ) then
        if( maxval(coordsys%node_ID(1:3)) > 10 .or. minval(coordsys%node_ID(1:3)) < 1 ) then
          write(*,*) "!ORIENTATION, DEFINITION=LOCAL_NODES needs three local node numbers (1-10)!"
          write(idbg,*) "!ORIENTATION, DEFINITION=LOCAL_NODES needs three local node numbers (1-10)!"
          return
        endif
        fstr_setup_orientation = 0
        return
      endif
      if( coordsys%node_ID(3) == 0 ) then
        nid = node_global_to_local( hecmesh, coordsys%node_ID(1:2), 2 )
        if( nid/=0 .and. nid/=2 ) then
          write(*,*) "We cannot define coordinate system using nodes in other CPU!"
          write(idbg,*) "We cannot define coordinate system using nodes in other CPU!"
          return
        endif
      else
        nid = node_global_to_local( hecmesh, coordsys%node_ID, 3 )
        if( nid/=0 .and. nid/=3 ) then
          write(*,*) "We cannot define coordinate system using nodes in other CPU!"
          write(idbg,*) "We cannot define coordinate system using nodes in other CPU!"
          return
        endif
      endif
    endif
 
    fstr_setup_orientation = 0
  end function fstr_setup_orientation
 
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_step( ctrl, counter, P )
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
    character(HECMW_NAME_LEN) :: amp
    integer(kind=kint) :: amp_id
 
    integer(kind=kint) :: rcode, iproc
 
    amp = ' '
    rcode = fstr_ctrl_get_step( ctrl, amp, iproc )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
    call amp_name_to_id( p%MESH, '!STEP', amp, amp_id )
    !    P%SOLID%NLSTATIC_ngrp_amp = amp_id;
 
  end subroutine fstr_setup_step
 
  subroutine fstr_setup_step_convergence( stepinfo, cpname, cnvparam, istep )
    implicit none
    type( step_info ), intent(inout)  :: stepinfo
    character(len=*), intent(in)      :: cpname
    type( tparamconverg ), intent(in) :: cnvparam(0:)
    integer(kind=kint), intent(in)    :: istep
 
    integer(kind=kint) :: i, id
 
    id = 0
    if( len_trim(cpname) > 0 ) then
      do i=1,ubound(cnvparam,1)
        if( hecmw_streqr( cnvparam(i)%name, cpname ) ) then
          id = i; exit
        endif
      enddo
      if( id == 0 ) then
        write(*,*) '### Error: CONVERG_PARAM referred to by STEP is not defined before it : ', trim(cpname)
        write(ilog,*) '### Error: CONVERG_PARAM referred to by STEP is not defined before it : ', trim(cpname)
        stop hecmw_exit_input
      endif
    endif
 
    call setup_stepinfo_converg( stepinfo, cnvparam(id) )
 
    if( .not. any( stepinfo%cnv_check(kcnvresidual:kcnvcorrection, kcnvtranslation:kcnvrotation, :) ) ) then
      write(*,*) '### Error: RESIDUAL and CORRECTION are all SKIP in step : ', istep
      write(ilog,*) '### Error: RESIDUAL and CORRECTION are all SKIP in step : ', istep
      stop hecmw_exit_input
    endif
 
  end subroutine fstr_setup_step_convergence
 
  integer(kind=kint) function fstr_setup_initial( ctrl, cond, hecMESH )
    implicit none
    integer(kind=kint)        :: ctrl
    type( tinitialcondition ) :: cond
    type(hecmwst_local_mesh)  :: hecmesh
    integer, pointer          :: grp_id(:), dof(:)
    real(kind=kreal), pointer :: temp(:)
    character(len=HECMW_NAME_LEN), pointer :: grp_id_name(:)
    character(len=HECMW_NAME_LEN) :: data_fmt, ss
    integer :: i,j,n, is, ie, gid, nid, rcode
 
    fstr_setup_initial = -1
 
    ss = 'TEMPERATURE,VELOCITY,ACCELERATION '
    rcode = fstr_ctrl_get_param_ex( ctrl, 'TYPE ', ss, 1, 'P', nid )
    if( nid==1 ) then
      cond%cond_name = "temperature"
      allocate( cond%intval(hecmesh%n_node) )
      allocate( cond%realval(hecmesh%n_node) )
    elseif( nid==2 ) then
      cond%cond_name = "velocity"
      allocate( cond%intval(hecmesh%n_node) )
      allocate( cond%realval(hecmesh%n_node) )
    elseif( nid==3 ) then
      cond%cond_name = "acceleration"
      allocate( cond%intval(hecmesh%n_node) )
      allocate( cond%realval(hecmesh%n_node) )
    else
      return
    endif
 
    cond%intval = -1
    cond%realval = 0.d0
 
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n<=0 ) return
    allocate( temp(n), grp_id_name(n), grp_id(n), dof(n) )
    dof = 0
    write(ss,*)  hecmw_name_len
    if( nid==1 ) then
      write(data_fmt,'(a,a,a)') 'S',trim(adjustl(ss)),'R '
      fstr_setup_initial = &
            fstr_ctrl_get_data_array_ex( ctrl, data_fmt, grp_id_name, temp )
    else
        write(data_fmt,'(a,a,a)') 'S',trim(adjustl(ss)),'IR '
        fstr_setup_initial = &
            fstr_ctrl_get_data_array_ex( ctrl, data_fmt, grp_id_name, dof, temp )
    endif
 
    if( fstr_setup_initial /= 0 ) then
        if( associated(grp_id) )    deallocate( grp_id )
        if( associated(temp) )    deallocate( temp )
        if( associated(dof) )       deallocate( dof )
        if( associated(grp_id_name) )    deallocate( grp_id_name )
        return
    end if
 
    call node_grp_name_to_id_ex( hecmesh, '!INITIAL CONDITION', n, grp_id_name, grp_id )
    do i=1,n
      gid = grp_id(i)
      is = hecmesh%node_group%grp_index(gid-1) + 1
      ie = hecmesh%node_group%grp_index(gid  )
      do j=is, ie
        nid = hecmesh%node_group%grp_item(j)
        cond%realval(nid) = temp(i)
        cond%intval(nid) = dof(i)
      enddo
    enddo
 
    if( associated(grp_id) )    deallocate( grp_id )
    if( associated(temp) )      deallocate( temp )
    if( associated(dof) )       deallocate( dof )
    if( associated(grp_id_name) )    deallocate( grp_id_name )
end function fstr_setup_initial
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_write( ctrl, counter, P )
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
    integer(kind=kint) :: res, visual, neutral
 
    integer(kind=kint) :: rcode
 
    rcode = fstr_ctrl_get_write( ctrl, res, visual, neutral )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
    if( res == 1 ) p%PARAM%fg_result = 1
    if( visual == 1 ) p%PARAM%fg_visual = 1
    if( neutral == 1 ) p%PARAM%fg_neutral = 1
 
  end subroutine fstr_setup_write
 
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
  subroutine fstr_setup_echo( ctrl, counter, P )
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
 
    rcode = fstr_ctrl_get_echo( ctrl,        &
      p%PARAM%fg_echo )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
  end subroutine fstr_setup_echo
 
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
  subroutine fstr_setup_restart( ctrl, nout, version )
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: nout
    integer(kind=kint) :: version
 
    integer(kind=kint) :: rcode
    nout = 0
    rcode = fstr_ctrl_get_param_ex( ctrl, 'FREQUENCY ', '# ', 0, 'I', nout )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
    rcode = fstr_ctrl_get_param_ex( ctrl, 'VERSION ', '# ', 0, 'I', version )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
  end subroutine fstr_setup_restart
 
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_couple( ctrl, counter, P )
    use hecmw_setup_util, only &
      : hecmw_expand_integer_array
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
    integer(kind=kint) :: rcode
    character(HECMW_NAME_LEN), pointer :: grp_id_name(:)
    integer(kind=kint) :: n, old_size, new_size
 
    if( p%SOLID%file_type /= kbcffstr ) return
 
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n == 0 ) return
    old_size = p%SOLID%COUPLE_ngrp_tot
    new_size = old_size + n
    p%SOLID%COUPLE_ngrp_tot = new_size
 
    call hecmw_expand_integer_array ( p%SOLID%COUPLE_ngrp_ID,  old_size, new_size )
 
    allocate( grp_id_name(n))
    rcode = fstr_ctrl_get_couple( ctrl,           &
      p%PARAM%fg_couple_type,       &
      p%PARAM%fg_couple_first,      &
      p%PARAM%fg_couple_window,     &
      grp_id_name, hecmw_name_len  )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
    call surf_grp_name_to_id_ex( p%MESH, '!COUPLE', &
      n, grp_id_name, p%SOLID%COUPLE_ngrp_ID(old_size+1:))
 
    deallocate( grp_id_name )
    p%PARAM%fg_couple = 1
 
  end subroutine fstr_setup_couple
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_amplitude( ctrl, P )
    implicit none
    integer(kind=kint) :: ctrl
    type(fstr_param_pack) :: P
    real(kind=kreal), pointer :: val(:), table(:)
    character(len=HECMW_NAME_LEN) :: name
    integer :: nline, n, type_def, type_time, type_val, rcode
 
    nline = fstr_ctrl_get_data_line_n( ctrl )
    if( nline<=0 ) return
    allocate( val(nline*4) )
    allocate( table(nline*4) )
    rcode = fstr_ctrl_get_amplitude( ctrl, nline, name, type_def, type_time, type_val, &
        n, val, table )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
    call append_new_amplitude( p%MESH%amp, name, type_def, type_time, type_val, n, val, table )
 
    if( associated(val) ) deallocate( val )
    if( associated(table) ) deallocate( table )
  end subroutine fstr_setup_amplitude
 
 
  subroutine fstr_setup_element_activation( ctrl, counter, P )
    use hecmw_setup_util, only &
      : hecmw_expand_integer_array &
      , hecmw_expand_real_array
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
    character(HECMW_NAME_LEN) :: amp
    integer(kind=kint) :: amp_id
    character(HECMW_NAME_LEN), pointer :: grp_id_name(:)
    integer(kind=kint) :: i, n, old_size, new_size
    integer(kind=kint) :: gid, mode, measure, state
    real(kind=kreal)   :: eps
    real(kind=kreal), pointer :: thlow(:), thup(:)
 
    gid = 1
    rcode = fstr_ctrl_get_param_ex( ctrl, 'GRPID ',  '# ',            0, 'I', gid  )
 
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n == 0 ) return
    old_size = p%SOLID%elemact%ELEMACT_egrp_tot
    new_size = old_size + n
    p%SOLID%elemact%ELEMACT_egrp_tot = new_size
 
    call hecmw_expand_integer_array ( p%SOLID%elemact%ELEMACT_egrp_GRPID, old_size, new_size )
    call hecmw_expand_integer_array ( p%SOLID%elemact%ELEMACT_egrp_ID,  old_size, new_size )
    call hecmw_expand_integer_array ( p%SOLID%elemact%ELEMACT_egrp_amp, old_size, new_size )
    call hecmw_expand_real_array ( p%SOLID%elemact%ELEMACT_egrp_eps, old_size, new_size )
    call hecmw_expand_integer_array ( p%SOLID%elemact%ELEMACT_egrp_depends, old_size, new_size )
    call hecmw_expand_real_array ( p%SOLID%elemact%ELEMACT_egrp_ts_lower, old_size, new_size )
    call hecmw_expand_real_array ( p%SOLID%elemact%ELEMACT_egrp_ts_upper, old_size, new_size )
    call hecmw_expand_integer_array ( p%SOLID%elemact%ELEMACT_egrp_state, old_size, new_size )
 
    allocate( grp_id_name(n), thlow(n), thup(n) )
    amp = ' '
    eps = 1.d-3
    rcode = fstr_ctrl_get_element_activation( ctrl, amp, eps, grp_id_name, mode, measure, state, thlow, thup )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
    call amp_name_to_id( p%MESH, '!ELEMENT_ACTIVATION', amp, amp_id )
    do i=1,n
      p%SOLID%elemact%ELEMACT_egrp_amp(old_size+i) = amp_id
      p%SOLID%elemact%ELEMACT_egrp_eps(old_size+i) = eps
    end do
    p%SOLID%elemact%ELEMACT_egrp_GRPID(old_size+1:new_size) = gid
    p%SOLID%elemact%ELEMACT_egrp_depends(old_size+1:new_size) = measure
    p%SOLID%elemact%ELEMACT_egrp_ts_lower(old_size+1:new_size) = thlow(1:n)
    p%SOLID%elemact%ELEMACT_egrp_ts_upper(old_size+1:new_size) = thup(1:n)
    p%SOLID%elemact%ELEMACT_egrp_state(old_size+1:new_size) = state
 
    call elem_grp_name_to_id_ex( p%MESH, '!ELEMENT_ACTIVATION', n, grp_id_name, p%SOLID%elemact%ELEMACT_egrp_ID(old_size+1:))
 
    deallocate( grp_id_name )
  end subroutine fstr_setup_element_activation
 
 
  !*****************************************************************************!
  !* HEADERS FOR STATIC ANALYSIS ***********************************************!
  !*****************************************************************************!
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_static( ctrl, counter, P )
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
    integer(kind=kint) :: rcode
 
    integer :: nout, nout_monit,node_monit_1 ,elem_monit_1 ,intg_monit_1
    integer :: ipt, idx_elpl, iout_list(6)
    real(kind=kreal) :: sig_y0, h_dash
 
    if( counter > 1 ) then
      write(*,*)
    endif
 
    ipt = 0
    if( fstr_ctrl_get_param_ex( ctrl, 'TYPE ', 'INFINITESIMAL,NLGEOM,INFINITE ', 0, 'P', ipt  )/=0 )  &
      return
    if( ipt == 2 ) p%PARAM%nlgeom = .true.
 
    ! for backward compatibility
    if( ipt == 3 ) then
      write(*,*) "Warning : !STATIC : parameter 'TYPE=INFINITE' is deprecated." &
           & //  " Please use the replacement parameter 'TYPE=INFINITESIMAL'"
    endif
 
    rcode = fstr_ctrl_get_static( ctrl, &
      dt, etime, itmax, eps, p%SOLID%restart_nout, &
      idx_elpl, &
      iout_list, &
      sig_y0, h_dash, &
      nout, nout_monit, node_monit_1, &
      elem_monit_1, intg_monit_1 )
 
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
  end subroutine fstr_setup_static
 
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_boundary( ctrl, counter, P )
    use hecmw_setup_util, only &
      : hecmw_expand_integer_array &
      , hecmw_expand_real_array
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
    character(HECMW_NAME_LEN) :: amp, rotc_name(1)
    integer(kind=kint) :: amp_id, rotc_id(1), n_rotc
    character(HECMW_NAME_LEN), pointer :: grp_id_name(:)
    integer(kind=kint),pointer :: dof_ids (:)
    integer(kind=kint),pointer :: dof_ide (:)
    real(kind=kreal),pointer :: val_ptr(:)
    integer(kind=kint) :: i, n, old_size, new_size
 
    integer(kind=kint) :: gid, istot
 
    gid = 1
    rcode = fstr_ctrl_get_param_ex( ctrl, 'GRPID ',  '# ',            0, 'I', gid  )
 
    istot = 0
    rcode = fstr_ctrl_get_param_ex( ctrl, 'TOTAL ', '# ', 0, 'E', istot )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
    ! get center of torque load
    rotc_name = ' '
    rotc_id = -1
    n_rotc = -1
    rcode = fstr_ctrl_get_param_ex( ctrl, 'ROT_CENTER ', '# ', 0, 'S', rotc_name )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
    if(  rotc_name(1) /= ' ' ) then
      if( istot /= 0 ) then
        write(*,*) 'fstr control file error : !BOUNDARY : rotational boundary cannot be specified with total value'
        write(ilog,*) 'fstr control file error : !BOUNDARY : rotational boundary cannot be specified with total value'
        call fstr_ctrl_err_stop
      endif
      p%SOLID%BOUNDARY_ngrp_rot = p%SOLID%BOUNDARY_ngrp_rot + 1
      n_rotc = p%SOLID%BOUNDARY_ngrp_rot
      call node_grp_name_to_id_ex( p%MESH, '!BOUNDARY,ROT_CENTER=', 1, rotc_name, rotc_id)
    endif
 
 
    ! ENTIRE -----------------------------------------------
    p%SOLID%file_type = kbcffstr
 
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n == 0 ) return
    old_size = p%SOLID%BOUNDARY_ngrp_tot
    new_size = old_size + n
    p%SOLID%BOUNDARY_ngrp_tot = new_size
    call hecmw_expand_integer_array (p%SOLID%BOUNDARY_ngrp_GRPID, old_size, new_size )
    call hecmw_expand_integer_array (p%SOLID%BOUNDARY_ngrp_ID,    old_size, new_size )
    call hecmw_expand_integer_array (p%SOLID%BOUNDARY_ngrp_type,  old_size, new_size )
    call hecmw_expand_real_array    (p%SOLID%BOUNDARY_ngrp_val,   old_size, new_size )
    call hecmw_expand_integer_array (p%SOLID%BOUNDARY_ngrp_amp,   old_size, new_size )
    call hecmw_expand_integer_array (p%SOLID%BOUNDARY_ngrp_istot, old_size, new_size )
    call hecmw_expand_integer_array (p%SOLID%BOUNDARY_ngrp_rotID, old_size, new_size )
    call hecmw_expand_integer_array (p%SOLID%BOUNDARY_ngrp_centerID, old_size, new_size )
 
    allocate( grp_id_name(n) )
    allocate( dof_ids(n) )
    allocate( dof_ide(n) )
    allocate( val_ptr(n) )
 
    amp = ' '
    val_ptr = 0.0d0
    rcode = fstr_ctrl_get_boundary( ctrl, amp, grp_id_name, hecmw_name_len, dof_ids, dof_ide, val_ptr)
    if( rcode /= 0 ) call fstr_ctrl_err_stop
    call amp_name_to_id( p%MESH, '!BOUNDARY', amp, amp_id )
    p%SOLID%BOUNDARY_ngrp_GRPID(old_size+1:new_size) = gid
    call node_grp_name_to_id_ex( p%MESH, '!BOUNDARY', n, grp_id_name, p%SOLID%BOUNDARY_ngrp_ID(old_size+1:))
    p%SOLID%BOUNDARY_ngrp_istot(old_size+1:new_size) = istot
 
    ! set up information about rotation ( default value is set if ROT_CENTER is not given.)
    p%SOLID%BOUNDARY_ngrp_rotID(old_size+1:) = n_rotc
    p%SOLID%BOUNDARY_ngrp_centerID(old_size+1:) = rotc_id(1)
 
    do i = 1, n
      if( (dof_ids(i) < 1).or.(6 < dof_ids(i)).or.(dof_ide(i) < 1).or.(6 < dof_ide(i)) ) then
        write(*,*) 'fstr control file error : !BOUNDARY : range of dof_ids and dof_ide is from 1 to 6'
        write(ilog,*) 'fstr control file error : !BOUNDARY : range of dof_ids and dof_ide is from 1 to 6'
        call fstr_ctrl_err_stop
      end if
      p%SOLID%BOUNDARY_ngrp_val(old_size+i) = val_ptr(i)
      p%SOLID%BOUNDARY_ngrp_type(old_size+i) = 10 * dof_ids(i) + dof_ide(i)
      p%SOLID%BOUNDARY_ngrp_amp(old_size+i) = amp_id
    end do
 
    deallocate( grp_id_name )
    deallocate( dof_ids )
    deallocate( dof_ide )
    deallocate( val_ptr )
    nullify( grp_id_name )
    nullify( dof_ids )
    nullify( dof_ide )
    nullify( val_ptr )
 
  end subroutine fstr_setup_boundary
 
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_cload( ctrl, counter, P )
    use hecmw_setup_util, only &
      : hecmw_expand_integer_array &
      , hecmw_expand_real_array
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
    character(HECMW_NAME_LEN) :: amp, rotc_name(1)
    integer(kind=kint) :: amp_id, rotc_id(1), n_rotc
    character(HECMW_NAME_LEN), pointer :: grp_id_name(:)
    real(kind=kreal),pointer :: val_ptr(:)
    integer(kind=kint),pointer :: id_ptr(:)
    integer(kind=kint) :: i, n, old_size, new_size
    integer(kind=kint) :: gid, ngrp_id, nnode_grp
    character(len=HECMW_MSG_LEN) :: msg(2)
 
    if( p%SOLID%file_type /= kbcffstr ) return
    gid = 1
    rcode = fstr_ctrl_get_param_ex( ctrl, 'GRPID ',  '# ',            0, 'I', gid  )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
    ! get center of torque load
    rotc_name = ' '
    rotc_id = -1
    n_rotc = -1
    rcode = fstr_ctrl_get_param_ex( ctrl, 'ROT_CENTER ', '# ', 0, 'S', rotc_name )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
    if(  rotc_name(1) /= ' ' ) then
      p%SOLID%CLOAD_ngrp_rot = p%SOLID%CLOAD_ngrp_rot + 1
      n_rotc = p%SOLID%CLOAD_ngrp_rot
      call node_grp_name_to_id_ex( p%MESH, '!CLOAD,ROT_CENTER=', 1, rotc_name, rotc_id)
    endif
 
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n == 0 ) return
    old_size = p%SOLID%CLOAD_ngrp_tot
    new_size = old_size + n
    p%SOLID%CLOAD_ngrp_tot = new_size
    ! Keiji Suemitsu (20140624) <
    call hecmw_expand_integer_array ( p%SOLID%CLOAD_ngrp_GRPID, old_size, new_size )
    call hecmw_expand_integer_array ( p%SOLID%CLOAD_ngrp_ID,  old_size, new_size )
    call hecmw_expand_integer_array ( p%SOLID%CLOAD_ngrp_DOF, old_size, new_size )
    call hecmw_expand_real_array    ( p%SOLID%CLOAD_ngrp_val, old_size, new_size )
    call hecmw_expand_integer_array ( p%SOLID%CLOAD_ngrp_amp, old_size, new_size )
    call hecmw_expand_integer_array ( p%SOLID%CLOAD_ngrp_rotID, old_size, new_size )
    call hecmw_expand_integer_array ( p%SOLID%CLOAD_ngrp_centerID, old_size, new_size )
    ! > Keiji Suemitsu (20140624)
 
    allocate( grp_id_name(n))
    allocate( id_ptr(n) )
    allocate( val_ptr(n) )
    amp = ' '
    id_ptr = 0
    val_ptr = 0.0d0
    rcode = fstr_ctrl_get_cload( ctrl, amp, grp_id_name, hecmw_name_len, id_ptr, val_ptr )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
    ! set up information about torque load ( default value is set if ROT_CENTER is not given.)
    p%SOLID%CLOAD_ngrp_rotID(old_size+1:) = n_rotc
    p%SOLID%CLOAD_ngrp_centerID(old_size+1:) = rotc_id(1)
 
    call amp_name_to_id( p%MESH, '!CLOAD', amp, amp_id )
    do i=1,n
      p%SOLID%CLOAD_ngrp_amp(old_size+i) = amp_id
      p%SOLID%CLOAD_ngrp_DOF(old_size+i) = id_ptr(i)
      p%SOLID%CLOAD_ngrp_val(old_size+i) = val_ptr(i)
    end do
    p%SOLID%CLOAD_ngrp_GRPID(old_size+1:new_size) = gid
    call node_grp_name_to_id_ex( p%MESH, '!CLOAD', n, grp_id_name, p%SOLID%CLOAD_ngrp_ID(old_size+1:))
 
    deallocate( grp_id_name )
    deallocate( id_ptr )
    deallocate( val_ptr )
    nullify( grp_id_name )
    nullify( id_ptr )
    nullify( val_ptr )
 
    ! the refinement adds nodes to the group, so the load given to each node of the group
    ! sums up to more than intended before the refinement
    if( p%MESH%n_refine > 0 ) then
      do i=1,n
        ngrp_id = p%SOLID%CLOAD_ngrp_ID(old_size+i)
        nnode_grp = hecmw_ngrp_get_number(p%MESH, ngrp_id)
        if( nnode_grp > 1 .and. p%MESH%my_rank == 0 ) then
          msg(1) = ' ### Warning: !CLOAD on a refined mesh'
          write(msg(2),'(a,a,a,i0,a)') '   node group "', trim(p%MESH%node_group%grp_name(ngrp_id)), &
            '" now has ', nnode_grp, ' nodes; the same load is applied to each of them'
          write(*,'(/a/a/)') trim(msg(1)), trim(msg(2))
          write(ilog,'(/a/a/)') trim(msg(1)), trim(msg(2))
        endif
      enddo
    endif
 
  end subroutine fstr_setup_cload
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
  subroutine fstr_setup_fload( ctrl, counter, P )
    use hecmw_setup_util, only &
      : hecmw_expand_integer_array &
      , hecmw_expand_real_array
    implicit none
  !---- args
    integer(kind=kint)   :: ctrl
    integer(kind=kint)   :: counter
    type(fstr_param_pack) :: P
  !---- vals
    integer(kind=kint)                  :: rcode
    character(HECMW_NAME_LEN), pointer :: grp_id_name(:)
    real(kind=kreal), pointer           :: val_ptr(:)
    integer(kind=kint), pointer        :: id_ptr(:)
    integer(kind=kint)                  :: i, n, old_size, new_size
    integer(kind=kint)                  :: gid, loadcase
  !---- body
 
    if( p%SOLID%file_type /= kbcffstr) return
 
    !read grpid
    gid = 1
    rcode = fstr_ctrl_get_param_ex( ctrl, 'GRPID ',  '# ',  0, 'I', gid )
    !read loadcase (real=1:default, img=2)
    loadcase = kfloadcase_re
    rcode = fstr_ctrl_get_param_ex( ctrl, 'LOAD CASE ', '# ', 0, 'I', loadcase)
    !write(*,*) "loadcase=", loadcase
    !pause
 
    !read the num of dataline
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n == 0 ) return
    old_size = p%FREQ%FLOAD_ngrp_tot
    new_size = old_size + n
 
    !expand data array
    p%FREQ%FLOAD_ngrp_tot = new_size
    call hecmw_expand_integer_array( p%FREQ%FLOAD_ngrp_GRPID, old_size, new_size )
    call hecmw_expand_integer_array( p%FREQ%FLOAD_ngrp_ID,    old_size, new_size )
    call hecmw_expand_integer_array( p%FREQ%FLOAD_ngrp_TYPE,  old_size, new_size )
    call hecmw_expand_integer_array( p%FREQ%FLOAD_ngrp_DOF,   old_size, new_size )
    call hecmw_expand_real_array   ( p%FREQ%FLOAD_ngrp_valre, old_size, new_size )
    call hecmw_expand_real_array   ( p%FREQ%FLOAD_ngrp_valim, old_size, new_size )
 
    !fill bc data
    allocate( grp_id_name(n) )
    allocate( id_ptr(n) )
    allocate( val_ptr(n) )
    id_ptr  = 0
    val_ptr = 0.0d0
    rcode = fstr_ctrl_get_fload( ctrl, grp_id_name, hecmw_name_len, id_ptr, val_ptr)
    if( rcode /= 0 ) call fstr_ctrl_err_stop
    if(loadcase == kfloadcase_re) then
      do i = 1, n
        p%FREQ%FLOAD_ngrp_DOF(old_size+i) = id_ptr(i)
        p%FREQ%FLOAD_ngrp_valre(old_size+i) = val_ptr(i)
      enddo
    else if(loadcase == kfloadcase_im) then
      do i = 1, n
        p%FREQ%FLOAD_ngrp_DOF(old_size+i) = id_ptr(i)
        p%FREQ%FLOAD_ngrp_valim(old_size+i) = val_ptr(i)
      enddo
    else
      !error
      write(*,*)    "Error this load set is not defined!"
      write(ilog,*) "Error this load set is not defined!"
      stop hecmw_exit_model
    end if
    p%FREQ%FLOAD_ngrp_GRPID(old_size+1:new_size) = gid
    call nodesurf_grp_name_to_id_ex( p%MESH, '!FLOAD', n, grp_id_name, &
         p%FREQ%FLOAD_ngrp_ID(old_size+1:), p%FREQ%FLOAD_ngrp_TYPE(old_size+1:))
 
    deallocate( grp_id_name )
    deallocate( id_ptr )
    deallocate( val_ptr )
    nullify( grp_id_name )
    nullify( id_ptr )
    nullify( val_ptr )
    return
 
    contains
 
    function fstr_ctrl_get_fload(ctrl, node_id, node_id_len, dof_id, value)
      integer(kind=kint)                    :: ctrl
      character(len=HECMW_NAME_LEN)        :: node_id(:)  !Node group name
      integer(kind=kint), pointer          :: dof_id(:)
      integer(kind=kint)                    :: node_id_len
      real(kind=kreal), pointer             :: value(:)
      integer(kind=kint)                    :: fstr_ctrl_get_fload !return value
      character(len=HECMW_NAME_LEN)        :: data_fmt, ss
 
      write(ss,*) node_id_len
      write(data_fmt, '(a,a,a)') 'S', trim(adjustl(ss)), 'IR '
 
      fstr_ctrl_get_fload = fstr_ctrl_get_data_array_ex(ctrl, data_fmt, node_id, dof_id, value)
    end function
 
  end subroutine
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
  subroutine fstr_setup_facceleration( ctrl, counter, P )
    use hecmw_setup_util, only &
      : hecmw_expand_integer_array &
      , hecmw_expand_real_array
    implicit none
  !---- args
    integer(kind=kint)    :: ctrl
    integer(kind=kint)    :: counter
    type(fstr_param_pack) :: P
  !---- vals
    integer(kind=kint)                 :: rcode
    character(HECMW_NAME_LEN), pointer :: grp_id_name(:)
    integer(kind=kint), pointer        :: dof_ids(:)
    integer(kind=kint), pointer        :: dof_ide(:)
    real(kind=kreal), pointer          :: val_ptr(:)
    integer(kind=kint)                 :: i, n, old_size, new_size
    integer(kind=kint)                 :: loadcase
  !---- body
 
    !only the FrontISTR format is supported
    if( p%SOLID%file_type /= kbcffstr) return
 
    !read loadcase (real=1:default, img=2)
    loadcase = kfloadcase_re
    rcode = fstr_ctrl_get_param_ex( ctrl, 'LOAD CASE ', '# ', 0, 'I', loadcase)
    if( loadcase /= kfloadcase_re .and. loadcase /= kfloadcase_im ) then
      write(*,*)    'fstr control file error : !FACCELERATION : LOAD CASE must be 1 or 2'
      write(ilog,*) 'fstr control file error : !FACCELERATION : LOAD CASE must be 1 or 2'
      stop hecmw_exit_input
    end if
 
    !append the data lines of this header to the entries read so far
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n == 0 ) return
    old_size = p%FREQ%FACC_ngrp_tot
    new_size = old_size + n
 
    !the expanded part is zero, so the part not given by loadcase stays zero
    p%FREQ%FACC_ngrp_tot = new_size
    call hecmw_expand_integer_array( p%FREQ%FACC_ngrp_ID,    old_size, new_size )
    call hecmw_expand_integer_array( p%FREQ%FACC_ngrp_type,  old_size, new_size )
    call hecmw_expand_real_array   ( p%FREQ%FACC_ngrp_valre, old_size, new_size )
    call hecmw_expand_real_array   ( p%FREQ%FACC_ngrp_valim, old_size, new_size )
 
    !read the data lines into temporary arrays
    allocate( grp_id_name(n) )
    allocate( dof_ids(n) )
    allocate( dof_ide(n) )
    allocate( val_ptr(n) )
    val_ptr = 0.0d0
    rcode = fstr_ctrl_get_facceleration( ctrl, grp_id_name, hecmw_name_len, dof_ids, dof_ide, val_ptr )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
    !node ID or node group name -> node group ID
    call node_grp_name_to_id_ex( p%MESH, '!FACCELERATION', n, grp_id_name, p%FREQ%FACC_ngrp_ID(old_size+1:) )
 
    do i = 1, n
      !check the range of the DOFs as !ACCELERATION does
      if( (dof_ids(i) < 1).or.(6 < dof_ids(i)).or.(dof_ide(i) < dof_ids(i)).or.(6 < dof_ide(i)) ) then
        write(*,*)    'fstr control file error : !FACCELERATION : range of dof_ids and dof_ide is from 1 to 6'
        write(ilog,*) 'fstr control file error : !FACCELERATION : range of dof_ids and dof_ide is from 1 to 6'
        stop hecmw_exit_input
      end if
      !encode the DOF range as !ACCELERATION does
      p%FREQ%FACC_ngrp_type(old_size+i) = 10 * dof_ids(i) + dof_ide(i)
      !store the amplitude in the part selected by loadcase
      if( loadcase == kfloadcase_re ) then
        p%FREQ%FACC_ngrp_valre(old_size+i) = val_ptr(i)
      else
        p%FREQ%FACC_ngrp_valim(old_size+i) = val_ptr(i)
      end if
    end do
 
    deallocate( grp_id_name )
    deallocate( dof_ids )
    deallocate( dof_ide )
    deallocate( val_ptr )
    return
 
    contains
 
    !read all data lines in the format of NODE_ID, DOF_idS, DOF_idE, Value
    function fstr_ctrl_get_facceleration(ctrl, node_id, node_id_len, dof_ids, dof_ide, value)
      integer(kind=kint)            :: ctrl
      character(len=HECMW_NAME_LEN) :: node_id(:)
      integer(kind=kint)            :: node_id_len
      integer(kind=kint), pointer   :: dof_ids(:)
      integer(kind=kint), pointer   :: dof_ide(:)
      real(kind=kreal), pointer     :: value(:)
      integer(kind=kint)            :: fstr_ctrl_get_facceleration
      character(len=HECMW_NAME_LEN) :: data_fmt, ss
 
      !format: string of node_id_len characters, two integers and a real
      write(ss,*) node_id_len
      write(data_fmt, '(a,a,a)') 'S', trim(adjustl(ss)), 'IIr '
 
      fstr_ctrl_get_facceleration = fstr_ctrl_get_data_array_ex(ctrl, data_fmt, node_id, dof_ids, dof_ide, value)
    end function
 
  end subroutine
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
  subroutine fstr_setup_eigenread( ctrl, counter, P )
  !---- args
    integer(kind=kint)    :: ctrl
    integer(kind=kint)    :: counter
    type(fstr_param_pack) :: P
  !---- vals
    integer(kind=kint)                :: filename_len
    character(len=HECMW_NAME_LEN) :: datafmt, ss
  !---- body
 
    filename_len = hecmw_filename_len
    write(ss,*) filename_len
    write(datafmt, '(a,a,a)') 'F', trim(adjustl(ss)), ' '
 
    if( fstr_ctrl_get_data_ex( ctrl, 1, datafmt, p%FREQ%eigenlog_filename ) /= 0) return
    if( fstr_ctrl_get_data_ex( ctrl, 2, 'ii ', p%FREQ%start_mode, p%FREQ%end_mode ) /= 0) return
 
    return
 
  end subroutine
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_expand_dload_array( array, old_size, new_size )
    implicit none
    real(kind=kreal), pointer :: array(:,:)
    integer(kind=kint) :: old_size, new_size, i, j
    real(kind=kreal), pointer :: temp(:,:)
 
    if( old_size >= new_size ) then
      return
    end if
 
    if( associated( array ) ) then
      allocate(temp(0:6, old_size))
      temp = array
      deallocate(array)
      allocate(array(0:6, new_size))
      array = 0
      do i=1,old_size
        do j=0,6
          array(j,i) = temp(j,i)
        end do
      end do
      deallocate(temp)
    else
      allocate(array(0:6, new_size))
      array = 0
    end if
  end subroutine fstr_expand_dload_array
 
  subroutine fstr_setup_dload( ctrl, counter, P )
    use hecmw_setup_util, only &
      : hecmw_expand_integer_array
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
    character(HECMW_NAME_LEN) :: amp
    integer(kind=kint) :: amp_id
    integer(kind=kint) :: follow
    character(HECMW_NAME_LEN), pointer :: grp_id_name(:)
    real(kind=kreal),pointer :: new_params(:,:)
    logical,pointer :: fg_surface(:)
    integer(kind=kint),pointer :: lid_ptr(:)
    integer(kind=kint) :: i, j, n, old_size, new_size
    integer(kind=kint) :: gid
 
    if( p%SOLID%file_type /= kbcffstr ) return
 
    gid = 1
    rcode = fstr_ctrl_get_param_ex( ctrl, 'GRPID ',  '# ',            0, 'I', gid  )
 
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n == 0 ) return
    old_size = p%SOLID%DLOAD_ngrp_tot
    new_size = old_size + n
    p%SOLID%DLOAD_ngrp_tot = new_size
    ! Keiji Suemitsu (20140624) <
    call hecmw_expand_integer_array ( p%SOLID%DLOAD_ngrp_GRPID, old_size, new_size )
    call hecmw_expand_integer_array ( p%SOLID%DLOAD_ngrp_ID,  old_size, new_size )
    call hecmw_expand_integer_array ( p%SOLID%DLOAD_ngrp_LID, old_size, new_size )
    call hecmw_expand_integer_array ( p%SOLID%DLOAD_ngrp_amp, old_size, new_size )
    call fstr_expand_dload_array ( p%SOLID%DLOAD_ngrp_params, old_size, new_size )
    ! > Keiji Suemitsu (20140624)
 
    allocate( grp_id_name(n))
    allocate( lid_ptr(n) )
    allocate( new_params(0:6,n))
    allocate( fg_surface(n))
    new_params = 0
    amp = ' '
    follow = p%SOLID%DLOAD_follow
    if( .not. p%PARAM%nlgeom ) follow = 0
    rcode = fstr_ctrl_get_dload( ctrl, amp, follow, &
      grp_id_name, hecmw_name_len,    &
      lid_ptr, new_params )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
    call amp_name_to_id( p%MESH, '!DLOAD', amp, amp_id )
    p%SOLID%DLOAD_follow = follow
    do i=1,n
      p%SOLID%DLOAD_ngrp_amp(old_size+i) = amp_id
      p%SOLID%DLOAD_ngrp_LID(old_size+i) = lid_ptr(i)
      do j=0, 6
        p%SOLID%DLOAD_ngrp_params(j,old_size+i) = new_params(j,i)
      end do
      fg_surface(i) =  ( lid_ptr(i) == 100 )
    end do
    p%SOLID%DLOAD_ngrp_GRPID(old_size+1:new_size) = gid
    call dload_grp_name_to_id_ex( p%MESH, n, grp_id_name, fg_surface, p%SOLID%DLOAD_ngrp_ID(old_size+1:))
    deallocate( grp_id_name )
    deallocate( lid_ptr )
    deallocate( new_params )
    deallocate( fg_surface )
    nullify( grp_id_name )
    nullify( lid_ptr )
    nullify( new_params )
    nullify( fg_surface )
  end subroutine fstr_setup_dload
 
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_temperature( ctrl, counter, P )
    use hecmw_setup_util, only &
      : hecmw_expand_integer_array &
      , hecmw_expand_real_array
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode, gid
    character(HECMW_NAME_LEN), pointer :: grp_id_name(:)
    real(kind=kreal),pointer :: val_ptr(:)
    integer(kind=kint) :: i, n, old_size, new_size
 
    if( p%SOLID%file_type /= kbcffstr ) return
 
    gid = 1
    rcode = fstr_ctrl_get_param_ex( ctrl, 'GRPID ',  '# ',            0, 'I', gid  )
 
    n = fstr_ctrl_get_data_line_n( ctrl )
    old_size = p%SOLID%TEMP_ngrp_tot
    if( n > 0 ) then
      new_size = old_size + n
    else
      new_size = old_size + 1
    endif
    call hecmw_expand_integer_array ( p%SOLID%TEMP_ngrp_GRPID, old_size, new_size )
    call hecmw_expand_integer_array ( p%SOLID%TEMP_ngrp_ID, old_size, new_size )
    call hecmw_expand_real_array    ( p%SOLID%TEMP_ngrp_val,old_size, new_size )
 
    allocate( grp_id_name(n))
    allocate( val_ptr(n) )
    val_ptr = 0.0d0
 
    rcode = fstr_ctrl_get_temperature( ctrl,      &
      p%SOLID%TEMP_irres,           &
      p%SOLID%TEMP_tstep,           &
      p%SOLID%TEMP_interval,        &
      p%SOLID%TEMP_rtype,           &
      grp_id_name, hecmw_name_len,  &
      val_ptr )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
    do i = 1, n
      p%SOLID%TEMP_ngrp_val(old_size+i) = val_ptr(i)
    enddo
    deallocate( val_ptr )
    nullify( val_ptr )
 
    p%SOLID%TEMP_ngrp_GRPID(old_size+1:new_size) = gid
    if( n > 0 ) then
      if( p%SOLID%TEMP_irres == 0 ) then
        p%SOLID%TEMP_ngrp_tot = new_size
        call node_grp_name_to_id_ex( p%MESH, '!TEMPERATURE', &
          n, grp_id_name, p%SOLID%TEMP_ngrp_ID(old_size+1:))
      endif
      deallocate( grp_id_name )
    endif
 
  end subroutine fstr_setup_temperature
 
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_spring( ctrl, counter, P )
    use hecmw_setup_util, only &
      : hecmw_expand_integer_array &
      , hecmw_expand_real_array
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode, incremental
    character(HECMW_NAME_LEN) :: amp
    integer(kind=kint) :: amp_id
    character(HECMW_NAME_LEN), pointer :: grp_id_name(:)
    real(kind=kreal),pointer :: val_ptr(:)
    integer(kind=kint),pointer :: id_ptr(:)
    integer(kind=kint) :: i, n, old_size, new_size
    integer(kind=kint) :: gid
 
    if( p%SOLID%file_type /= kbcffstr ) return
    gid = 1
    rcode = fstr_ctrl_get_param_ex( ctrl, 'GRPID ',  '# ',            0, 'I', gid  )
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n == 0 ) return
    old_size = p%SOLID%SPRING_ngrp_tot
    new_size = old_size + n
    p%SOLID%SPRING_ngrp_tot = new_size
    call hecmw_expand_integer_array ( p%SOLID%SPRING_ngrp_GRPID, old_size, new_size )
    call hecmw_expand_integer_array ( p%SOLID%SPRING_ngrp_ID,  old_size, new_size )
    call hecmw_expand_integer_array ( p%SOLID%SPRING_ngrp_DOF, old_size, new_size )
    call hecmw_expand_real_array    ( p%SOLID%SPRING_ngrp_val, old_size, new_size )
    call hecmw_expand_integer_array ( p%SOLID%SPRING_ngrp_amp, old_size, new_size )
    call hecmw_expand_integer_array ( p%SOLID%SPRING_incremental, old_size, new_size )
 
    allocate( grp_id_name(n))
    allocate( id_ptr(n) )
    allocate( val_ptr(n) )
    amp = ' '
    incremental = 0
    id_ptr = 0
    val_ptr = 0.0d0
    rcode = fstr_ctrl_get_spring( ctrl, amp, grp_id_name, hecmw_name_len, id_ptr, val_ptr, incremental )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
    call amp_name_to_id( p%MESH, '!SPRING', amp, amp_id )
    do i=1,n
      p%SOLID%SPRING_ngrp_amp(old_size+i) = amp_id
      p%SOLID%SPRING_ngrp_DOF(old_size+i) = id_ptr(i)
      p%SOLID%SPRING_ngrp_val(old_size+i) = val_ptr(i)
    end do
    p%SOLID%SPRING_ngrp_GRPID(old_size+1:new_size) = gid
    call node_grp_name_to_id_ex( p%MESH, '!SPRING', n, grp_id_name, p%SOLID%SPRING_ngrp_ID(old_size+1:))
    p%SOLID%SPRING_incremental(old_size+1:new_size) = incremental
 
    deallocate( grp_id_name )
    deallocate( id_ptr )
    deallocate( val_ptr )
    nullify( grp_id_name )
    nullify( id_ptr )
    nullify( val_ptr )
 
  end subroutine fstr_setup_spring
 
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_reftemp( ctrl, counter, P )
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
 
    rcode = fstr_ctrl_get_reftemp( ctrl, p%PARAM%ref_temp )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
  end subroutine fstr_setup_reftemp
 
 
  !*****************************************************************************!
  !* HEADERS FOR HEAT ANALYSIS *************************************************!
  !*****************************************************************************!
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_heat( ctrl, counter, P )
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
    integer(kind=kint) :: n
    character(len=HECMW_NAME_LEN) :: mName
    integer(kind=kint) :: i
 
    n = fstr_ctrl_get_data_line_n( ctrl )
 
    if( n == 0 ) return
 
    call reallocate_real( p%PARAM%dtime, n)
    call reallocate_real( p%PARAM%etime, n)
    call reallocate_real( p%PARAM%dtmin, n)
    call reallocate_real( p%PARAM%delmax,n)
    call reallocate_integer( p%PARAM%itmax, n)
    call reallocate_real( p%PARAM%eps,   n)
    p%PARAM%analysis_n = n
 
    p%PARAM%dtime = 0
    p%PARAM%etime = 0
    p%PARAM%dtmin = 0
    p%PARAM%delmax = 0
    p%PARAM%itmax = 20
    p%PARAM%eps = 1.0e-6
    p%PARAM%timepoint_id = 0
 
    rcode = fstr_ctrl_get_heat(   ctrl,        &
      p%PARAM%dtime,     &
      p%PARAM%etime,     &
      p%PARAM%dtmin,     &
      p%PARAM%delmax,    &
      p%PARAM%itmax,     &
      p%PARAM%eps,       &
      mname,             &
      p%HEAT%beta)
    if( rcode /= 0 ) then
      call fstr_ctrl_err_stop
    end if
 
    if( associated(p%PARAM%timepoints) ) then
      do i=1,size(p%PARAM%timepoints)
        if( hecmw_streqr( p%PARAM%timepoints(i)%name, mname ) ) then
          p%PARAM%timepoint_id = i; exit
        endif
      enddo
    endif
 
    call reallocate_real( p%HEAT%STEP_DLTIME, n)
    call reallocate_real( p%HEAT%STEP_EETIME, n)
    call reallocate_real( p%HEAT%STEP_DELMIN, n)
    call reallocate_real( p%HEAT%STEP_DELMAX, n)
    p%HEAT%STEPtot = n
 
    p%HEAT%STEP_DLTIME = p%PARAM%dtime
    p%HEAT%STEP_EETIME = p%PARAM%etime
    p%HEAT%STEP_DELMIN = p%PARAM%dtmin
    p%HEAT%STEP_DELMAX = p%PARAM%delmax
    p%HEAT%timepoint_id = p%PARAM%timepoint_id
 
  end subroutine fstr_setup_heat
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_fixtemp( ctrl, counter, P )
    use hecmw_setup_util, only &
      : get_grp_member &
      , get_grp_member_n &
      , hecmw_expand_integer_array &
      , hecmw_expand_real_array
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack),target :: P
 
    integer(kind=kint) :: rcode
    character(HECMW_NAME_LEN) :: amp
    integer(kind=kint) :: amp_id
    character(HECMW_NAME_LEN), pointer :: grp_id_name(:)
    real(kind=kreal),pointer :: value(:)
    integer(kind=kint) :: i, j, n, m, head, id, member_n, old_size, new_size
    integer(kind=kint),pointer :: member(:)
    integer(kind=kint) :: local_id, rtc
    ! ------------------------------------------------
 
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n == 0 ) return
 
    allocate( grp_id_name(n))
    allocate( value(n))
 
    amp = ' '
    rcode = fstr_ctrl_get_fixtemp( ctrl, amp, &
      grp_id_name, hecmw_name_len, value )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
    call amp_name_to_id( p%MESH, '!FIXTEMP', amp, amp_id )
 
    m = 0
    do i = 1, n
      !rtc = get_local_member_index( P%MESH, 'node', grp_id_name(i), local_id )
      rtc = get_sorted_local_member_index( p%MESH, p%PARAM, 'node', grp_id_name(i), local_id )
      if( rtc > 0 ) then
        m = m + 1
      else if( rtc < 0 ) then
        m = m + get_grp_member_n( p%MESH, 'node_grp', grp_id_name(i) )
      end if
    end do
 
    if (m == 0) then
      deallocate( grp_id_name )
      deallocate( value )
      return
    endif
 
    ! JP-8
    old_size = p%HEAT%T_FIX_tot
    new_size = old_size + m
    call hecmw_expand_integer_array( p%HEAT%T_FIX_node, old_size, new_size )
    call hecmw_expand_integer_array( p%HEAT%T_FIX_ampl, old_size, new_size )
    call hecmw_expand_real_array(    p%HEAT%T_FIX_val,  old_size, new_size )
    p%HEAT%T_FIX_tot = new_size
 
    head = old_size + 1
    member => p%HEAT%T_FIX_node(head:)
    id = head
    do i = 1, n
      !rtc = get_local_member_index( P%MESH, 'node', grp_id_name(i), local_id )
      rtc = get_sorted_local_member_index( p%MESH, p%PARAM, 'node', grp_id_name(i), local_id )
      if( rtc > 0 ) then
        member(1) = local_id
        member_n = 1
      else if( rtc < 0 ) then
        member_n = get_grp_member( p%MESH, 'node_grp', grp_id_name(i), member )
      else
        cycle
      end if
      if( i<n ) then
        member => member( member_n+1 : )
      endif
      do j = 1, member_n
        p%HEAT%T_FIX_val  (id) = value(i)
        p%HEAT%T_FIX_ampl (id) = amp_id
        id = id + 1
      end do
    end do
 
    deallocate( grp_id_name )
    deallocate( value )
  end subroutine fstr_setup_fixtemp
 
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_cflux( ctrl, counter, P )
    use hecmw_setup_util, only &
      : get_grp_member &
      , get_grp_member_n &
      , hecmw_expand_integer_array &
      , hecmw_expand_real_array
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
    character(HECMW_NAME_LEN) :: amp
    integer(kind=kint) :: amp_id
    character(HECMW_NAME_LEN), pointer :: grp_id_name(:)
    real(kind=kreal),pointer :: value(:)
    integer(kind=kint) :: i, j, n, m, head, id, member_n, old_size, new_size
    integer(kind=kint),pointer :: member(:)
    integer(kind=kint) :: local_id, rtc
    ! ------------------------------------------------
 
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n == 0 ) return
 
    allocate( grp_id_name(n))
    allocate( value(n))
 
    amp = ' '
    rcode = fstr_ctrl_get_cflux( ctrl, amp, &
      grp_id_name, hecmw_name_len, value )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
    call amp_name_to_id( p%MESH, '!CFLUX', amp, amp_id )
 
    m = 0
 
    do i = 1, n
      rtc = get_local_member_index( p%MESH, 'node', grp_id_name(i), local_id )
      if( rtc > 0 ) then
        m = m + 1
      else if( rtc < 0 ) then
        m = m + get_grp_member_n( p%MESH, 'node_grp', grp_id_name(i) )
      end if
    end do
 
    if (m == 0) then
      deallocate( grp_id_name )
      deallocate( value )
      return
    endif
 
    ! JP-9
    old_size = p%HEAT%Q_NOD_tot
    new_size = old_size + m
    call hecmw_expand_integer_array( p%HEAT%Q_NOD_node, old_size, new_size )
    call hecmw_expand_integer_array( p%HEAT%Q_NOD_ampl, old_size, new_size )
    call hecmw_expand_real_array(    p%HEAT%Q_NOD_val,  old_size, new_size )
    p%HEAT%Q_NOD_tot = new_size
 
    head = old_size + 1
    member => p%HEAT%Q_NOD_node(head:)
    id = head
    do i = 1, n
      rtc = get_local_member_index( p%MESH, 'node', grp_id_name(i), local_id )
      if( rtc > 0 ) then
        member(1) = local_id
        member_n = 1
      else if( rtc < 0 ) then
        member_n = get_grp_member( p%MESH, 'node_grp', grp_id_name(i), member )
      else
        cycle
      end if
      if( i<n ) member => member( member_n+1 : )
      do j = 1, member_n
        p%HEAT%Q_NOD_val  (id) = value(i)
        p%HEAT%Q_NOD_ampl (id) = amp_id
        id = id + 1
      end do
    end do
 
    deallocate( grp_id_name )
    deallocate( value )
  end subroutine fstr_setup_cflux
 
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
 
  subroutine fstr_setup_dflux( ctrl, counter, P )
    use hecmw_setup_util, only &
      : get_grp_member &
      , get_grp_member_n &
      , hecmw_expand_integer_array &
      , hecmw_expand_real_array
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
    character(HECMW_NAME_LEN) :: amp
    integer(kind=kint) :: amp_id
    character(HECMW_NAME_LEN), pointer :: grp_id_name(:)
    integer(kind=kint),pointer :: load_type(:)
    real(kind=kreal),pointer :: value(:)
    integer(kind=kint) :: i, j, n, m, head, id, member_n, old_size, new_size
    integer(kind=kint),pointer :: member(:)
    integer(kind=kint) :: local_id, rtc
    ! ------------------------------------------------
 
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n == 0 ) return
 
    allocate( grp_id_name(n))
    allocate( load_type(n))
    allocate( value(n))
 
    amp = ' '
    rcode = fstr_ctrl_get_dflux( ctrl, amp, &
      grp_id_name, hecmw_name_len, load_type, value )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
    call amp_name_to_id( p%MESH, '!DFLUX', amp, amp_id )
 
    m = 0
    do i = 1, n
      rtc = get_local_member_index( p%MESH, 'element', grp_id_name(i), local_id )
      if( rtc > 0 ) then
        m = m + 1
      else if( rtc < 0 ) then
        m = m + get_grp_member_n( p%MESH, 'elem_grp', grp_id_name(i) )
      end if
    end do
 
    if (m == 0) then
      deallocate( grp_id_name )
      deallocate( load_type )
      deallocate( value )
      return
    endif
 
    ! JP-10
    old_size = p%HEAT%Q_SUF_tot
    new_size = old_size + m
    call hecmw_expand_integer_array( p%HEAT%Q_SUF_elem, old_size, new_size )
    call hecmw_expand_integer_array( p%HEAT%Q_SUF_ampl, old_size, new_size )
    call hecmw_expand_integer_array( p%HEAT%Q_SUF_surf, old_size, new_size )
    call hecmw_expand_real_array(    p%HEAT%Q_SUF_val,  old_size, new_size )
    p%HEAT%Q_SUF_tot = new_size
 
    head = old_size + 1
    member => p%HEAT%Q_SUF_elem(head:)
    id = head
    do i = 1, n
      rtc = get_local_member_index( p%MESH, 'element', grp_id_name(i), local_id )
      if( rtc > 0 ) then
        member(1) = local_id
        member_n = 1
      else if( rtc < 0 ) then
        member_n = get_grp_member( p%MESH, 'elem_grp', grp_id_name(i), member )
      else
        cycle
      end if
      if( i<n ) member => member( member_n+1 : )
      do j = 1, member_n
        p%HEAT%Q_SUF_surf (id) = load_type(i)
        p%HEAT%Q_SUF_val  (id) = value(i)
        p%HEAT%Q_SUF_ampl (id) = amp_id
        id = id + 1
      end do
    end do
 
    deallocate( grp_id_name )
    deallocate( load_type )
    deallocate( value )
  end subroutine fstr_setup_dflux
 
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
 
  subroutine fstr_setup_sflux( ctrl, counter, P )
    use hecmw_setup_util, only &
      : get_grp_member &
      , get_grp_member_n &
      , hecmw_expand_integer_array &
      , hecmw_expand_real_array
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
    character(HECMW_NAME_LEN) :: amp
    integer(kind=kint) :: amp_id
    character(HECMW_NAME_LEN), pointer :: grp_id_name(:)
    real(kind=kreal),pointer :: value(:)
    integer(kind=kint) :: i, j, n, m, head, id, member_n, old_size, new_size
    integer(kind=kint),pointer :: member1(:), member2(:)
    ! ------------------------------------------------
 
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n == 0 ) return
 
    allocate( grp_id_name(n))
    allocate( value(n))
 
    amp = ' '
    rcode = fstr_ctrl_get_sflux( ctrl, amp, &
      grp_id_name, hecmw_name_len, value )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
    call amp_name_to_id( p%MESH, '!SFLUX', amp, amp_id )
 
    m = 0
    do i = 1, n
      m = m + get_grp_member_n( p%MESH, 'surf_grp', grp_id_name(i) )
    end do
 
    if (m == 0) then
      deallocate( grp_id_name )
      deallocate( value )
      return
    endif
 
    ! JP-11
    old_size = p%HEAT%Q_SUF_tot
    new_size = old_size + m
    call hecmw_expand_integer_array( p%HEAT%Q_SUF_elem, old_size, new_size )
    call hecmw_expand_integer_array( p%HEAT%Q_SUF_ampl, old_size, new_size )
    call hecmw_expand_integer_array( p%HEAT%Q_SUF_surf, old_size, new_size )
    call hecmw_expand_real_array(    p%HEAT%Q_SUF_val,  old_size, new_size )
    p%HEAT%Q_SUF_tot = new_size
 
    head = old_size + 1
    member1 => p%HEAT%Q_SUF_elem(head:)
    member2 => p%HEAT%Q_SUF_surf(head:)
    id = head
    do i = 1, n
      member_n = get_grp_member( p%MESH, 'surf_grp', grp_id_name(i), member1, member2 )
      if( i<n ) then
        member1 => member1( member_n+1 : )
        member2 => member2( member_n+1 : )
      end if
      do j = 1, member_n
        p%HEAT%Q_SUF_val  (id) = value(i)
        p%HEAT%Q_SUF_ampl (id) = amp_id
        id = id + 1
      end do
    end do
 
    deallocate( grp_id_name )
    deallocate( value )
  end subroutine fstr_setup_sflux
 
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
 
  subroutine fstr_setup_film( ctrl, counter, P )
    use hecmw_setup_util, only &
      : get_grp_member &
      , get_grp_member_n &
      , hecmw_expand_integer_array &
      , hecmw_expand_integer_array2 &
      , hecmw_expand_real_array2
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
    character(HECMW_NAME_LEN) :: amp1, amp2
    integer(kind=kint) :: amp_id1, amp_id2
    character(HECMW_NAME_LEN), pointer :: grp_id_name(:)
    integer(kind=kint),pointer :: load_type(:)
    real(kind=kreal),pointer :: value(:)
    real(kind=kreal),pointer :: shink(:)
    integer(kind=kint) :: i, j, n, m, head, id, member_n, old_size, new_size
    integer(kind=kint),pointer :: member(:)
    integer(kind=kint) :: local_id, rtc
    ! ------------------------------------------------
 
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n == 0 ) return
 
    allocate( grp_id_name(n))
    allocate( load_type(n))
    allocate( value(n))
    allocate( shink(n))
 
    amp1 = ' '
    amp2 = ' '
 
    rcode = fstr_ctrl_get_film( ctrl, amp1, amp2, &
      grp_id_name, hecmw_name_len, load_type, value, shink )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
    call amp_name_to_id( p%MESH, '!FILM', amp1, amp_id1 )
    call amp_name_to_id( p%MESH, '!FILM', amp2, amp_id2 )
 
    m = 0
    do i = 1, n
      rtc = get_local_member_index( p%MESH, 'element', grp_id_name(i), local_id )
      if( rtc > 0 ) then
        m = m + 1
      else if( rtc < 0 ) then
        m = m + get_grp_member_n( p%MESH, 'elem_grp', grp_id_name(i) )
      end if
    end do
 
    if (m == 0) then
      deallocate( grp_id_name )
      deallocate( load_type )
      deallocate( value )
      deallocate( shink )
      return
    endif
 
    ! JP-12
    old_size = p%HEAT%H_SUF_tot
    new_size = old_size + m
    call hecmw_expand_integer_array(  p%HEAT%H_SUF_elem,    old_size, new_size )
    call hecmw_expand_integer_array2( p%HEAT%H_SUF_ampl, 2, old_size, new_size )
    call hecmw_expand_integer_array(  p%HEAT%H_SUF_surf,    old_size, new_size )
    call hecmw_expand_real_array2(    p%HEAT%H_SUF_val,  2, old_size, new_size )
    p%HEAT%H_SUF_tot = new_size
 
    head = old_size + 1
    member => p%HEAT%H_SUF_elem(head:)
    id = head
    do i = 1, n
      rtc = get_local_member_index( p%MESH, 'element', grp_id_name(i), local_id )
      if( rtc > 0 ) then
        member(1) = local_id
        member_n = 1
      else if( rtc < 0 ) then
        member_n = get_grp_member( p%MESH, 'elem_grp', grp_id_name(i), member )
      else
        cycle
      end if
      if( i<n ) member => member( member_n+1 : )
      do j = 1, member_n
        p%HEAT%H_SUF_surf (id)   = load_type(i)
        p%HEAT%H_SUF_val  (id,1) = value(i)
        p%HEAT%H_SUF_val  (id,2) = shink(i)
        p%HEAT%H_SUF_ampl (id,1) = amp_id1
        p%HEAT%H_SUF_ampl (id,2) = amp_id2
        id= id + 1
      end do
    end do
 
    deallocate( grp_id_name )
    deallocate( load_type )
    deallocate( value )
    deallocate( shink )
  end subroutine fstr_setup_film
 
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
 
  subroutine fstr_setup_sfilm( ctrl, counter, P )
    use hecmw_setup_util, only &
      : get_grp_member &
      , get_grp_member_n &
      , hecmw_expand_integer_array &
      , hecmw_expand_integer_array2 &
      , hecmw_expand_real_array2
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
    character(HECMW_NAME_LEN) :: amp1, amp2
    integer(kind=kint) :: amp_id1, amp_id2
    character(HECMW_NAME_LEN), pointer :: grp_id_name(:)
    real(kind=kreal),pointer :: value(:)
    real(kind=kreal),pointer :: shink(:)
    integer(kind=kint) :: i, j, n, m, head, id, member_n, old_size, new_size
    integer(kind=kint),pointer :: member1(:), member2(:)
    ! ------------------------------------------------
 
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n == 0 ) return
 
    allocate( grp_id_name(n))
    allocate( value(n))
    allocate( shink(n))
 
    amp1 = ' '
    amp2 = ' '
    rcode = fstr_ctrl_get_sfilm( ctrl, amp1, amp2, &
      grp_id_name, hecmw_name_len, value, shink )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
    call amp_name_to_id( p%MESH, '!SFILM', amp1, amp_id1 )
    call amp_name_to_id( p%MESH, '!SFILM', amp2, amp_id2 )
 
    m = 0
    do i = 1, n
      m = m + get_grp_member_n( p%MESH, 'surf_grp', grp_id_name(i) )
    end do
 
    if (m == 0) then
      deallocate( grp_id_name )
      deallocate( value )
      deallocate( shink )
      return
    endif
 
    ! JP-13
    old_size = p%HEAT%H_SUF_tot
    new_size = old_size + m
    call hecmw_expand_integer_array(  p%HEAT%H_SUF_elem,    old_size, new_size )
    call hecmw_expand_integer_array2( p%HEAT%H_SUF_ampl, 2, old_size, new_size )
    call hecmw_expand_integer_array(  p%HEAT%H_SUF_surf,    old_size, new_size )
    call hecmw_expand_real_array2(    p%HEAT%H_SUF_val,  2, old_size, new_size )
    p%HEAT%H_SUF_tot = new_size
 
    head = old_size + 1
    member1 => p%HEAT%H_SUF_elem(head:)
    member2 => p%HEAT%H_SUF_surf(head:)
    id = head
    do i = 1, n
      member_n = get_grp_member( p%MESH, 'surf_grp', grp_id_name(i), member1, member2 )
      if( i<n ) then
        member1 => member1( member_n+1 : )
        member2 => member2( member_n+1 : )
      end if
      do j = 1, member_n
        p%HEAT%H_SUF_val  (id,1) = value(i)
        p%HEAT%H_SUF_val  (id,2) = shink(i)
        p%HEAT%H_SUF_ampl (id,1) = amp_id1
        p%HEAT%H_SUF_ampl (id,2) = amp_id2
        id = id + 1
      end do
    end do
 
    deallocate( grp_id_name )
    deallocate( value )
    deallocate( shink )
  end subroutine fstr_setup_sfilm
 
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
 
  subroutine fstr_setup_radiate( ctrl, counter, P )
    use hecmw_setup_util, only &
      : get_grp_member &
      , get_grp_member_n &
      , hecmw_expand_integer_array &
      , hecmw_expand_integer_array2 &
      , hecmw_expand_real_array2
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
    character(HECMW_NAME_LEN) :: amp1, amp2
    integer(kind=kint) :: amp_id1, amp_id2
    character(HECMW_NAME_LEN), pointer :: grp_id_name(:)
    integer(kind=kint),pointer :: load_type(:)
    real(kind=kreal),pointer :: value(:)
    real(kind=kreal),pointer :: shink(:)
    integer(kind=kint) :: i, j, n, m, head, id, member_n, old_size, new_size
    integer(kind=kint),pointer :: member(:)
    integer(kind=kint) :: local_id, rtc
    ! ------------------------------------------------
 
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n == 0 ) return
 
    allocate( grp_id_name(n))
    allocate( load_type(n))
    allocate( value(n))
    allocate( shink(n))
 
    amp1 = ' '
    amp2 = ' '
    rcode = fstr_ctrl_get_radiate( ctrl, amp1, amp2, &
      grp_id_name, hecmw_name_len, load_type, value, shink )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
    call amp_name_to_id( p%MESH, '!RADIATE', amp1, amp_id1 )
    call amp_name_to_id( p%MESH, '!RADIATE', amp2, amp_id2 )
 
    m = 0
    do i = 1, n
      rtc = get_local_member_index( p%MESH, 'element', grp_id_name(i), local_id )
      if( rtc > 0 ) then
        m = m + 1
      else if( rtc < 0 ) then
        m = m + get_grp_member_n( p%MESH, 'elem_grp', grp_id_name(i) )
      end if
    end do
 
    if (m == 0) then
      deallocate( grp_id_name )
      deallocate( load_type )
      deallocate( value )
      deallocate( shink )
      return
    endif
 
    ! JP-14
    old_size = p%HEAT%R_SUF_tot
    new_size = old_size + m
    call hecmw_expand_integer_array(  p%HEAT%R_SUF_elem,    old_size, new_size )
    call hecmw_expand_integer_array2( p%HEAT%R_SUF_ampl, 2, old_size, new_size )
    call hecmw_expand_integer_array(  p%HEAT%R_SUF_surf,    old_size, new_size )
    call hecmw_expand_real_array2(    p%HEAT%R_SUF_val,  2, old_size, new_size )
    p%HEAT%R_SUF_tot = new_size
 
    head = old_size + 1
    member => p%HEAT%R_SUF_elem(head:)
    id = head
    do i = 1, n
      rtc = get_local_member_index( p%MESH, 'element', grp_id_name(i), local_id )
      if( rtc > 0 ) then
        member(1) = local_id
        member_n = 1
      else if( rtc < 0 ) then
        member_n = get_grp_member( p%MESH, 'elem_grp', grp_id_name(i), member )
      else
        cycle
      end if
      if( i<n ) member => member( member_n+1 : )
      do j = 1, member_n
        p%HEAT%R_SUF_surf (id)   = load_type(i)
        p%HEAT%R_SUF_val  (id,1) = value(i)
        p%HEAT%R_SUF_val  (id,2) = shink(i)
        p%HEAT%R_SUF_ampl (id,1) = amp_id1
        p%HEAT%R_SUF_ampl (id,2) = amp_id2
        id = id + 1
      end do
    end do
 
    deallocate( grp_id_name )
    deallocate( load_type )
    deallocate( value )
    deallocate( shink )
  end subroutine fstr_setup_radiate
 
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
 
  subroutine fstr_setup_sradiate( ctrl, counter, P )
    use hecmw_setup_util, only &
      : get_grp_member &
      , get_grp_member_n &
      , hecmw_expand_integer_array &
      , hecmw_expand_integer_array2 &
      , hecmw_expand_real_array2
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
    character(HECMW_NAME_LEN) :: amp1, amp2
    integer(kind=kint) :: amp_id1, amp_id2
    character(HECMW_NAME_LEN), pointer :: grp_id_name(:)
    real(kind=kreal),pointer :: value(:)
    real(kind=kreal),pointer :: shink(:)
    integer(kind=kint) :: i, j, n, m, head, id, member_n, old_size, new_size
    integer(kind=kint),pointer :: member1(:), member2(:)
    ! ------------------------------------------------
 
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n == 0 ) return
 
    allocate( grp_id_name(n))
    allocate( value(n))
    allocate( shink(n))
 
    amp1 = ' '
    amp2 = ' '
    rcode = fstr_ctrl_get_sradiate( ctrl, amp1, amp2, grp_id_name, hecmw_name_len, value, shink )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
    call amp_name_to_id( p%MESH, '!SRADIATE', amp1, amp_id1 )
    call amp_name_to_id( p%MESH, '!SRADIATE', amp2, amp_id2 )
 
    m = 0
    do i = 1, n
      m = m + get_grp_member_n( p%MESH, 'surf_grp', grp_id_name(i) )
    end do
 
    if (m == 0) then
      deallocate( grp_id_name )
      deallocate( value )
      deallocate( shink )
      return
    endif
 
    ! JP-15
    old_size = p%HEAT%R_SUF_tot
    new_size = old_size + m
    call hecmw_expand_integer_array(  p%HEAT%R_SUF_elem,    old_size, new_size )
    call hecmw_expand_integer_array2( p%HEAT%R_SUF_ampl, 2, old_size, new_size )
    call hecmw_expand_integer_array(  p%HEAT%R_SUF_surf,    old_size, new_size )
    call hecmw_expand_real_array2(    p%HEAT%R_SUF_val,  2, old_size, new_size )
    p%HEAT%R_SUF_tot = new_size
 
    head = old_size + 1
    member1 => p%HEAT%R_SUF_elem(head:)
    member2 => p%HEAT%R_SUF_surf(head:)
    id = head
    do i = 1, n
      member_n = get_grp_member( p%MESH, 'surf_grp', grp_id_name(i), member1, member2 )
      if( i<n ) then
        member1 => member1( member_n+1 : )
        member2 => member2( member_n+1 : )
      end if
      do j = 1, member_n
        p%HEAT%R_SUF_val  (id,1) = value(i)
        p%HEAT%R_SUF_val  (id,2) = shink(i)
        p%HEAT%R_SUF_ampl (id,1) = amp_id1
        p%HEAT%R_SUF_ampl (id,2) = amp_id2
        id = id + 1
      end do
    end do
 
    deallocate( grp_id_name )
    deallocate( value )
    deallocate( shink )
  end subroutine fstr_setup_sradiate
 
 
  !*****************************************************************************!
  !* HEADERS FOR EIGEN ANALYSIS ************************************************!
  !*****************************************************************************!
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_eigen( ctrl, counter, P )
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
 
    rcode = fstr_ctrl_get_eigen( ctrl, p%EIGEN%nget, p%EIGEN%tolerance, p%EIGEN%maxiter, p%EIGEN%sigma)
    if( rcode /= 0) call fstr_ctrl_err_stop
 
  end subroutine fstr_setup_eigen
 
 
  !*****************************************************************************!
  !* HEADERS FOR DYNAMIC ANALYSIS **********************************************!
  !*****************************************************************************!
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_dynamic( ctrl, counter, P )
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
    integer(kind=kint) :: rcode
    character(HECMW_NAME_LEN) :: grp_id_name(1)
    integer(kind=kint) :: grp_id(1)
 
    rcode = fstr_ctrl_get_dynamic( ctrl, &
      p%PARAM%nlgeom,  &
      p%DYN%idx_eqa, &
      p%DYN%idx_resp,&
      p%DYN%n_step,  &
      p%DYN%t_start, &
      p%DYN%t_end,   &
      p%DYN%t_delta, &
      p%DYN%gamma,   &
      p%DYN%beta,    &
      p%DYN%idx_mas, &
      p%DYN%idx_dmp, &
      p%DYN%ray_m,   &
      p%DYN%ray_k,   &
      p%DYN%nout,    &
      grp_id_name(1), hecmw_name_len,  &
      p%DYN%nout_monit,  &
      p%DYN%iout_list )
 
    if( rcode /= 0) call fstr_ctrl_err_stop
 
    if( p%DYN%idx_mas /= kmasslumped .and. p%DYN%idx_mas /= kmassconsistent ) then
      write(*,*) '### Error: !DYNAMIC mass matrix type must be 1 (lumped) or 2 (consistent)'
      write(ilog,*) '### Error: !DYNAMIC mass matrix type must be 1 (lumped) or 2 (consistent)'
      call flush(6)
      call flush(ilog)
      call fstr_ctrl_err_stop
    endif
    if( p%DYN%idx_mas == kmassconsistent .and. &
        (p%DYN%idx_eqa /= 1 .or. p%DYN%idx_resp /= 1) ) then
      write(*,*) '### Error: consistent mass matrix is available only for implicit transient analysis'
      write(ilog,*) '### Error: consistent mass matrix is available only for implicit transient analysis'
      call flush(6)
      call flush(ilog)
      call fstr_ctrl_err_stop
    endif
 
    if (p%DYN%idx_resp == 1) then
      call node_grp_name_to_id_ex( p%MESH, '!DYNAMIC', 1, grp_id_name, grp_id)
      p%DYN%ngrp_monit = grp_id(1)
    else
      read(grp_id_name,*) p%DYN%ngrp_monit
    endif
 
  end subroutine fstr_setup_dynamic
 
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_velocity( ctrl, counter, P )
    use hecmw_setup_util, only: hecmw_expand_integer_array, hecmw_expand_real_array
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
    integer(kind=kint) :: vType
    character(HECMW_NAME_LEN) :: amp, rotc_name(1)
    integer(kind=kint) :: amp_id, rotc_id(1), n_rotc
    character(HECMW_NAME_LEN), pointer :: grp_id_name(:)
    integer(kind=kint),pointer :: dof_ids (:)
    integer(kind=kint),pointer :: dof_ide (:)
    real(kind=kreal),pointer :: val_ptr(:)
    integer(kind=kint) :: i, j, n, old_size, new_size
    integer(kind=kint) :: gid
 
    gid = 1
    rcode = fstr_ctrl_get_param_ex( ctrl, 'GRPID ',  '# ',  0, 'I', gid  )
 
    ! get center of rotation ( same syntax as !BOUNDARY,ROT_CENTER= )
    rotc_name = ' '
    rotc_id = -1
    n_rotc = -1
    rcode = fstr_ctrl_get_param_ex( ctrl, 'ROT_CENTER ', '# ', 0, 'S', rotc_name )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
    if( rotc_name(1) /= ' ' ) then
      p%SOLID%VELOCITY_ngrp_rot = p%SOLID%VELOCITY_ngrp_rot + 1
      n_rotc = p%SOLID%VELOCITY_ngrp_rot
      call node_grp_name_to_id_ex( p%MESH, '!VELOCITY,ROT_CENTER=', 1, rotc_name, rotc_id)
    endif
 
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n == 0 ) return
    old_size = p%SOLID%VELOCITY_ngrp_tot
    new_size = old_size + n
    p%SOLID%VELOCITY_ngrp_tot = new_size
 
    call hecmw_expand_integer_array (p%SOLID%VELOCITY_ngrp_GRPID, old_size, new_size )
    call hecmw_expand_integer_array (p%SOLID%VELOCITY_ngrp_ID  , old_size, new_size )
    call hecmw_expand_integer_array (p%SOLID%VELOCITY_ngrp_type, old_size, new_size )
    call hecmw_expand_real_array    (p%SOLID%VELOCITY_ngrp_val , old_size, new_size )
    call hecmw_expand_integer_array (p%SOLID%VELOCITY_ngrp_amp , old_size, new_size )
    call hecmw_expand_integer_array (p%SOLID%VELOCITY_ngrp_rotID, old_size, new_size )
    call hecmw_expand_integer_array (p%SOLID%VELOCITY_ngrp_centerID, old_size, new_size )
 
    allocate( grp_id_name(n))
    allocate( dof_ids(n))
    allocate( dof_ide(n))
    allocate( val_ptr(n) )
 
    amp = ''
    val_ptr = 0.0d0
    rcode = fstr_ctrl_get_velocity( ctrl,  vtype, amp,   &
      grp_id_name, hecmw_name_len,  &
      dof_ids, dof_ide, val_ptr )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
    p%SOLID%VELOCITY_type = vtype
    if( vtype == kbcinitial ) p%DYN%VarInitialize = .true.
    call amp_name_to_id( p%MESH, '!VELOCITY', amp, amp_id )
    call node_grp_name_to_id_ex( p%MESH, '!VELOCITY', &
      n, grp_id_name, p%SOLID%VELOCITY_ngrp_ID(old_size+1:))
    p%SOLID%VELOCITY_ngrp_GRPID(old_size+1:new_size) = gid
 
    ! set up information about rotation ( default value is set if ROT_CENTER is not given.)
    p%SOLID%VELOCITY_ngrp_rotID(old_size+1:) = n_rotc
    p%SOLID%VELOCITY_ngrp_centerID(old_size+1:) = rotc_id(1)
 
    j = old_size+1
    do i = 1, n
      if( (dof_ids(i) < 1).or.(6 < dof_ids(i)).or.(dof_ide(i) < 1).or.(6 < dof_ide(i)) ) then
        write(ilog,*) 'fstr control file error : !VELOCITY : range of dof_ids and dof_ide is from 1 to 6'
        stop hecmw_exit_input
      end if
      p%SOLID%VELOCITY_ngrp_type(j) = 10 * dof_ids(i) + dof_ide(i)
      p%SOLID%VELOCITY_ngrp_amp(j) = amp_id
      p%SOLID%VELOCITY_ngrp_val(old_size+i) = val_ptr(i)
      j = j+1
    end do
 
    deallocate( grp_id_name )
    deallocate( dof_ids )
    deallocate( dof_ide )
    deallocate( val_ptr )
    nullify( grp_id_name )
    nullify( dof_ids )
    nullify( dof_ide )
    nullify( val_ptr )
 
  end subroutine fstr_setup_velocity
 
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_acceleration( ctrl, counter, P )
    use hecmw_setup_util, only: hecmw_expand_integer_array, hecmw_expand_real_array
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
    integer(kind=kint) :: aType
    character(HECMW_NAME_LEN) :: amp
    integer(kind=kint) :: amp_id
    character(HECMW_NAME_LEN), pointer :: grp_id_name(:)
    integer(kind=kint),pointer :: dof_ids (:)
    integer(kind=kint),pointer :: dof_ide (:)
    real(kind=kreal),pointer :: val_ptr(:)
    integer(kind=kint) :: i, j, n, old_size, new_size
    integer(kind=kint) :: gid
 
    gid = 1
    rcode = fstr_ctrl_get_param_ex( ctrl, 'GRPID ',  '# ',  0, 'I', gid  )
 
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n == 0 ) return
    old_size = p%SOLID%ACCELERATION_ngrp_tot
    new_size = old_size + n
    p%SOLID%ACCELERATION_ngrp_tot = new_size
 
    call hecmw_expand_integer_array (p%SOLID%ACCELERATION_ngrp_GRPID, old_size, new_size )
    call hecmw_expand_integer_array (p%SOLID%ACCELERATION_ngrp_ID  , old_size, new_size )
    call hecmw_expand_integer_array (p%SOLID%ACCELERATION_ngrp_type, old_size, new_size )
    call hecmw_expand_real_array    (p%SOLID%ACCELERATION_ngrp_val , old_size, new_size )
    call hecmw_expand_integer_array (p%SOLID%ACCELERATION_ngrp_amp , old_size, new_size )
 
    allocate( grp_id_name(n))
    allocate( dof_ids(n))
    allocate( dof_ide(n))
    allocate( val_ptr(n))
 
    amp = ' '
    val_ptr = 0.0d0
    rcode = fstr_ctrl_get_acceleration( ctrl,  atype, amp,   &
      grp_id_name, hecmw_name_len,  &
      dof_ids, dof_ide,  val_ptr)
    if( rcode /= 0 ) call fstr_ctrl_err_stop
    p%SOLID%ACCELERATION_type = atype
    if( atype == kbcinitial )p%DYN%VarInitialize = .true.
    call amp_name_to_id( p%MESH, '!ACCELERATION', amp, amp_id )
    call node_grp_name_to_id_ex( p%MESH, '!ACCELERATION', &
      n, grp_id_name, p%SOLID%ACCELERATION_ngrp_ID(old_size+1:))
    p%SOLID%ACCELERATION_ngrp_GRPID(old_size+1:new_size) = gid
 
    j = old_size+1
    do i = 1, n
      if( (dof_ids(i) < 1).or.(6 < dof_ids(i)).or.(dof_ide(i) < 1).or.(6 < dof_ide(i)) ) then
        write(ilog,*) 'fstr control file error : !ACCELERATION : range of dof_ids and dof_ide is from 1 to 6'
        stop hecmw_exit_input
      end if
      p%SOLID%ACCELERATION_ngrp_type(j) = 10 * dof_ids(i) + dof_ide(i)
      p%SOLID%ACCELERATION_ngrp_amp(j) = amp_id
      p%SOLID%ACCELERATION_ngrp_val(old_size+i) = val_ptr(i)
      j = j+1
    end do
 
    deallocate( grp_id_name )
    deallocate( dof_ids )
    deallocate( dof_ide )
    deallocate( val_ptr )
    nullify( grp_id_name )
    nullify( dof_ids )
    nullify( dof_ide )
    nullify( val_ptr )
  end subroutine fstr_setup_acceleration
 
 
  !*****************************************************************************!
  !* MPC ***********************************************************************!
  !*****************************************************************************!
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_mpc( ctrl, counter, P )
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack), target :: P
 
    integer(kind=kint) :: rcode
    !        integer(kind=kint) :: type
    !        integer(kind=kint),pointer :: node1_ptr(:)
    !        integer(kind=kint),pointer :: node2_ptr(:)
    !        integer(kind=kint),pointer :: dof_ptr(:)
    !        integer(kind=kint) :: n, old_size, new_size
    !
    !        rcode = fstr_ctrl_get_param_ex( ctrl, 'TYPE ', 'RIGID ', 1, 'P', type )
    !        if( rcode < 0 ) call fstr_ctrl_err_stop
    !
    !        n = fstr_ctrl_get_data_line_n( ctrl )
    !        if( n == 0 ) return
    !        old_size = P%MPC_RD%nmpc
    !        new_size = old_size + n
    !        P%MPC_RD%nmpc = new_size
    !
    !        call hecmw_expand_integer_array ( P%MPC_RD%node1,  old_size, new_size )
    !        call hecmw_expand_integer_array ( P%MPC_RD%node2,  old_size, new_size )
    !        call hecmw_expand_integer_array ( P%MPC_RD%dof,    old_size, new_size )
    !
    !        node1_ptr => P%MPC_RD%node1(old_size+1:)
    !        node2_ptr => P%MPC_RD%node2(old_size+1:)
    !        dof_ptr   => P%MPC_RD%dof(old_size+1:)
    !
    !        rcode = fstr_ctrl_get_MPC( ctrl, type, node1_ptr, node2_ptr, dof_ptr )
    !        if( rcode /= 0 ) call fstr_ctrl_err_stop
    !
    !        if( node_global_to_local( P%MESH, node1_ptr, n ) /= n ) then
    !                call fstr_setup_util_err_stop( '### Error : not exist node (!MPC)' )
    !        endif
    !        if( node_global_to_local( P%MESH, node2_ptr, n ) /= n ) then
    !                call fstr_setup_util_err_stop( '### Error : not exist node (!MPC)' )
    !        endif
 
    !   penalty => svRarray(11)
    rcode = fstr_ctrl_get_mpc( ctrl, svrarray(11))
    if( rcode /= 0) call fstr_ctrl_err_stop
  end subroutine fstr_setup_mpc
 
 
  !*****************************************************************************!
  !* IMPORTING NASTRAN BOUNDARY CONDITIONS *************************************!
  !*****************************************************************************!
 
  subroutine fstr_setup_solid_nastran( ctrl, hecMESH, fstrSOLID )
    implicit none
    integer(kind=kint) :: ctrl
    type (hecmwST_local_mesh) :: hecMESH
    type (fstr_solid        ) :: fstrSOLID
    write(ilog,*) '### Error : In !BOUNDARY, TYPE=NASTRAN is not supported.'
    call fstr_abort( hecmw_exit_model )
  end subroutine fstr_setup_solid_nastran
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_contactalgo( ctrl, P )
    implicit none
    integer(kind=kint) :: ctrl
    !        integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
 
 
    rcode = fstr_ctrl_get_contactalgo( ctrl, p%PARAM%contact_algo, p%PARAM%augiter, &
      p%PARAM%fric_cone_follow, p%PARAM%eps_fric_band )
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
  end subroutine fstr_setup_contactalgo
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_output_sstype( ctrl, P )
    implicit none
    integer(kind=kint) :: ctrl
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode, nid
    character(len=HECMW_NAME_LEN) :: data_fmt
 
    data_fmt = 'SOLUTION,MATERIAL '
    rcode = fstr_ctrl_get_param_ex( ctrl, 'TYPE ', data_fmt, 0, 'P', nid )
    opsstype = nid
    if( rcode /= 0 ) call fstr_ctrl_err_stop
 
  end subroutine fstr_setup_output_sstype
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_convert_contact_type( hecMESH )
    implicit none
    type(hecmwst_local_mesh), pointer :: hecMESH
    integer(kind=kint) :: n, i, sgrp_id, ngrp_id
    ! convert SURF_SURF to NODE_SURF
    n = hecmesh%contact_pair%n_pair
    do i = 1,n
      if( hecmesh%contact_pair%type(i) /= hecmw_contact_type_surf_surf ) cycle
      sgrp_id = hecmesh%contact_pair%slave_grp_id(i)
      call append_node_grp_from_surf_grp( hecmesh, sgrp_id, ngrp_id )
      ! change type of contact and slave group ID
      hecmesh%contact_pair%type(i) = hecmw_contact_type_node_surf
      hecmesh%contact_pair%slave_grp_id(i) = ngrp_id
    enddo
  end subroutine fstr_convert_contact_type
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_convert_rotation_dof( P )
    use hecmw_setup_util, only : append_new_groups, get_grp_id
    implicit none
    type(fstr_param_pack) :: P
    type(hecmwst_local_mesh), pointer :: hecMESH
    integer(kind=kint), allocatable :: rot_node(:), n_rot(:), n_rot_all(:), n_explicit(:), rot_grp(:), src(:)
    integer(kind=kint), allocatable :: new_index(:), new_list(:)
    character(len=HECMW_NAME_LEN), allocatable :: new_name(:)
    logical, allocatable :: is_rot(:), in_all(:), is_explicit(:)
    integer(kind=kint) :: icel, is, ie, nn, i, j, k, m, n, ig, n_new, first_id, ids, ide
 
    hecmesh => p%MESH
    if( hecmesh%n_dof /= 3 ) return
 
    ! explicit dummy nodes are in ALL, and the generated ones are not
    allocate( in_all(hecmesh%n_node) )
    in_all = .false.
    ig = get_grp_id( hecmesh, 'node_grp', 'ALL' )
    if( ig > 0 ) then
      is = hecmesh%node_group%grp_index(ig-1)
      in_all(hecmesh%node_group%grp_item(is+1:hecmesh%node_group%grp_index(ig))) = .true.
    endif
 
    allocate( rot_node(hecmesh%n_node), is_explicit(hecmesh%n_node) )
    rot_node = 0
    is_explicit = .false.
    do icel = 1, hecmesh%n_elem
      if( .not. hecmw_is_etype_33struct(hecmesh%elem_type(icel)) ) cycle
      nn = hecmw_get_max_node(hecmesh%elem_type(icel))/2
      is = hecmesh%elem_node_index(icel-1)
      do j = 1, nn
        rot_node(hecmesh%elem_node_item(is+j)) = hecmesh%elem_node_item(is+nn+j)
        if( in_all(hecmesh%elem_node_item(is+nn+j)) ) is_explicit(hecmesh%elem_node_item(is+j)) = .true.
      enddo
    enddo
    ! a rank may hold an external translational node without its 641/761/781 element,
    ! so the conditions are converted on every rank by the counts over all the ranks
    n = count( rot_node > 0 )
    call hecmw_allreduce_i1( hecmesh, n, hecmw_max )
    if( n == 0 ) return
 
    allocate( n_rot(hecmesh%node_group%n_grp), n_rot_all(hecmesh%node_group%n_grp) )
    allocate( n_explicit(hecmesh%node_group%n_grp) )
    do ig = 1, hecmesh%node_group%n_grp
      is = hecmesh%node_group%grp_index(ig-1)
      n_rot(ig) = count( rot_node(hecmesh%node_group%grp_item(is+1:hecmesh%node_group%grp_index(ig))) > 0 )
      n_explicit(ig) = count( is_explicit(hecmesh%node_group%grp_item(is+1:hecmesh%node_group%grp_index(ig))) )
    enddo
    n_rot_all = n_rot
    call hecmw_allreduce_i( hecmesh, n_rot_all, hecmesh%node_group%n_grp, hecmw_max )
 
    do i = 1, p%SOLID%SPRING_ngrp_tot
      if( p%SOLID%SPRING_ngrp_DOF(i) > 3 .and. n_rot_all(p%SOLID%SPRING_ngrp_ID(i)) > 0 ) &
        call rotation_dof_err_stop( '!SPRING' )
    enddo
    do i = 1, p%SOLID%VELOCITY_ngrp_tot
      if( p%SOLID%VELOCITY_ngrp_rotID(i) > 0 ) cycle
      if( mod(p%SOLID%VELOCITY_ngrp_type(i), 10) > 3 .and. n_rot_all(p%SOLID%VELOCITY_ngrp_ID(i)) > 0 ) &
        call rotation_dof_err_stop( '!VELOCITY' )
    enddo
    do i = 1, p%SOLID%ACCELERATION_ngrp_tot
      if( mod(p%SOLID%ACCELERATION_ngrp_type(i), 10) > 3 .and. n_rot_all(p%SOLID%ACCELERATION_ngrp_ID(i)) > 0 ) &
        call rotation_dof_err_stop( '!ACCELERATION' )
    enddo
    do i = 1, p%FREQ%FLOAD_ngrp_tot
      if( p%FREQ%FLOAD_ngrp_TYPE(i) /= kfloadtype_node ) cycle
      if( p%FREQ%FLOAD_ngrp_DOF(i) > 3 .and. n_rot_all(p%FREQ%FLOAD_ngrp_ID(i)) > 0 ) &
        call rotation_dof_err_stop( '!FLOAD' )
    enddo
    !the rotational DOFs of !FACCELERATION are also not allowed (the end DOF is mod(type,10))
    do i = 1, p%FREQ%FACC_ngrp_tot
      if( mod(p%FREQ%FACC_ngrp_type(i), 10) > 3 .and. n_rot_all(p%FREQ%FACC_ngrp_ID(i)) > 0 ) &
        call rotation_dof_err_stop( '!FACCELERATION' )
    enddo
    do i = 1, hecmesh%mpc%n_mpc
      do j = hecmesh%mpc%mpc_index(i-1)+1, hecmesh%mpc%mpc_index(i)
        if( hecmesh%mpc%mpc_dof(j) > 3 .and. rot_node(hecmesh%mpc%mpc_item(j)) > 0 ) &
          call rotation_dof_err_stop( '!EQUATION' )
      enddo
    enddo
 
    ! node groups whose rotation is moved: rot_grp(ig) = -1, then the id of the dummy node group
    allocate( rot_grp(hecmesh%node_group%n_grp) )
    rot_grp = 0
    do i = 1, p%SOLID%BOUNDARY_ngrp_tot
      if( p%SOLID%BOUNDARY_ngrp_rotID(i) > 0 ) cycle
      ig = p%SOLID%BOUNDARY_ngrp_ID(i)
      if( mod(p%SOLID%BOUNDARY_ngrp_type(i), 10) <= 3 .or. n_rot_all(ig) == 0 ) cycle
      if( n_explicit(ig) > 0 ) call explicit_dummy_node_err_stop( '!BOUNDARY' )
      rot_grp(ig) = -1
    enddo
    do i = 1, p%SOLID%CLOAD_ngrp_tot
      if( p%SOLID%CLOAD_ngrp_rotID(i) > 0 .or. p%SOLID%CLOAD_ngrp_DOF(i) <= 3 ) cycle
      ig = p%SOLID%CLOAD_ngrp_ID(i)
      is = hecmesh%node_group%grp_index(ig-1)
      ie = hecmesh%node_group%grp_index(ig)
      if( any( rot_node(hecmesh%node_group%grp_item(is+1:ie)) == 0 .and. &
        hecmesh%node_group%grp_item(is+1:ie) <= hecmesh%nn_internal ) ) call fstr_setup_util_err_stop( &
        '### Error: !CLOAD : dof 4-6 is given to nodes without rotational dof' )
      if( n_explicit(ig) > 0 ) call explicit_dummy_node_err_stop( '!CLOAD' )
      rot_grp(ig) = -1
    enddo
 
    n_new = count( rot_grp == -1 )
    if( n_new == 0 ) return
    allocate( new_name(n_new), new_index(0:n_new), new_list(sum(n_rot, mask=rot_grp == -1)) )
    new_index(0) = 0
    k = 0
    m = 0
    do ig = 1, hecmesh%node_group%n_grp
      if( rot_grp(ig) /= -1 ) cycle
      k = k + 1
      write( new_name(k), '(a,i0)' ) 'FSTR_ROT_', ig
      do j = hecmesh%node_group%grp_index(ig-1)+1, hecmesh%node_group%grp_index(ig)
        if( rot_node(hecmesh%node_group%grp_item(j)) == 0 ) cycle
        m = m + 1
        new_list(m) = rot_node(hecmesh%node_group%grp_item(j))
      enddo
      new_index(k) = m
      rot_grp(ig) = k
    enddo
    call append_new_groups( hecmesh, 'node_grp', n_new, new_name, new_index, new_list, first_id )
    where( rot_grp > 0 ) rot_grp = rot_grp + first_id - 1
 
    ! split each !BOUNDARY entry into dof 1-3 of its nodes and dof 1-3 of their dummy nodes
    n = p%SOLID%BOUNDARY_ngrp_tot
    if( n > 0 ) then
      allocate( src(2*n), is_rot(2*n) )
      m = 0
      do i = 1, n
        ide = mod(p%SOLID%BOUNDARY_ngrp_type(i), 10)
        if( p%SOLID%BOUNDARY_ngrp_rotID(i) <= 0 .and. ide > 3 .and. rot_grp(p%SOLID%BOUNDARY_ngrp_ID(i)) > 0 ) then
          if( p%SOLID%BOUNDARY_ngrp_type(i)/10 <= 3 ) then
            m = m + 1
            src(m) = i
            is_rot(m) = .false.
          endif
          m = m + 1
          src(m) = i
          is_rot(m) = .true.
        else
          m = m + 1
          src(m) = i
          is_rot(m) = .false.
        endif
      enddo
      call remap_integer_array( p%SOLID%BOUNDARY_ngrp_GRPID, src, m )
      call remap_integer_array( p%SOLID%BOUNDARY_ngrp_ID, src, m )
      call remap_integer_array( p%SOLID%BOUNDARY_ngrp_type, src, m )
      call remap_integer_array( p%SOLID%BOUNDARY_ngrp_amp, src, m )
      call remap_real_array( p%SOLID%BOUNDARY_ngrp_val, src, m )
      call remap_integer_array( p%SOLID%BOUNDARY_ngrp_istot, src, m )
      call remap_integer_array( p%SOLID%BOUNDARY_ngrp_rotID, src, m )
      call remap_integer_array( p%SOLID%BOUNDARY_ngrp_centerID, src, m )
      p%SOLID%BOUNDARY_ngrp_tot = m
      do k = 1, m
        if( p%SOLID%BOUNDARY_ngrp_rotID(k) > 0 ) cycle
        ig = p%SOLID%BOUNDARY_ngrp_ID(k)
        ids = p%SOLID%BOUNDARY_ngrp_type(k)/10
        ide = mod(p%SOLID%BOUNDARY_ngrp_type(k), 10)
        if( ide <= 3 .or. rot_grp(ig) <= 0 ) cycle
        if( is_rot(k) ) then
          p%SOLID%BOUNDARY_ngrp_ID(k) = rot_grp(ig)
          p%SOLID%BOUNDARY_ngrp_type(k) = 10 * (max(ids, 4) - 3) + (ide - 3)
        else
          p%SOLID%BOUNDARY_ngrp_type(k) = 10 * ids + 3
        endif
      enddo
    endif
 
    do i = 1, p%SOLID%CLOAD_ngrp_tot
      if( p%SOLID%CLOAD_ngrp_rotID(i) > 0 .or. p%SOLID%CLOAD_ngrp_DOF(i) <= 3 ) cycle
      p%SOLID%CLOAD_ngrp_ID(i) = rot_grp(p%SOLID%CLOAD_ngrp_ID(i))
      p%SOLID%CLOAD_ngrp_DOF(i) = p%SOLID%CLOAD_ngrp_DOF(i) - 3
    enddo
  end subroutine fstr_convert_rotation_dof
 
  subroutine rotation_dof_err_stop( header_name )
    implicit none
    character(len=*) :: header_name
 
    call fstr_setup_util_err_stop( '### Error: '//header_name// &
      ' : dof 4-6 of beam/shell nodes in a mesh with solid elements is not supported' )
  end subroutine rotation_dof_err_stop
 
  subroutine explicit_dummy_node_err_stop( header_name )
    implicit none
    character(len=*) :: header_name
 
    call fstr_setup_util_err_stop( '### Error: '//header_name// &
      ' : dof 4-6 is given to beam/shell nodes with explicit dummy nodes; give it as dof 1-3 of the dummy nodes' )
  end subroutine explicit_dummy_node_err_stop
 
  subroutine remap_integer_array( array, src, n )
    implicit none
    integer(kind=kint), pointer :: array(:)
    integer(kind=kint), intent(in) :: src(:), n
    integer(kind=kint), pointer :: new_array(:)
 
    allocate( new_array(n) )
    new_array(1:n) = array(src(1:n))
    deallocate( array )
    array => new_array
  end subroutine remap_integer_array
 
  subroutine remap_real_array( array, src, n )
    implicit none
    real(kind=kreal), pointer :: array(:)
    integer(kind=kint), intent(in) :: src(:), n
    real(kind=kreal), pointer :: new_array(:)
 
    allocate( new_array(n) )
    new_array(1:n) = array(src(1:n))
    deallocate( array )
    array => new_array
  end subroutine remap_real_array
 
end module m_fstr_setup