! Copyright 2011 Max-Planck-Institut für Eisenforschung GmbH
!
! This file is part of DAMASK,
! the Düsseldorf Advanced MAterial Simulation Kit.
!
! DAMASK is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! DAMASK 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 General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with DAMASK. If not, see <http://www.gnu.org/licenses/>.
!
!##############################################################
!* $Id$
!##############################################################
MODULE numerics
!##############################################################

use prec, only: pInt, pReal
use IO, only: IO_warning
implicit none

character(len=64), parameter :: numerics_configFile = 'numerics.config' ! name of configuration file
integer(pInt) ::                iJacoStiffness, &                       ! frequency of stiffness update
                                iJacoLpresiduum, &                      ! frequency of Jacobian update of residuum in Lp
                                nHomog, &                               ! homogenization loop limit (only for debugging info, loop limit is determined by "subStepMinHomog")
                                nMPstate, &                             ! materialpoint state loop limit
                                nCryst, &                               ! crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst")
                                nState, &                               ! state loop limit
                                nStress, &                              ! stress loop limit
                                pert_method, &                          ! method used in perturbation technique for tangent
                                numerics_integrationMode                ! integration mode 1 = central solution ; integration mode 2 = perturbation
integer(pInt), dimension(2) ::  numerics_integrator                     ! method used for state integration (central & perturbed state)
real(pReal) ::                  relevantStrain, &                       ! strain increment considered significant (used by crystallite to determine whether strain inc is considered significant)
                                defgradTolerance, &                     ! deviation of deformation gradient that is still allowed (used by CPFEM to determine outdated ffn1)
                                pert_Fg, &                              ! strain perturbation for FEM Jacobi
                                subStepMinCryst, &                      ! minimum (relative) size of sub-step allowed during cutback in crystallite
                                subStepMinHomog, &                      ! minimum (relative) size of sub-step allowed during cutback in homogenization
                                subStepSizeCryst, &                     ! size of first substep when cutback in crystallite
                                subStepSizeHomog, &                     ! size of first substep when cutback in homogenization
                                stepIncreaseCryst, &                    ! increase of next substep size when previous substep converged in crystallite
                                stepIncreaseHomog, &                    ! increase of next substep size when previous substep converged in homogenization
                                rTol_crystalliteState, &                ! relative tolerance in crystallite state loop 
                                rTol_crystalliteTemperature, &          ! relative tolerance in crystallite temperature loop 
                                rTol_crystalliteStress, &               ! relative tolerance in crystallite stress loop
                                aTol_crystalliteStress, &               ! absolute tolerance in crystallite stress loop

!* RGC parameters: added <<<updated 17.12.2009>>>
                                absTol_RGC, &                           ! absolute tolerance of RGC residuum
                                relTol_RGC, &                           ! relative tolerance of RGC residuum
                                absMax_RGC, &                           ! absolute maximum of RGC residuum
                                relMax_RGC, &                           ! relative maximum of RGC residuum
                                pPert_RGC, &                            ! perturbation for computing RGC penalty tangent
                                xSmoo_RGC, &                            ! RGC penalty smoothing parameter (hyperbolic tangent)
                                viscPower_RGC, &                        ! power (sensitivity rate) of numerical viscosity in RGC scheme
                                viscModus_RGC, &                        ! stress modulus of RGC numerical viscosity
                                refRelaxRate_RGC, &                     ! reference relaxation rate in RGC viscosity
                                maxdRelax_RGC, &                        ! threshold of maximum relaxation vector increment (if exceed this then cutback)
                                maxVolDiscr_RGC, &                      ! threshold of maximum volume discrepancy allowed
                                volDiscrMod_RGC, &                      ! stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint)
                                volDiscrPow_RGC, &                      ! powerlaw penalty for volume discrepancy
!* spectral parameters:
                                err_div_tol, &                          ! error of divergence in fourier space
                                err_stress_tolrel, &                    ! factor to multiply with highest stress to get err_stress_tol
                                fftw_timelimit, &                       ! sets the timelimit of plan creation for FFTW, see manual on www.fftw.org
                                rotation_tol                            ! tolerance of rotation specified in loadcase
