READ_INPUT_FILES.f90

Go to the documentation of this file.

!    Copyright (C) 2012 The SPEED FOUNDATION
!    Author: Ilario Mazzieri
!
!    This file is part of SPEED.
!
!    SPEED is free software; you can redistribute it and/or modify it
!    under the terms of the GNU Affero General Public License as
!    published by the Free Software Foundation, either version 3 of the
!    License, or (at your option) any later version.
!
!    SPEED is distributed in the hope that it will be useful, but
!    WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!    Affero General Public License for more details.
!
!    You should have received a copy of the GNU Affero General Public License
!    along with SPEED.  If not, see <http://www.gnu.org/licenses/>.
 
 
 
 
     subroutine read_input_files()
     
      use max_var
      use speed_par
      
 
      implicit none
 
      include 'SPEED.MPI'
 
      head_file = 'SPEED.input'
      
      if (mpi_id.eq.0) then
         write(*,'(A)') 
         write(*,'(A)')'------------------Reading Header File------------------'
         write(*,'(A,A50)') 'Header File: ',head_file
      endif
      
      inquire(file=head_file,exist=filefound);
      if(filefound .eqv. .false.) then
        write(*,*) 'File ', head_file, ' is missing!'
        call exit(exit_missing_file)
      endif
       
      opt_out_var= 1;  opt_out_form = 1; opt_out_data = 1
      
      restart = 0;
      trestart = 0;
      
      call read_dime_header(head_file,restart,ntime_err,debug) 
 
      if(ntime_err .ne. 0) allocate(time_error(ntime_err))
            
      tstart = 0.0
      tstop = 0.0
 
      
      ! Set default variables
      b_setuponly = b_setuponly_default
      b_failoncoeffs = b_failoncoeffs_default
      b_instabilitycontrol = b_instabilitycontrol_default
      b_failcfl = b_failcfl_default
      instability_maxval = instability_maxval_default
      dg_c = -100.d0;
      pen_c = -100.d0;
 
      call read_header(head_file,grid_file,mat_file, mpi_file, monitor_file, bkp_file, &
                       sdof_file, deltat,tstart,tstop,&
                       opt_out_var,& 
                       trestart,&
                       ndt_mon_lst,&        
                       deltat_fixed,&        
                       depth_search_mon_pgm,ndt_mon_pgm,num_pgm,&
                       rotation_angle_mon_pgm,file_mon_pgm,&        
                       depth_search_mon_lst,num_lst,file_mon_lst,&        
                       depth_search_sys_lst,sys_lst,file_sys_lst,&
                       dg_c,pen_c,&
                       rk_scheme, rk_order, rk_stages, testmode, ntime_err, time_error,&
                       damping_type, &
                       num_testcase, label_testcase, &
                       b_failoncoeffs, b_setuponly, b_failcfl, b_instabilitycontrol, instability_maxval)
 
      if (mpi_id.eq.0) write(*,'(A)')'Read.'
 
      ! Instability control
      if (mpi_id .eq. 0) then
          if (b_instabilitycontrol) then
              write(*,'(A,L)') 'Instability control  :   ', b_instabilitycontrol
              write(*,'(A,E12.4,A,E12.4,A)') 'Instability threshold: ', instability_maxval,' (default: ', instability_maxval_default, ')'
          endif
      endif
      
      
      if(restart .eq. 1 .and. mpi_id .eq. 0) then 
         write(*,*) 'Restart Active ...' 
         write(*,*) 'Backup every ', trestart, ' time steps'          
      endif
 
      
      
      if (mpi_id.eq.0) then
         if(dg_c .eq. 1.d0)   write(*,'(A)') 'DG METHOD    : NIPG'
         if(dg_c .eq. 0.d0)   write(*,'(A)') 'DG METHOD    : IIPG'
         if(dg_c .eq. -1.d0)  write(*,'(A)') 'DG METHOD    : SIPG'
         if(pen_c .ge. 0.d0)  write(*,'(A,E12.4)') 'PENALIZATION : ', pen_c
         if(dg_c .eq. -100)   write(*,'(A)') 'SEM METHOD    '
      endif
 
 
      mat_file = mat_file(1:len_trim(mat_file)) // '.mate'
      
      if (mpi_id.eq.0) then
         write(*,'(A)')    
         write(*,'(A,A20)')'-----------------Reading Material File-----------------'
         write(*,'(A,A20)')'Material File    : ',mat_file
      endif
      
      inquire(file=mat_file,exist=filefound);
      if(filefound .eqv. .false.) call exit(exit_missing_file)
      
      call read_dime_filemate(mat_file,nmat, &
                            nmat_nle, nmat_rnd, &                                                                        
                            nload_dirx_el,nload_diry_el,nload_dirz_el, &
                            nload_neux_el,nload_neuy_el,nload_neuz_el, &
                            nload_neun_el, &                                                            
                            nload_poix_el,nload_poiy_el,nload_poiz_el, &
                            nload_trax_el,nload_tray_el,nload_traz_el, &
                            nload_plax_el,nload_play_el,nload_plaz_el, &                                
                            nload_forx_el,nload_fory_el,nload_forz_el, &
                            nload_forc_el,nload_pres_el,nload_shea_el, &
                            nload_abc_el,nload_dg_el,nfunc,nfunc_data,&
                            nload_sism_el,&                                                                 
                            nload_expl_el,&                                                                
                            n_case, nmat_nhe, n_test,n_frac,srcmodflag)                   
      
      sdofflag = 0                      
      if(sys_lst.eq.1) sdofflag = 1                      
      
      if(n_test.gt.0 .and. mpi_id .eq. 0)  write(*,'(A)')'*********TEST MODE*********'
 
      if (mpi_id.eq.0) then
         write(*,'(A,I8)') 'Materials        : ',nmat
         if(nmat_nle.gt.0) write(*,'(A,I8)')     'Materials NL-El. : ',nmat_nle        
         if(nmat_rnd.gt.0) write(*,'(A,I8)')     'Materials Random : ',nmat_rnd        
         if(nload_dirx_el.gt.0) write(*,'(A,I8)')'Dirichlet X B.C. : ',nload_dirx_el
         if(nload_diry_el.gt.0) write(*,'(A,I8)')'Dirichlet Y B.C. : ',nload_diry_el
         if(nload_dirz_el.gt.0) write(*,'(A,I8)')'Dirichlet Z B.C. : ',nload_dirz_el
         if(nload_neux_el.gt.0) write(*,'(A,I8)')'Neumann X B.C.   : ',nload_neux_el
         if(nload_neuy_el.gt.0) write(*,'(A,I8)')'Neumann Y B.C.   : ',nload_neuy_el
         if(nload_neuz_el.gt.0) write(*,'(A,I8)')'Neumann Z B.C.   : ',nload_neuz_el
         if(nload_neun_el.gt.0) write(*,'(A,I8)')'Neumann N B.C.   : ',nload_neun_el              
         if(nload_poix_el.gt.0) write(*,'(A,I8)')'Point Loads X    : ',nload_poix_el
         if(nload_poiy_el.gt.0) write(*,'(A,I8)')'Point Loads Y    : ',nload_poiy_el
         if(nload_poiz_el.gt.0) write(*,'(A,I8)')'Point Loads Z    : ',nload_poiz_el
