added sanity check for positive rho,C_p,K in thermal

src/phase.f90

@@ -427,12 +427,12 @@ subroutine phase_init
     phase => phases%get_dict(ph)
     refs = config_listReferences(phase,indent=3)
     if (len(refs) > 0) print'(/,1x,a)', refs
+    phase_rho(ph) = phase%get_asReal('rho',defaultVal=0.0_pREAL)
     phase_lattice(ph) = phase%get_asStr('lattice')
     if (all(phase_lattice(ph) /= ['cF','cI','hP','tI'])) &
       call IO_error(130,ext_msg='phase_init: '//phase%get_asStr('lattice'))
     if (any(phase_lattice(ph) == ['hP','tI'])) &
       phase_cOverA(ph) = phase%get_asReal('c/a')
-    phase_rho(ph) = phase%get_asReal('rho',defaultVal=0.0_pREAL)
     allocate(phase_O_0(ph)%data(count(material_ID_phase==ph)))
   end do
 

src/phase_thermal.f90

@@ -76,7 +76,9 @@ module subroutine thermal_init(phases)
     thermal
   type(tList), pointer :: &
     sources
-  character(len=:), allocatable :: refs
+  character(len=:), allocatable :: &
+    refs, &
+    extmsg
   integer :: &
     ph, so, &
     Nmembers
@@ -88,6 +90,7 @@ module subroutine thermal_init(phases)
   allocate(thermalState(phases%length))
   allocate(thermal_Nsources(phases%length),source = 0)
   allocate(param(phases%length))
+  extmsg = ''
 
   do ph = 1, phases%length
     Nmembers = count(material_ID_phase == ph)
@@ -106,6 +109,15 @@ module subroutine thermal_init(phases)
       if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%K(3,3) = thermal%get_asReal('K_33')
       param(ph)%K = crystal_symmetrize_33(param(ph)%K,phase_lattice(ph))
 
+      ! sanity checks
+      if (    param(ph)%C_p <= 0.0_pREAL )  extmsg = trim(extmsg)//' C_p'
+      if (any(param(ph)%K   <  0.0_pREAL))  extmsg = trim(extmsg)//' K'
+      if (    phase_rho(ph) <= 0.0_pREAL )  extmsg = trim(extmsg)//' rho'
+
+!--------------------------------------------------------------------------------------------------
+!  exit if any parameter is out of range
+      if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg))
+
 #if defined(__GFORTRAN__)
       param(ph)%output = output_as1dStr(thermal)
 #else