character(len=64) ::            fftw_planner_string                     ! reads the planing-rigor flag, see manual on www.fftw.org
integer(pInt) ::                fftw_planner_flag                       ! conversion of fftw_planner_string to integer, basically what is usually done in the include file of fftw
logical ::                      memory_efficient,&                      ! for fast execution (pre calculation of gamma_hat)
                                divergence_correction,&                 ! correct divergence calculation in fourier space
                                update_gamma,&                          ! update gamma operator with current stiffness  
                                simplified_algorithm                    ! use short algorithm without fluctuation field
real(pReal) ::                  cut_off_value                           ! percentage of frequencies to cut away
integer(pInt) ::                itmax , &                               ! maximum number of iterations


!* Random seeding parameters
                                fixedSeed                               ! fixed seeding for pseudo-random number generator
!* OpenMP variable
integer(pInt)                   DAMASK_NumThreadsInt                    ! value stored in environment variable DAMASK_NUM_THREADS


CONTAINS
 
!*******************************************
!    initialization subroutine
!*******************************************
subroutine numerics_init()
  
  !*** variables and functions from other modules ***!
  use prec, only:                             pInt, & 
                                              pReal  
  use IO, only:                               IO_error, &
                                              IO_open_file, &
                                              IO_isBlank, &
                                              IO_stringPos, &
                                              IO_stringValue, &
                                              IO_lc, &
                                              IO_floatValue, &
                                              IO_intValue
!$ use OMP_LIB                                                                ! the openMP function library
  
  implicit none

  !*** input variables ***!
  
  !*** output variables ***!
  
  !*** local variables ***!
  integer(pInt), parameter ::                 fileunit = 300  
  integer(pInt), parameter ::                 maxNchunks = 2
  integer(pInt), dimension(1+2*maxNchunks) :: positions
  character(len=64)                           tag
  character(len=1024)                         line
  
! OpenMP variable
!$ character(len=4) DAMASK_NumThreadsString                               !environment variable DAMASK_NUM_THREADS

  !$OMP CRITICAL (write2out)
    write(6,*)
    write(6,*) '<<<+-  numerics init  -+>>>'
    write(6,*) '$Id$'
    write(6,*)
  !$OMP END CRITICAL (write2out)
  
  ! initialize all parameters with standard values
  relevantStrain          = 1.0e-7_pReal
  defgradTolerance        = 1.0e-7_pReal
  iJacoStiffness          = 1_pInt
  iJacoLpresiduum         = 1_pInt
  pert_Fg                 = 1.0e-7_pReal
  pert_method             = 1
  nHomog                  = 20_pInt
  subStepMinHomog         = 1.0e-3_pReal
  subStepSizeHomog        = 0.25
  stepIncreaseHomog       = 1.5
  nMPstate                = 10_pInt
  nCryst                  = 20_pInt
  subStepMinCryst         = 1.0e-3_pReal
  subStepsizeCryst        = 0.25
  stepIncreaseCryst       = 1.5
  nState                  = 10_pInt
  nStress                 = 40_pInt
  rTol_crystalliteState   = 1.0e-6_pReal
  rTol_crystalliteTemperature = 1.0e-6_pReal
  rTol_crystalliteStress  = 1.0e-6_pReal
  aTol_crystalliteStress  = 1.0e-8_pReal            ! residuum is in Lp (hence strain on the order of 1e-8 here)
  numerics_integrator(1)  = 1                       ! fix-point iteration
  numerics_integrator(2)  = 1                       ! fix-point iteration
  
!* RGC parameters: added <<<updated 17.12.2009>>> with moderate setting
  absTol_RGC              = 1.0e+4
  relTol_RGC              = 1.0e-3
  absMax_RGC              = 1.0e+10
  relMax_RGC              = 1.0e+2
  pPert_RGC               = 1.0e-7
  xSmoo_RGC               = 1.0e-5
  viscPower_RGC           = 1.0e+0  ! Newton viscosity (linear model)
  viscModus_RGC           = 0.0e+0  ! No viscosity is applied
  refRelaxRate_RGC        = 1.0e-3
  maxdRelax_RGC           = 1.0e+0
  maxVolDiscr_RGC         = 1.0e-5  ! tolerance for volume discrepancy allowed
  volDiscrMod_RGC         = 1.0e+12
  volDiscrPow_RGC         = 5.0

