anchor.diff:         alpha = unit_field(nu_dpdt,nunit) !dPaw/dt, updated each iter
          Qinit = elem_field(ne_Vdot0,ne) !terminal element flow, updated each dt
          ! beta is rate of change of 'external' pressure, incl muscle and entrance
  -       beta = dp_external/dt ! == dPmus/dt (-ve for insp), updated each dt
  +       beta = dp_external/dt ! == dP_muscle/dt (-ve for insp), updated each dt
   
   !!!    Q = C*(alpha-beta)+(Qinit-C*(alpha-beta))*exp(-dt/(C*R))
          Q = unit_field(nu_comp,nunit)*(alpha-beta)+ &
  
positive.code:     write(*,'('' Chest wall RV = '',F8.3,'' L'')') chestwall_restvol/1.0e+6_dp
          
      call write_flow_step_results(init_vol, &
           current_vol,ppl_current,pptrans,Pcw,p_mus,0.0_dp,0.0_dp)
      
      continue = .true.
      do while (continue)
         n = n + 1 ! increment the breath number
         ttime = 0.0_dp ! each breath starts with ttime=0
         endtime = ventilation_values%time_breath * dble(n) - 0.5_dp * dt ! the end time of this breath
         p_mus = 0.0_dp 
         ptrans_frc = SUM(unit_field(nu_pe,1:num_units))/num_units !ptrans at frc
  
         if(n.gt.1)then !write out 'end of breath' information
  
negative.code:     real(dp), intent(in) :: undef
      ! Local variables
      integer :: ne,nunit
      real(dp) :: exp_term,lambda,ratio
      character(len=60) :: sub_name
  
  