Commits

Anonymous committed b388f3b

Updates for general phase

Comments (0)

Files changed (10)

src/pflotran/condition_control.F90

   use MFD_module, only : MFDInitializeMassMatrices
 #endif
 
+  use Global_Aux_module
   use General_Aux_module
   
   implicit none
   PetscReal, pointer :: xx_p(:), iphase_loc_p(:), xx_faces_p(:)
   PetscErrorCode :: ierr
   
-  PetscReal :: temperature, p_sat
   character(len=MAXSTRINGLENGTH) :: string
   
   type(option_type), pointer :: option
   type(patch_type), pointer :: cur_patch
   type(flow_general_condition_type), pointer :: general
   type(dataset_type), pointer :: dataset
+  type(global_auxvar_type) :: global_aux
+  type(general_auxvar_type) :: general_aux
   PetscBool :: use_dataset
   PetscBool :: dataset_flag(realization%option%nflowdof)
   PetscInt :: num_connections
   PetscInt, pointer :: conn_id_ptr(:)
   PetscInt :: ghosted_offset
+  PetscReal :: x(realization%option%nflowdof)
+  PetscReal :: temperature, p_sat
 
   option => realization%option
   discretization => realization%discretization
   field => realization%field
   patch => realization%patch
 
-
+  if (option%iflowmode == G_MODE) then
+    call GlobalAuxVarInit(global_aux,option)
+    call GeneralAuxVarInit(general_aux,option)
+  endif
+  
   cur_level => realization%level_list%first
   do 
     if (.not.associated(cur_level)) exit
                   iphase_loc_p(ghosted_id) = 0
                   cycle
                 endif
+                ! decrement ibegin to give a local offset of 0
+                ibegin = ibegin - 1
                 select case(initial_condition%flow_condition%iphase)
                   case(TWO_PHASE_STATE)
-                    xx_p(ibegin+GENERAL_GAS_PRESSURE_DOF-1) = &
+                    xx_p(ibegin+GENERAL_GAS_PRESSURE_DOF) = &
                       general%gas_pressure%flow_dataset%time_series%cur_value(1)
-                    xx_p(ibegin+GENERAL_GAS_SATURATION_DOF-1) = &
+                    xx_p(ibegin+GENERAL_GAS_SATURATION_DOF) = &
                       general%gas_saturation%flow_dataset%time_series%cur_value(1)
                     temperature = general%temperature%flow_dataset%time_series%cur_value(1)
                     call psat(temperature,p_sat,ierr)
                     ! p_a = p_g - p_s(T)
-                    xx_p(ibegin+GENERAL_AIR_PRESSURE_DOF-1) = &
+                    xx_p(ibegin+GENERAL_AIR_PRESSURE_DOF) = &
                       general%gas_pressure%flow_dataset%time_series%cur_value(1) - &
                       p_sat
                   case(LIQUID_STATE)
-                    xx_p(ibegin+GENERAL_LIQUID_PRESSURE_DOF-1) = &
+                    xx_p(ibegin+GENERAL_LIQUID_PRESSURE_DOF) = &
                       general%liquid_pressure%flow_dataset%time_series%cur_value(1)
-                    xx_p(ibegin+GENERAL_MOLE_FRACTION_DOF-1) = &
+                    xx_p(ibegin+GENERAL_LIQUID_STATE_MOLE_FRACTION_DOF) = &
                       general%mole_fraction%flow_dataset%time_series%cur_value(1)
-                    xx_p(ibegin+GENERAL_TEMPERATURE_DOF-1) = &
+                    xx_p(ibegin+GENERAL_LIQUID_STATE_TEMPERATURE_DOF) = &
                       general%temperature%flow_dataset%time_series%cur_value(1)
                   case(GAS_STATE)
-                    xx_p(ibegin+GENERAL_GAS_PRESSURE_DOF-1) = &
+                    xx_p(ibegin+GENERAL_GAS_PRESSURE_DOF) = &
                       general%gas_pressure%flow_dataset%time_series%cur_value(1)
-                    xx_p(ibegin+GENERAL_AIR_PRESSURE_DOF-1) = &
+                    xx_p(ibegin+GENERAL_AIR_PRESSURE_DOF) = &
                       general%gas_pressure%flow_dataset%time_series%cur_value(1) * &
                       general%mole_fraction%flow_dataset%time_series%cur_value(1)
-                    xx_p(ibegin+GENERAL_TEMPERATURE_DOF-1) = &
+                    xx_p(ibegin+GENERAL_GAS_STATE_TEMPERATURE_DOF) = &
                       general%temperature%flow_dataset%time_series%cur_value(1)
                 end select
                 iphase_loc_p(ghosted_id) = initial_condition%flow_condition%iphase
                   cycle
                 endif
                 xx_p(ibegin:iend) = &
-                      initial_condition%flow_aux_real_var(1:option%nflowdof,iconn)
-                      iphase_loc_p(ghosted_id)=initial_condition%flow_aux_int_var(1,iconn)
+                  initial_condition%flow_aux_real_var(1:option%nflowdof,iconn)
+                iphase_loc_p(ghosted_id) = initial_condition%flow_condition%iphase
                 cur_patch%aux%Global%aux_vars(ghosted_id)%istate = &
-                  int(iphase_loc_p(ghosted_id))
+                  initial_condition%flow_condition%iphase
               enddo
             endif
             initial_condition => initial_condition%next
     enddo
     cur_level => cur_level%next
   enddo
+  
+  if (option%iflowmode == G_MODE) then
+    call GlobalAuxVarStrip(global_aux)
+    call GeneralAuxVarStrip(general_aux)
+  endif  
    
   ! update dependent vectors
   call DiscretizationGlobalToLocal(discretization,field%flow_xx,field%flow_xx_loc,NFLOWDOF)  

src/pflotran/definitions.h

 PetscInt, parameter :: LIQUID_PHASE = 1
 PetscInt, parameter :: GAS_PHASE = 2
 
-! thermodynamic state of fluid ids
-PetscInt, parameter :: LIQUID_STATE = 1
-PetscInt, parameter :: GAS_STATE = 2
-PetscInt, parameter :: TWO_PHASE_STATE = 3
-
 ! approaches to coupling reactive transport
 PetscInt, parameter :: GLOBAL_IMPLICIT = 0
 PetscInt, parameter :: OPERATOR_SPLIT = 1