!* spectral parameters:  
  err_div_tol             = 1.0e-4       ! 1.0e-4 proposed by Suquet
  err_stress_tolrel       = 0.01         ! relative tolerance for fullfillment of stress BC (1% of maximum stress)
  itmax                   = 20_pInt      ! Maximum iteration number
  memory_efficient        = .true.       ! Precalculate Gamma-operator (81 double per point)
  fftw_timelimit          = -1.0_pReal   ! no timelimit of plan creation for FFTW
  fftw_planner_string     ='FFTW_PATIENT'
  rotation_tol            = 1.0e-12 
  divergence_correction   = .true.       ! correct divergence by empirical factor
  simplified_algorithm    = .true.       ! use algorithm without fluctuation field
  update_gamma            = .false.      ! do not update gamma operator with current stiffness
  cut_off_value           = 0.0_pReal    ! use all frequencies

!* Random seeding parameters
  fixedSeed               = 0_pInt


!* determin number of threads from environment variable DAMASK_NUM_THREADS
  DAMASK_NumThreadsInt = 0_pInt
!$ call GetEnv('DAMASK_NUM_THREADS',DAMASK_NumThreadsString)                   ! get environment variable DAMASK_NUM_THREADS...
!$ read(DAMASK_NumThreadsString,'(i4)') DAMASK_NumThreadsInt                   ! ...convert it to integer...
!$ if (DAMASK_NumThreadsInt < 1) DAMASK_NumThreadsInt = 1                      ! ...ensure that its at least one...
!$ call omp_set_num_threads(DAMASK_NumThreadsInt)                              ! ...and use it as number of threads for parallel execution

  ! try to open the config file
  if(IO_open_file(fileunit,numerics_configFile)) then 
  
    !$OMP CRITICAL (write2out)
      write(6,*) '   ... using values from config file'
      write(6,*)
    !$OMP END CRITICAL (write2out)
    
    line = ''
    ! read variables from config file and overwrite parameters
    do
      read(fileunit,'(a1024)',END=100) line
      if (IO_isBlank(line)) cycle                           ! skip empty lines
      positions = IO_stringPos(line,maxNchunks)
      tag = IO_lc(IO_stringValue(line,positions,1))         ! extract key
      select case(tag)
        case ('relevantstrain')
              relevantStrain = IO_floatValue(line,positions,2)
        case ('defgradtolerance')
              defgradTolerance = IO_floatValue(line,positions,2)
        case ('ijacostiffness')
              iJacoStiffness = IO_intValue(line,positions,2)
        case ('ijacolpresiduum')
              iJacoLpresiduum = IO_intValue(line,positions,2)
        case ('pert_fg')
              pert_Fg = IO_floatValue(line,positions,2)
        case ('pert_method')
              pert_method = IO_intValue(line,positions,2)
        case ('nhomog')
              nHomog = IO_intValue(line,positions,2)
        case ('nmpstate')
              nMPstate = IO_intValue(line,positions,2)
        case ('ncryst')
              nCryst = IO_intValue(line,positions,2)
        case ('nstate')
              nState = IO_intValue(line,positions,2)
        case ('nstress')
              nStress = IO_intValue(line,positions,2)
        case ('substepmincryst')
              subStepMinCryst = IO_floatValue(line,positions,2)
        case ('substepsizecryst')
              subStepSizeCryst = IO_floatValue(line,positions,2)
        case ('stepincreasecryst')
              stepIncreaseCryst = IO_floatValue(line,positions,2)
        case ('substepminhomog')
              subStepMinHomog = IO_floatValue(line,positions,2)
        case ('substepsizehomog')
              subStepSizeHomog = IO_floatValue(line,positions,2)
        case ('stepincreasehomog')
              stepIncreaseHomog = IO_floatValue(line,positions,2)
        case ('rtol_crystallitestate')
              rTol_crystalliteState = IO_floatValue(line,positions,2)
        case ('rtol_crystallitetemperature')
              rTol_crystalliteTemperature = IO_floatValue(line,positions,2)
        case ('rtol_crystallitestress')
              rTol_crystalliteStress = IO_floatValue(line,positions,2)
        case ('atol_crystallitestress')
              aTol_crystalliteStress = IO_floatValue(line,positions,2)
        case ('integrator')
              numerics_integrator(1) = IO_intValue(line,positions,2)
        case ('integratorstiffness')
              numerics_integrator(2) = IO_intValue(line,positions,2)