!         if(nload_traX_el.gt.0) write(*,'(A,I8)')'Travel. Load X   : ',nload_traX_el
!         if(nload_traY_el.gt.0) write(*,'(A,I8)')'Travel. Load Y   : ',nload_traY_el
         if(nload_traz_el.gt.0) write(*,'(A,I8)')'Travel. Load Z   : ',nload_traz_el         
         if(nload_plax_el.gt.0) write(*,'(A,I8)')'Plane Loads X    : ',nload_plax_el              
         if(nload_play_el.gt.0) write(*,'(A,I8)')'Plane Loads Y    : ',nload_play_el              
         if(nload_plaz_el.gt.0) write(*,'(A,I8)')'Plane Loads Z    : ',nload_plaz_el              
         if(nload_forx_el.gt.0) write(*,'(A,I8)')'Force X          : ',nload_forx_el
         if(nload_fory_el.gt.0) write(*,'(A,I8)')'Force Y          : ',nload_fory_el
         if(nload_forz_el.gt.0) write(*,'(A,I8)')'Force Z          : ',nload_forz_el
         if(nload_forc_el.gt.0) write(*,'(A,I8)')'Force load       : ',nload_forc_el
         if(nload_pres_el.gt.0) write(*,'(A,I8)')'Pressure Load    : ',nload_pres_el
         if(nload_shea_el.gt.0) write(*,'(A,I8)')'Shear Load       : ',nload_shea_el
         if(nload_abc_el.gt.0)  write(*,'(A,I8)')'ABSO Boundaries  : ',nload_abc_el
         if(nload_dg_el.gt.0)   write(*,'(A,I8)')'DG Interfaces    : ',nload_dg_el
         if(nfunc.gt.0)         write(*,'(A,I8)')'Functions        : ',nfunc
         if(srcmodflag.eq.1) write(*,*)'Using Not Honoring Fault Method For Siesmic Sources'
         if(sdofflag.eq.1) write(*,*)'Coupling SDOF response :: True'
         if(nload_sism_el.gt.0) write(*,'(A,I8)')'Seis. Mom. Load  : ',nload_sism_el        
         if(nload_expl_el.gt.0) write(*,'(A,I8)')'Explosive Load   : ',nload_expl_el            
         if (n_case.gt.1) then
            write(*,'(A)')'CASE WARNING: More than one case defined,'        
            write(*,'(A)')'              only the 1st case will be adopted'  
         endif
         write(*,'(A,I8)') 'CASE             : ',n_case                               
 
         if(nmat_nhe.gt.0) then
            write(*,'(A,I8)')     'Not_Honoring Enhanced Blocks : ',nmat_nhe        
          endif
      endif
      
      if (nmat.le.0) then
         write(*,*)'Error ! nmat = 0'
         call exit(exit_no_materials)
      endif
      
      allocate (type_mat(nmat), tref_mat(nmat), prop_mat(nmat,4), sdeg_mat(nmat), tag_mat(nmat))
 
      if (nmat_nle.gt.0) allocate(type_mat_nle(nmat_nle), prop_mat_nle(nmat_nle,1), val_mat_nle(nmat_nle,1), tag_mat_nle(nmat_nle))
      if (nmat_rnd.gt.0) allocate(rand_mat(nmat_rnd))
      
      if (nload_dirx_el.gt.0) allocate (val_dirx_el(nload_dirx_el,4), fun_dirx_el(nload_dirx_el), tag_dirx_el(nload_dirx_el))      
      if (nload_diry_el.gt.0) allocate (val_diry_el(nload_diry_el,4), fun_diry_el(nload_diry_el), tag_diry_el(nload_diry_el))      
      if (nload_dirz_el.gt.0) allocate (val_dirz_el(nload_dirz_el,4), fun_dirz_el(nload_dirz_el), tag_dirz_el(nload_dirz_el))
      
      if (nload_neux_el.gt.0) allocate (val_neux_el(nload_neux_el,4), fun_neux_el(nload_neux_el), tag_neux_el(nload_neux_el))      
      if (nload_neuy_el.gt.0) allocate (val_neuy_el(nload_neuy_el,4), fun_neuy_el(nload_neuy_el), tag_neuy_el(nload_neuy_el))      
      if (nload_neuz_el.gt.0) allocate (val_neuz_el(nload_neuz_el,4), fun_neuz_el(nload_neuz_el), tag_neuz_el(nload_neuz_el))
      if (nload_neun_el.gt.0) allocate (val_neun_el(nload_neun_el,4), fun_neun_el(nload_neun_el), tag_neun_el(nload_neun_el))
      
      if (nload_poix_el.gt.0) allocate (val_poix_el(nload_poix_el,4), fun_poix_el(nload_poix_el))
      if (nload_poiy_el.gt.0) allocate (val_poiy_el(nload_poiy_el,4), fun_poiy_el(nload_poiy_el))   
      if (nload_poiz_el.gt.0) allocate (val_poiz_el(nload_poiz_el,4), fun_poiz_el(nload_poiz_el))
 