src/pflotran/general.F90

          GeneralInitializeTimestep, GeneralUpdateAuxVars, &
          GeneralMaxChange, GeneralUpdateSolution, &
          GeneralGetTecplotHeader, GeneralComputeMassBalance, &
-         GeneralDestroy
+         GeneralDestroy, GeneralSetPlotVariables
 
 contains
 
     enddo
     cur_level => cur_level%next
   enddo
+  
+  call GeneralSetPlotVariables(realization)  
 
 end subroutine GeneralSetup
 
 ! date: 03/10/11
 !
 ! ************************************************************************** !
-subroutine GeneralUpdateAuxVars(realization)
+subroutine GeneralUpdateAuxVars(realization,update_state)
 
   use Realization_class
   use Level_module
   use Patch_module
 
   type(realization_type) :: realization
+  PetscBool :: update_state
   
   type(level_type), pointer :: cur_level
   type(patch_type), pointer :: cur_patch
     cur_patch => cur_level%patch_list%first
     do
       if (.not.associated(cur_patch)) exit
-      realization%patch => cur_patch             ! do not update state
-      call GeneralUpdateAuxVarsPatch(realization,PETSC_FALSE)
+      realization%patch => cur_patch             
+      call GeneralUpdateAuxVarsPatch(realization,update_state)
       cur_patch => cur_patch%next
     enddo
     cur_level => cur_level%next
                        porosity_loc_p(ghosted_id),perm_xx_loc_p(ghosted_id), &                       
                        option)
     if (update_state) then
-      call GeneralUpdateState(xx_loc_p(ghosted_start:ghosted_end), &
-                              gen_aux_vars(ZERO_INTEGER,ghosted_id), &
-                              global_aux_vars(ghosted_id), &
-                              patch%saturation_function_array(patch%sat_func_id(ghosted_id))%ptr, &
-                              porosity_loc_p(ghosted_id),perm_xx_loc_p(ghosted_id), &
-                              ghosted_id, &  ! for debugging
-                              option)
+      call GeneralAuxVarUpdateState(xx_loc_p(ghosted_start:ghosted_end), &
+                                    gen_aux_vars(ZERO_INTEGER,ghosted_id), &
+                                    global_aux_vars(ghosted_id), &
+                                    patch%saturation_function_array( &
+                                      patch%sat_func_id(ghosted_id))%ptr, &
+                                    porosity_loc_p(ghosted_id), &
+                                    perm_xx_loc_p(ghosted_id), &
+                                    ghosted_id, &  ! for debugging
+                                    option)
     endif
   enddo
 
       ! set this based on data given 
       global_aux_vars_bc(sum_connection)%istate = &
         boundary_condition%flow_condition%iphase
+      ! update state and update aux var; this could result in two update to 
+      ! the aux var as update state updates if the state changes
+      call GeneralAuxVarUpdateState(xxbc,gen_aux_vars_bc(sum_connection), &
+                         global_aux_vars_bc(sum_connection), &
+                         patch%saturation_function_array(patch%sat_func_id(ghosted_id))%ptr, &
+                         porosity_loc_p(ghosted_id),perm_xx_loc_p(ghosted_id), &                         
+                         ghosted_id,option)
       call GeneralAuxVarCompute(xxbc,gen_aux_vars_bc(sum_connection), &
                          global_aux_vars_bc(sum_connection), &
                          patch%saturation_function_array(patch%sat_func_id(ghosted_id))%ptr, &
   select case(global_aux_var%istate)
     case(LIQUID_STATE)
        x(GENERAL_LIQUID_PRESSURE_DOF) = gen_aux_var(ZERO_INTEGER)%pres(option%liquid_phase)
-       x(GENERAL_MOLE_FRACTION_DOF) = gen_aux_var(ZERO_INTEGER)%xmol(option%air_id,option%liquid_phase)
-       x(GENERAL_TEMPERATURE_DOF) = gen_aux_var(ZERO_INTEGER)%temp
+       x(GENERAL_LIQUID_STATE_MOLE_FRACTION_DOF) = gen_aux_var(ZERO_INTEGER)%xmol(option%air_id,option%liquid_phase)
+       x(GENERAL_LIQUID_STATE_TEMPERATURE_DOF) = gen_aux_var(ZERO_INTEGER)%temp
        pert(GENERAL_LIQUID_PRESSURE_DOF) = 1.d0
-       pert(GENERAL_MOLE_FRACTION_DOF) = -1.d0*perturbation_tolerance*x(GENERAL_MOLE_FRACTION_DOF)
-       pert(GENERAL_TEMPERATURE_DOF) = -1.d0*perturbation_tolerance*x(GENERAL_TEMPERATURE_DOF)
+       pert(GENERAL_LIQUID_STATE_MOLE_FRACTION_DOF) = -1.d0*perturbation_tolerance*x(GENERAL_LIQUID_STATE_MOLE_FRACTION_DOF)
+       pert(GENERAL_LIQUID_STATE_TEMPERATURE_DOF) = -1.d0*perturbation_tolerance*x(GENERAL_LIQUID_STATE_TEMPERATURE_DOF)
     case(GAS_STATE)
        x(GENERAL_GAS_PRESSURE_DOF) = gen_aux_var(ZERO_INTEGER)%pres(option%gas_phase)
        x(GENERAL_AIR_PRESSURE_DOF) = gen_aux_var(ZERO_INTEGER)%pres(option%air_pressure_id)
-       x(GENERAL_TEMPERATURE_DOF) = gen_aux_var(ZERO_INTEGER)%temp
+       x(GENERAL_GAS_STATE_TEMPERATURE_DOF) = gen_aux_var(ZERO_INTEGER)%temp
        pert(GENERAL_GAS_PRESSURE_DOF) = 1.d0
        if (x(GENERAL_GAS_PRESSURE_DOF) - x(GENERAL_AIR_PRESSURE_DOF) > 1.d0) then 
          pert(GENERAL_AIR_PRESSURE_DOF) = 1.d0
        else
          pert(GENERAL_AIR_PRESSURE_DOF) = -1.d0
        endif
-       pert(GENERAL_TEMPERATURE_DOF) = perturbation_tolerance*x(GENERAL_TEMPERATURE_DOF)
+       pert(GENERAL_GAS_STATE_TEMPERATURE_DOF) = perturbation_tolerance*x(GENERAL_GAS_STATE_TEMPERATURE_DOF)
     case(TWO_PHASE_STATE)
        x(GENERAL_GAS_PRESSURE_DOF) = gen_aux_var(ZERO_INTEGER)%pres(option%gas_phase)
        x(GENERAL_AIR_PRESSURE_DOF) = gen_aux_var(ZERO_INTEGER)%pres(option%air_pressure_id)
 #ifdef DEBUG_GENERAL
     call GlobalAuxVarCopy(global_aux_var,global_aux_var_debug,option)
     call GeneralAuxVarCopy(gen_aux_var(idof),general_aux_var_debug,option)
-    call GeneralUpdateState(x_pert,general_aux_var_debug,global_aux_var,&
-                            saturation_function,0.d0,0.d0,ghosted_id,option)
+    call GeneralAuxVarUpdateState(x_pert,general_aux_var_debug, &
+                                  global_aux_var,saturation_function,0.d0, &
+                                  0.d0,ghosted_id,option)
     if (global_aux_var%istate /= global_aux_var_debug%istate) then
       write(option%io_buffer,'(''Change in state: '',i3,'' -> '',i3)') &
         global_aux_var%istate, global_aux_var_debug%istate
               bc_type == CONDUCTANCE_BC) then
                 ! flow in         ! boundary cell is <= pref
             if (delta_pressure > 0.d0 .and. &
-                global_aux_var_up%pres(iphase)-option%reference_pressure < eps) then
+                gen_aux_var_up%pres(iphase)-option%reference_pressure < eps) then
               delta_pressure = 0.d0
             endif
           endif
         if (dabs(aux_vars(idof)) > floweps) then
           v_darcy(iphase) = aux_vars(idof)
           if (v_darcy(iphase) > 0.d0) then 
