Merge remote-tracking branch 'origin/development' into long-YAML-files

VERSION

@@ -1 +1 @@
-3.0.0-alpha7-395-gb4e706f09
+3.0.0-alpha7-414-g23f896631

python/damask/grid_filters.py

@@ -564,17 +564,10 @@ def unravel_index(idx: _np.ndarray) -> _np.ndarray:
     >>> seq = np.arange(6).reshape((3,2,1),order='F')
     >>> (coord_idx := damask.grid_filters.unravel_index(seq))
     array([[[[0, 0, 0]],
-
             [[0, 1, 0]]],
-
-
            [[[1, 0, 0]],
-
             [[1, 1, 0]]],
-
-
            [[[2, 0, 0]],
-
             [[2, 1, 0]]]])
     >>> coord_idx[1,1,0]
     array([1, 1, 0])
@@ -608,17 +601,12 @@ def ravel_index(idx: _np.ndarray) -> _np.ndarray:
     >>> import damask
     >>> (rev := np.array([[1,1,0],[0,1,0],[1,0,0],[0,0,0]]).reshape((2,2,1,3)))
     array([[[[1, 1, 0]],
-
             [[0, 1, 0]]],
-
-
            [[[1, 0, 0]],
-
             [[0, 0, 0]]]])
     >>> (flat_idx := damask.grid_filters.ravel_index(rev))
     array([[[3],
             [2]],
-
            [[1],
             [0]]])
 

python/damask/mechanics.py

@@ -14,7 +14,7 @@ from . import _rotation
 
 
 def deformation_Cauchy_Green_left(F: _np.ndarray) -> _np.ndarray:
-    """
+    r"""
     Calculate left Cauchy-Green deformation tensor (Finger deformation tensor).
 
     Parameters
@@ -27,12 +27,18 @@ def deformation_Cauchy_Green_left(F: _np.ndarray) -> _np.ndarray:
     B : numpy.ndarray, shape (...,3,3)
         Left Cauchy-Green deformation tensor.
 
+    Notes
+    -----
+    .. math::
+
+       \vb{B} = \vb{F} \vb{F}^\text{T}
+
     """
     return _np.matmul(F,_tensor.transpose(F))
 
 
 def deformation_Cauchy_Green_right(F: _np.ndarray) -> _np.ndarray:
-    """
+    r"""
     Calculate right Cauchy-Green deformation tensor.
 
     Parameters
@@ -45,12 +51,18 @@ def deformation_Cauchy_Green_right(F: _np.ndarray) -> _np.ndarray:
     C : numpy.ndarray, shape (...,3,3)
         Right Cauchy-Green deformation tensor.
 
+    Notes
+    -----
+    .. math::
+
+       \vb{C} = \vb{F}^\text{T} \vb{F}
+
     """
     return _np.matmul(_tensor.transpose(F),F)
 
 
 def equivalent_strain_Mises(epsilon: _np.ndarray) -> _np.ndarray:
-    """
+    r"""
     Calculate the Mises equivalent of a strain tensor.
 
     Parameters
@@ -63,12 +75,23 @@ def equivalent_strain_Mises(epsilon: _np.ndarray) -> _np.ndarray:
     epsilon_vM : numpy.ndarray, shape (...)
         Von Mises equivalent strain of epsilon.
 
+    Notes
+    -----
+    The von Mises equivalent of a strain tensor is defined as:
+
+    .. math::
+
+       \epsilon_\text{vM} = \sqrt{2/3 \epsilon^\prime_{ij} \epsilon^\prime_{ij}}
+
+    where :math:`\vb*{\epsilon}^\prime` is the deviatoric part
+    of the strain tensor.
+
     """
     return _equivalent_Mises(epsilon,2.0/3.0)
 
 
 def equivalent_stress_Mises(sigma: _np.ndarray) -> _np.ndarray:
-    """
+    r"""
     Calculate the Mises equivalent of a stress tensor.
 
     Parameters
@@ -81,6 +104,17 @@ def equivalent_stress_Mises(sigma: _np.ndarray) -> _np.ndarray:
     sigma_vM : numpy.ndarray, shape (...)
         Von Mises equivalent stress of sigma.
 
+    Notes
+    -----
+    The von Mises equivalent of a stress tensor is defined as:
+
+    .. math::
+
+       \sigma_\text{vM} = \sqrt{3/2 \sigma^\prime_{ij} \sigma^\prime_{ij}}
+
+    where :math:`\vb*{\sigma}^\prime` is the deviatoric part
+    of the stress tensor.
+
     """
     return _equivalent_Mises(sigma,3.0/2.0)
 
@@ -105,7 +139,7 @@ def maximum_shear(T_sym: _np.ndarray) -> _np.ndarray:
 
 
 def rotation(T: _np.ndarray) -> _rotation.Rotation:
-    """
+    r"""
     Calculate the rotational part of a tensor.
 
     Parameters
@@ -118,6 +152,17 @@ def rotation(T: _np.ndarray) -> _rotation.Rotation:
     R : damask.Rotation, shape (...)
         Rotational part of the vector.
 
+    Notes
+    -----
+    The rotational part is calculated from the polar decomposition:
+
+    .. math::
+
+       \vb{R} = \vb{T} \vb{U}^{-1} = \vb{V}^{-1} \vb{T}
+
+    where :math:`\vb{V}` and :math:`\vb{U}` are the left
+    and right stretch tensor, respectively.
+
     """
     return _rotation.Rotation.from_matrix(_polar_decomposition(T,'R')[0])
 
@@ -212,7 +257,7 @@ def stress_second_Piola_Kirchhoff(P: _np.ndarray,
 
 
 def stretch_left(T: _np.ndarray) -> _np.ndarray:
-    """
+    r"""
     Calculate left stretch of a tensor.
 
     Parameters
@@ -225,12 +270,23 @@ def stretch_left(T: _np.ndarray) -> _np.ndarray:
     V : numpy.ndarray, shape (...,3,3)
         Left stretch tensor from Polar decomposition of T.
 
+    Notes
+    -----
+    The left stretch tensor is calculated from the
+    polar decomposition:
+
+    .. math::
+
+       \vb{V} = \vb{T} \vb{R}^\text{T}
+
+    where :math:`\vb{R}` is a rotation.
+
     """
     return _polar_decomposition(T,'V')[0]
 
 
 def stretch_right(T: _np.ndarray) -> _np.ndarray:
-    """
+    r"""
     Calculate right stretch of a tensor.
 
     Parameters
@@ -243,6 +299,17 @@ def stretch_right(T: _np.ndarray) -> _np.ndarray:
     U : numpy.ndarray, shape (...,3,3)
         Left stretch tensor from Polar decomposition of T.
 
+    Notes
+    -----
+    The right stretch tensor is calculated from the
+    polar decomposition:
+
+    .. math::
+
+       \vb{U} = \vb{R}^\text{T} \vb{T}
+
+    where :math:`\vb{R}` is a rotation.
+
     """
     return _polar_decomposition(T,'U')[0]
 
@@ -260,6 +327,11 @@ def _polar_decomposition(T: _np.ndarray,
         Requested outputs: ‘R’ for the rotation tensor,
         ‘V’ for left stretch tensor, and ‘U’ for right stretch tensor.
 
+    Returns
+    -------
+    VRU : tuple of numpy.ndarray, shape (...,3,3)
+       Requested components of the singular value decomposition.
+
     """
     u, _, vh = _np.linalg.svd(T)
     R = _np.einsum('...ij,...jk',u,vh)