!      if (nload_traX_el.gt.0) allocate (val_traX_el(nload_traX_el,4), fun_traX_el(nload_traX_el))    
!      if (nload_traY_el.gt.0) allocate (val_traY_el(nload_traY_el,4), fun_traY_el(nload_traY_el))   
      if (nload_traz_el.gt.0) allocate (val_traz_el(nload_traz_el,2), fun_traz_el(nload_traz_el))
 
 
      if (nload_plax_el.gt.0) allocate (val_plax_el(nload_plax_el,1), fun_plax_el(nload_plax_el), tag_plax_el(nload_plax_el))
      if (nload_play_el.gt.0) allocate (val_play_el(nload_play_el,1), fun_play_el(nload_play_el), tag_play_el(nload_play_el))
      if (nload_plaz_el.gt.0) allocate (val_plaz_el(nload_plaz_el,1), fun_plaz_el(nload_plaz_el), tag_plaz_el(nload_plaz_el))
      
      if (nload_forx_el.gt.0) allocate (val_forx_el(nload_forx_el,4), fun_forx_el(nload_forx_el))
      if (nload_fory_el.gt.0) allocate (val_fory_el(nload_fory_el,4), fun_fory_el(nload_fory_el))
      if (nload_forz_el.gt.0) allocate (val_forz_el(nload_forz_el,4), fun_forz_el(nload_forz_el))
     
     
      if (n_test.gt.0) allocate (fun_test(n_test)) !val_fun_test(n_test,10),      
                                       
      if (nload_forc_el.gt.0) allocate (val_forc_el(nload_forc_el,10), fun_forc_el(nload_forc_el))     
      if (nload_pres_el.gt.0) allocate (val_pres_el(nload_pres_el,10), fun_pres_el(nload_pres_el))      
      if (nload_shea_el.gt.0) allocate (val_shea_el(nload_shea_el,10), fun_shea_el(nload_shea_el))
      
      if (nload_abc_el.gt.0) allocate (tag_abc_el(nload_abc_el))
      if (nload_dg_el .gt. 0) allocate (tag_dg_el(nload_dg_el), tag_dg_yn(nload_dg_el),tag_dg_link(nload_dg_el), tag_dg_frc(nload_dg_el), val_dg_frc(nload_dg_el,2))
      if (nload_dg_el .gt. 0) then
              tag_dg_frc=0; val_dg_frc=0; tag_dg_link =0;
      endif
 
      if (srcmodflag.eq.0) then
          szsism = 21
          if (nload_sism_el.gt.0) allocate (val_sism_el(nload_sism_el,21), &
                                            fun_sism_el(nload_sism_el), tag_sism_el(nload_sism_el))
      elseif (srcmodflag.eq.1) then
          szsism = 15
          if (nload_sism_el.gt.0) allocate (val_sism_el(nload_sism_el,15), &
                                            fun_sism_el(nload_sism_el), tag_sism_el(nload_sism_el))
      endif
 
 
      if (nload_expl_el.gt.0) allocate (val_expl_el(nload_expl_el,20), &
                                        fun_expl_el(nload_expl_el), tag_expl_el(nload_expl_el))
      
      if (n_case.gt.0) allocate (val_case(n_case), tag_case(n_case), tol_case(n_case))
      if (n_case .eq. 0) allocate(tag_case(1)); tag_case(1)=0;
 
      if (nmat_nhe.gt.0) allocate (val_nhe(nmat_nhe), tol_nhe(nmat_nhe))
      
      if (nfunc.gt.0) allocate (tag_func(nfunc), func_type(nfunc), func_indx(nfunc +1), func_data(nfunc_data))
      
      
      allocate(qs(nmat), qp(nmat)); qs = 0.d0; qp = 0.d0;
 
      call read_filemate(mat_file,nmat,prop_mat,type_mat,tref_mat,tag_mat, qs, qp,&
                nmat_nle,prop_mat_nle,val_mat_nle,type_mat_nle,tag_mat_nle, &   
                nmat_rnd, rand_mat, &     
                nload_dirx_el,val_dirx_el,fun_dirx_el,tag_dirx_el, &
                nload_diry_el,val_diry_el,fun_diry_el,tag_diry_el, &
                nload_dirz_el,val_dirz_el,fun_dirz_el,tag_dirz_el, &
                nload_neux_el,val_neux_el,fun_neux_el,tag_neux_el, &
                nload_neuy_el,val_neuy_el,fun_neuy_el,tag_neuy_el, &
                nload_neuz_el,val_neuz_el,fun_neuz_el,tag_neuz_el, &
                nload_neun_el,val_neun_el,fun_neun_el,tag_neun_el, &            
                nload_poix_el,val_poix_el,fun_poix_el, &
                nload_poiy_el,val_poiy_el,fun_poiy_el, &
                nload_poiz_el,val_poiz_el,fun_poiz_el, &
                nload_trax_el,val_trax_el,fun_trax_el,&
                nload_tray_el,val_tray_el,fun_tray_el,&
                nload_traz_el,val_traz_el,fun_traz_el,&               
                nload_plax_el,val_plax_el,fun_plax_el,tag_plax_el, &            
                nload_play_el,val_play_el,fun_play_el,tag_play_el, &                 
                nload_plaz_el,val_plaz_el,fun_plaz_el,tag_plaz_el, &                 
                nload_forx_el,val_forx_el,fun_forx_el,&
                nload_fory_el,val_fory_el,fun_fory_el,&
                nload_forz_el,val_forz_el,fun_forz_el,&
                nload_forc_el,val_forc_el,fun_forc_el,&
                nload_pres_el,val_pres_el,fun_pres_el,&
                nload_shea_el,val_shea_el,fun_shea_el,&
                n_test,fun_test,& !val_fun_test,&
                nload_abc_el,tag_abc_el,&
                nload_dg_el,tag_dg_el,tag_dg_yn, tag_dg_link, &
                tag_dg_frc, val_dg_frc, n_frac, &
                srcmodflag, szsism, &
                nload_sism_el,val_sism_el,fun_sism_el,tag_sism_el, &                 
                nload_expl_el,val_expl_el,fun_expl_el,tag_expl_el, &                 
                n_case,val_case,tag_case,tol_case, &
                nmat_nhe,val_nhe,tol_nhe, &                                 
                nfunc,func_type,func_indx,func_data,tag_func,nfunc_data, &
                fmax,fpeak,damping_type)
                 
 
      if(fmax .eq. 0.d0) then
         fmax = 3.d0 
         if (mpi_id .eq. 0) write(*,*) 'ATTENTION: FMAX not defined!'
         if (mpi_id .eq. 0) write(*,*) 'FMAX ASSUMED = 3'
      endif   
 
 
 
 
      if (damping_type .eq. 1) then 
 
         do im = 1, nmat
            if(qs(im) .eq. 0.d0) then 
               prop_mat(im,4) = 0.d0;
            else
               prop_mat(im,4) = 4.d0*datan(1.d0)*(fmax)/qs(im)
            endif   
 
         enddo
         
         
       else 
       
            prop_mat(:,4) = 0.d0;
 
       endif   
      
      
      