-            density_ave = global_aux_var_up%den(iphase)
+            density_ave = gen_aux_var_up%den(iphase)
           else 
-            density_ave = global_aux_var_dn%den(iphase)
+            density_ave = gen_aux_var_dn%den(iphase)
           endif 
         endif
     end select
 
   call VecWAXPY(field%flow_dxx,-1.d0,field%flow_xx,field%flow_yy,ierr)
   call VecStrideNorm(field%flow_dxx,ZERO_INTEGER,NORM_INFINITY,option%dpmax,ierr)
-  call VecWAXPY(field%flow_dxx,-1.d0,field%flow_xx,field%flow_yy,ierr)
-  call VecStrideNorm(field%flow_dxx,ONE_INTEGER,NORM_INFINITY,option%dcmax,ierr)
-  call VecWAXPY(field%flow_dxx,-1.d0,field%flow_xx,field%flow_yy,ierr)
-  call VecStrideNorm(field%flow_dxx,TWO_INTEGER,NORM_INFINITY,option%dtmpmax,ierr)
+!  call VecWAXPY(field%flow_dxx,-1.d0,field%flow_xx,field%flow_yy,ierr)
+!  call VecStrideNorm(field%flow_dxx,ONE_INTEGER,NORM_INFINITY,option%dcmax,ierr)
+!  call VecWAXPY(field%flow_dxx,-1.d0,field%flow_xx,field%flow_yy,ierr)
+!  call VecStrideNorm(field%flow_dxx,TWO_INTEGER,NORM_INFINITY,option%dtmpmax,ierr)
 
 end subroutine GeneralMaxChange
 
 ! ************************************************************************** !
 !
-! GeneralUpdateState: Updates the state and swaps primary variables
-! author: Glenn Hammond
-! date: 05/25/11
-!
-! ************************************************************************** !
-subroutine GeneralUpdateState(x,gen_aux_var,global_aux_var, &
-                              saturation_function,por,perm,ghosted_id,option)
-
-  use Option_module
-  use Global_Aux_module
-  use water_eos_module
-  use Gas_Eos_module
-  use Saturation_Function_module
-  
-  implicit none
-
-  type(option_type) :: option
-  PetscInt :: ghosted_id
-  type(saturation_function_type) :: saturation_function
-  type(general_auxvar_type) :: gen_aux_var
-  type(global_auxvar_type) :: global_aux_var
-
-  PetscReal, parameter :: epsilon = 1.d-6
-  PetscReal :: x(option%nflowdof)
-  PetscReal :: por, perm
-  PetscInt :: apid, cpid, vpid
-  PetscInt :: gid, lid, acid, wid, eid
-  PetscReal :: dummy, guess
-  PetscReal :: Ps
-  PetscBool :: flag
-  PetscErrorCode :: ierr
-
-  lid = option%liquid_phase
-  gid = option%gas_phase
-  apid = option%air_pressure_id
-  cpid = option%capillary_pressure_id
-  vpid = option%vapor_pressure_id
-
-  acid = option%air_id ! air component id
-  wid = option%water_id
-  eid = option%energy_id
-
-  flag = PETSC_FALSE
-  
-  select case(global_aux_var%istate)
-    case(LIQUID_STATE)
-      call psat(gen_aux_var%temp,Ps,ierr)
-      if (gen_aux_var%pres(vpid) <= Ps) then
-        global_aux_var%istate = TWO_PHASE_STATE
-!geh        x(GENERAL_GAS_PRESSURE_DOF) = gen_aux_var%pres(vpid)
-!geh        x(GENERAL_AIR_PRESSURE_DOF) = epsilon
-        ! vapor pressure needs to be >= epsilon
-        x(GENERAL_GAS_PRESSURE_DOF) = gen_aux_var%pres(lid)
-        x(GENERAL_AIR_PRESSURE_DOF) = gen_aux_var%pres(lid) - Ps
-        x(GENERAL_GAS_SATURATION_DOF) = epsilon
-        flag = PETSC_TRUE
-#ifdef DEBUG_GENERAL
-        write(option%io_buffer,'(''Liquid -> 2 Phase at Cell '',i11)') ghosted_id
-        call printMsg(option)
-#endif        
-      endif
-    case(GAS_STATE)
-      call psat(gen_aux_var%temp,Ps,ierr)
-      if (gen_aux_var%pres(vpid) >= Ps) then
-        global_aux_var%istate = TWO_PHASE_STATE
-!geh        x(GENERAL_GAS_PRESSURE_DOF) = gen_aux_var%pres(vpid)
-        x(GENERAL_GAS_PRESSURE_DOF) = gen_aux_var%pres(vpid)+gen_aux_var%pres(apid)
-        x(GENERAL_AIR_PRESSURE_DOF) = gen_aux_var%pres(apid)
-        x(GENERAL_GAS_SATURATION_DOF) = 1.d0 - epsilon
-        flag = PETSC_TRUE
-#ifdef DEBUG_GENERAL
-        write(option%io_buffer,'(''Gas -> 2 Phase at Cell '',i11)') ghosted_id
-        call printMsg(option)
-#endif        
-      endif
-    case(TWO_PHASE_STATE)
-      if (gen_aux_var%sat(gid) < 0.d0) then
-        ! convert to liquid state
-        global_aux_var%istate = LIQUID_STATE
-        x(GENERAL_LIQUID_PRESSURE_DOF) = (1.d0+epsilon)* &
-                                         gen_aux_var%pres(vpid)
-        x(GENERAL_MOLE_FRACTION_DOF) = gen_aux_var%xmol(acid,lid)
-        x(GENERAL_TEMPERATURE_DOF) = gen_aux_var%temp
-        flag = PETSC_TRUE
-#ifdef DEBUG_GENERAL
-        write(option%io_buffer,'(''2 Phase -> Liquid at Cell '',i11)') ghosted_id
-        call printMsg(option)
-#endif        
-      else if (gen_aux_var%sat(gid) > 1.d0) then
-        ! convert to gas state
-        global_aux_var%istate = GAS_STATE
-        x(GENERAL_GAS_PRESSURE_DOF) = (1.d0-epsilon)* &
-                                      gen_aux_var%pres(vpid)
-        x(GENERAL_AIR_PRESSURE_DOF) = gen_aux_var%pres(apid)
-        x(GENERAL_TEMPERATURE_DOF) = gen_aux_var%temp
-        flag = PETSC_TRUE
-#ifdef DEBUG_GENERAL
-        write(option%io_buffer,'(''2 Phase -> Gas at Cell '',i11)') ghosted_id
-        call printMsg(option)
-#endif        
-      endif
-  end select
-  
-  if (flag) then
-    call GeneralAuxVarCompute(x,gen_aux_var, global_aux_var,&
-                              saturation_function,por,perm,option)
-    option%variables_swapped = PETSC_TRUE
-  endif
-
-end subroutine GeneralUpdateState
-
-! ************************************************************************** !
-!
 ! GeneralGetTecplotHeader: Returns General Lite contribution to 
 !                               Tecplot file header
 ! author: Glenn Hammond
 
 ! ************************************************************************** !
 !