!* RGC parameters: 
        case ('atol_rgc')
              absTol_RGC = IO_floatValue(line,positions,2)
        case ('rtol_rgc')
              relTol_RGC = IO_floatValue(line,positions,2)
        case ('amax_rgc')
              absMax_RGC = IO_floatValue(line,positions,2)
        case ('rmax_rgc')
              relMax_RGC = IO_floatValue(line,positions,2)
        case ('perturbpenalty_rgc')
              pPert_RGC = IO_floatValue(line,positions,2)
        case ('relevantmismatch_rgc')
              xSmoo_RGC = IO_floatValue(line,positions,2)
        case ('viscositypower_rgc')
              viscPower_RGC = IO_floatValue(line,positions,2)
        case ('viscositymodulus_rgc')
              viscModus_RGC = IO_floatValue(line,positions,2)
        case ('refrelaxationrate_rgc')
              refRelaxRate_RGC = IO_floatValue(line,positions,2)
        case ('maxrelaxation_rgc')
              maxdRelax_RGC = IO_floatValue(line,positions,2)
        case ('maxvoldiscrepancy_rgc')
              maxVolDiscr_RGC = IO_floatValue(line,positions,2)
        case ('voldiscrepancymod_rgc')
              volDiscrMod_RGC = IO_floatValue(line,positions,2)
        case ('discrepancypower_rgc')
              volDiscrPow_RGC = IO_floatValue(line,positions,2)

!* spectral parameters
        case ('err_div_tol')
              err_div_tol = IO_floatValue(line,positions,2)
        case ('err_stress_tolrel')
              err_stress_tolrel = IO_floatValue(line,positions,2)
        case ('itmax')
              itmax = IO_intValue(line,positions,2)
        case ('memory_efficient')
              memory_efficient = IO_intValue(line,positions,2)  > 0_pInt
        case ('fftw_timelimit')
              fftw_timelimit = IO_floatValue(line,positions,2)
        case ('fftw_planner_string')
              fftw_planner_string = IO_stringValue(line,positions,2)
        case ('rotation_tol')
              rotation_tol = IO_floatValue(line,positions,2)
        case ('divergence_correction')
              divergence_correction = IO_intValue(line,positions,2)  > 0_pInt
        case ('update_gamma')
              update_gamma = IO_intValue(line,positions,2)  > 0_pInt
        case ('simplified_algorithm')
              simplified_algorithm = IO_intValue(line,positions,2)  > 0_pInt
        case ('cut_off_value')
              cut_off_value = IO_floatValue(line,positions,2)