!**************************************************
! Materials can have different spectral degree
 
      do i = 1,nmat
         sdeg_mat(i) = type_mat(i)
      enddo
 
!************************************************** 
           
      if (mpi_id.eq.0) write(*,'(A)')'Read.'
 
      if (mpi_id.eq.0) then                        
         write(*,'(A)')                                
         if (n_case.gt.0) then
           do i = 1, n_case
                               
              write(*,'(A)') '------------------Not-Honoring case-------------------'
              write(*,'(A,I8)') 'CASE :', tag_case(i)    !tag_case(1)
            select case (tag_case(i))  !(tag_case(1))
              case(1)
                write(*,'(A)')'GRENOBLE H.'                        
              case(2)
                write(*,'(A)')'GRENOBLE'                        
              case(3)
                write(*,'(A)')'GUBBIO'                        
              case(4)
                write(*,'(A)')'SULMONA'                        
              case(5)
                write(*,'(A)')'VOLVI'                        
              case(6)
                write(*,'(A)')'FRIULI'                        
              case(7)
                write(*,'(A)')'AQUILA'   
              case(8)
                write(*,'(A)')'SANTIAGO'   
              case(11)
                write(*,'(A)')'CHRISTCHURCH NEW TOPO'        
              case(12)
                write(*,'(A)')'PO PLAIN (new model)'    
              case(13)
                write(*,'(A)')'PO PLAIN-BEDROCK (new-model)'                                 
              case(14)
                write(*,'(A)')'WELLINGTON (Benites)'                                 
              case(15)
                write(*,'(A)')'MARSICA-FUCINO'   
              case(16)
                write(*,'(A)')'ISTANBUL'                     
              case(18)
                write(*,'(A)')'BEIJING-TUTORIAL'                     
              case(19)
                write(*,'(A)')'THESSALONIKI'                     
              case(20)
                write(*,'(A)')'ATHENS'                     
              case(21)
                write(*,'(A)')'BEIJING'
              case(22)
                write(*,'(A)')'NORCIA'      
              case(27)
                write(*,'(A)')'AQUILA-OB'        
              case(28)
                write(*,'(A)')'NORCIA-OB' 
              case(29)
                write(*,'(A)')'THESS-OB'                             
              case(30) 
                write(*,'(A)')'ATHENS-Parthenon'
              case(31) 
                write(*,'(A)')'GRONINGEN'                     
              case(32) 
                write(*,'(A)')'GRONINGEN-LAYERED MODEL' 
              case(33) 
                write(*,'(A)')'GRONINGEN-ZE' 
              case(35) 
                write(*,'(A)')'THESS+MYGD-FINAL'                  
              case(38) 
                write(*,'(A)')'MONTELIMAR'                  
              case(40)
                write(*,'(A)')'KUTCH BASIN, INDIA'                                            
              case(46)
                write(*,'(A)')'KUMAMOTO, JAPAN'
              case(60)
                write(*,'(A)')'JAKARTA, INDONESIA'
              case(70)
                write(*,'(A)')'AQUILA MULTI-BASIN'
              case(91)
                write(*,'(A)')'IRPINIA'  
              case(98,99,100)
                write(*,'(A)')'TEST MODE'                                
              case default
                 write(*,'(A)')'ERROR.. this case was not implemented!'
             end select 
                                   
               write(*,'(A,I8)')   'MATERIAL TYPE    : ',val_case(i)   !(val_case(1,1))        
              write(*,'(A,E12.4)') 'MATERIAL TOL.    : ',tol_case(i)   !(tol_case(1))
            write(*,*)
            enddo
         endif
      endif
      
      if (mpi_id.eq.0) then
         write(*,'(A)')
         do im = 1,nmat    
            write(*,'(A,I8)')    'MATERIAL : ',tag_mat(im)
            !write(*,'(A,3E12.4)') 'PROPERTIES : ',(prop_mat(im,i),i=1,3)
            write(*,'(A,I8)')    'DEGREE   : ',type_mat(im)
            write(*,'(A,E12.4)') 'rho      : ',prop_mat(im,1)
            write(*,'(A,E12.4)') 'Vp       : ',((prop_mat(im,2) + 2*prop_mat(im,3))/prop_mat(im,1))**0.5
            write(*,'(A,E12.4)') 'Vs       : ',(prop_mat(im,3)/prop_mat(im,1))**0.5
            !write(*,'(A,E12.4)') 'gamma    : ',prop_mat(im,4)
            write(*,'(A,E12.4)') 'Qs       : ',qs(im)
            write(*,'(A,E12.4)') 'Qp       : ',qp(im)
            write(*,*)
         enddo
      endif
 
      if (mpi_id.eq.0) then                                                                                
          write(*,'(A)')                                                                                
         do im_nle = 1,nmat_nle                                                                                
            write(*,'(A,I8)') 'NON LINEAR ELASTIC - MATERIAL : ',tag_mat_nle(im_nle)                        
            !write(*,'(A,I8)') 'PROPERTIES : ',(prop_mat_nle(im_nle,1))                                        
            write(*,'(A,E12.4)') 'NLE DEPTH : ',(val_mat_nle(im_nle,1))                                        
            write(*,'(A,I8)')    'DEGREE    : ',type_mat_nle(im_nle)                                                
            write(*,*)                                                                                        
         enddo                                                                                                
         if ((fpeak.eq.0).and.(nmat_nle.gt.0)) then                                                        
             write(*,'(A)')'ERROR: Peak frequency for damping not defined!'                                
             call exit(exit_damping_peak)
         endif                                                                                                   
       endif                                                                                                
      
 
      grid_file = grid_file(1:len_trim(grid_file)) // '.mesh'    
      if (mpi_id.eq.0)  write(*,'(A)') '-------------------Reading Grid File-------------------'
      if (mpi_id.eq.0)  write(*,'(A,A20)') 'Grid File : ',grid_file
      
      inquire(file=grid_file,exist=filefound);
      if(filefound .eqv. .false.) call exit(exit_missing_file)
      
 