+! GeneralSetPlotVariables: Adds variables to be printed to list
+! author: Glenn Hammond
+! date: 02/15/13
+!
+! ************************************************************************** !
+subroutine GeneralSetPlotVariables(realization)
+  
+  use Realization_class
+  use Output_Aux_module
+  use Variables_module
+    
+  implicit none
+  
+  type(realization_type) :: realization
+  
+  character(len=MAXWORDLENGTH) :: name, units
+  type(output_variable_list_type), pointer :: list
+  type(output_variable_type), pointer :: output_variable
+  
+  list => realization%output_option%output_variable_list
+
+  if (associated(list%first)) then
+    return
+  endif
+  
+  name = 'Temperature'
+  units = 'C'
+  call OutputVariableAddToList(list,name,OUTPUT_GENERIC,units, &
+                               TEMPERATURE)
+
+  name = 'Liquid Pressure'
+  units = 'Pa'
+  call OutputVariableAddToList(list,name,OUTPUT_PRESSURE,units, &
+                               LIQUID_PRESSURE)
+
+  name = 'Gas Pressure'
+  units = 'Pa'
+  call OutputVariableAddToList(list,name,OUTPUT_PRESSURE,units, &
+                               GAS_PRESSURE)
+
+  name = 'Liquid Saturation'
+  units = ''
+  call OutputVariableAddToList(list,name,OUTPUT_SATURATION,units, &
+                               LIQUID_SATURATION)
+  
+  name = 'Gas Saturation'
+  units = ''
+  call OutputVariableAddToList(list,name,OUTPUT_SATURATION,units, &
+                               GAS_SATURATION)
+  
+  name = 'Liquid Density'
+  units = ''
+  call OutputVariableAddToList(list,name,OUTPUT_GENERIC,units, &
+                               LIQUID_DENSITY)
+  
+  name = 'Gas Density'
+  units = ''
+  call OutputVariableAddToList(list,name,OUTPUT_GENERIC,units, &
+                               GAS_DENSITY)
+  
+  name = 'X_g^l'
+  units = ''
+  call OutputVariableAddToList(list,name,OUTPUT_GENERIC,units, &
+                               LIQUID_MOLE_FRACTION, &
+                               realization%option%air_id)
+  
+  name = 'X_l^l'
+  units = ''
+  call OutputVariableAddToList(list,name,OUTPUT_GENERIC,units, &
+                               LIQUID_MOLE_FRACTION, &
+                               realization%option%water_id)
+  
+  name = 'X_g^g'
+  units = ''
+  call OutputVariableAddToList(list,name,OUTPUT_GENERIC,units, &
+                               GAS_MOLE_FRACTION, &
+                               realization%option%air_id)
+  
+  name = 'X_l^g'
+  units = ''
+  call OutputVariableAddToList(list,name,OUTPUT_GENERIC,units, &
+                               GAS_MOLE_FRACTION, &
+                               realization%option%water_id)
+  
+  name = 'Thermodynamic State'
+  units = ''
+  output_variable => OutputVariableCreate(name,OUTPUT_DISCRETE,units,STATE)
+  output_variable%plot_only = PETSC_TRUE ! toggle output off for observation
+  output_variable%iformat = 1 ! integer
+  call OutputVariableAddToList( &
+         realization%output_option%output_variable_list,output_variable)   
+  
+end subroutine GeneralSetPlotVariables
+
+! ************************************************************************** !
+!
 ! GeneralDestroy: Deallocates variables associated with Richard
 ! author: Glenn Hammond
 ! date: 03/09/11

src/pflotran/general_aux.F90

 
 #include "definitions.h"
 
+  ! thermodynamic state of fluid ids
+  PetscInt, parameter, public :: LIQUID_STATE = 1
+  PetscInt, parameter, public :: GAS_STATE = 2
+  PetscInt, parameter, public :: TWO_PHASE_STATE = 3
+
   PetscInt, parameter, public :: GENERAL_LIQUID_PRESSURE_DOF = 1
   PetscInt, parameter, public :: GENERAL_GAS_PRESSURE_DOF = 1
   PetscInt, parameter, public :: GENERAL_AIR_PRESSURE_DOF = 2
   PetscInt, parameter, public :: GENERAL_GAS_SATURATION_DOF = 3
   PetscInt, parameter, public :: GENERAL_LIQUID_FLUX_DOF = 1
   PetscInt, parameter, public :: GENERAL_GAS_FLUX_DOF = 1