@@ -290,6 +362,11 @@ def _equivalent_Mises(T_sym: _np.ndarray,
     s : float
         Scaling factor (2/3 for strain, 3/2 for stress).
 
+    Returns
+    -------
+    eq : numpy.ndarray, shape (...)
+        Scaled second invariant of the deviatoric part of T_sym.
+
     """
     d = _tensor.deviatoric(T_sym)
     return _np.sqrt(s*_np.sum(d**2.0,axis=(-1,-2)))

python/setup.cfg

@@ -26,4 +26,3 @@ install_requires =
     vtk>=8.1
     matplotlib>=3.0                                                                                 # requires numpy, pillow
     pyyaml>=3.12
-setup_requires = setuptools

src/HDF5_utilities.f90

@@ -3,6 +3,7 @@
 !> @author Yi-Chin Yang, Max-Planck-Institut für Eisenforschung GmbH
 !> @author Jennifer Nastola, Max-Planck-Institut für Eisenforschung GmbH
 !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Philip Eisenlohr, Michigan State University
 !--------------------------------------------------------------------------------------------------
 module HDF5_utilities
   use HDF5
@@ -17,6 +18,7 @@ module HDF5_utilities
 
   use prec
   use parallelization
+  use misc
 
 #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
   implicit none(type,external)
@@ -26,7 +28,7 @@ module HDF5_utilities
   private
 
 !--------------------------------------------------------------------------------------------------
-!> @brief reads integer or float data of defined shape from file
+!> @brief Read integer or float data of defined shape from file.
 !> @details for parallel IO, all dimension except for the last need to match
 !--------------------------------------------------------------------------------------------------
   interface HDF5_read
@@ -48,7 +50,7 @@ module HDF5_utilities
   end interface HDF5_read
 
 !--------------------------------------------------------------------------------------------------
-!> @brief writes integer or real data of defined shape to file
+!> @brief Write integer or real data of defined shape to file.
 !> @details for parallel IO, all dimension except for the last need to match
 !--------------------------------------------------------------------------------------------------
   interface HDF5_write
@@ -74,7 +76,7 @@ module HDF5_utilities
   end interface HDF5_write
 
 !--------------------------------------------------------------------------------------------------
-!> @brief attached attributes of type char, integer or real to a file/dataset/group
+!> @brief Attach attribute of type char, integer, or real to a file/dataset/group.
 !--------------------------------------------------------------------------------------------------
   interface HDF5_addAttribute
     module procedure HDF5_addAttribute_str
@@ -110,7 +112,7 @@ contains
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief initialize HDF5 libary and do sanity checks
+!> @brief Initialize HDF5 libary and perform sanity checks.
 !--------------------------------------------------------------------------------------------------
 subroutine HDF5_utilities_init
 
@@ -176,31 +178,22 @@ integer(HID_T) function HDF5_openFile(fileName,mode,parallel)
   integer :: hdferr
 
 
-  if (present(mode)) then
+  m = misc_optional(mode,'r')
-    m = mode
-  else
-    m = 'r'
-  end if
 
   call H5Pcreate_f(H5P_FILE_ACCESS_F, plist_id, hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
 
 #ifdef PETSC
-  if (present(parallel)) then
+  if (misc_optional(parallel,.true.)) &
 #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
-    if (parallel) call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL_F90, hdferr)
-  else
     call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL_F90, hdferr)
 #else
-    if (parallel) call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr)
+    call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL,     hdferr)
-  else
-    call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr)
 #endif
-  end if
   if (hdferr < 0) error stop 'HDF5 error'
 #endif
 
-  if    (m == 'w') then
+  if     (m == 'w') then
     call H5Fcreate_f(fileName,H5F_ACC_TRUNC_F,HDF5_openFile,hdferr,access_prp = plist_id)
     if (hdferr < 0) error stop 'HDF5 error'
   elseif (m == 'a') then
@@ -220,7 +213,7 @@ end function HDF5_openFile
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief close the opened HDF5 output file
+!> @brief Close an open HDF5 file.
 !--------------------------------------------------------------------------------------------------
 subroutine HDF5_closeFile(fileHandle)
 
@@ -235,7 +228,7 @@ end subroutine HDF5_closeFile
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief adds a new group to the fileHandle
+!> @brief Add a new group to the file.
 !--------------------------------------------------------------------------------------------------
 integer(HID_T) function HDF5_addGroup(fileHandle,groupName)
 
@@ -268,7 +261,7 @@ end function HDF5_addGroup
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief open an existing group of a file
+!> @brief Open an existing group in a file.
 !--------------------------------------------------------------------------------------------------
 integer(HID_T) function HDF5_openGroup(fileHandle,groupName)
 
@@ -304,7 +297,7 @@ end function HDF5_openGroup
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief close a group
+!> @brief Close a group.
 !--------------------------------------------------------------------------------------------------
 subroutine HDF5_closeGroup(group_id)
 
@@ -330,11 +323,7 @@ logical function HDF5_objectExists(loc_id,path)
   character(len=:), allocatable :: p
 
 
-  if (present(path)) then
+  p = trim(misc_optional(path,'.'))
-    p = trim(path)
-  else
-    p = '.'
-  end if
 
   call H5Lexists_f(loc_id, p, HDF5_objectExists, hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
@@ -364,11 +353,7 @@ subroutine HDF5_addAttribute_str(loc_id,attrLabel,attrValue,path)
   type(C_PTR), target, dimension(1) :: ptr
 
 
-  if (present(path)) then
+  p = trim(misc_optional(path,'.'))
-    p = trim(path)
-  else
-    p = '.'
-  end if
 
   attrValue_(1) = trim(attrValue)//C_NULL_CHAR
   ptr(1) = c_loc(attrValue_(1))
@@ -376,14 +361,14 @@ subroutine HDF5_addAttribute_str(loc_id,attrLabel,attrValue,path)
   call H5Screate_f(H5S_SCALAR_F,space_id,hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
 
-  call H5Aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr)
+  call H5Aexists_by_name_f(loc_id,p,attrLabel,attrExists,hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
   if (attrExists) then
-    call H5Adelete_by_name_f(loc_id, trim(p), attrLabel, hdferr)
+    call H5Adelete_by_name_f(loc_id,p,attrLabel,hdferr)
     if (hdferr < 0) error stop 'HDF5 error'
   end if
 
-  call H5Acreate_by_name_f(loc_id,trim(p),trim(attrLabel),H5T_STRING,space_id,attr_id,hdferr)
+  call H5Acreate_by_name_f(loc_id,p,trim(attrLabel),H5T_STRING,space_id,attr_id,hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
   call H5Awrite_f(attr_id, H5T_STRING, c_loc(ptr), hdferr)                                          ! ptr instead of c_loc(ptr) works on gfortran, not on ifort
   if (hdferr < 0) error stop 'HDF5 error'
@@ -412,23 +397,19 @@ subroutine HDF5_addAttribute_int(loc_id,attrLabel,attrValue,path)
   character(len=:), allocatable :: p
 
 
-  if (present(path)) then
+  p = trim(misc_optional(path,'.'))
-    p = trim(path)
-  else
-    p = '.'
-  end if
 
   call H5Screate_f(H5S_SCALAR_F,space_id,hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
 
-  call H5Aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr)
+  call H5Aexists_by_name_f(loc_id,p,attrLabel,attrExists,hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
   if (attrExists) then
-    call H5Adelete_by_name_f(loc_id, trim(p), attrLabel, hdferr)
+    call H5Adelete_by_name_f(loc_id,p,attrLabel,hdferr)
     if (hdferr < 0) error stop 'HDF5 error'
   end if
 
-  call H5Acreate_by_name_f(loc_id,trim(p),trim(attrLabel),H5T_NATIVE_INTEGER,space_id,attr_id,hdferr)
+  call H5Acreate_by_name_f(loc_id,p,trim(attrLabel),H5T_NATIVE_INTEGER,space_id,attr_id,hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
   call H5Awrite_f(attr_id, H5T_NATIVE_INTEGER, attrValue, int([1],HSIZE_T), hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
@@ -457,23 +438,19 @@ subroutine HDF5_addAttribute_real(loc_id,attrLabel,attrValue,path)
   character(len=:), allocatable :: p
 
 
-  if (present(path)) then
+  p = trim(misc_optional(path,'.'))
-    p = trim(path)
-  else
-    p = '.'
-  end if
 
   call H5Screate_f(H5S_SCALAR_F,space_id,hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
 
-  call H5Aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr)
+  call H5Aexists_by_name_f(loc_id,p,attrLabel,attrExists,hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
   if (attrExists) then
-    call H5Adelete_by_name_f(loc_id, trim(p), attrLabel, hdferr)
+    call H5Adelete_by_name_f(loc_id,p,attrLabel,hdferr)
     if (hdferr < 0) error stop 'HDF5 error'
   end if
 
-  call H5Acreate_by_name_f(loc_id,trim(p),trim(attrLabel),H5T_NATIVE_DOUBLE,space_id,attr_id,hdferr)
+  call H5Acreate_by_name_f(loc_id,p,trim(attrLabel),H5T_NATIVE_DOUBLE,space_id,attr_id,hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
   call H5Awrite_f(attr_id, H5T_NATIVE_DOUBLE, attrValue, int([1],HSIZE_T), hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
@@ -504,11 +481,7 @@ subroutine HDF5_addAttribute_str_array(loc_id,attrLabel,attrValue,path)
   type(C_PTR), target, dimension(size(attrValue))  :: ptr
 
 
-  if (present(path)) then
+  p = trim(misc_optional(path,'.'))
-    p = trim(path)
-  else
-    p = '.'
-  end if
 
   do i=1,size(attrValue)
     attrValue_(i) = attrValue(i)//C_NULL_CHAR
@@ -518,14 +491,14 @@ subroutine HDF5_addAttribute_str_array(loc_id,attrLabel,attrValue,path)
   call H5Screate_simple_f(1,shape(attrValue_,kind=HSIZE_T),space_id,hdferr,shape(attrValue_,kind=HSIZE_T))
   if (hdferr < 0) error stop 'HDF5 error'
 
-  call H5Aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr)
+  call H5Aexists_by_name_f(loc_id,p,attrLabel,attrExists,hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
   if (attrExists) then
-    call H5Adelete_by_name_f(loc_id, trim(p), attrLabel, hdferr)
+    call H5Adelete_by_name_f(loc_id,p,attrLabel,hdferr)
     if (hdferr < 0) error stop 'HDF5 error'
   end if
 
-  call H5Acreate_by_name_f(loc_id,trim(p),trim(attrLabel),H5T_STRING,space_id,attr_id,hdferr)
+  call H5Acreate_by_name_f(loc_id,p,trim(attrLabel),H5T_STRING,space_id,attr_id,hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
   call H5Awrite_f(attr_id, H5T_STRING, c_loc(ptr), hdferr)                                          ! ptr instead of c_loc(ptr) works on gfortran, not on ifort
   if (hdferr < 0) error stop 'HDF5 error'
@@ -555,25 +528,21 @@ subroutine HDF5_addAttribute_int_array(loc_id,attrLabel,attrValue,path)
   character(len=:), allocatable :: p
 
 
-  if (present(path)) then
+  p = trim(misc_optional(path,'.'))
-    p = trim(path)
-  else
-    p = '.'
-  end if
 
   array_size = size(attrValue,kind=HSIZE_T)
 
   call H5Screate_simple_f(1, array_size, space_id, hdferr, array_size)
   if (hdferr < 0) error stop 'HDF5 error'
 
-  call H5Aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr)
+  call H5Aexists_by_name_f(loc_id,p,attrLabel,attrExists,hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
   if (attrExists) then
-    call H5Adelete_by_name_f(loc_id, trim(p), attrLabel, hdferr)
+    call H5Adelete_by_name_f(loc_id,p,attrLabel,hdferr)
     if (hdferr < 0) error stop 'HDF5 error'
   end if
 
-  call H5Acreate_by_name_f(loc_id,trim(p),trim(attrLabel),H5T_NATIVE_INTEGER,space_id,attr_id,hdferr)
+  call H5Acreate_by_name_f(loc_id,p,trim(attrLabel),H5T_NATIVE_INTEGER,space_id,attr_id,hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
   call H5Awrite_f(attr_id, H5T_NATIVE_INTEGER, attrValue, array_size, hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
@@ -603,25 +572,21 @@ subroutine HDF5_addAttribute_real_array(loc_id,attrLabel,attrValue,path)
   character(len=:), allocatable :: p
 
 
-  if (present(path)) then
+  p = trim(misc_optional(path,'.'))
-    p = trim(path)
-  else
-    p = '.'
-  end if
 
   array_size = size(attrValue,kind=HSIZE_T)
 
   call H5Screate_simple_f(1, array_size, space_id, hdferr, array_size)
   if (hdferr < 0) error stop 'HDF5 error'
 
-  call H5Aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr)
+  call H5Aexists_by_name_f(loc_id,p,attrLabel,attrExists,hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
   if (attrExists) then
-    call H5Adelete_by_name_f(loc_id, trim(p), attrLabel, hdferr)
+    call H5Adelete_by_name_f(loc_id,p,attrLabel,hdferr)
     if (hdferr < 0) error stop 'HDF5 error'
   end if
 
-  call H5Acreate_by_name_f(loc_id,trim(p),trim(attrLabel),H5T_NATIVE_DOUBLE,space_id,attr_id,hdferr)
+  call H5Acreate_by_name_f(loc_id,p,trim(attrLabel),H5T_NATIVE_DOUBLE,space_id,attr_id,hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
   call H5Awrite_f(attr_id, H5T_NATIVE_DOUBLE, attrValue, array_size, hdferr)
   if (hdferr < 0) error stop 'HDF5 error'
@@ -676,13 +641,10 @@ subroutine HDF5_read_real1(dataset,loc_id,datasetName,parallel)
 
   myShape = int(shape(dataset),HSIZE_T)
 
-  if (present(parallel)) then
+  call initialize_read(dset_id,filespace_id,memspace_id,plist_id, aplist_id, &
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
+                       myStart,totalShape,loc_id,myShape,datasetName, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel)
+                       misc_optional(parallel,parallel_default))
-  else
+
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
-  end if
   if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
@@ -713,13 +675,10 @@ subroutine HDF5_read_real2(dataset,loc_id,datasetName,parallel)
 
   myShape = int(shape(dataset),HSIZE_T)
 
-  if (present(parallel)) then
+  call initialize_read(dset_id,filespace_id,memspace_id,plist_id,aplist_id, &
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
+                       myStart,totalShape,loc_id,myShape,datasetName, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel)
+                       misc_optional(parallel,parallel_default))
-  else
+
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
-  end if
   if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
@@ -750,13 +709,10 @@ subroutine HDF5_read_real3(dataset,loc_id,datasetName,parallel)
 
   myShape = int(shape(dataset),HSIZE_T)
 
-  if (present(parallel)) then
+  call initialize_read(dset_id,filespace_id,memspace_id,plist_id,aplist_id, &
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
+                       myStart,totalShape,loc_id,myShape,datasetName, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel)
+                       misc_optional(parallel,parallel_default))
-  else
+
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
-  end if
   if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
@@ -788,13 +744,10 @@ subroutine HDF5_read_real4(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_read(dset_id,filespace_id,memspace_id,plist_id,aplist_id, &
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
+                       myStart,totalShape,loc_id,myShape,datasetName, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel)
+                       misc_optional(parallel,parallel_default))
-  else
+
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
-  end if
   if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
@@ -826,13 +779,10 @@ subroutine HDF5_read_real5(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_read(dset_id,filespace_id,memspace_id,plist_id,aplist_id, &
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
+                       myStart,totalShape,loc_id,myShape,datasetName, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel)
+                       misc_optional(parallel,parallel_default))
-  else
+
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
-  end if
   if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
@@ -864,14 +814,11 @@ subroutine HDF5_read_real6(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_read(dset_id,filespace_id,memspace_id,plist_id,aplist_id, &
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
+                       myStart,totalShape,loc_id,myShape,datasetName, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel)
+                       misc_optional(parallel,parallel_default))
-  else
+
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
+if (any(totalShape == 0)) return
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
-  end if
-  if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
                  file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
@@ -902,13 +849,10 @@ subroutine HDF5_read_real7(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_read(dset_id,filespace_id,memspace_id,plist_id,aplist_id, &
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
+                       myStart,totalShape,loc_id,myShape,datasetName, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel)
+                       misc_optional(parallel,parallel_default))
-  else
+
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
-  end if
   if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
@@ -940,14 +884,11 @@ subroutine HDF5_read_int1(dataset,loc_id,datasetName,parallel)
 
   myShape = int(shape(dataset),HSIZE_T)
 
-  if (present(parallel)) then
+  call initialize_read(dset_id,filespace_id,memspace_id,plist_id,aplist_id, &
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
+                       myStart,totalShape,loc_id,myShape,datasetName, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel)
+                       misc_optional(parallel,parallel_default))
-  else
+
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
+if (any(totalShape == 0)) return
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
-  end if
-  if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
                  file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
@@ -978,13 +919,10 @@ subroutine HDF5_read_int2(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_read(dset_id,filespace_id,memspace_id,plist_id,aplist_id, &
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
+                       myStart,totalShape,loc_id,myShape,datasetName, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel)
+                       misc_optional(parallel,parallel_default))
-  else
+
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
-  end if
   if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
@@ -1015,13 +953,10 @@ subroutine HDF5_read_int3(dataset,loc_id,datasetName,parallel)
 
   myShape = int(shape(dataset),HSIZE_T)
 
-  if (present(parallel)) then
+  call initialize_read(dset_id,filespace_id,memspace_id,plist_id,aplist_id, &
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
+                         myStart,totalShape,loc_id,myShape,datasetName, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel)
+                         misc_optional(parallel,parallel_default))
-  else
+
-   call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
-                        myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
-  end if
   if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
@@ -1052,13 +987,10 @@ subroutine HDF5_read_int4(dataset,loc_id,datasetName,parallel)
 
   myShape = int(shape(dataset),HSIZE_T)
 
-  if (present(parallel)) then
+  call initialize_read(dset_id,filespace_id,memspace_id,plist_id,aplist_id, &
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
+                         myStart,totalShape,loc_id,myShape,datasetName, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel)
+                         misc_optional(parallel,parallel_default))
-  else
+
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
-  end if
   if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
@@ -1089,13 +1021,10 @@ subroutine HDF5_read_int5(dataset,loc_id,datasetName,parallel)
 
   myShape = int(shape(dataset),HSIZE_T)
 
-  if (present(parallel)) then
+  call initialize_read(dset_id,filespace_id,memspace_id,plist_id,aplist_id, &
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
+                         myStart,totalShape,loc_id,myShape,datasetName, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel)
+                         misc_optional(parallel,parallel_default))
-  else
+
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
-  end if
   if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
@@ -1127,13 +1056,10 @@ subroutine HDF5_read_int6(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_read(dset_id,filespace_id,memspace_id,plist_id,aplist_id, &
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
+                       myStart,totalShape,loc_id,myShape,datasetName, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel)
+                       misc_optional(parallel,parallel_default))
-  else
+
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
-  end if
   if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
@@ -1165,13 +1091,10 @@ subroutine HDF5_read_int7(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_read(dset_id,filespace_id,memspace_id,plist_id,aplist_id, &
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
+                       myStart,totalShape,loc_id,myShape,datasetName, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel)
+                       misc_optional(parallel,parallel_default))
-  else
+
-    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
-                         myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
-  end if
   if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
@@ -1207,13 +1130,9 @@ subroutine HDF5_write_real1(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_write(dset_id,filespace_id,memspace_id,plist_id, &
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
+                        myStart,totalShape,loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE, &
-                          myStart, totalShape,loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel)
+                        misc_optional(parallel,parallel_default))
-  else
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
-                          myStart, totalShape,loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel_default)
-  end if
 
   if (product(totalShape) /= 0) then
     call H5Dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
@@ -1248,13 +1167,9 @@ subroutine HDF5_write_real2(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_write(dset_id,filespace_id,memspace_id,plist_id, &
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
+                        myStart,totalShape,loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel)
+                        misc_optional(parallel,parallel_default))
-  else
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel_default)
-  end if
 
   if (product(totalShape) /= 0) then
     call H5Dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
@@ -1289,13 +1204,9 @@ subroutine HDF5_write_real3(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_write(dset_id,filespace_id,memspace_id,plist_id, &
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
+                        myStart,totalShape,loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel)
+                        misc_optional(parallel,parallel_default))
-  else
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel_default)
-  end if
 
   if (product(totalShape) /= 0) then
     call H5Dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
@@ -1330,13 +1241,9 @@ subroutine HDF5_write_real4(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_write(dset_id,filespace_id,memspace_id,plist_id, &
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
+                        myStart,totalShape,loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel)
+                        misc_optional(parallel,parallel_default))
-  else
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel_default)
-  end if
 
   if (product(totalShape) /= 0) then
     call H5Dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
@@ -1372,13 +1279,9 @@ subroutine HDF5_write_real5(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_write(dset_id,filespace_id,memspace_id,plist_id, &
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
+                        myStart,totalShape,loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel)
+                        misc_optional(parallel,parallel_default))
-  else
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel_default)
-  end if
 
   if (product(totalShape) /= 0) then
     call H5Dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
@@ -1413,13 +1316,9 @@ subroutine HDF5_write_real6(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_write(dset_id,filespace_id,memspace_id,plist_id, &
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
+                        myStart,totalShape,loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel)
+                        misc_optional(parallel,parallel_default))
-  else
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel_default)
-  end if
 
   if (product(totalShape) /= 0) then
     call H5Dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
@@ -1454,13 +1353,9 @@ subroutine HDF5_write_real7(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_write(dset_id,filespace_id,memspace_id,plist_id, &
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
+                        myStart,totalShape,loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel)
+                        misc_optional(parallel,parallel_default))
-  else
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel_default)
-  end if
 
   if (product(totalShape) /= 0) then
     call H5Dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
@@ -1497,13 +1392,9 @@ subroutine HDF5_write_real(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_write(dset_id,filespace_id,memspace_id,plist_id, &
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
+                        myStart,totalShape,loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE, &
-                          myStart, totalShape,loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel)
+                        misc_optional(parallel,parallel_default))
-  else
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
-                          myStart, totalShape,loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel_default)
-  end if
 
   if (product(totalShape) /= 0) then
     select rank(dataset)
@@ -1624,13 +1515,9 @@ subroutine HDF5_write_int1(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_write(dset_id,filespace_id,memspace_id,plist_id, &
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
+                        myStart,totalShape,loc_id,myShape,datasetName,H5T_NATIVE_INTEGER, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel)
+                        misc_optional(parallel,parallel_default))
-  else
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel_default)
-  end if
 
   if (product(totalShape) /= 0) then
     call H5Dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
@@ -1665,13 +1552,9 @@ subroutine HDF5_write_int2(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_write(dset_id,filespace_id,memspace_id,plist_id, &
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
+                        myStart,totalShape,loc_id,myShape,datasetName,H5T_NATIVE_INTEGER, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel)
+                        misc_optional(parallel,parallel_default))
-  else
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel_default)
-  end if
 
   if (product(totalShape) /= 0) then
     call H5Dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
@@ -1706,13 +1589,9 @@ subroutine HDF5_write_int3(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_write(dset_id,filespace_id,memspace_id,plist_id, &
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
+                        myStart,totalShape,loc_id,myShape,datasetName,H5T_NATIVE_INTEGER, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel)
+                        misc_optional(parallel,parallel_default))
-  else
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel_default)
-  end if
 
   if (product(totalShape) /= 0) then
     call H5Dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
@@ -1747,13 +1626,9 @@ subroutine HDF5_write_int4(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_write(dset_id,filespace_id,memspace_id,plist_id, &
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
+                        myStart,totalShape,loc_id,myShape,datasetName,H5T_NATIVE_INTEGER, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel)
+                        misc_optional(parallel,parallel_default))
-  else
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel_default)
-  end if
 
   if (product(totalShape) /= 0) then
     call H5Dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
@@ -1788,18 +1663,14 @@ subroutine HDF5_write_int5(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_write(dset_id,filespace_id,memspace_id,plist_id, &
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
+                        myStart,totalShape,loc_id,myShape,datasetName,H5T_NATIVE_INTEGER, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel)
+                        misc_optional(parallel,parallel_default))
-  else
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel_default)
-  end if
 
   if (product(totalShape) /= 0) then
     call H5Dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
                    file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
-   if (hdferr < 0) error stop 'HDF5 error'
+    if (hdferr < 0) error stop 'HDF5 error'
   end if
 
   call finalize_write(plist_id, dset_id, filespace_id, memspace_id)
@@ -1829,13 +1700,9 @@ subroutine HDF5_write_int6(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_write(dset_id,filespace_id,memspace_id,plist_id, &
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
+                        myStart,totalShape,loc_id,myShape,datasetName,H5T_NATIVE_INTEGER, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel)
+                        misc_optional(parallel,parallel_default))
-  else
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel_default)
-  end if
 
   if (product(totalShape) /= 0) then
     call H5Dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
@@ -1870,13 +1737,9 @@ subroutine HDF5_write_int7(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_write(dset_id,filespace_id,memspace_id,plist_id, &
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
+                        myStart,totalShape,loc_id,myShape,datasetName,H5T_NATIVE_INTEGER, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel)
+                        misc_optional(parallel,parallel_default))
-  else
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel_default)
-  end if
 
   if (product(totalShape) /= 0) then
     call H5Dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
@@ -1913,13 +1776,9 @@ subroutine HDF5_write_int(dataset,loc_id,datasetName,parallel)
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-  if (present(parallel)) then
+  call initialize_write(dset_id,filespace_id,memspace_id,plist_id, &
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
+                        myStart,totalShape,loc_id,myShape,datasetName,H5T_NATIVE_INTEGER, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel)
+                        misc_optional(parallel,parallel_default))
-  else
-    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
-                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel_default)
-  end if
 
   if (product(totalShape) /= 0) then
     select rank(dataset)

src/IO.f90

@@ -727,6 +727,7 @@ subroutine selfTest()
   if (CRLF2LF('A'//CR//LF//'B') /= 'A'//LF//'B')     error stop 'CRLF2LF/3'
   if (CRLF2LF('A'//CR//LF//'B'//CR//LF) /= &
               'A'//LF//'B'//LF)                      error stop 'CRLF2LF/4'
+  if (CRLF2LF('A'//LF//CR//'B') /= 'A'//LF//CR//'B') error stop 'CRLF2LF/5'
 
   str='  ';        if (.not. IO_isBlank(str))        error stop 'IO_isBlank/1'
   str='  #isBlank';if (.not. IO_isBlank(str))        error stop 'IO_isBlank/2'

src/Marc/DAMASK_Marc.f90

@@ -142,6 +142,7 @@ end function solverIsSymmetric
 end module DAMASK_interface
 
 #include "../parallelization.f90"
+#include "../misc.f90"
 #include "../constants.f90"
 #include "../IO.f90"
 #include "../YAML_types.f90"

src/YAML_parse.f90

@@ -1,11 +1,12 @@
 !----------------------------------------------------------------------------------------------------
 !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
 !> @author Sharan Roongta, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief Parser for YAML files
+!> @brief Parser for YAML files.
-!> @details module converts a YAML input file to an equivalent YAML flow style which is then parsed.
+!> @details Module converts a YAML input file to an equivalent YAML flow style which is then parsed.
 !----------------------------------------------------------------------------------------------------
 module YAML_parse
   use prec
+  use misc
   use IO
   use YAML_types
 #ifdef FYAML
@@ -54,7 +55,7 @@ end subroutine YAML_parse_init
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Parse a YAML string with list as root into a a structure of nodes.
+!> @brief Parse a YAML string with list at root into a structure of nodes.
 !> @details The string needs to end with a newline (unless using libfyaml).
 !--------------------------------------------------------------------------------------------------
 function YAML_parse_str_asList(str) result(list)
@@ -72,7 +73,7 @@ end function YAML_parse_str_asList
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Parse a YAML string with dict as root into a a structure of nodes.
+!> @brief Parse a YAML string with dict at root into a structure of nodes.
 !> @details The string needs to end with a newline (unless using libfyaml).
 !--------------------------------------------------------------------------------------------------
 function YAML_parse_str_asDict(str) result(dict)
@@ -90,12 +91,12 @@ end function YAML_parse_str_asDict
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Read the flow style string and store it in the form of dictionaries, lists and scalars.
+!> @brief Read a string in flow style and store it in the form of dictionaries, lists, and scalars.
-!> @details A node type pointer can either point to a dictionary, list or scalar type entities.
+!> @details A node-type pointer can either point to a dictionary, list, or scalar type entities.
 !--------------------------------------------------------------------------------------------------
 recursive function parse_flow(YAML_flow) result(node)
 
-  character(len=*), intent(in) :: YAML_flow                                                      !< YAML file in flow style
+  character(len=*), intent(in) :: YAML_flow                                                         !< YAML file in flow style
   class(tNode), pointer        :: node
 
   class(tNode), pointer :: &
@@ -158,7 +159,7 @@ end function parse_flow
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Find location of chunk end: ',' or '}' or  ']'.
+!> @brief Find location of chunk end: ',' '}', or ']'.
 !> @details leaves nested lists ( '[...]' and dicts '{...}') intact
 !--------------------------------------------------------------------------------------------------
 integer function find_end(str,e_char)
@@ -168,7 +169,7 @@ integer function find_end(str,e_char)
 
   integer(pI64) :: N_sq, &                                                                          !< number of open square brackets
                    N_cu, &                                                                          !< number of open curly brackets
-                   i
+                                  i
 
   N_sq = 0_pI64
   N_cu = 0_pI64
@@ -209,7 +210,7 @@ end function quotedString
 
 #ifdef FYAML
 !--------------------------------------------------------------------------------------------------
-! @brief Convert all block style YAML parts to flow style.
+! @brief Convert all block-style YAML parts to flow style.
 !--------------------------------------------------------------------------------------------------
 function to_flow(mixed) result(flow)
 
@@ -237,41 +238,43 @@ end function to_flow
 
 #else
 !--------------------------------------------------------------------------------------------------
-! @brief Determine Indentation.
+! @brief Determine indentation depth.
-! @details It determines the indentation level for a given block/line.
+! @details Indentation level is determined for a given block/line.
-! In cases for nested lists, an offset is added to determine the indent of the item block (skip
+! In case of nested lists, an offset is added to determine the indent of the item block (skip
-! leading dashes)
+! leading dashes).
 !--------------------------------------------------------------------------------------------------
 integer function indentDepth(line,offset)
 
   character(len=*), intent(in) :: line
   integer, optional,intent(in) :: offset
 
-  indentDepth = verify(line,IO_WHITESPACE) -1
+
-  if (present(offset)) indentDepth = indentDepth + offset
+  indentDepth = verify(line,IO_WHITESPACE) - 1 + misc_optional(offset,0)
 
 end function indentDepth
 
 
 !--------------------------------------------------------------------------------------------------
-! @brief Check whether a string is in flow style, i.e. starts with '{' or '['
+! @brief Check whether a string is in flow style, i.e. starts with '{' or '['.
 !--------------------------------------------------------------------------------------------------
 logical function isFlow(line)
 
   character(len=*), intent(in) :: line
 
+
   isFlow = index(adjustl(line),'[') == 1 .or. index(adjustl(line),'{') == 1
 
 end function isFlow
 
 
 !--------------------------------------------------------------------------------------------------
-! @brief check whether a string is a scalar item, i.e. starts without any special symbols
+! @brief Check whether a string is a scalar item, i.e. starts without any special symbols.
 !--------------------------------------------------------------------------------------------------
 logical function isScalar(line)
 
   character(len=*), intent(in) :: line
 
+
   isScalar = (.not. isKeyValue(line) .and. &
               .not. isKey(line) .and. &
               .not. isListItem(line) .and. &
@@ -281,12 +284,13 @@ end function isScalar
 
 
 !--------------------------------------------------------------------------------------------------
-! @brief check whether a string is a list item, i.e. starts with '-'
+! @brief Check whether a string is a list item, i.e. starts with '-'.
 !--------------------------------------------------------------------------------------------------
 logical function isListItem(line)
 
   character(len=*), intent(in) :: line
 
+
   isListItem = .false.
   if (len_trim(adjustl(line))> 2 .and. index(trim(adjustl(line)), '-') == 1) then
     isListItem = scan(trim(adjustl(line)),' ') == 2
@@ -298,13 +302,14 @@ end function isListItem
 
 
 !--------------------------------------------------------------------------------------------------
-! @brief check whether a string contains a key value pair of the for '<key>: <value>'
+! @brief Check whether a string contains a key-value pair of the form '<key>: <value>'.
 !--------------------------------------------------------------------------------------------------
 logical function isKeyValue(line)
 
   character(len=*), intent(in) :: line
   isKeyValue = .false.
 
+
   if ( .not. isKey(line) .and. index(IO_rmComment(line),':') > 0 .and. .not. isFlow(line)) then
     if (index(IO_rmComment(line),': ') > 0) isKeyValue = .true.
   end if
@@ -313,13 +318,14 @@ end function isKeyValue
 
 
 !--------------------------------------------------------------------------------------------------
-! @brief check whether a string contains a key without a value, i.e. it ends with ':'
+! @brief Check whether a string contains a key without a value, i.e. it ends in ':'.
-! ToDo: check whether this is safe for trailing spaces followed by a new line character
+! ToDo: check whether this is safe for trailing spaces followed by a newline character
 !--------------------------------------------------------------------------------------------------
 logical function isKey(line)
 
   character(len=*), intent(in) :: line
 
+
   if (len(IO_rmComment(line)) == 0) then
     isKey = .false.
   else
@@ -332,20 +338,21 @@ end function isKey
 
 
 !--------------------------------------------------------------------------------------------------
-! @brief check whether a string is a list in flow style
+! @brief Check whether a string is a list in flow style.
 !--------------------------------------------------------------------------------------------------
 logical function isFlowList(line)
 
   character(len=*), intent(in) :: line
 
+
   isFlowList = index(adjustl(line),'[') == 1
 
 end function isFlowList
 
 
 !--------------------------------------------------------------------------------------------------
-! @brief skip empty lines
+! @brief Skip empty lines.
-! @details update start position in the block by skipping empty lines if present.
+! @details Update start position in the block by skipping empty lines if present.
 !--------------------------------------------------------------------------------------------------
 subroutine skip_empty_lines(blck,s_blck)
 
@@ -354,6 +361,7 @@ subroutine skip_empty_lines(blck,s_blck)
 
   logical :: empty
 
+
   empty = .true.
   do while(empty .and. len_trim(blck(s_blck:)) /= 0)
     empty = len_trim(IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))) == 0
@@ -364,8 +372,8 @@ end subroutine skip_empty_lines
 
 
 !--------------------------------------------------------------------------------------------------
-! @brief skip file header
+! @brief Skip file header.
-! @details update start position in the block by skipping file header if present.
+! @details Update start position in the block by skipping file header if present.
 !--------------------------------------------------------------------------------------------------
 subroutine skip_file_header(blck,s_blck)
 
@@ -374,6 +382,7 @@ subroutine skip_file_header(blck,s_blck)
 
   character(len=:), allocatable    :: line
 
+
   line = IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))
   if (index(adjustl(line),'%YAML') == 1) then
     s_blck = s_blck + index(blck(s_blck:),IO_EOL)
@@ -389,7 +398,7 @@ end subroutine skip_file_header
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief check if a line in flow YAML starts and ends in the same line
+!> @brief Check whether a line in flow style starts and ends on the same line.
 !--------------------------------------------------------------------------------------------------
 logical function flow_is_closed(str,e_char)
 
@@ -400,6 +409,7 @@ logical function flow_is_closed(str,e_char)
                                   i
   character(len=:), allocatable:: line
 
+
   flow_is_closed = .false.
   N_sq = 0
   N_cu = 0
@@ -418,7 +428,7 @@ end function flow_is_closed
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief return the flow YAML line without line break
+!> @brief Return a flow-style line without line break.
 !--------------------------------------------------------------------------------------------------
 subroutine remove_line_break(blck,s_blck,e_char,flow_line)
 
@@ -428,6 +438,7 @@ subroutine remove_line_break(blck,s_blck,e_char,flow_line)
   character(len=:), allocatable, intent(out) :: flow_line
   logical :: line_end
 
+
   line_end =.false.
   flow_line = ''
 
@@ -441,7 +452,7 @@ end subroutine remove_line_break
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief return the scalar list item without line break
+!> @brief Return a scalar list item without line break.
 !--------------------------------------------------------------------------------------------------
 subroutine list_item_inline(blck,s_blck,inline,offset)
 
@@ -453,6 +464,7 @@ subroutine list_item_inline(blck,s_blck,inline,offset)
   character(len=:), allocatable :: line
   integer :: indent,indent_next
 
+
   indent = indentDepth(blck(s_blck:),offset)
   line   = IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))
   inline = line(indent-offset+3:)
@@ -472,8 +484,8 @@ end subroutine list_item_inline
 
 
 !--------------------------------------------------------------------------------------------------
-! @brief reads a line of YAML block which is already in flow style
+! @brief Read a line of YAML block that is already in flow style.
-! @details A dict should be enclosed within '{}' for it to be consistent with the DAMASK YAML parser
+! @details A dict should be enclosed within '{}' for it to be consistent with the DAMASK YAML parser.
 !--------------------------------------------------------------------------------------------------
 recursive subroutine line_isFlow(flow,s_flow,line)
 
@@ -486,6 +498,7 @@ recursive subroutine line_isFlow(flow,s_flow,line)
     list_chunk, &
     dict_chunk
 
+
   if (index(adjustl(line),'[') == 1) then
     s = index(line,'[')
     flow(s_flow:s_flow) = '['
@@ -536,8 +549,8 @@ end subroutine line_isFlow
 
 
 !-------------------------------------------------------------------------------------------------
-! @brief reads a line of YAML block of type <key>: <value> and places it in the YAML flow style structure
+! @brief Transform a line of YAML of type <key>: <value> to flow style.
-! @details Makes sure that the <value> is consistent with the input required in DAMASK YAML parser
+! @details Ensures that the <value> is consistent with the input required in the DAMASK YAML parser.
 !-------------------------------------------------------------------------------------------------
 recursive subroutine keyValue_toFlow(flow,s_flow,line)
 
@@ -551,6 +564,7 @@ recursive subroutine keyValue_toFlow(flow,s_flow,line)
     col_pos, &
     offset_value
 
+
   col_pos = index(line,':')
   if (line(col_pos+1:col_pos+1) /= ' ') call IO_error(704,ext_msg=line)
   if (isFlow(line(col_pos+1:))) then
@@ -568,7 +582,7 @@ end subroutine keyValue_toFlow
 
 
 !-------------------------------------------------------------------------------------------------
-! @brief reads a line of YAML block and places it in the YAML flow style structure
+! @brief Transform a line of YAML to flow style.
 !-------------------------------------------------------------------------------------------------
 subroutine line_toFlow(flow,s_flow,line)
 
@@ -587,7 +601,7 @@ end subroutine line_toFlow
 
 
 !-------------------------------------------------------------------------------------------------
-! @brief convert a yaml list in block style to a yaml list in flow style
+! @brief Transform a block-style list to flow style.
 ! @details enters the function when encountered with the list indicator '- '
 ! reads each scalar list item and separates each other with a ','
 ! If list item is non scalar, it stores the offset for that list item block
@@ -680,7 +694,7 @@ end subroutine lst
 
 
 !--------------------------------------------------------------------------------------------------
-! @brief convert a yaml dict in block style to a yaml dict in flow style
+! @brief Transform a block-style dict to flow style.
 ! @details enters the function when encountered with the dictionary indicator ':'
 ! parses each line in the block and compares indentation of a line with the preceding line
 ! upon increase in indentation level -> 'decide' function decides if the line is a list or dict
@@ -818,7 +832,7 @@ end subroutine decide
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Convert all block style YAML parts to flow style.
+!> @brief Convert all block-style parts to flow style.
 !> @details The input needs to end with a newline.
 !--------------------------------------------------------------------------------------------------
 function to_flow(blck)

src/lattice.f90

@@ -8,6 +8,7 @@
 !--------------------------------------------------------------------------------------------------
 module lattice
   use prec
+  use misc
   use IO
   use config
   use math
@@ -2159,12 +2160,8 @@ pure function lattice_isotropic_nu(C,assumption,lattice) result(nu)
   real(pReal) :: K, mu
   logical                       :: error
   real(pReal), dimension(6,6)   :: S
-  character(len=:), allocatable :: lattice_
 
 
-  lattice_ = IO_WHITESPACE
-  if (present(lattice)) lattice_ = lattice
-
   if     (IO_lc(assumption) == 'isostrain') then
     K = sum(C(1:3,1:3)) / 9.0_pReal
   elseif (IO_lc(assumption) == 'isostress') then
@@ -2175,7 +2172,7 @@ pure function lattice_isotropic_nu(C,assumption,lattice) result(nu)
     error stop 'invalid assumption'
   end if
 
-  mu = lattice_isotropic_mu(C,assumption,lattice_)
+  mu = lattice_isotropic_mu(C,assumption,misc_optional(lattice,''))
   nu = (1.5_pReal*K-mu)/(3.0_pReal*K+mu)
 
 end function lattice_isotropic_nu
@@ -2195,14 +2192,10 @@ pure function lattice_isotropic_mu(C,assumption,lattice) result(mu)
 
   logical                       :: error
   real(pReal), dimension(6,6)   :: S
-  character(len=:), allocatable :: lattice_
 
 
-  lattice_ = IO_WHITESPACE
-  if (present(lattice)) lattice_ = lattice
-
   if     (IO_lc(assumption) == 'isostrain') then
-      select case(lattice_)
+      select case(misc_optional(lattice,''))
         case('cF','cI')
           mu = ( C(1,1) - C(1,2) + C(4,4)*3.0_pReal) / 5.0_pReal
         case default
@@ -2213,7 +2206,7 @@ pure function lattice_isotropic_mu(C,assumption,lattice) result(mu)
       end select
 
   elseif (IO_lc(assumption) == 'isostress') then
-      select case(lattice_)
+      select case(misc_optional(lattice,''))
         case('cF','cI')
           mu = 5.0_pReal &
                / (4.0_pReal/(C(1,1)-C(1,2)) + 3.0_pReal/C(4,4))

src/math.f90

@@ -7,6 +7,7 @@
 !--------------------------------------------------------------------------------------------------
 module math
   use prec
+  use misc
   use IO
   use config
   use YAML_types
@@ -140,23 +141,9 @@ pure recursive subroutine math_sort(a, istart, iend, sortDim)
   integer :: ipivot,s,e,d
 
 
-  if (present(istart)) then
+  s = misc_optional(istart,lbound(a,2))
-    s = istart
+  e = misc_optional(iend,ubound(a,2))
-  else
+  d = misc_optional(sortDim,1)
-    s = lbound(a,2)
-  end if
-
-  if (present(iend)) then
-    e = iend
-  else
-    e = ubound(a,2)
-  end if
-
-  if (present(sortDim)) then
-    d = sortDim
-  else
-    d = 1
-  end if
 
   if (s < e) then
     call qsort_partition(a,ipivot, s,e, d)
@@ -448,20 +435,14 @@ pure function math_exp33(A,n)
   real(pReal), dimension(3,3) :: B, math_exp33
 
   real(pReal) :: invFac
-  integer     :: n_,i
+  integer     :: i
 
 
-  if (present(n)) then
-    n_ = n
-  else
-    n_ = 5
-  end if
-
   invFac     = 1.0_pReal                                                                            ! 0!
   B          = math_I3
   math_exp33 = math_I3                                                                              ! A^0 = I
 
-  do i = 1, n_
+  do i = 1, misc_optional(n,5)
     invFac = invFac/real(i,pReal)                                                                   ! invfac = 1/(i!)
     B = matmul(B,A)
     math_exp33 = math_exp33 + invFac*B                                                              ! exp = SUM (A^i)/(i!)
@@ -729,12 +710,7 @@ pure function math_sym33to6(m33,weighted)
   real(pReal), dimension(6) :: w
   integer :: i
 
-
+  w = merge(NRMMANDEL,1.0_pReal,misc_optional(weighted,.true.))
-  if (present(weighted)) then
-    w = merge(NRMMANDEL,1.0_pReal,weighted)
-  else
-    w = NRMMANDEL
-  end if
 
   math_sym33to6 = [(w(i)*m33(MAPNYE(1,i),MAPNYE(2,i)),i=1,6)]
 
@@ -757,11 +733,7 @@ pure function math_6toSym33(v6,weighted)
   integer :: i
 
 
-  if (present(weighted)) then
+  w = merge(INVNRMMANDEL,1.0_pReal,misc_optional(weighted,.true.))
-    w = merge(INVNRMMANDEL,1.0_pReal,weighted)
-  else
-    w = INVNRMMANDEL
-  end if
 
   do i=1,6
     math_6toSym33(MAPNYE(1,i),MAPNYE(2,i)) = w(i)*v6(i)
@@ -831,11 +803,7 @@ pure function math_sym3333to66(m3333,weighted)
   integer :: i,j
 
 
-  if (present(weighted)) then
+  w = merge(NRMMANDEL,1.0_pReal,misc_optional(weighted,.true.))
-    w = merge(NRMMANDEL,1.0_pReal,weighted)
-  else
-    w = NRMMANDEL
-  end if
 
 #ifndef __INTEL_COMPILER
   do concurrent(i=1:6, j=1:6)
@@ -864,11 +832,7 @@ pure function math_66toSym3333(m66,weighted)
   integer :: i,j
 
 
-  if (present(weighted)) then
+  w = merge(INVNRMMANDEL,1.0_pReal,misc_optional(weighted,.true.))
-    w = merge(INVNRMMANDEL,1.0_pReal,weighted)
-  else
-    w = INVNRMMANDEL
-  end if
 
   do i=1,6; do j=1,6
     math_66toSym3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j)) = w(i)*w(j)*m66(i,j)
@@ -996,24 +960,12 @@ impure elemental subroutine math_normal(x,mu,sigma)
   real(pReal), intent(out) :: x
   real(pReal), intent(in), optional :: mu, sigma
 
-  real(pReal) :: sigma_, mu_
   real(pReal), dimension(2) :: rnd
 
 
-  if (present(mu)) then
-    mu_ = mu
-  else
-    mu_ = 0.0_pReal
-  end if
-
-  if (present(sigma)) then
-    sigma_ = sigma
-  else
-    sigma_ = 1.0_pReal
-  end if
-
   call random_number(rnd)
-  x = mu_ + sigma_ * sqrt(-2.0_pReal*log(1.0_pReal-rnd(1)))*cos(TAU*(1.0_pReal - rnd(2)))
+  x = misc_optional(mu,0.0_pReal) &
+    + misc_optional(sigma,1.0_pReal) * sqrt(-2.0_pReal*log(1.0_pReal-rnd(1)))*cos(TAU*(1.0_pReal - rnd(2)))
 
 end subroutine math_normal
 

src/misc.f90

@@ -0,0 +1,100 @@
+!--------------------------------------------------------------------------------------------------
+!> @author Martin Diehl, KU Leuven
+!> @author Philip Eisenlohr, Michigan State University
+!> @brief Miscellaneous tools.
+!--------------------------------------------------------------------------------------------------
+module misc
+  use prec
+
+  implicit none(type,external)
+  private
+
+  interface misc_optional
+    module procedure misc_optional_bool
+    module procedure misc_optional_integer
+    module procedure misc_optional_real
+    module procedure misc_optional_string
+  end interface misc_optional
+
+  public :: &
+    misc_optional
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Return bool value if given, otherwise default.
+!--------------------------------------------------------------------------------------------------
+pure function misc_optional_bool(given,default) result(var)
+
+  logical, intent(in), optional :: given
+  logical, intent(in)           :: default
+  logical                       :: var
+
+
+  if (present(given)) then
+    var = given
+  else
+    var = default
+  end if
+
+end function misc_optional_bool
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Return integer value if given, otherwise default.
+!--------------------------------------------------------------------------------------------------
+pure function misc_optional_integer(given,default) result(var)
+
+  integer, intent(in), optional :: given
+  integer, intent(in)           :: default
+  integer                       :: var
+
+
+  if (present(given)) then
+    var = given
+  else
+    var = default
+  end if
+
+end function misc_optional_integer
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Return real value if given, otherwise default.
+!--------------------------------------------------------------------------------------------------
+pure function misc_optional_real(given,default) result(var)
+
+  real(pReal), intent(in), optional :: given
+  real(pReal), intent(in)           :: default
+  real(pReal)                       :: var
+
+
+  if (present(given)) then
+    var = given
+  else
+    var = default
+  end if
+
+end function misc_optional_real
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Return string value if given, otherwise default.
+!--------------------------------------------------------------------------------------------------
+pure function misc_optional_string(given,default) result(var)
+
+  character(len=*), intent(in), optional :: given
+  character(len=*), intent(in)           :: default
+  character(len=:), allocatable          :: var
+
+
+  if (present(given)) then
+    var = given
+  else
+    var = default
+  end if
+
+end function misc_optional_string
+
+end module misc

src/result.f90

@@ -6,6 +6,7 @@
 !--------------------------------------------------------------------------------------------------
 module result
   use prec
+  use misc
   use parallelization
   use IO
   use HDF5_utilities
@@ -224,11 +225,7 @@ subroutine result_addAttribute_str(attrLabel,attrValue,path)
   character(len=*), intent(in), optional :: path
 
 
-  if (present(path)) then
+  call HDF5_addAttribute(resultFile,attrLabel, attrValue, path)
-    call HDF5_addAttribute(resultFile,attrLabel, attrValue, path)
-  else
-    call HDF5_addAttribute(resultFile,attrLabel, attrValue)
-  end if
 
 end subroutine result_addAttribute_str
 
@@ -243,11 +240,7 @@ subroutine result_addAttribute_int(attrLabel,attrValue,path)
   character(len=*), intent(in), optional :: path
 
 
-  if (present(path)) then
+  call HDF5_addAttribute(resultFile,attrLabel, attrValue, path)
-    call HDF5_addAttribute(resultFile,attrLabel, attrValue, path)
-  else
-    call HDF5_addAttribute(resultFile,attrLabel, attrValue)
-  end if
 
 end subroutine result_addAttribute_int
 
@@ -262,11 +255,7 @@ subroutine result_addAttribute_real(attrLabel,attrValue,path)
   character(len=*), intent(in), optional :: path
 
 
-  if (present(path)) then
+  call HDF5_addAttribute(resultFile,attrLabel, attrValue, path)
-    call HDF5_addAttribute(resultFile,attrLabel, attrValue, path)
-  else
-    call HDF5_addAttribute(resultFile,attrLabel, attrValue)
-  end if
 
 end subroutine result_addAttribute_real
 
@@ -281,11 +270,7 @@ subroutine result_addAttribute_str_array(attrLabel,attrValue,path)
   character(len=*), intent(in), optional     :: path
 
 
-  if (present(path)) then
+  call HDF5_addAttribute(resultFile,attrLabel, attrValue, path)
-    call HDF5_addAttribute(resultFile,attrLabel, attrValue, path)
-  else
-    call HDF5_addAttribute(resultFile,attrLabel, attrValue)
-  end if
 
 end subroutine result_addAttribute_str_array
 
@@ -300,11 +285,7 @@ subroutine result_addAttribute_int_array(attrLabel,attrValue,path)
   character(len=*), intent(in), optional     :: path
 
 
-  if (present(path)) then
+  call HDF5_addAttribute(resultFile,attrLabel, attrValue, path)
-    call HDF5_addAttribute(resultFile,attrLabel, attrValue, path)
-  else
-    call HDF5_addAttribute(resultFile,attrLabel, attrValue)
-  end if
 
 end subroutine result_addAttribute_int_array
 
@@ -319,11 +300,7 @@ subroutine result_addAttribute_real_array(attrLabel,attrValue,path)
   character(len=*), intent(in), optional     :: path
 
 
-  if (present(path)) then
+  call HDF5_addAttribute(resultFile,attrLabel, attrValue, path)
-    call HDF5_addAttribute(resultFile,attrLabel, attrValue, path)
-  else
-    call HDF5_addAttribute(resultFile,attrLabel, attrValue)
-  end if
 
 end subroutine result_addAttribute_real_array
 
@@ -416,19 +393,12 @@ subroutine result_writeTensorDataset_real(dataset,group,label,description,SIunit
   real(pReal),      intent(in), dimension(:,:,:) :: dataset
 
   integer :: i
-  logical :: transposed_
   integer(HID_T) :: groupHandle
   real(pReal), dimension(:,:,:), allocatable :: dataset_transposed
 
 
-  if (present(transposed)) then
-    transposed_ = transposed
-  else
-    transposed_ = .true.
-  end if
-
   groupHandle = result_openGroup(group)
-  if (transposed_) then
+  if (misc_optional(transposed,.true.)) then
     if (size(dataset,1) /= size(dataset,2)) error stop 'transpose non-symmetric tensor'
     allocate(dataset_transposed,mold=dataset)
     do i=1,size(dataset_transposed,3)

src/rotations.f90

@@ -178,11 +178,7 @@ subroutine fromEulers(self,eu,degrees)
   real(pReal), dimension(3)             :: Eulers
 
 
-  if (.not. present(degrees)) then
+  Eulers = merge(eu*INRAD,eu,misc_optional(degrees,.false.))
-    Eulers = eu
-  else
-    Eulers = merge(eu*INRAD,eu,degrees)
-  end if
 
   if (any(Eulers<0.0_pReal) .or. any(Eulers>TAU) .or. Eulers(2) > PI) &
     call IO_error(402,ext_msg='fromEulers')
@@ -202,18 +198,10 @@ subroutine fromAxisAngle(self,ax,degrees,P)
   real(pReal),dimension(3)              :: axis
 
 
-  if (.not. present(degrees)) then
+  angle = merge(ax(4)*INRAD,ax(4),misc_optional(degrees,.false.))
-    angle = ax(4)
-  else
-    angle = merge(ax(4)*INRAD,ax(4),degrees)
-  end if
 
-  if (.not. present(P)) then
+  axis = ax(1:3) * merge(-1.0_pReal,1.0_pReal,misc_optional(P,-1) == 1)
-    axis = ax(1:3)
+  if (abs(misc_optional(P,-1)) /= 1) call IO_error(402,ext_msg='fromAxisAngle (P)')
-  else
-    axis = ax(1:3) * merge(-1.0_pReal,1.0_pReal,P == 1)
-    if (abs(P) /= 1) call IO_error(402,ext_msg='fromAxisAngle (P)')
-  end if
 
   if (dNeq(norm2(axis),1.0_pReal) .or. angle < 0.0_pReal .or. angle > PI) &
     call IO_error(402,ext_msg='fromAxisAngle')
@@ -277,22 +265,15 @@ pure function rotVector(self,v,active) result(vRot)
   logical,         intent(in), optional     :: active
 
   real(pReal), dimension(4) :: v_normed, q
-  logical                   :: passive
 
 
-  if (present(active)) then
-    passive = .not. active
-  else
-    passive = .true.
-  end if
-
   if (dEq0(norm2(v))) then
     vRot = v
   else
     v_normed = [0.0_pReal,v]/norm2(v)
-    q = merge(multiplyQuaternion(self%q, multiplyQuaternion(v_normed, conjugateQuaternion(self%q))), &
+    q = merge(multiplyQuaternion(conjugateQuaternion(self%q), multiplyQuaternion(v_normed, self%q)), &
-              multiplyQuaternion(conjugateQuaternion(self%q), multiplyQuaternion(v_normed, self%q)), &
+              multiplyQuaternion(self%q, multiplyQuaternion(v_normed, conjugateQuaternion(self%q))), &
-              passive)
+              misc_optional(active,.false.))
     vRot = q(2:4)*norm2(v)
   end if
 
@@ -311,18 +292,10 @@ pure function rotTensor2(self,T,active) result(tRot)
   real(pReal),     intent(in), dimension(3,3) :: T
   logical,         intent(in), optional       :: active
 
-  logical           :: passive
 
-
+  tRot = merge(matmul(matmul(transpose(self%asMatrix()),T),self%asMatrix()), &
-  if (present(active)) then
+               matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix())), &
-    passive = .not. active
+               misc_optional(active,.false.))
-  else
-    passive = .true.
-  end if
-
-  tRot = merge(matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix())), &
-               matmul(matmul(transpose(self%asMatrix()),T),self%asMatrix()), &
-               passive)
 
 end function rotTensor2
 
@@ -342,12 +315,7 @@ pure function rotTensor4(self,T,active) result(tRot)
   real(pReal), dimension(3,3) :: R
   integer :: i,j,k,l,m,n,o,p
 
-
+  R = merge(transpose(self%asMatrix()),self%asMatrix(),misc_optional(active,.false.))
-  if (present(active)) then
-    R = merge(transpose(self%asMatrix()),self%asMatrix(),active)
-  else
-    R = self%asMatrix()
-  end if
 
   tRot = 0.0_pReal
   do i = 1,3;do j = 1,3;do k = 1,3;do l = 1,3
@@ -375,11 +343,7 @@ pure function rotStiffness(self,C,active) result(cRot)
   real(pReal), dimension(6,6) :: M
 
 
-  if (present(active)) then
+  R = merge(transpose(self%asMatrix()),self%asMatrix(),misc_optional(active,.false.))
-    R = merge(transpose(self%asMatrix()),self%asMatrix(),active)
-  else
-    R = self%asMatrix()
-  end if
 
   M = reshape([R(1,1)**2,                   R(2,1)**2,                   R(3,1)**2, &
                R(2,1)*R(3,1),               R(1,1)*R(3,1),               R(1,1)*R(2,1), &