fix number of bugs regarding setting c_11

src/Physics/PDE_Coefficients.F90

@@ -1063,6 +1063,9 @@ End Subroutine Initialize_Reference_Heating
     Subroutine Augment_Reference()
         Implicit None
         Real*8, Allocatable :: temp_functions(:,:), temp_constants(:)
+        Character(len=2) :: intstring
+        Character*12 :: dstring
+        Character*8 :: dofmt = '(ES12.5)'
 
         If (my_rank .eq. 0) Then
             Call stdout%print('Reference state will be augmented.')
@@ -1111,10 +1114,29 @@ Subroutine Augment_Reference()
             temp_constants(2) = ra_constants(2)
         Endif
 
-        If (use_custom_function(14)) Then
+        If (with_custom_reference .and. use_custom_function(14)) Then
+            ! Set c_11 = 1 by default
+            If (.not. use_custom_constant(11)) Then
+                If (my_rank .eq. 0) Then
+                    Call stdout%print('User didn''t set c_11; now setting c_11 to 1.')
+                Endif
+                ra_constants(11) = 1.0d0
+            Else
+                If (my_rank .eq. 0) Then
+                    Write(dstring,dofmt) ra_constants(11)
+                    Call stdout%print('User set c_11 to '//Trim(dstring))
+                Endif
+            Endif
             If (my_rank .eq. 0) Then
                 Call stdout%print('Background thermal gradient is set to:')
                 Call stdout%print('c_11*f_14')
+                ! Make a warning if the user is trying to change a dsdr that physically was assumed = 0
+                If (any(reference_type .eq. (/1, 2/))) Then 
+                    Write(intstring, '(I2)') reference_type
+                    Call stdout%print('WARNING: reference_type = '//Adjustl(intstring)//' assumes dsdr = 0 by default.')
+                    Call stdout%print('make sure the combination of') 
+                    Call stdout%print('reference_type = '//Adjustl(intstring)//' and nonzero dsdr was intended.')
+                Endif
                 Call stdout%print(' ')
             Endif        
             ref%dsdr(:) = ra_constants(11)*ra_functions(:,14)
@@ -2017,19 +2039,21 @@ Subroutine Set_Reference_Equation_Coefficients
         ! the combination c_11 and f_14 cannot be backed out from the "ref" structure
         ! (which only has ref%dsdr = c_11*f_14)
         ! and will depend on the convention of the different reference_types
-        Select Case(reference_type)
-            Case(1,2)
-                ra_constants(11) = 0.0d0
-                ra_functions(:,14) = 0.0d0
-            Case(3)
-                ra_constants(11) = 1.0d0
-                ra_functions(:,14) = ref%dsdr
-            ! For reference_type = 4, c_11 and f_14 have been set already
-            Case(5)
-                ra_constants(11) = Prandtl_Number*Buoyancy_Number_Visc/Rayleigh_Number
-                ra_functions(:,14) = ref%dsdr/ra_constants(11)
-        End Select
-
+        If (.not. use_custom_function(14)) Then
+            ! If user set dsdr via "with custom", then c_11 and f_14 have already been set
+            Select Case(reference_type)
+                Case(1,2)
+                    ra_constants(11) = 0.0d0
+                    ra_functions(:,14) = 0.0d0
+                Case(3)
+                    ra_constants(11) = 1.0d0
+                    ra_functions(:,14) = ref%dsdr
+                ! For reference_type = 4, c_11 and f_14 have been set already
+                Case(5)
+                    ra_constants(11) = Prandtl_Number*Buoyancy_Number_Visc/Rayleigh_Number
+                    ra_functions(:,14) = ref%dsdr/ra_constants(11)
+            End Select
+        Endif
 
         ra_functions(:,1) = ref%density
         ra_functions(:,4) = ref%temperature