-  PetscInt, parameter, public :: GENERAL_TEMPERATURE_DOF = 3
-  PetscInt, parameter, public :: GENERAL_MOLE_FRACTION_DOF = 2
+  
+  PetscInt, parameter, public :: GENERAL_LIQUID_STATE_TEMPERATURE_DOF = 3
+  PetscInt, parameter, public :: GENERAL_GAS_STATE_TEMPERATURE_DOF = 3
+  PetscInt, parameter, public :: GENERAL_2PHASE_STATE_TEMPERATURE_DOF = 2
+  
+  PetscInt, parameter, public :: GENERAL_LIQUID_STATE_MOLE_FRACTION_DOF = 2
+  
   PetscInt, parameter, public :: GENERAL_LIQUID_CONDUCTANCE_DOF = -1
   PetscInt, parameter, public :: GENERAL_GAS_CONDUCTANCE_DOF = -2
   PetscInt, parameter, public :: GENERAL_FLUX_DOF = 4
   public :: GeneralAuxCreate, GeneralAuxDestroy, &
             GeneralAuxVarCompute, GeneralAuxVarInit, &
             GeneralAuxVarCopy, GeneralAuxVarDestroy, &
-            GeneralAuxVarStrip
+            GeneralAuxVarStrip, GeneralAuxVarUpdateState
 
 contains
 
   aux_var%den = 0.d0
   allocate(aux_var%den_kg(option%nphase))
   aux_var%den_kg = 0.d0
+  ! keep at 25 C.
+  aux_var%temp = 25.d0
   allocate(aux_var%xmol(option%nflowspec,option%nphase))
   aux_var%xmol = 0.d0
   allocate(aux_var%H(option%nphase))
   select case(global_aux_var%istate)
     case(LIQUID_STATE)
       gen_aux_var%pres(lid) = x(GENERAL_LIQUID_PRESSURE_DOF)
-      gen_aux_var%xmol(acid,lid) = x(GENERAL_MOLE_FRACTION_DOF)
-      gen_aux_var%temp = x(GENERAL_TEMPERATURE_DOF)
+      gen_aux_var%xmol(acid,lid) = x(GENERAL_LIQUID_STATE_MOLE_FRACTION_DOF)
+      gen_aux_var%temp = x(GENERAL_LIQUID_STATE_TEMPERATURE_DOF)
 
       gen_aux_var%xmol(wid,lid) = 1.d0 - gen_aux_var%xmol(acid,lid)
       gen_aux_var%sat(lid) = 1.d0
     case(GAS_STATE)
       gen_aux_var%pres(gid) = x(GENERAL_GAS_PRESSURE_DOF)
       gen_aux_var%pres(apid) = x(GENERAL_AIR_PRESSURE_DOF)
-      gen_aux_var%temp = x(GENERAL_TEMPERATURE_DOF)
+      gen_aux_var%temp = x(GENERAL_GAS_STATE_TEMPERATURE_DOF)
 
       gen_aux_var%sat(lid) = 0.d0
       gen_aux_var%sat(gid) = 1.d0
 
 end subroutine GeneralAuxVarCompute
 
+
+! ************************************************************************** !
+!
+! GeneralUpdateState: Updates the state and swaps primary variables
+! author: Glenn Hammond
+! date: 05/25/11
+!
+! ************************************************************************** !
+subroutine GeneralAuxVarUpdateState(x,gen_aux_var,global_aux_var, &
+                                    saturation_function,por,perm,ghosted_id, &
+                                    option)
+
+  use Option_module
+  use Global_Aux_module
+  use water_eos_module
+  use Gas_Eos_module
+  use Saturation_Function_module
+  
+  implicit none
+
+  type(option_type) :: option
+  PetscInt :: ghosted_id
+  type(saturation_function_type) :: saturation_function
+  type(general_auxvar_type) :: gen_aux_var
+  type(global_auxvar_type) :: global_aux_var
+
+  PetscReal, parameter :: epsilon = 1.d-6
+  PetscReal :: x(option%nflowdof)
+  PetscReal :: por, perm
+  PetscInt :: apid, cpid, vpid
+  PetscInt :: gid, lid, acid, wid, eid
+  PetscReal :: dummy, guess
+  PetscReal :: Ps
+  PetscBool :: flag
+  PetscErrorCode :: ierr
+
+  lid = option%liquid_phase
+  gid = option%gas_phase
+  apid = option%air_pressure_id
+  cpid = option%capillary_pressure_id
+  vpid = option%vapor_pressure_id
+
+  acid = option%air_id ! air component id
+  wid = option%water_id
+  eid = option%energy_id
+
+  flag = PETSC_FALSE
+  
+  select case(global_aux_var%istate)
+    case(LIQUID_STATE)
+      call psat(gen_aux_var%temp,Ps,ierr)
+      if (gen_aux_var%pres(vpid) <= Ps) then
+        global_aux_var%istate = TWO_PHASE_STATE
+!geh        x(GENERAL_GAS_PRESSURE_DOF) = gen_aux_var%pres(vpid)
+!geh        x(GENERAL_AIR_PRESSURE_DOF) = epsilon
+        ! vapor pressure needs to be >= epsilon
+        x(GENERAL_GAS_PRESSURE_DOF) = gen_aux_var%pres(lid)
+        x(GENERAL_AIR_PRESSURE_DOF) = gen_aux_var%pres(lid) - Ps
+        x(GENERAL_GAS_SATURATION_DOF) = epsilon
+        flag = PETSC_TRUE
+#ifdef DEBUG_GENERAL
+        write(option%io_buffer,'(''Liquid -> 2 Phase at Cell '',i11)') ghosted_id
+        call printMsg(option)
+#endif        
+      endif
+    case(GAS_STATE)
+      call psat(gen_aux_var%temp,Ps,ierr)
+      if (gen_aux_var%pres(vpid) >= Ps) then
+        global_aux_var%istate = TWO_PHASE_STATE
+!geh        x(GENERAL_GAS_PRESSURE_DOF) = gen_aux_var%pres(vpid)
+        x(GENERAL_GAS_PRESSURE_DOF) = gen_aux_var%pres(vpid)+gen_aux_var%pres(apid)
+        x(GENERAL_AIR_PRESSURE_DOF) = gen_aux_var%pres(apid)
+        x(GENERAL_GAS_SATURATION_DOF) = 1.d0 - epsilon
+        flag = PETSC_TRUE
+#ifdef DEBUG_GENERAL
+        write(option%io_buffer,'(''Gas -> 2 Phase at Cell '',i11)') ghosted_id
+        call printMsg(option)
+#endif        
+      endif
+    case(TWO_PHASE_STATE)
+      if (gen_aux_var%sat(gid) < 0.d0) then
+        ! convert to liquid state
+        global_aux_var%istate = LIQUID_STATE
+        x(GENERAL_LIQUID_PRESSURE_DOF) = (1.d0+epsilon)* &
+                                         gen_aux_var%pres(vpid)
+        x(GENERAL_LIQUID_STATE_MOLE_FRACTION_DOF) = &
+          gen_aux_var%xmol(acid,lid)
+        x(GENERAL_LIQUID_STATE_TEMPERATURE_DOF) = gen_aux_var%temp
+        flag = PETSC_TRUE
+#ifdef DEBUG_GENERAL
+        write(option%io_buffer,'(''2 Phase -> Liquid at Cell '',i11)') ghosted_id
+        call printMsg(option)
+#endif        
+      else if (gen_aux_var%sat(gid) > 1.d0) then
+        ! convert to gas state
+        global_aux_var%istate = GAS_STATE
+        x(GENERAL_GAS_PRESSURE_DOF) = (1.d0-epsilon)* &
+                                      gen_aux_var%pres(vpid)
+        x(GENERAL_AIR_PRESSURE_DOF) = gen_aux_var%pres(apid)
+        x(GENERAL_GAS_STATE_TEMPERATURE_DOF) = gen_aux_var%temp
+        flag = PETSC_TRUE
+#ifdef DEBUG_GENERAL
+        write(option%io_buffer,'(''2 Phase -> Gas at Cell '',i11)') ghosted_id
+        call printMsg(option)
+#endif        
+      endif
+  end select
+  
+  if (flag) then
+    call GeneralAuxVarCompute(x,gen_aux_var, global_aux_var,&
+                              saturation_function,por,perm,option)
+    option%variables_swapped = PETSC_TRUE
+  endif
+
+end subroutine GeneralAuxVarUpdateState
+
 ! ************************************************************************** !
 !
 ! GeneralAuxVarSingleDestroy: Deallocates a mode auxiliary object
   if (associated(aux_vars)) then
     do iaux = 1, size(aux_vars,2)
       do idof = 1, size(aux_vars,1)
