new vgrad terms

Author: Loren

src/Diagnostics/Diagnostics_Base.F90

@@ -94,6 +94,7 @@ Module Diagnostics_Base
     ! actual induction (not multiplied by B)
     Integer, Parameter :: ialt_offset3 = 3000
     Integer, Parameter :: ialt_offset4 = ialt_offset3+100 ! 3100
+    Integer, Parameter :: ialt_offset5 = ialt_offset4+200 ! 3300
 
     ! total
     Integer, Parameter :: induct_work_r     = custom_offset + 1
@@ -300,6 +301,24 @@ Module Diagnostics_Base
     Integer, Parameter :: ialtcomp_t2     = ialt_offset4 + 5 + comp_offset
     Integer, Parameter :: ialtcomp_p2     = ialt_offset4 + 6 + comp_offset
 
+    ! bit more breakup, start at 3300
+    ! comp
+    ! vgrad
+    Integer, Parameter :: ialtcomp_t3     = ialt_offset5 + 1
+    Integer, Parameter :: ialtcomp_t4     = ialt_offset5 + 2
+    ! induct
+    Integer, Parameter :: ialtcomp_t5     = ialt_offset5 + 3
+    Integer, Parameter :: ialtcomp_t6     = ialt_offset5 + 4
+    ! shear
+    ! vgrad
+    Integer, Parameter :: ialtshear_t3   = ialt_offset5 + 5
+
+    ! mag torque things
+    Integer, Parameter :: alttorque1   = ialt_offset5 + 6
+    Integer, Parameter :: alttorque2   = ialt_offset5 + 7
+    Integer, Parameter :: alttorque3   = ialt_offset5 + 8
+    Integer, Parameter :: alttorque4   = ialt_offset5 + 9
+
 
     include "turbKE_codes.F"
     include "axial_field_codes.F"

src/Diagnostics/Diagnostics_Custom.F90

@@ -1334,6 +1334,73 @@ Subroutine Custom_MHD_Diagnostics(buffer)
             Endif
         Enddo
 
+        ! terms bit more broken up
+        ! comp, vgrad
+        If (compute_quantity(ialtcomp_t3)) Then
+            DO_PSI
+                qty(PSI) = -buff1(PSI, dvpdp)/radius(r)/sintheta(t)
+            END_DO
+            Call Add_Quantity(qty)
+        Endif
+        If (compute_quantity(ialtcomp_t4)) Then
+            DO_PSI
+                qty(PSI) = -buff1(PSI, vtheta)*cottheta(t)/radius(r)
+            END_DO
+            Call Add_Quantity(qty)
+        Endif
+        ! comp, induct
+        If (compute_quantity(ialtcomp_t5)) Then
+            DO_PSI
+                qty(PSI) = -buff2(PSI,btheta)*buff1(PSI, dvpdp)/radius(r)/sintheta(t)
+            END_DO
+            Call Add_Quantity(qty)
+        Endif
+        If (compute_quantity(ialtcomp_t6)) Then
+            DO_PSI
+                qty(PSI) = -buff2(PSI,btheta)*buff1(PSI, vtheta)*cottheta(t)/radius(r)
+            END_DO
+            Call Add_Quantity(qty)
+        Endif
+        ! shear, vgrad
+        If (compute_quantity(ialtshear_t3)) Then
+            DO_PSI
+                qty(PSI) = buff1(PSI, dvtdr) - buff1(PSI, vtheta)/radius(r)
+            END_DO
+            Call Add_Quantity(qty)
+        Endif
+
+
+        ! more TORQUE stuff
+        ! torque, vgrad
+        If (compute_quantity(alttorque1)) Then
+            DO_PSI
+                qty(PSI) = radius(r)*sintheta(t)/four_pi*(buffer(PSI, dbpdr) +&
+                    &buffer(PSI, bphi)/radius(r))
+            END_DO
+            Call Add_Quantity(qty)
+        Endif
+        If (compute_quantity(alttorque2)) Then
+            DO_PSI
+                qty(PSI) = buffer(PSI,br)*radius(r)*sintheta(t)/four_pi*(buffer(PSI, dbpdr) +&
+                    &buffer(PSI, bphi)/radius(r))
+            END_DO
+            Call Add_Quantity(qty)
+        Endif
+        If (compute_quantity(alttorque3)) Then
+            DO_PSI
+                qty(PSI) = (sintheta(t)*buffer(PSI, dbpdt) +&
+                    &costheta(t)*buffer(PSI, bphi))/four_pi
+            END_DO
+            Call Add_Quantity(qty)
+        Endif
+        If (compute_quantity(alttorque4)) Then
+            DO_PSI
+                qty(PSI) = buffer(PSI,btheta)/four_pi*(sintheta(t)*buffer(PSI, dbpdt) +&
+                    &costheta(t)*buffer(PSI, bphi))
+            END_DO
+            Call Add_Quantity(qty)
+        Endif
+
         DeAllocate(ind_work_r)
         DeAllocate(ind_work_t)
         DeAllocate(ind_work_p)