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 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 mcontact
  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, j, cid, nout, nin, ierror, cparam_id
    character(len=HECMW_NAME_LEN) :: header_name, fname(MAXOUTFILE)
    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
 
    ! counters
    integer(kind=kint) :: c_solution, c_solver, c_nlsolver, c_step, c_write, 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_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
    integer(kind=kint) :: c_output, islog
    integer(kind=kint) :: k
    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_eigenread = 0
    c_elemopt  = 0;
    c_aincparam= 0; c_timepoints = 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
    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 == '!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 == '!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, 'S', 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
        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
    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
    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
    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%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_exapnsion)=alpha
        fstrsolid%materials(cid)%variables(m_thick)=thick
        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
 
    ! 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
    enddo
 
    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
    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
        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
        endif
        cparam_id = 0
        do i=1,size(fstrparam%contactparam)-1
          if( fstr_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
          else
            if(paracontactflag) then
              isok = fstr_contact_init( fstrsolid%contacts(c_contact+i), p%MESH, fstrparam%contactparam(cparam_id), myrank)
            else
              isok = fstr_contact_init( fstrsolid%contacts(c_contact+i), p%MESH, fstrparam%contactparam(cparam_id))
            endif
            !       call fstr_write_contact( 6, fstrSOLID%contacts(c_contact+i) )
          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
        endif
        cparam_id = 0
        do i=1,size(fstrparam%contactparam)-1
          if( fstr_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
          else
            if(paracontactflag) then
              isok = fstr_embed_init( fstrsolid%embeds(c_embed+i), p%MESH, fstrparam%contactparam(cparam_id), myrank)
            else
              isok = fstr_embed_init( fstrsolid%embeds(c_embed+i), p%MESH, fstrparam%contactparam(cparam_id))
            endif
          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 ) ) then
          write(*,*) '### Error: Fail in read in step definition : ' , c_istep
          write(ilog,*) '### Error: Fail in read in step definition : ', c_istep
          stop
        endif
        if( associated(fstrparam%timepoints) ) then
          do i=1,size(fstrparam%timepoints)
            if( fstr_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( fstr_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 ) ) then
          write(*,*) '### Error: Fail in read in step definition : ' , c_istep
          write(ilog,*) '### Error: Fail in read in step definition : ', c_istep
          stop
        endif
        ! 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( fstr_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( fstr_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
        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
        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
        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
        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
        endif
        cid = 0
        if(cache < hecmesh%material%n_mat) then
          if(fstr_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( fstr_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
        endif
        fstrsolid%materials(cid)%name = mname
        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
          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
          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
          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
          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
            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
          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
          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
          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
          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
          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
          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
          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
          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
          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
          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
        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
        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
        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
        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
        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:n+1) ) /= 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
        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
          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
        endif
 
      else if( header_name == '!INCLUDE' ) then
        ctrl_list(ictrl) = ctrl
        input_filename   = ""
        ierror = fstr_ctrl_get_param_ex( ctrl, 'INPUT ', '# ', 0, 'S', 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
        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 )
      !call fstr_print_steps( 6, 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 hecmw_abort( hecmw_comm_get_comm())
      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 )
    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_solid_init( hecMESH, fstrSOLID )
    use m_fstr
    type(hecmwst_local_mesh),target :: hecMESH
    type(fstr_solid)                :: fstrSOLID
 
    integer :: ndof, ntotal, ierror, ic_type
 
    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%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, ic_type
 
    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
    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, isect, nodlocal(fstrSOLID%max_ncon), iiS, nn, con_stf
 
    if( fstrsolid%max_ncon_stf > 20 ) fstrsolid%max_ncon_stf = 20   
 
    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 ) cycle
      iis = hecmesh%elem_node_index(i-1)
      nn = hecmesh%elem_node_index(i-1) - 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
    integer :: ncon_stf
 
    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) ) ) &
        nthick = fstr_shell_num_thickness_points( fstrsolid%elements(i)%etype )
      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
 
  subroutine fstr_solid_finalize( fstrSOLID )
    use mmechgauss, only: fstr_finalize_shell_layer_gausses
    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
 
    ! 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
 
  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  = 1
    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"
    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
    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)
    !   dumptype   => svIarray(31)
    !   dumpexit   => svIarray(32)
    !   usejad     => svIarray(33)
    !   ncolor_in  => svIarray(34)
    !   mpc_method => svIarray(13)
    !   estcond    => svIarray(14)
    !   contact_elim=> svIarray(15)
    !   method2    => svIarray(8)
    !   recyclepre => svIarray(35)
    !   solver_opt => svIarray(41:50)
    !   nBFGS      => svIarray(60)
 
    !   resid      => svRarray(1)
    !   sigma_diag => svRarray(2)
    !   sigma      => svRarray(3)
    !   thresh     => svRarray(4)
    !   filter     => svRarray(5)
 
    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) )
    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 )
    implicit none
    integer(kind=kint)         :: ctrl
    type( hecmwst_local_mesh ) :: hecmesh
    integer                    :: cnt
    type( tlocalcoordsys )     :: coordsys
 
    integer                   :: j, is, ie, grp_id(1)
    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 '
    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)
    call fstr_strupr( coordsys%sys_name )
 
    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
      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
 
  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 )
    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) :: i, 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 fstr_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 )
    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 fstr_expand_integer_array ( p%SOLID%elemact%ELEMACT_egrp_GRPID, old_size, new_size )
    call fstr_expand_integer_array ( p%SOLID%elemact%ELEMACT_egrp_ID,  old_size, new_size )
    call fstr_expand_integer_array ( p%SOLID%elemact%ELEMACT_egrp_amp, old_size, new_size )
    call fstr_expand_real_array ( p%SOLID%elemact%ELEMACT_egrp_eps, old_size, new_size )
    call fstr_expand_integer_array ( p%SOLID%elemact%ELEMACT_egrp_depends, old_size, new_size )
    call fstr_expand_real_array ( p%SOLID%elemact%ELEMACT_egrp_ts_lower, old_size, new_size )
    call fstr_expand_real_array ( p%SOLID%elemact%ELEMACT_egrp_ts_upper, old_size, new_size )
    call fstr_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 )
    implicit none
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: counter
    type(fstr_param_pack) :: P
 
    integer(kind=kint) :: rcode
    integer(kind=kint) :: type = 0
    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  )
    !  rcode = fstr_ctrl_get_param_ex( ctrl, 'TYPE ', 'FSTR,NASTRAN ', 0, 'P', type )
    !  if( rcode < 0 ) call fstr_ctrl_err_stop
    !  if( rcode == 1 ) type = 0 ! PARAM_NOTHING
 
    !  if( type == 0 ) then
 
    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 fstr_expand_integer_array (p%SOLID%BOUNDARY_ngrp_GRPID, old_size, new_size )
    call fstr_expand_integer_array (p%SOLID%BOUNDARY_ngrp_ID,    old_size, new_size )
    call fstr_expand_integer_array (p%SOLID%BOUNDARY_ngrp_type,  old_size, new_size )
    call fstr_expand_real_array    (p%SOLID%BOUNDARY_ngrp_val,   old_size, new_size )
    call fstr_expand_integer_array (p%SOLID%BOUNDARY_ngrp_amp,   old_size, new_size )
    call fstr_expand_integer_array (p%SOLID%BOUNDARY_ngrp_istot, old_size, new_size )
    call fstr_expand_integer_array (p%SOLID%BOUNDARY_ngrp_rotID, old_size, new_size )
    call fstr_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 )
    !  else
    !   ! NASTRAN ---------------------------------------------
    !
    !     P%SOLID%file_type = kbcfNASTRAN
    !     call fstr_setup_solid_nastran( ctrl, P%MESH, P%SOLID )
    !  end if
 
  end subroutine fstr_setup_boundary
 
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_cload( 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) :: 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
 
    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 fstr_expand_integer_array ( p%SOLID%CLOAD_ngrp_GRPID, old_size, new_size )
    call fstr_expand_integer_array ( p%SOLID%CLOAD_ngrp_ID,  old_size, new_size )
    call fstr_expand_integer_array ( p%SOLID%CLOAD_ngrp_DOF, old_size, new_size )
    call fstr_expand_real_array    ( p%SOLID%CLOAD_ngrp_val, old_size, new_size )
    call fstr_expand_integer_array ( p%SOLID%CLOAD_ngrp_amp, old_size, new_size )
    call fstr_expand_integer_array ( p%SOLID%CLOAD_ngrp_rotID, old_size, new_size )
    call fstr_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 )
 
    if( p%MESH%n_refine > 0 ) then
      do i=1,n
        if( hecmw_ngrp_get_number(p%MESH, p%SOLID%CLOAD_NGRP_ID(old_size+i)) > 1 ) then
          write(*,*) 'fstr control file error : !CLOAD : cannot be used with NGRP when mesh is refined'
          write(ilog,*) 'fstr control file error : !CLOAD : cannot be used with NGRP when mesh is refined'
          call fstr_ctrl_err_stop
        endif
      enddo
    endif
 
  end subroutine fstr_setup_cload
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
  subroutine fstr_setup_fload( ctrl, counter, P )
  !---- args
    integer(kind=kint)   :: ctrl
    integer(kind=kint)   :: counter
    type(fstr_param_pack) :: P
  !---- vals
    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           :: 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 fstr_expand_integer_array( p%FREQ%FLOAD_ngrp_GRPID, old_size, new_size )
    call fstr_expand_integer_array( p%FREQ%FLOAD_ngrp_ID,    old_size, new_size )
    call fstr_expand_integer_array( p%FREQ%FLOAD_ngrp_TYPE,  old_size, new_size )
    call fstr_expand_integer_array( p%FREQ%FLOAD_ngrp_DOF,   old_size, new_size )
    call fstr_expand_real_array   ( p%FREQ%FLOAD_ngrp_valre, old_size, new_size )
    call fstr_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
    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_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)') 'S', 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 )
    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 fstr_expand_integer_array ( p%SOLID%DLOAD_ngrp_GRPID, old_size, new_size )
    call fstr_expand_integer_array ( p%SOLID%DLOAD_ngrp_ID,  old_size, new_size )
    call fstr_expand_integer_array ( p%SOLID%DLOAD_ngrp_LID, old_size, new_size )
    call fstr_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 )
    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 fstr_expand_integer_array ( p%SOLID%TEMP_ngrp_GRPID, old_size, new_size )
    call fstr_expand_integer_array ( p%SOLID%TEMP_ngrp_ID, old_size, new_size )
    call fstr_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 )
    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 :: 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 fstr_expand_integer_array ( p%SOLID%SPRING_ngrp_GRPID, old_size, new_size )
    call fstr_expand_integer_array ( p%SOLID%SPRING_ngrp_ID,  old_size, new_size )
    call fstr_expand_integer_array ( p%SOLID%SPRING_ngrp_DOF, old_size, new_size )
    call fstr_expand_real_array    ( p%SOLID%SPRING_ngrp_val, old_size, new_size )
    call fstr_expand_integer_array ( p%SOLID%SPRING_ngrp_amp, old_size, new_size )
 
    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_spring( ctrl, amp, grp_id_name, hecmw_name_len, id_ptr, val_ptr )
    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:))
 
    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( fstr_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 )
    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 fstr_expand_integer_array( p%HEAT%T_FIX_node, old_size, new_size )
    call fstr_expand_integer_array( p%HEAT%T_FIX_ampl, old_size, new_size )
    call fstr_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 )
    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 fstr_expand_integer_array( p%HEAT%Q_NOD_node, old_size, new_size )
    call fstr_expand_integer_array( p%HEAT%Q_NOD_ampl, old_size, new_size )
    call fstr_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 )
    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 fstr_expand_integer_array( p%HEAT%Q_SUF_elem, old_size, new_size )
    call fstr_expand_integer_array( p%HEAT%Q_SUF_ampl, old_size, new_size )
    call fstr_expand_integer_array( p%HEAT%Q_SUF_surf, old_size, new_size )
    call fstr_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 )
    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 fstr_expand_integer_array( p%HEAT%Q_SUF_elem, old_size, new_size )
    call fstr_expand_integer_array( p%HEAT%Q_SUF_ampl, old_size, new_size )
    call fstr_expand_integer_array( p%HEAT%Q_SUF_surf, old_size, new_size )
    call fstr_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 )
    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 fstr_expand_integer_array(  p%HEAT%H_SUF_elem,    old_size, new_size )
    call fstr_expand_integer_array2( p%HEAT%H_SUF_ampl, 2, old_size, new_size )
    call fstr_expand_integer_array(  p%HEAT%H_SUF_surf,    old_size, new_size )
    call fstr_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 )
    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 fstr_expand_integer_array(  p%HEAT%H_SUF_elem,    old_size, new_size )
    call fstr_expand_integer_array2( p%HEAT%H_SUF_ampl, 2, old_size, new_size )
    call fstr_expand_integer_array(  p%HEAT%H_SUF_surf,    old_size, new_size )
    call fstr_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 )
    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 fstr_expand_integer_array(  p%HEAT%R_SUF_elem,    old_size, new_size )
    call fstr_expand_integer_array2( p%HEAT%R_SUF_ampl, 2, old_size, new_size )
    call fstr_expand_integer_array(  p%HEAT%R_SUF_surf,    old_size, new_size )
    call fstr_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 )
    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 fstr_expand_integer_array(  p%HEAT%R_SUF_elem,    old_size, new_size )
    call fstr_expand_integer_array2( p%HEAT%R_SUF_ampl, 2, old_size, new_size )
    call fstr_expand_integer_array(  p%HEAT%R_SUF_surf,    old_size, new_size )
    call fstr_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_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 )
    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
    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%VELOCITY_ngrp_tot
    new_size = old_size + n
    p%SOLID%VELOCITY_ngrp_tot = new_size
 
    call fstr_expand_integer_array (p%SOLID%VELOCITY_ngrp_GRPID, old_size, new_size )
    call fstr_expand_integer_array (p%SOLID%VELOCITY_ngrp_ID  , old_size, new_size )
    call fstr_expand_integer_array (p%SOLID%VELOCITY_ngrp_type, old_size, new_size )
    call fstr_expand_real_array    (p%SOLID%VELOCITY_ngrp_val , old_size, new_size )
    call fstr_expand_integer_array (p%SOLID%VELOCITY_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_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
 
    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
      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 )
    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 fstr_expand_integer_array (p%SOLID%ACCELERATION_ngrp_GRPID, old_size, new_size )
    call fstr_expand_integer_array (p%SOLID%ACCELERATION_ngrp_ID  , old_size, new_size )
    call fstr_expand_integer_array (p%SOLID%ACCELERATION_ngrp_type, old_size, new_size )
    call fstr_expand_real_array    (p%SOLID%ACCELERATION_ngrp_val , old_size, new_size )
    call fstr_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
      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 fstr_expand_integer_array ( P%MPC_RD%node1,  old_size, new_size )
    !        call fstr_expand_integer_array ( P%MPC_RD%node2,  old_size, new_size )
    !        call fstr_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 hecmw_abort( hecmw_comm_get_comm())
  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 )
    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, ngrp_id2
    ! 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
      ! ! for DEBUG
      ! sgrp_id = hecMESH%contact_pair%master_grp_id(i)
      ! call append_node_grp_from_surf_grp( hecMESH, sgrp_id, ngrp_id2 )
      ! ! intersection node group of slave and master
      ! call append_intersection_node_grp( hecMESH, ngrp_id, ngrp_id2 )
      ! ! intersection node_group of original slave and patch-slave
      ! ngrp_id=get_grp_id( hecMESH, 'node_grp', 'SLAVE' )
      ! ngrp_id2=get_grp_id( hecMESH, 'node_grp', '_PT_SLAVE_S' )
      ! call append_intersection_node_grp( hecMESH, ngrp_id, ngrp_id2 )
    enddo
  end subroutine fstr_convert_contact_type
 
end module m_fstr_setup