-        call GeneralAuxVarStrip(aux_vars(idof,iaux))
+        call GeneralAuxVarStrip(aux_vars(idof-1,iaux))
       enddo
     enddo  
     deallocate(aux_vars)

src/pflotran/hydrostatic.F90

                     temperature_gradient(X_DIRECTION)*dist_x + & ! gradient in K/m
                     temperature_gradient(Y_DIRECTION)*dist_y + &
                     temperature_gradient(Z_DIRECTION)*dist_z 
-        coupler%flow_aux_real_var(GENERAL_TEMPERATURE_DOF,iconn) = temperature
-        coupler%flow_aux_real_var(GENERAL_MOLE_FRACTION_DOF,iconn) = concentration_at_datum
+        coupler%flow_aux_real_var(GENERAL_LIQUID_STATE_TEMPERATURE_DOF,iconn) = &
+          temperature
+        coupler%flow_aux_real_var(GENERAL_LIQUID_STATE_MOLE_FRACTION_DOF,iconn) = &
+          concentration_at_datum
 
         coupler%flow_aux_int_var(GENERAL_LIQUID_PRESSURE_DOF,iconn) = condition%iphase
       case default
-        coupler%flow_aux_int_var(1,iconn) = 1
+        coupler%flow_aux_int_var(COUPLER_IPHASE_INDEX,iconn) = 1
     end select
 
   enddo

src/pflotran/init.F90

       case(FLASH2_MODE)
         call Flash2UpdateAuxVars(realization)
       case(G_MODE)
-        call GeneralUpdateAuxVars(realization)
+        call GeneralUpdateAuxVars(realization,PETSC_TRUE)
     end select
   else ! no flow mode specified
     if (len_trim(realization%nonuniform_velocity_filename) > 0) then

src/pflotran/level.F90

   enddo
 
 end subroutine LevelConvertListToArray
-#if 0
-! ************************************************************************** !
-!
-! LevelProcessCouplers: Deallocates a realization
-! author: Glenn Hammond
-! date: 02/22/08
-!
-! ************************************************************************** !
-subroutine LevelProcessCouplers(level,flow_conditions,transport_conditions, &
-                                materials,option)
 
