Adding some extra scalar field diagnostics.

These are offset to allow for the future addition of other codes in approximately the order in which they appear for the main
fields.  The codes for these anticipated outputs are added but commented out and not implemented yet.

The docs have also been updated.

Author: cianwilson

doc/source/User_Guide/arbitrary_scalar_fields.txt

@@ -255,6 +255,14 @@ Active scalar fields
  :math:`\frac{\partial\chi_{a_i}}{\partial \phi}`               10010+200*(i-1)     longitudinal derivative of active scalar field i
  :math:`\frac{\partial\chi_{a_i}^\prime}{\partial \phi}`        10011+200*(i-1)     longitudinal derivative of the active scalar field i perturbation
  :math:`\frac{\partial\overline{\chi_{a_i}}}{\partial \phi}`    10012+200*(i-1)     longitudinal derivative of the active scalar field i mean
+ :math:`d_{7_i}\mathrm{g}_{7_i}\chi`                            10071+200*(i-1)     chi\_buoyancy\_force for active scalar field i
+ :math:`d_{7_i}\mathrm{g}_{7_i}\chi'`                           10072+200*(i-1)     chi\_buoyancy\_pforce for active scalar field i perturbation
+ :math:`d_{7_i}\mathrm{g}_{7_i}\overline{\chi}`                 10073+200*(i-1)     chi\_buoyancy\_mforce for active scalar field i mean
+ :math:`d_{7_i}\mathrm{g}_{7_i}\chi_{00}`                       10074+200*(i-1)     chi\_buoyancy\_force\_ell0 for active scalar field i l0
+ :math:`d_{7_i}v_r\mathrm{g}_{7_i}\chi`                         10171+200*(i-1)     buoy\_work for active scalar field i
+ :math:`d_{7_i}v'_r\mathrm{g}_{7_i}\chi'`                       10072+200*(i-1)     buoy\_work\_pp  for active scalar field i perturbation
+ :math:`d_{7_i}\overline{v_r}\mathrm{g}_{7_i}\overline{\chi}`   10073+200*(i-1)     buoy\_work\_mm for active scalar field i mean
+ :math:`d_{7_i}\mathrm{g}_{7_i}{\chi^\prime}{v^\prime}_r`       10181+200*(i-1)     chi\_production\_buoyant\_pKE for active scalar field i
 =============================================================== =================== ================================================================
 
 Passive scalar fields

src/Diagnostics/Diagnostics_Scalars.F90

@@ -139,6 +139,75 @@ Subroutine Compute_Scalars(buffer)
                 END_DO
                 Call Add_Quantity(qty)
             Endif
+
+            ! only valid for active scalars
+            if (ii .le. n_active_scalars) then
+                ! NOTE:  The ell=0 component of the r-momentum equation is entirely described by
+                !        hydrostatic balance between the pressure and entropy and composition perturbations
+                !        (the reference state is assumed to also be in hydrostatic balance).
+                !        As such, the ell=0 buoancy force is uninteresting from the point of
+                !        view of the flow.  We explicitly separate the ell=0 component for this
+                !        term (as with the entropy and pressure terms).
+                ! -- full buoyancy
+                If (compute_quantity(chi_buoyancy_force+scoff) .or. compute_quantity(chi_buoy_work+scoff)) Then
+                    DO_PSI
+                        qty(PSI) = ref%chi_buoyancy_coeff(ii,r)*(buffer(PSI,chivar)-&
+                                & ell0_values(r,chivar))
+                    END_DO
+                    If (compute_quantity(chi_buoyancy_force+scoff))Call Add_Quantity(qty)
+                    If (compute_quantity(chi_buoy_work+scoff)) Then
+                        DO_PSI
+                            qty(PSI)=buffer(PSI,vr)*qty(PSI)
+                        END_DO
+                        Call Add_Quantity(qty)
+                    Endif
+                Endif
+
+                ! -- fluctuating buoyancy (ell = 0, m =0 already subtracted)
+                If (compute_quantity(chi_buoyancy_pforce+scoff) .or. compute_quantity(chi_buoy_work_pp+scoff)) Then
+                    DO_PSI
+                        qty(PSI) = ref%chi_buoyancy_coeff(ii,r)*fbuffer(PSI,chivar)
+                    END_DO
+                    If (compute_quantity(chi_buoyancy_pforce+scoff)) Call Add_Quantity(qty)
+                    If (compute_quantity(chi_buoy_work_pp+scoff)) Then
+                        DO_PSI
+                            qty(PSI)=buffer(PSI,vr)*qty(PSI)
+                        END_DO
+                        Call Add_Quantity(qty)
+                    Endif
+                Endif
+
+                ! -- mean buoyancy
+                If (compute_quantity(chi_buoyancy_mforce+scoff) .or. compute_quantity(chi_buoy_work_mm+scoff)) Then
+                    DO_PSI
+                        qty(PSI) = ref%chi_buoyancy_coeff(ii,r)*(m0_values(PSI2,chivar)-&
+                                & ell0_values(r,chivar))
+                    END_DO
+                    If (compute_quantity(chi_buoyancy_mforce+scoff)) Call Add_Quantity(qty)
+                    If (compute_quantity(chi_buoy_work_mm+scoff)) Then
+                        DO_PSI
+                            qty(PSI)=buffer(PSI,vr)*qty(PSI)
+                        END_DO
+                        Call Add_Quantity(qty)
+                    Endif
+                Endif
+
+                If (compute_quantity(chi_buoyancy_force_ell0+scoff)) Then
+                    DO_PSI
+                        qty(PSI) = ref%chi_buoyancy_coeff(ii,r)*ell0_values(r,chivar)
+                    END_DO
+                    Call Add_Quantity(qty)
+                Endif
+
+                ! Chi Buoyant Production of turbulent kinetic energy.
+                If (compute_quantity(chi_production_buoyant_pKE+scoff)) Then
+                    DO_PSI
+                        qty(PSI) = ref%chi_buoyancy_coeff(ii,r)*fbuffer(PSI,chivar)*fbuffer(PSI,vr)
+                    END_DO
+                    Call Add_Quantity(qty)
+                Endif
+            endif
+            
         Enddo
     End Subroutine Compute_Scalars