(dust) apply save limits and log to all dust species in interactive mode

src/calc_quantities.f90

@@ -554,6 +554,37 @@ subroutine append_grain_size_label(string,idust,tags,vals,ierr)
 
 end subroutine append_grain_size_label
 
+!-----------------------------------------------------------------
+!
+!  utility (private) to get mass of a particular grain species
+!  from the header tags
+!
+!-----------------------------------------------------------------
+real function get_mass_of_species(string,idust,tags,vals,ierr)
+ use labels, only:count_non_blank
+ character(len=*), intent(inout) :: string
+ integer, intent(in)  :: idust
+ character(len=*), intent(in) :: tags(:)
+ real,    intent(in) :: vals(:)
+ integer, intent(out) :: ierr
+ integer :: ntags,nd,i
+
+ get_mass_of_species = 0.
+ ntags = count_non_blank(tags)
+ nd = 0
+ ierr = 1
+ do i=1,ntags
+    if (index(tags(i),'mdust_in') > 0) then
+       nd = nd + 1
+       if (nd==idust) then
+          ierr = 0
+          get_mass_of_species = vals(i)
+       endif
+    endif
+ enddo
+
+end function get_mass_of_species
+
 !-----------------------------------------------------------------
 !
 !  utility (private) to either print the example quantity or
@@ -599,7 +630,7 @@ end subroutine print_or_prefill
 subroutine check_calculated_quantities(ncalcok,ncalctot,incolumn,verbose)
  use settings_data,  only:ncolumns,iRescale
  use fparser,        only:checkf
- use labels,         only:label,unitslabel,shortstring
+ use labels,         only:label,unitslabel,shortstring,irhodust_start,irhodust_end
  integer, intent(out) :: ncalcok,ncalctot
  integer, dimension(maxcalc), intent(out), optional :: incolumn
  logical, intent(in), optional :: verbose
@@ -649,6 +680,11 @@ subroutine check_calculated_quantities(ncalcok,ncalctot,incolumn,verbose)
              print "(1x,i2,') ',a50,' [OK]')",ncolumns+ncalcok,trim(calclabel(i))//' = '//calcstring(i)
           endif
        endif
+       !
+       !--recognise the dust density in the list of calculated quantities
+       !
+       if (trim(calcstring(i))=='density*dustfrac1') irhodust_start = ncolumns+ncalcok
+       if (calcstring(i)(1:16)=='density*dustfrac')    irhodust_end = ncolumns+ncalcok ! overwrite until last density*dustfrac
     else
        indexinactive = indexinactive - 1
        if (isverbose) then
@@ -940,29 +976,35 @@ subroutine identify_calculated_quantity(labelcol,ncolumns,icolumn)
  !
  select case(label_synonym(labelcol))
  case('r','radius','rad')
-    if (irad <= 0 .or. irad > ncolumns) then
-       irad = icolumn
-       if (debugmode) print "(1x,a,i2,a)",'identifying column ',icolumn,' as the radius'
-    endif
+    call assign_column(irad,icolumn,ncolumns,debugmode,'radius')
  case('kinetic energy','ke','1/2 v^2','v^2/2')
-    if (ike <= 0 .or. irad > ncolumns) then
-       ike = icolumn
-       if (debugmode) print "(1x,a,i2,a)",'identifying column ',icolumn,' as the kinetic energy'
-    endif
+    call assign_column(ike,icolumn,ncolumns,debugmode,'kinetic energy')
  case('pressure','pr','p')
-    if (ipr <= 0 .or. ipr > ncolumns) then
-       ipr = icolumn
-       if (debugmode) print "(1x,a,i2,a)",'identifying column ',icolumn,' as the pressure'
-    endif
+    call assign_column(ipr,icolumn,ncolumns,debugmode,'pressure')
  case('kappa','opacity')
-    if (ikappa <= 0 .or. ikappa > ncolumns) then
-       ikappa = icolumn
-       if (debugmode) print "(1x,a,i2,a)",'identifying column ',icolumn,' as the opacity'
-    endif
+    call assign_column(ikappa,icolumn,ncolumns,debugmode,'opacity')
  end select
 
 end subroutine identify_calculated_quantity
 