-  use Option_module
-  use Material_module
-  use Condition_module
-  
-  implicit none
-  
-  type(level_type) :: level
-  type(material_property_type), pointer :: materials(:)
-  type(condition_list_type) :: flow_conditions
-  type(condition_list_type) :: transport_conditions
-  type(option_type) :: option  
-  
-  type(patch_type), pointer :: cur_patch
-
-  cur_patch => level%patch_list%first
-  do 
-    if (.not.associated(cur_patch)) exit
-    call PatchProcessCouplers(cur_patch,flow_conditions,transport_conditions, &
-                              materials,option)
-    cur_patch => cur_patch%next
-  enddo
- 
-end subroutine LevelProcessCouplers
-
-! ************************************************************************** !
-!
-! LevelInitAllCouplerAuxVars: Initializes coupler auxillary variables 
-!                                within list
-! author: Glenn Hammond
-! date: 02/22/08
-!
-! ************************************************************************** !
-subroutine LevelInitAllCouplerAuxVars(level,option)
-
-  use Option_module
-  
-  implicit none
-  
-  type(level_type) :: level
-  type(option_type) :: option
-  
-  type(patch_type), pointer :: cur_patch
-
-  cur_patch => level%patch_list%first
-  do 
-    if (.not.associated(cur_patch)) exit
-    call PatchInitAllCouplerAuxVars(cur_patch,option)
-    cur_patch => cur_patch%next
-  enddo
- 
-end subroutine LevelInitAllCouplerAuxVars
-
-! ************************************************************************** !
-!
-! LevelInitCouplerAuxVars: Initializes coupler auxillary variables 
-!                                within list
-! author: Glenn Hammond
-! date: 02/22/08
-!
-! ************************************************************************** !
-subroutine LevelInitCouplerAuxVars(level,option)
-
-  use Option_module
-  
-  implicit none
-  
-  type(level_type) :: level
-  type(option_type) :: option
-  
-  type(patch_type), pointer :: cur_patch
-
-  cur_patch => level%patch_list%first
-  do 
-    if (.not.associated(cur_patch)) exit
-    call PatchInitCouplerAuxVars(cur_patch,option)
-    cur_patch => cur_patch%next
-  enddo
- 
-end subroutine LevelInitCouplerAuxVars
-
-! ************************************************************************** !
-!
-! LevelUpdateAllCouplerAuxVars: Updates auxiliary variables associated 
-!                                  with couplers in lis
-! author: Glenn Hammond
-! date: 02/22/08
-!
-! ************************************************************************** !
-subroutine LevelUpdateAllCouplerAuxVars(level,force_update_flag,option)
-
-  use Option_module
-
-  implicit none
-  
-  type(level_type) :: level
-  PetscBool :: force_update_flag
-  type(option_type) :: option
-    
-  type(patch_type), pointer :: cur_patch
-
-  cur_patch => level%patch_list%first
-  do 
-    if (.not.associated(cur_patch)) exit
-    call PatchUpdateAllCouplerAuxVars(cur_patch,force_update_flag,option)
-    cur_patch => cur_patch%next
-  enddo
- 
-end subroutine LevelUpdateAllCouplerAuxVars
-
-! ************************************************************************** !
-!
-! LevelUpdateCouplerAuxVars: Updates auxiliary variables associated 
-!                                  with couplers in lis
-! author: Glenn Hammond
-! date: 02/22/08
-!
-! ************************************************************************** !
-subroutine LevelUpdateCouplerAuxVars(level,force_update_flag,field,option)
-
-  use Option_module
-
-  implicit none
-  
-  type(level_type) :: level
-  PetscBool :: force_update_flag
-  type(field_type) :: field
-  type(option_type) :: option
-    
-  type(patch_type), pointer :: cur_patch
-
-  cur_patch => level%patch_list%first
-  do 
-    if (.not.associated(cur_patch)) exit
-    call PatchUpdateCouplerAuxVars(cur_patch,force_update_flag,field,option)
-    cur_patch => cur_patch%next
-  enddo
- 
-end subroutine LevelUpdateCouplerAuxVars
-#endif
 ! ************************************************************************** !
 !
 ! LevelDestroyList: Deallocates a level list and array of leveles

src/pflotran/patch.F90

   use Dataset_Aux_module
   
   use General_Aux_module
+  use Grid_module
 
   implicit none
   
   PetscBool :: update
   PetscBool :: dof1, dof2, dof3
   PetscReal :: temperature, p_sat
+  PetscReal :: x(option%nflowdof)
   character(len=MAXSTRINGLENGTH) :: string, string2
   PetscErrorCode :: ierr
   
   PetscInt :: idof, num_connections,sum_connection
+  PetscInt :: iconn, local_id, ghosted_id
   
 
   if (.not.associated(coupler_list)) return
                     call printErrMsg(option)
                   endif
                   call HydrostaticUpdateCoupler(coupler,option,patch%grid)
+#if 0                  
+                  do iconn=1,coupler%connection_set%num_connections
+                    local_id = coupler%connection_set%id_dn(iconn)
+                    ghosted_id = patch%grid%nL2G(local_id)
+                    x(1:option%nflowdof) = &
+                      coupler%flow_aux_real_var(1:option%nflowdof,iconn)
+                    patch%aux%Global%aux_vars(ghosted_id)%istate = LIQUID_STATE
+                    call GeneralAuxVarUpdateState(x, &
+                           patch%aux%General%aux_vars(ZERO_INTEGER,ghosted_id), &
+                           patch%aux%Global%aux_vars(ghosted_id), &
+                           patch%saturation_function_array( &
+                             patch%sat_func_id(ghosted_id))%ptr, &
+                           0.d0,0.d0,ghosted_id,option)
+                    coupler%flow_aux_real_var(1:option%nflowdof,iconn) = &
+                      x(1:option%nflowdof)
+                    coupler%flow_aux_int_var(ONE_INTEGER,iconn) = &
+                      patch%aux%Global%aux_vars(ghosted_id)%istate  
+                    patch%aux%Global%aux_vars(ghosted_id)%istate = &
+                      patch%aux%Global%aux_vars(ghosted_id)%istate    
+                  enddo
+#endif                  
                 else
                   select case(general%liquid_pressure%itype)
                     case(DIRICHLET_BC)
                   end select
                   select case(general%mole_fraction%itype)
                     case(DIRICHLET_BC)
-                      coupler%flow_aux_real_var(GENERAL_MOLE_FRACTION_DOF,1:num_connections) = &
+                      coupler%flow_aux_real_var(GENERAL_LIQUID_STATE_MOLE_FRACTION_DOF,1:num_connections) = &
                         general%mole_fraction%flow_dataset%time_series%cur_value(1)
                       dof2 = PETSC_TRUE
                   end select
                   select case(general%temperature%itype)
                     case(DIRICHLET_BC)
-                      coupler%flow_aux_real_var(GENERAL_TEMPERATURE_DOF,1:num_connections) = &
+                      coupler%flow_aux_real_var(GENERAL_LIQUID_STATE_TEMPERATURE_DOF,1:num_connections) = &
                         general%temperature%flow_dataset%time_series%cur_value(1)
                       dof3 = PETSC_TRUE
                   end select
                 end select                
                 select case(general%temperature%itype)
                   case(DIRICHLET_BC)
-                    coupler%flow_aux_real_var(GENERAL_TEMPERATURE_DOF,1:num_connections) = &
+                    coupler%flow_aux_real_var(GENERAL_GAS_STATE_TEMPERATURE_DOF,1:num_connections) = &
                       general%temperature%flow_dataset%time_series%cur_value(1)
                     dof3 = PETSC_TRUE
                 end select
             end select  
-            if (.not.dof1 .and. .not.dof2 .and. .not.dof3 .and. &
-                general%liquid_pressure%itype /= HYDROSTATIC_BC) then
+            !geh: is this really correct, or should it be .or.
+            if (.not.dof1 .or. .not.dof2 .or. .not.dof3) then
               option%io_buffer = 'Error with general phase boundary condition'