!* Random seeding parameters
        case ('fixed_seed')
              fixedSeed = IO_intValue(line,positions,2)
      endselect
    enddo
    100 close(fileunit)
  
  ! no config file, so we use standard values
  else 
  
    !$OMP CRITICAL (write2out)
      write(6,*) '   ... using standard values'
      write(6,*)
    !$OMP END CRITICAL (write2out)
    
  endif
  select case(IO_lc(fftw_planner_string))                        ! setting parameters for the plan creation of FFTW. Basically a translation from fftw3.f
    case('estimate','fftw_estimate')                             ! ordered from slow execution (but fast plan creation) to fast execution
       fftw_planner_flag = 64
     case('measure','fftw_measure')
       fftw_planner_flag = 0
     case('patient','fftw_patient')
       fftw_planner_flag= 32
     case('exhaustive','fftw_exhaustive')
       fftw_planner_flag = 8 
     case default
       call IO_warning(warning_ID=47_pInt,ext_msg=trim(IO_lc(fftw_planner_string)))
       fftw_planner_flag = 32
  end select

  ! writing parameters to output file
  !$OMP CRITICAL (write2out)
    write(6,'(a24,x,e8.1)') ' relevantStrain:         ',relevantStrain
    write(6,'(a24,x,e8.1)') ' defgradTolerance:       ',defgradTolerance
    write(6,'(a24,x,i8)')   ' iJacoStiffness:         ',iJacoStiffness
    write(6,'(a24,x,i8)')   ' iJacoLpresiduum:        ',iJacoLpresiduum
    write(6,'(a24,x,e8.1)') ' pert_Fg:                ',pert_Fg
    write(6,'(a24,x,i8)')   ' pert_method:            ',pert_method
    write(6,'(a24,x,i8)')   ' nCryst:                 ',nCryst
    write(6,'(a24,x,e8.1)') ' subStepMinCryst:        ',subStepMinCryst
    write(6,'(a24,x,e8.1)') ' subStepSizeCryst:       ',subStepSizeCryst
    write(6,'(a24,x,e8.1)') ' stepIncreaseCryst:      ',stepIncreaseCryst
    write(6,'(a24,x,i8)')   ' nState:                 ',nState
    write(6,'(a24,x,i8)')   ' nStress:                ',nStress
    write(6,'(a24,x,e8.1)') ' rTol_crystalliteState:  ',rTol_crystalliteState
    write(6,'(a24,x,e8.1)') ' rTol_crystalliteTemp:   ',rTol_crystalliteTemperature
    write(6,'(a24,x,e8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress
    write(6,'(a24,x,e8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress
    write(6,'(a24,2(x,i8),/)')' integrator:             ',numerics_integrator
  
    write(6,'(a24,x,i8)')   ' nHomog:                 ',nHomog
    write(6,'(a24,x,e8.1)') ' subStepMinHomog:        ',subStepMinHomog
    write(6,'(a24,x,e8.1)') ' subStepSizeHomog:       ',subStepSizeHomog
    write(6,'(a24,x,e8.1)') ' stepIncreaseHomog:      ',stepIncreaseHomog
    write(6,'(a24,x,i8,/)') ' nMPstate:               ',nMPstate

!* RGC parameters
    write(6,'(a24,x,e8.1)') ' aTol_RGC:               ',absTol_RGC
    write(6,'(a24,x,e8.1)') ' rTol_RGC:               ',relTol_RGC
    write(6,'(a24,x,e8.1)') ' aMax_RGC:               ',absMax_RGC
    write(6,'(a24,x,e8.1)') ' rMax_RGC:               ',relMax_RGC
    write(6,'(a24,x,e8.1)') ' perturbPenalty_RGC:     ',pPert_RGC
    write(6,'(a24,x,e8.1)') ' relevantMismatch_RGC:   ',xSmoo_RGC
    write(6,'(a24,x,e8.1)') ' viscosityrate_RGC:      ',viscPower_RGC
    write(6,'(a24,x,e8.1)') ' viscositymodulus_RGC:   ',viscModus_RGC
    write(6,'(a24,x,e8.1)') ' maxrelaxation_RGC:      ',maxdRelax_RGC
    write(6,'(a24,x,e8.1)') ' maxVolDiscrepancy_RGC:  ',maxVolDiscr_RGC
    write(6,'(a24,x,e8.1)') ' volDiscrepancyMod_RGC:  ',volDiscrMod_RGC
    write(6,'(a24,x,e8.1,/)') ' discrepancyPower_RGC:   ',volDiscrPow_RGC

!* spectral parameters
    write(6,'(a24,x,e8.1)')   ' err_div_tol:            ',err_div_tol
    write(6,'(a24,x,e8.1)')   ' err_stress_tolrel:      ',err_stress_tolrel
    write(6,'(a24,x,i8)')     ' itmax:                  ',itmax
    write(6,'(a24,x,L8)')     ' memory_efficient:       ',memory_efficient
    if(fftw_timelimit<0.0_pReal) then
      write(6,'(a24,x,L8)')   ' fftw_timelimit:         ',.false.
    else    
      write(6,'(a24,x,e8.1)') ' fftw_timelimit:         ',fftw_timelimit
    endif
    write(6,'(a24,x,a)')      ' fftw_planner_string:    ',trim(fftw_planner_string)
    write(6,'(a24,x,i8)')     ' fftw_planner_flag:      ',fftw_planner_flag
    write(6,'(a24,x,e8.1)')   ' rotation_tol:           ',rotation_tol
    write(6,'(a24,x,L8,/)')   ' divergence_correction:  ',divergence_correction
    write(6,'(a24,x,L8,/)')   ' update_gamma:           ',update_gamma
    write(6,'(a24,x,L8,/)')   ' simplified_algorithm:   ',simplified_algorithm
    write(6,'(a24,x,e8.1)')   ' cut_off_value:          ',cut_off_value
!* Random seeding parameters
    write(6,'(a24,x,i16,/)')   ' fixed_seed:             ',fixedSeed
  !$OMP END CRITICAL (write2out)

!* openMP parameter
!$  write(6,'(a24,x,i8,/)')   ' number of threads:      ',DAMASK_NumThreadsInt
  
  ! sanity check  
  if (relevantStrain <= 0.0_pReal)          call IO_error(260)
  if (defgradTolerance <= 0.0_pReal)        call IO_error(294)
  if (iJacoStiffness < 1_pInt)              call IO_error(261)
  if (iJacoLpresiduum < 1_pInt)             call IO_error(262)
  if (pert_Fg <= 0.0_pReal)                 call IO_error(263)
  if (pert_method <= 0_pInt .or. pert_method >= 4_pInt) &
                                            call IO_error(299)
  if (nHomog < 1_pInt)                      call IO_error(264)
  if (nMPstate < 1_pInt)                    call IO_error(279)  !! missing in IO !!
  if (nCryst < 1_pInt)                      call IO_error(265)
  if (nState < 1_pInt)                      call IO_error(266)
  if (nStress < 1_pInt)                     call IO_error(267)
  if (subStepMinCryst <= 0.0_pReal)         call IO_error(268)
  if (subStepSizeCryst <= 0.0_pReal)        call IO_error(268)
  if (stepIncreaseCryst <= 0.0_pReal)       call IO_error(268)
  if (subStepMinHomog <= 0.0_pReal)         call IO_error(268)
  if (subStepSizeHomog <= 0.0_pReal)        call IO_error(268)
  if (stepIncreaseHomog <= 0.0_pReal)       call IO_error(268)
  if (rTol_crystalliteState <= 0.0_pReal)   call IO_error(269)
  if (rTol_crystalliteTemperature <= 0.0_pReal) call IO_error(276) !! oops !!
  if (rTol_crystalliteStress <= 0.0_pReal)  call IO_error(270)
  if (aTol_crystalliteStress <= 0.0_pReal)  call IO_error(271)
  if (any(numerics_integrator <= 0_pInt) .or. any(numerics_integrator >= 6_pInt)) &
                                            call IO_error(298)

!* RGC parameters: added <<<updated 17.11.2009>>>
  if (absTol_RGC <= 0.0_pReal)              call IO_error(272)
  if (relTol_RGC <= 0.0_pReal)              call IO_error(273)
  if (absMax_RGC <= 0.0_pReal)              call IO_error(274)
  if (relMax_RGC <= 0.0_pReal)              call IO_error(275)
  if (pPert_RGC <= 0.0_pReal)               call IO_error(276)   !! oops !!
  if (xSmoo_RGC <= 0.0_pReal)               call IO_error(277)
  if (viscPower_RGC < 0.0_pReal)            call IO_error(278)
  if (viscModus_RGC < 0.0_pReal)            call IO_error(278)
  if (refRelaxRate_RGC <= 0.0_pReal)        call IO_error(278)
  if (maxdRelax_RGC <= 0.0_pReal)           call IO_error(288)
  if (maxVolDiscr_RGC <= 0.0_pReal)         call IO_error(289)
  if (volDiscrMod_RGC < 0.0_pReal)          call IO_error(289)
  if (volDiscrPow_RGC <= 0.0_pReal)         call IO_error(289)

!* spectral parameters
  if (err_div_tol <= 0.0_pReal)             call IO_error(49)
  if (err_stress_tolrel <= 0.0_pReal)       call IO_error(49)
  if (itmax <= 1.0_pInt)                    call IO_error(49)
  
  if (fixedSeed <= 0_pInt) then
    !$OMP CRITICAL (write2out)
      write(6,'(a)') 'Random is random!'
    !$OMP END CRITICAL (write2out)
  endif
endsubroutine

END MODULE numerics