!counting hexahedras and squares
       call read_dime_filemesh(grid_file,nmat,tag_mat,&
                     nload_dirx_el,tag_dirx_el,nload_diry_el,tag_diry_el,&
                     nload_dirz_el,tag_dirz_el,&
                     nload_neux_el,tag_neux_el,nload_neuy_el,tag_neuy_el,&
                     nload_neuz_el,tag_neuz_el,&
                     nload_neun_el,tag_neun_el,& 
                     nload_abc_el,tag_abc_el,&
                     nload_dg_el,tag_dg_el,&
                     nnod_macro,nelem,nface)
      
      if (mpi_id.eq.0) then
         write(*,'(A,I8)')'Nodes     : ', nnod_macro
         write(*,'(A,I8)')'Elements  : ', nelem
         write(*,'(A,I8)')'Faces     : ', nface
      endif
      
      
      if (nelem.gt.0) then
         allocate (con(nelem,9))
      else
         write(*,*)'Error ! nelem = 0'
         call exit(exit_no_elements)
      endif
      
      if (nface.gt.0) allocate (con_bc(nface,5))
      
            
      call read_filemesh(grid_file,nmat,tag_mat,&
                        nload_dirx_el,tag_dirx_el,nload_diry_el,tag_diry_el,&
                        nload_dirz_el,tag_dirz_el,&
                        nload_neux_el,tag_neux_el,nload_neuy_el,tag_neuy_el,&
                        nload_neuz_el,tag_neuz_el,&
                        nload_neun_el,tag_neun_el,& 
                        nload_abc_el,tag_abc_el,&
                        nload_dg_el,tag_dg_el,&
                        nnod_macro,xx_macro, yy_macro, zz_macro,  &
                        nelem,con,nface,con_bc)
 
      if(mpi_id .eq. 0) write(*,'(A)') 'Read.'
      if(mpi_id .eq. 0) write(*,'(A)')
     
     end subroutine read_input_files