+!-----------------------------------------------------------------
+!
+!  helper routine for above
+!
+!-----------------------------------------------------------------
+subroutine assign_column(i,icolumn,ncolumns,debugmode,string)
+ integer, intent(inout) :: i
+ integer, intent(in)    :: icolumn,ncolumns
+ logical, intent(in)    :: debugmode
+ character(len=*), intent(in) :: string
+
+ if (i <= 0 .or. i > ncolumns) then
+    i = icolumn
+    if (debugmode) print "(1x,a,i2,a)",'identifying column ',icolumn,' as the '//trim(string)
+ endif
+
+end subroutine assign_column
+
 !-----------------------------------------------------------------
 !
 !  utility (private) to fill the array of variable names

src/interactive.f90

@@ -2821,6 +2821,44 @@ end subroutine mvtitle
 !--saves current plot limits
 !
 subroutine save_limits(iplot,xmin,xmax,setlim2)
+ use labels, only:irhodust_start,irhodust_end
+ integer, intent(in) :: iplot
+ real,    intent(in) :: xmin,xmax
+ logical, intent(in), optional :: setlim2
+ real :: xmintemp,xmaxtemp
+ integer :: i
+
+ if (iplot > 0 .and. iplot >= irhodust_start .and. iplot <= irhodust_end) then
+    !
+    !--if we save the limits for one dust density, apply to whole grid
+    ! 
+    do i=irhodust_start,irhodust_end
+       xmintemp = xmin
+       xmaxtemp = xmax
+       if (present(setlim2)) then
+          call save_limits_i(i,xmintemp,xmaxtemp,setlim2)
+       else
+          call save_limits_i(i,xmintemp,xmaxtemp)
+       endif
+    enddo
+    print*,'> applying saved limits to all dust species <'
+ else
+    !
+    !--otherwise just pass options through to save_limits_i
+    !
+    if (present(setlim2)) then
+       call save_limits_i(iplot,xmin,xmax,setlim2)
+    else
+       call save_limits_i(iplot,xmin,xmax)
+    endif
+ endif
+
+end subroutine save_limits
+
+!
+!--save plot limits for one column
+!
+subroutine save_limits_i(iplot,xmin,xmax,setlim2)
  use limits,          only:lim,lim2
  use labels,          only:is_coord
  use multiplot,       only:itrans
@@ -2866,7 +2904,7 @@ subroutine save_limits(iplot,xmin,xmax,setlim2)
  endif
 
  return
-end subroutine save_limits
+end subroutine save_limits_i
 
 !
 !--implements parameter range restriction
@@ -2975,9 +3013,10 @@ end subroutine save_itrackpart_recalcradius
 !  note this only changes a pure log transform: will not change combinations
 !
 subroutine change_itrans(iplot,xmin,xmax)
- use multiplot, only:itrans
+ use multiplot,     only:itrans
  use settings_data, only:numplot
- use transforms, only:transform_limits,transform_limits_inverse
+ use transforms,    only:transform_limits,transform_limits_inverse
+ use labels,        only:irhodust_start,irhodust_end
  integer, intent(in) :: iplot
  real, intent(inout) :: xmin, xmax
 
@@ -2993,6 +3032,10 @@ subroutine change_itrans(iplot,xmin,xmax)
        xmin = max(xmax-4.,xmin) ! no more than 4 dex by default
     endif
  endif
+ if (iplot > 0 .and. iplot >= irhodust_start .and. iplot <= irhodust_end) then
+    print*,'>> applying transform to all dust densities <<'
+    itrans(irhodust_start:irhodust_end) = itrans(iplot)
+ endif
 
 end subroutine change_itrans
 

src/labels.f90

@@ -36,7 +36,7 @@ module labels
  character(len=20), dimension(maxparttypes) :: labeltype
  character(len=6), parameter :: labeldefault = 'column'
  character(len=lenunitslabel), dimension(0:maxplot), public :: unitslabel,unitslabel_default
- character(len=lenunitslabel), public :: labelzintegration,labelzintegration_default
+ character(len=lenunitslabel) :: labelzintegration,labelzintegration_default
  integer, dimension(3)       :: ix
  integer, dimension(maxplot) :: iamvec
  integer :: ivx,irho,iutherm,ipr,ih,irad,iBfirst,iBpol,iBtor,iax
@@ -49,6 +49,7 @@ module labels
  integer :: irhorestframe,idustfrac,ideltav
  integer :: idustfracsum,ideltavsum
  integer :: igrainsize,igraindens,ivrel
+ integer :: irhodust_start,irhodust_end
  integer :: nreq
 
  public
@@ -99,6 +100,8 @@ subroutine reset_columnids
  ideltav = 0
  ideltavsum = 0
  ipmomx = 0
+ irhodust_start = 0
+ irhodust_end = 0
  headertags = ''
 
 end subroutine reset_columnids