-            else if (general%liquid_pressure%itype == HYDROSTATIC_BC) then
-              call HydrostaticUpdateCoupler(coupler,option,patch%grid)
             endif
             
           case(MPH_MODE,IMS_MODE,FLASH2_MODE,THC_MODE) ! updated 10/17/11
   use Reaction_module
   use Reactive_Transport_Aux_module  
   use Surface_Complexation_Aux_module
+  use General_Aux_module
   use Output_Aux_module
   use Variables_module
   
             enddo
           case(LIQUID_PRESSURE)
             do local_id=1,grid%nlmax
-              ghosted_id = grid%nL2G(local_id)
-              istate = patch%aux%Global%aux_vars(ghosted_id)%istate
-              select case(istate)
-                case(LIQUID_STATE,TWO_PHASE_STATE)
-                  vec_ptr(local_id) = patch%aux%General%aux_vars(ZERO_INTEGER,ghosted_id)%pres(option%liquid_phase)
-                case(GAS_STATE)
-                  vec_ptr(local_id) = patch%aux%General%aux_vars(ZERO_INTEGER,ghosted_id)%pres(option%gas_phase)
-              end select
+              vec_ptr(local_id) = patch%aux%General%aux_vars(ZERO_INTEGER,grid%nL2G(local_id))%pres(option%liquid_phase)
+            enddo
+          case(GAS_PRESSURE)
+            do local_id=1,grid%nlmax
+              vec_ptr(local_id) = patch%aux%General%aux_vars(ZERO_INTEGER,grid%nL2G(local_id))%pres(option%gas_phase)
             enddo
           case(STATE)
             do local_id=1,grid%nlmax
   use Surface_Complexation_Aux_module
   use Output_Aux_module
   use Variables_module
+  use General_Aux_module  
 
   implicit none
 
           case(TEMPERATURE)
             value = patch%aux%General%aux_vars(ZERO_INTEGER,ghosted_id)%temp
           case(LIQUID_PRESSURE)
-            istate = patch%aux%Global%aux_vars(ghosted_id)%istate
-            select case(istate)
-              case(LIQUID_STATE,TWO_PHASE_STATE)
-                value = patch%aux%General%aux_vars(ZERO_INTEGER,ghosted_id)%pres(option%liquid_phase)
-              case(GAS_STATE)
-                value = patch%aux%General%aux_vars(ZERO_INTEGER,ghosted_id)%pres(option%gas_phase)
-            end select
+            value = patch%aux%General%aux_vars(ZERO_INTEGER,ghosted_id)%pres(option%liquid_phase)
+          case(GAS_PRESSURE)
+            value = patch%aux%General%aux_vars(ZERO_INTEGER,ghosted_id)%pres(option%gas_phase)
           case(STATE)
             value = patch%aux%Global%aux_vars(ghosted_id)%istate
           case(LIQUID_SATURATION)
   use Option_module
   use Field_module
   use Variables_module
+  use General_Aux_module  
 
   implicit none
 
             enddo
           case(LIQUID_PRESSURE)
             do local_id=1,grid%nlmax
-              ghosted_id = grid%nL2G(local_id)
-              istate = patch%aux%Global%aux_vars(ghosted_id)%istate
-              select case(istate)
-                case(LIQUID_STATE,TWO_PHASE_STATE)
-                  patch%aux%General%aux_vars(ZERO_INTEGER,ghosted_id)%pres(option%liquid_phase) = vec_ptr(local_id)
-                case(GAS_STATE)
-                  patch%aux%General%aux_vars(ZERO_INTEGER,ghosted_id)%pres(option%gas_phase) = vec_ptr(local_id)
-              end select
+              patch%aux%General%aux_vars(ZERO_INTEGER,grid%nL2G(local_id))%pres(option%liquid_phase) = vec_ptr(local_id)
+            enddo
+          case(GAS_PRESSURE)
+            do local_id=1,grid%nlmax
+              patch%aux%General%aux_vars(ZERO_INTEGER,grid%nL2G(local_id))%pres(option%gas_phase) = vec_ptr(local_id)
             enddo
           case(STATE)
             do local_id=1,grid%nlmax

src/pflotran/saturation_function.F90

   interface SaturationFunctionCompute
     module procedure SaturationFunctionCompute1
     module procedure SaturationFunctionCompute2
+    module procedure SaturationFunctionCompute3
   end interface
   
   public :: SaturationFunctionCreate, &
 
 end subroutine SaturationFunctionCompute1
 
-
 ! ************************************************************************** !
 !
 ! SaturationFunctionCompute2: Computes the saturation and relative permeability
 
 ! ************************************************************************** !
 !
+! SaturationFunctionCompute3: Computes just saturation as a function of 
+!                             capillary pressure
+! author: Glenn Hammond
+! date: 2/9/12
+!
+! ************************************************************************** !
+subroutine SaturationFunctionCompute3(pressure,saturation, &
+                                      saturation_function, &
+                                      option)
+  use Option_module
+  
+  implicit none
+
+  PetscReal :: pressure, saturation
+  PetscReal :: relative_perm_dummy, dsat_pres_dummy, dkr_pres_dummy
+  type(saturation_function_type) :: saturation_function
+  PetscReal :: auxvar1_dummy, auxvar2_dummy
+  type(option_type) :: option
+
+  PetscBool :: switch_to_saturated_dummy
+  
+  call SaturationFunctionCompute2(pressure,saturation,relative_perm_dummy, &
+                                  dsat_pres_dummy,dkr_pres_dummy, &
+                                  saturation_function, &
+                                  auxvar1_dummy,auxvar2_dummy, &
+                                  switch_to_saturated_dummy,option)
+
+end subroutine SaturationFunctionCompute3
+
+! ************************************************************************** !
+!
 ! SaturationFunctionComputeIce:Computes the saturation of ice, water vapor 
 !                              and liquid water given the saturation function
 !                              temperature, water vapor pressure and liquid

src/pflotran/timestepper.F90

     case(RICHARDS_MODE)
       call RichardsUpdateAuxVars(realization)
     case(G_MODE)
-      call GeneralUpdateAuxVars(realization)
+                                            ! do not update state
+      call GeneralUpdateAuxVars(realization,PETSC_FALSE)
   end select    
 
 end subroutine StepperUpdateFlowAuxVars