1. Satish Karra
  2. karra-pflotran-dev

Commits

Satish Karra  committed e520af4 Merge

Merge with pflotran-dev.

  • Participants
  • Parent commits 0014766, 11b8254
  • Branches default

Comments (0)

Files changed (35)

File docs/user_manual/governing_eqns.tex

View file
  • Ignore whitespace
 \EQ
 \bD_\a \eq %\varphi s 
 %\Big(
-\tau D_m + a_T \bI + \big(a_L-a_T\big)\frac{\bv\bv}{v},
+\tau D_m + a_T v\bI + \big(a_L-a_T\big)\frac{\bv\bv}{v},
 %\Big),
 \EN
 with longitudinal and transverse dispersivity coefficients $a_L$, $a_T$, respectively, $\tau$ refers to tortuosity, and $D_m$ to the molecular diffusion coefficient. Currently, only longitudinal dispersion is implemented in PFLOTRAN.

File docs/user_manual/installation.tex

View file
  • Ignore whitespace
 \end{Verbatim}
 where {\tt simn.in} is the usual PFLOTRAN input file. The names may be arbitrarily chosen. Then, launch the run as:
 
-\verb|mpirun -n XXX pflotran -pflotranin filename.in -multisimulation| 
+\verb|mpirun -n XXX pflotran -pflotranin filenames.in -multisimulation| 
 
 \noindent  
-Note that all simulations must run at once. The same logic used to allow a processor group to run multiple simulations with {\tt multirealization} is not implemented. Chose the number of cores to be a multiple of the number of input files. All output files have {\tt Gn} appended to the file name, e.g. {\tt sim10G10-001.tec}.
+Note that all simulations run at once. The same logic used to allow a processor group to run multiple simulations with {\tt multirealization} is {\em not} implemented. {\bf Choose the number of cores to be a multiple of the number of input filenames listed in the input file (i.e. filenames.in).} All output files have {\tt Gn} appended to the file name, e.g. {\tt sim10G10-001.tec}.
 
 \begin{comment}
 \subsection{Building SAMRAI Version 2.4.4}

File docs/user_manual/keywords.tex

View file
  • Ignore whitespace
 %\midrule
 \bf MPHASE (MPH) &Multiphase supercritical CO$_2$-brine-energy based on variable switching for phase changes\\
 %\midrule
-\bf FLASH2 &Multiphase supercritical CO$_2$-brine-energy based on the flash method for phase changes with a persistent set of unknowns---required for AMR\\
+\bf FLASH2 &Multiphase supercritical CO$_2$-brine-energy based on the flash method for phase changes with a persistent set of unknowns\\
 %\midrule
 \bf THC &Thermo-Hydro-Chemical coupled groundwater flow, thermal and solute transport\\
 %\midrule
 \item[NO\_PRINT\_INITIAL] \ the initial state of the system will not be printed to the output file if this card is activated
 \item[NO\_PRINT\_FINAL] \ the final state of the system will not be printed to the output file if this card is activated
 \item[PRINT\_COLUMN\_IDS] print column numbers in observation and mass balance output files
-\item[FORMAT] \ <file format>: \ specify the snapshot in time file type. Options available are: 
+\item[FORMAT] \ <file format>: \ specify the snapshot in time file type. File formats available are: 
 \begin{deflist}{0000000000000000000000000}
 \item[TECPLOT BLOCK] -TecPlot block format
 \item[TECPLOT POINT] -TecPlot point format (requires a single processor)

File regression_tests/default/543/543_hanford_srfcplx_param-np8.regression.gold

View file
  • Ignore whitespace
        56:   3.3390000000000E-01
        54:   3.5780000000000E-01
        59:   3.5940000000000E-01
--- GENERIC: Permeability X --
+-- GENERIC: Permeability --
       Max:   1.0000000000000E-09
       Min:   1.0000000000000E-12
      Mean:   9.6412762669168E-11

File regression_tests/default/543/543_hanford_srfcplx_param.regression.gold

View file
  • Ignore whitespace
        29:   3.2550000000000E-01
         1:   3.9410000000000E-01
        31:   2.9380000000000E-01
--- GENERIC: Permeability X --
+-- GENERIC: Permeability --
       Max:   1.0000000000000E-09
       Min:   1.0000000000000E-12
      Mean:   9.6412762669168E-11

File regression_tests/ngee/mass_transfer_1d_in_z.in

View file
  • Ignore whitespace
 END
 
 :=========================== mass transfer ====================================
-MASS_TRANSFER A(im)
+RT_MASS_TRANSFER A(im)
   IDOF 2
   FILENAME mass_transfer_1d_in_z_data.h5
   DATASET A(im) 
 END
 
-MASS_TRANSFER B(im)
+RT_MASS_TRANSFER B(im)
   IDOF 3
   FILENAME mass_transfer_1d_in_z_data.h5
   DATASET B(im) 
 END
 
-MASS_TRANSFER C(im)
+RT_MASS_TRANSFER C(im)
   IDOF 4
   FILENAME mass_transfer_1d_in_z_data.h5
   DATASET C(im) 

File regression_tests_refactor/default/543/543_hanford_srfcplx_param-np8.regression.gold

View file
  • Ignore whitespace
        56:   3.3390000000000E-01
        54:   3.5780000000000E-01
        59:   3.5940000000000E-01
--- GENERIC: Permeability X --
+-- GENERIC: Permeability --
       Max:   1.0000000000000E-09
       Min:   1.0000000000000E-12
      Mean:   9.6412762669168E-11

File regression_tests_refactor/default/543/543_hanford_srfcplx_param.regression.gold

View file
  • Ignore whitespace
        29:   3.2550000000000E-01
         1:   3.9410000000000E-01
        31:   2.9380000000000E-01
--- GENERIC: Permeability X --
+-- GENERIC: Permeability --
       Max:   1.0000000000000E-09
       Min:   1.0000000000000E-12
      Mean:   9.6412762669168E-11

File regression_tests_refactor/ngee/mass_transfer_1d_in_z.in

View file
  • Ignore whitespace
 END
 
 :=========================== mass transfer ====================================
-MASS_TRANSFER A(im)
+RT_MASS_TRANSFER A(im)
   IDOF 2
   FILENAME mass_transfer_1d_in_z_data.h5
   DATASET A(im) 
 END
 
-MASS_TRANSFER B(im)
+RT_MASS_TRANSFER B(im)
   IDOF 3
   FILENAME mass_transfer_1d_in_z_data.h5
   DATASET B(im) 
 END
 
-MASS_TRANSFER C(im)
+RT_MASS_TRANSFER C(im)
   IDOF 4
   FILENAME mass_transfer_1d_in_z_data.h5
   DATASET C(im) 

File src/pflotran/hydrogeophysics_factory.F90

View file
  • Ignore whitespace
   use Hydrogeophysics_Wrapper_module
   use Input_module
   use Simulation_Base_class 
+  use String_module
   
   implicit none
 
 
   class(hydrogeophysics_simulation_type), pointer :: simulation
   PetscMPIInt :: mycolor_mpi, mykey_mpi
-  PetscInt :: i, num_subsurface_processes, offset
+  PetscInt :: i, num_subsurface_processes, offset, num_slaves
   PetscInt :: local_size
+  PetscBool :: option_found
   PetscMPIInt :: mpi_int, process_range(3)
+  character(len=MAXSTRINGLENGTH) :: string
   PetscErrorCode :: ierr
   IS :: is
   Vec :: pflotran_solution_vec_mpi, pflotran_solution_vec_seq
   VecScatter :: pflotran_scatter
-
-  ! NOTE: PETSc must already have been initialized here!
-  if (option%global_commsize < 3) then
-    option%io_buffer = 'Must of at least processes allocates to ' // &
-      'simulation in order to run hydrogeophysics.'
-    call printErrMsg(option)
-  endif
-
+  
 #ifndef E4D
   option%io_buffer = 'Must compile with E4D defined during preprocessing ' // &
     'step in order to use HYDROGEOPHYSICS.'
   call printErrMsg(option)
 #endif
 
+  string = '-num_slaves'
+  num_slaves = -999
+  call InputGetCommandLineInt(string,i,option_found,option)
+  if (option_found) num_slaves = i
+
+  ! NOTE: PETSc must already have been initialized here!
+  if (option%global_commsize < 3) then
+    option%io_buffer = 'At least 3 processes must be allocated to ' // &
+      'simulation in order to run hydrogeophysics.'
+    call printErrMsg(option)
+  endif
+
   simulation => HydrogeophysicsCreate(option)
   
-  num_subsurface_processes = 1
+  if (num_slaves < -998) then
+    num_subsurface_processes = option%global_commsize / 2
+    num_slaves = option%global_commsize - num_subsurface_processes - 1
+  else if (num_slaves <= 0) then
+    option%io_buffer = 'Number of slaves must be greater than zero. ' // &
+      'Currently set to ' // StringFormatInt(num_slaves) // '.'
+    call printErrMsg(option)
+  else
+    if (num_slaves > option%global_commsize - 2) then
+      option%io_buffer = 'Too many slave processes allocated to ' // &
+        'simulation: ' // StringFormatInt(num_slaves)
+      call printErrMsg(option)
+    endif
+    num_subsurface_processes = option%global_commsize - num_slaves - 1
+  endif
+  
+  write(option%io_buffer,*) 'Number of E4D processes: ', &
+    StringFormatInt(num_slaves+1)
+  call printMsg(option)
+  write(option%io_buffer,*) 'Number of PFLOTRAN processes: ', &
+    StringFormatInt(num_subsurface_processes)
+  call printMsg(option)
   
   ! split the communicator
   option%mygroup_id = 0

File src/pflotran/hydrogeophysics_wrapper.F90

View file
  • Ignore whitespace
 !  call printMsg(option,'HydrogeophysicsWrapperStep()')
   
   ! Bcasting a 1 to E4D tells it to run a forward simulation 
-  call MPI_Bcast(ONE_INTEGER_MPI,ONE_INTEGER_MPI,MPI_INTEGER, &
-                 ZERO_INTEGER_MPI,comm,ierr)
+  if (comm /= MPI_COMM_NULL) then
+    call MPI_Bcast(ONE_INTEGER_MPI,ONE_INTEGER_MPI,MPI_INTEGER, &
+                   ZERO_INTEGER_MPI,comm,ierr)
+  endif
   call VecScatterBegin(scatter,solution_mpi,solution_seq, &
                        INSERT_VALUES,SCATTER_FORWARD,ierr)
   call VecScatterEnd(scatter,solution_mpi,solution_seq, &

File src/pflotran/init.F90

View file
  • Ignore whitespace
   endif
   if (realization%output_option%print_permeability) then
     ! add permeability to header
-    call OutputVariableAddToList( &
-           realization%output_option%output_variable_list, &
-           'Permeability X',OUTPUT_GENERIC,'m^2',PERMEABILITY)
+    if (MaterialAnisotropyExists(realization%material_properties)) then
+      call OutputVariableAddToList( &
+             realization%output_option%output_variable_list, &
+             'Permeability X',OUTPUT_GENERIC,'m^2',PERMEABILITY)
+      call OutputVariableAddToList( &
+             realization%output_option%output_variable_list, &
+             'Permeability Y',OUTPUT_GENERIC,'m^2',PERMEABILITY_Y)
+      call OutputVariableAddToList( &
+             realization%output_option%output_variable_list, &
+             'Permeability Z',OUTPUT_GENERIC,'m^2',PERMEABILITY_Z)
+#ifdef DASVYAT
+      if(realization%discretization%itype == STRUCTURED_GRID_MIMETIC .or. &
+         realization%discretization%itype == UNSTRUCTURED_GRID_MIMETIC) then
+        call OutputVariableAddToList( &
+               realization%output_option%output_variable_list, &
+               'Permeability XY',OUTPUT_GENERIC,'m^2',PERMEABILITY_XY)
+        call OutputVariableAddToList( &
+               realization%output_option%output_variable_list, &
+               'Permeability XZ',OUTPUT_GENERIC,'m^2',PERMEABILITY_XZ)
+        call OutputVariableAddToList( &
+               realization%output_option%output_variable_list, &
+               'Permeability YZ',OUTPUT_GENERIC,'m^2',PERMEABILITY_YZ)
+      endif
+#endif
+    else
+      call OutputVariableAddToList( &
+             realization%output_option%output_variable_list, &
+             'Permeability',OUTPUT_GENERIC,'m^2',PERMEABILITY)
+    endif
   endif
   if (realization%output_option%print_iproc) then
     output_variable => OutputVariableCreate('Processor ID',OUTPUT_DISCRETE,'', &
         call printErrMsgByRank(option)
     end select
   endif ! option%nsurfflowdof > 0
+
+  if (simulation%surf_realization%output_option%print_iproc) then
+    output_variable => OutputVariableCreate('Processor ID',OUTPUT_DISCRETE,'', &
+                                            PROCESSOR_ID)
+    output_variable%plot_only = PETSC_TRUE ! toggle output off for observation
+    output_variable%iformat = 1 ! integer
+    call OutputVariableAddToList( &
+           simulation%surf_realization%output_option%output_variable_list,output_variable)
+  endif
+
 #endif
 
   call printMsg(option," ")
   call printMsg(option,"  Finished Initialization")
   call PetscLogEventEnd(logging%event_init,ierr)
 
-#if 0
-  !----------------------------------------------------------------------------!
-  ! This section for setting up new process model approach
-  !----------------------------------------------------------------------------!
-  simulation%output_option => realization%output_option
-  simulation%option => realization%option
-!  simulation%synchronizer => SynchronizerCreate()
-!  simulation%synchronizer%option => realization%option
-!  simulation%synchronizer%output_option => realization%output_option
-!  simulation%synchronizer%waypoints => WaypointListCopy(realization%waypoints)
-!  simulation%synchronizer%waypoints => realization%waypoints
-  nullify(cur_process_model)
-
-  nullify(surf_flow_process_model_coupler)
-  nullify(sub_flow_process_model_coupler)
-  nullify(sub_tran_process_model_coupler)
-
-#ifdef SURFACE_FLOW
-  ! Create Surface-flow ProcessModel & ProcessModelCoupler
-  if (option%nsurfflowdof > 0) then
-    select case(option%iflowmode)
-      case(RICHARDS_MODE)
-        cur_process_model => PMSurfaceFlowCreate()
-      case(TH_MODE)
-        !cur_process_model => PMSurfaceTHCreate()
-    end select
-    cur_process_model%option => surf_realization%option
-    cur_process_model%output_option => surf_realization%output_option
-
-    surf_flow_process_model_coupler => ProcessModelCouplerCreate()
-    surf_flow_process_model_coupler%option => option
-    surf_flow_process_model_coupler%process_model_list => cur_process_model
-    surf_flow_process_model_coupler%pm_ptr%ptr => cur_process_model
-    surf_flow_process_model_coupler%depth = 0
-    nullify(cur_process_model)
-  endif
-#endif
-
-  ! Create Subsurface-flow ProcessModel & ProcessModelCoupler
-  if (option%nflowdof > 0) then
-    select case(option%iflowmode)
-      case(RICHARDS_MODE)
-        cur_process_model => PMRichardsCreate()
-      case(TH_MODE)
-        cur_process_model => PMTHCreate()
-      case(THC_MODE)
-        cur_process_model => PMTHCCreate()
-    end select
-    cur_process_model%option => realization%option
-    cur_process_model%output_option => realization%output_option
-
-    sub_flow_process_model_coupler => ProcessModelCouplerCreate()
-    sub_flow_process_model_coupler%option => option
-    sub_flow_process_model_coupler%process_model_list => cur_process_model
-    sub_flow_process_model_coupler%pm_ptr%ptr => cur_process_model
-    nullify(cur_process_model)
-  endif
-
-  ! Create Subsurface transport ProcessModel & ProcessModelCoupler
-  if (option%ntrandof > 0) then
-    cur_process_model => PMRTCreate()
-    cur_process_model%output_option => realization%output_option
-    cur_process_model%option => realization%option
-   
-    sub_tran_process_model_coupler => ProcessModelCouplerCreate()
-    sub_tran_process_model_coupler%option => option
-    sub_tran_process_model_coupler%process_model_list => cur_process_model
-    sub_tran_process_model_coupler%pm_ptr%ptr => cur_process_model
-    nullify(cur_process_model)
-  endif
-
-  ! Add the ProcessModelCouplers in a list
-  if (associated(surf_flow_process_model_coupler)) then
-    simulation%process_model_coupler_list => surf_flow_process_model_coupler
-    if (associated(sub_flow_process_model_coupler)) then
-      surf_flow_process_model_coupler%next => sub_flow_process_model_coupler
-      if (associated(sub_tran_process_model_coupler)) then
-        sub_flow_process_model_coupler%below => sub_tran_process_model_coupler
-      endif
-    else if (associated(sub_tran_process_model_coupler)) then
-      ! this will likely never be used (i.e. tranport without flow coupled
-      ! to surface)
-      surf_flow_process_model_coupler%next => sub_tran_process_model_coupler
-    endif
-  else if (associated(sub_flow_process_model_coupler)) then
-    simulation%process_model_coupler_list => sub_flow_process_model_coupler
-    if (associated(sub_tran_process_model_coupler)) then
-      sub_flow_process_model_coupler%below => sub_tran_process_model_coupler
-    endif
-  else
-    simulation%process_model_coupler_list => sub_tran_process_model_coupler
-  endif
-
-  ! For each ProcessModel, set:
-  ! - realization (subsurface or surface),
-  ! - stepper (flow/trans/surf_flow),
-  ! - SNES functions (Residual/Jacobain), or TS function (RHSFunction)
-  cur_process_model_coupler_top => simulation%process_model_coupler_list
-  do
-    if (.not.associated(cur_process_model_coupler_top)) exit
-    cur_process_model_coupler => cur_process_model_coupler_top
-    do
-      if (.not.associated(cur_process_model_coupler)) exit
-      cur_process_model => cur_process_model_coupler%process_model_list
-      do
-        if (.not.associated(cur_process_model)) exit
-        realization_class_ptr => realization
-        select type(cur_process_model)
-          class is (pm_richards_type)
-            call cur_process_model%PMRichardsSetRealization( &
-                                                         realization_class_ptr)
-            call cur_process_model_coupler%SetTimestepper(flow_stepper)
-            flow_stepper%dt = option%flow_dt
-          class is (pm_rt_type)
-            call cur_process_model%PMRTSetRealization(realization_class_ptr)
-            call cur_process_model_coupler%SetTimestepper(tran_stepper)
-            tran_stepper%dt = option%tran_dt
-          class is (pm_th_type)
-            call cur_process_model%PMTHSetRealization(realization_class_ptr)
-            call cur_process_model_coupler%SetTimestepper(flow_stepper)
-            flow_stepper%dt = option%flow_dt
-          class is (pm_thc_type)
-            call cur_process_model%PMTHCSetRealization(realization_class_ptr)
-            call cur_process_model_coupler%SetTimestepper(flow_stepper)
-            flow_stepper%dt = option%flow_dt
-#ifdef SURFACE_FLOW
-          class is (pm_surface_flow_type)
-            surf_realization_class_ptr => surf_realization
-            call cur_process_model%PMSurfaceFlowSetRealization(surf_realization_class_ptr)
-            call cur_process_model_coupler%SetTimestepper(surf_flow_stepper)
-            surf_flow_stepper%dt = option%surf_flow_dt
-#endif
-        end select
-
-        call cur_process_model%Init()
-        select type(cur_process_model)
-#ifdef SURFACE_FLOW
-          class is (pm_surface_flow_type)
-            call TSSetRHSFunction( &
-                            cur_process_model_coupler%timestepper%solver%ts, &
-                            cur_process_model%residual_vec, &
-                            PMRHSFunction, &
-                            cur_process_model_coupler%pm_ptr, &
-                            ierr)
-#endif
-          class default
-            call SNESSetFunction( &
-                           cur_process_model_coupler%timestepper%solver%snes, &
-                           cur_process_model%residual_vec, &
-                           PMResidual, &
-                           cur_process_model_coupler%pm_ptr,ierr)
-            call SNESSetJacobian( &
-                           cur_process_model_coupler%timestepper%solver%snes, &
-                           cur_process_model_coupler%timestepper%solver%J, &
-                           cur_process_model_coupler%timestepper%solver%Jpre, &
-                           PMJacobian, &
-                           cur_process_model_coupler%pm_ptr,ierr)
-        end select
-        cur_process_model => cur_process_model%next
-      enddo
-      cur_process_model_coupler => cur_process_model_coupler%below
-    enddo
-    cur_process_model_coupler_top => cur_process_model_coupler_top%next
-  enddo
-
-  ! If running with surface flow, add additional waypoints in waypoint-list
-  ! of synchronizer.
-#if 0  
-#ifdef SURFACE_FLOW
-  if (option%nsurfflowdof > 0) then
-    waypoint => realization%waypoints%first
-    dt_max = waypoint%dt_max
-    if(dt_max<1.d-40) then
-      option%io_buffer='waypoint%dt_max = 0.'
-      call printErrMsg(option)
-    endif
-    time = 0.d0
-    final_time = realization%waypoints%last%time
-    do
-      time = time + surf_realization%dt_coupling
-      if (time > final_time) exit
-      waypoint => WaypointCreate()
-      waypoint%time = time
-      waypoint%dt_max = dt_max
-      call WaypointInsertInList(waypoint,simulation%synchronizer%waypoints)
-   enddo
-  endif
-#endif
-#endif
-#endif
-
 end subroutine Init
 
 ! ************************************************************************** !
   type(output_option_type), pointer :: output_option
   type(uniform_velocity_dataset_type), pointer :: uniform_velocity_dataset
   class(dataset_base_type), pointer :: dataset
-  type(mass_transfer_type), pointer :: mass_transfer
+  type(mass_transfer_type), pointer :: flow_mass_transfer
+  type(mass_transfer_type), pointer :: rt_mass_transfer
   type(input_type), pointer :: input
 
   nullify(flow_stepper)
         nullify(coupler)        
       
 !....................
-      case ('MASS_TRANSFER')
-        mass_transfer => MassTransferCreate()
-        call InputReadWord(input,option,mass_transfer%name,PETSC_TRUE)
-        call InputDefaultMsg(input,option,'Mass Transfer name') 
-        call MassTransferRead(mass_transfer,input,option)
-        call MassTransferAddToList(mass_transfer, &
-                                   realization%mass_transfer_list)
-        nullify(mass_transfer)        
+      case ('FLOW_MASS_TRANSFER')
+        flow_mass_transfer => MassTransferCreate()
+        call InputReadWord(input,option,flow_mass_transfer%name,PETSC_TRUE)
+        call InputDefaultMsg(input,option,'Flow Mass Transfer name') 
+        call MassTransferRead(flow_mass_transfer,input,option)
+        call MassTransferAddToList(flow_mass_transfer, &
+                                   realization%flow_mass_transfer_list)
+        nullify(flow_mass_transfer)
+      
+!....................
+      case ('RT_MASS_TRANSFER')
+        rt_mass_transfer => MassTransferCreate()
+        call InputReadWord(input,option,rt_mass_transfer%name,PETSC_TRUE)
+        call InputDefaultMsg(input,option,'RT Mass Transfer name')
+        call MassTransferRead(rt_mass_transfer,input,option)
+        call MassTransferAddToList(rt_mass_transfer, &
+                                   realization%rt_mass_transfer_list)
+        nullify(rt_mass_transfer)
       
 !....................
       case ('STRATIGRAPHY','STRATA')

File src/pflotran/material.F90

View file
  • Ignore whitespace
             MaterialPropGetPtrFromList, &
             MaterialPropGetPtrFromArray, &
             MaterialPropConvertListToArray, &
+            MaterialAnisotropyExists, &
             MaterialPropertyRead
   
 contains
 
 ! ************************************************************************** !
 !
+! MaterialAnisotropyExists: Determines whether any of the material 
+!                           properties are anisotropic
+! author: Glenn Hammond
+! date: 07/11/13
+!
+! ************************************************************************** !
+function MaterialAnisotropyExists(material_property_list)
+
+  implicit none
+  
+  type(material_property_type), pointer :: material_property_list
+
+  PetscBool :: MaterialAnisotropyExists
+  
+  type(material_property_type), pointer :: cur_material_property
+    
+  MaterialAnisotropyExists = PETSC_FALSE
+  
+  cur_material_property => material_property_list
+  do 
+    if (.not.associated(cur_material_property)) exit
+    if (.not. cur_material_property%isotropic_permeability) then
+      MaterialAnisotropyExists = PETSC_TRUE
+      return
+    endif
+    cur_material_property => cur_material_property%next
+  enddo
+  
+end function MaterialAnisotropyExists
+
+! ************************************************************************** !
+!
 ! MaterialPropertyDestroy: Destroys a material_property
 ! author: Glenn Hammond
 ! date: 11/02/07

File src/pflotran/output_common.F90

View file
  • Ignore whitespace
             OutputGetVarFromArrayAtCoord, &
             OutputGetCellCenteredVelocities, &
             ConvertArrayToNatural, &
-            OutputFormatInt, &
-            OutputFormatDouble, &
             GetCellCoordinates, &
             GetVertexCoordinates, &
             OutputFilenameID, &
 
 ! ************************************************************************** !
 !
-! OutputFormatInt: Writes a integer to a string
-! author: Glenn Hammond
-! date: 01/13/12
-!
-! ************************************************************************** !  
-function OutputFormatInt(int_value)
-
-  implicit none
-  
-  PetscInt :: int_value
-  
-  character(len=MAXWORDLENGTH) :: OutputFormatInt
-
-  write(OutputFormatInt,'(1i12)') int_value
-  
-  OutputFormatInt = adjustl(OutputFormatInt)
-  
-end function OutputFormatInt
-
-! ************************************************************************** !
-!
-! OutputFormatDouble: Writes a double or real to a string
-! author: Glenn Hammond
-! date: 01/13/12
-!
-! ************************************************************************** !  
-function OutputFormatDouble(real_value)
-
-  implicit none
-  
-  PetscReal :: real_value
-  
-  character(len=MAXWORDLENGTH) :: OutputFormatDouble
-
-  write(OutputFormatDouble,'(1es13.5)') real_value
-  
-  OutputFormatDouble = adjustl(OutputFormatDouble)
-  
-end function OutputFormatDouble
-
-! ************************************************************************** !
-!
 ! GetCellCoordinates: Extracts coordinates of cells into a PetscVec
 ! author: Glenn Hammond
 ! date: 10/25/07

File src/pflotran/output_observation.F90

View file
  • Ignore whitespace
 ! ************************************************************************** !  
 subroutine OutputMassBalance(realization_base)
 
-  use Realization_class
+  use Realization_class, only : realization_type
   use Realization_Base_class, only : realization_base_type
   use Patch_module
   use Grid_module
   use Coupler_module
   use Utility_module
   
-  use Richards_module
-  use Mphase_module
-  use Immis_module
-  use Miscible_module
-  use TH_module
-  use THC_module
-  use THMC_module
-  use Reactive_Transport_module
-  use General_module
+  use Richards_module, only : RichardsComputeMassBalance
+  use Mphase_module, only : MphaseComputeMassBalance
+  use Immis_module, only : ImmisComputeMassBalance
+  use Miscible_module, only : MiscibleComputeMassBalance
+  use TH_module, only : THComputeMassBalance
+  use THC_module, only : THCComputeMassBalance
+  use THMC_module, only : THMCComputeMassBalance
+  use Reactive_Transport_module, only : RTComputeMassBalance
+  use General_module, only : GeneralComputeMassBalance
   
   use Global_Aux_module
   use Reactive_Transport_Aux_module

File src/pflotran/output_surface.F90

View file
  • Ignore whitespace
   use Surface_Realization_class
   use Grid_module
   use Option_module
+  use String_module
   
   implicit none
 
   
   
   string = 'ZONE T="' // &
-           trim(OutputFormatDouble(option%time/output_option%tconv)) // &
+           trim(StringFormatDouble(option%time/output_option%tconv)) // &
            '"'
   string2 = ''
   select case(tecplot_format)
       if ((surf_realization%discretization%itype == STRUCTURED_GRID).or. &
           (surf_realization%discretization%itype == STRUCTURED_GRID_MIMETIC)) then
         string2 = ', I=' // &
-                  trim(OutputFormatInt(grid%structured_grid%nx)) // &
+                  trim(StringFormatInt(grid%structured_grid%nx)) // &
                   ', J=' // &
-                  trim(OutputFormatInt(grid%structured_grid%ny)) // &
+                  trim(StringFormatInt(grid%structured_grid%ny)) // &
                   ', K=' // &
-                  trim(OutputFormatInt(grid%structured_grid%nz))
+                  trim(StringFormatInt(grid%structured_grid%nz))
       else
         string2 = 'POINT format currently not supported for unstructured'
       endif
       if ((surf_realization%discretization%itype == STRUCTURED_GRID).or. &
           (surf_realization%discretization%itype == STRUCTURED_GRID_MIMETIC)) then
         string2 = ', I=' // &
-                  trim(OutputFormatInt(grid%structured_grid%nx+1)) // &
+                  trim(StringFormatInt(grid%structured_grid%nx+1)) // &
                   ', J=' // &
-                  trim(OutputFormatInt(grid%structured_grid%ny+1)) // &
+                  trim(StringFormatInt(grid%structured_grid%ny+1)) // &
                   ', K=' // &
-                  trim(OutputFormatInt(grid%structured_grid%nz+1))
+                  trim(StringFormatInt(grid%structured_grid%nz+1))
       else
         string2 = ', N=' // &
-                  trim(OutputFormatInt(grid%unstructured_grid%num_vertices_global)) // &
+                  trim(StringFormatInt(grid%unstructured_grid%num_vertices_global)) // &
                   ', ELEMENTS=' // &
-                  trim(OutputFormatInt(grid%unstructured_grid%nmax))
+                  trim(StringFormatInt(grid%unstructured_grid%nmax))
         string2 = trim(string2) // ', ZONETYPE=FEQUADRILATERAL'
       endif
   
       if (variable_count > 4) then
         string3 = ', VARLOCATION=([4-' // &
-                  trim(OutputFormatInt(variable_count)) // &
+                  trim(StringFormatInt(variable_count)) // &
                   ']=CELLCENTERED)'
       else
         string3 = ', VARLOCATION=([4]=CELLCENTERED)'

File src/pflotran/output_tecplot.F90

View file
  • Ignore whitespace
   use Grid_module
   use Unstructured_Grid_Aux_module
   use Option_module
+  use String_module
   
   implicit none
 
   output_option => realization_base%output_option
 
   string = 'ZONE T="' // &
-           trim(OutputFormatDouble(option%time/output_option%tconv)) // &
+           trim(StringFormatDouble(option%time/output_option%tconv)) // &
            '"'
   string2 = ''
   select case(tecplot_format)
       if ((realization_base%discretization%itype == STRUCTURED_GRID).or. &
           (realization_base%discretization%itype == STRUCTURED_GRID_MIMETIC)) then
         string2 = ', I=' // &
-                  trim(OutputFormatInt(grid%structured_grid%nx)) // &
+                  trim(StringFormatInt(grid%structured_grid%nx)) // &
                   ', J=' // &
-                  trim(OutputFormatInt(grid%structured_grid%ny)) // &
+                  trim(StringFormatInt(grid%structured_grid%ny)) // &
                   ', K=' // &
-                  trim(OutputFormatInt(grid%structured_grid%nz))
+                  trim(StringFormatInt(grid%structured_grid%nz))
       else
         string2 = 'POINT format currently not supported for unstructured'
       endif  
       if ((realization_base%discretization%itype == STRUCTURED_GRID).or. &
           (realization_base%discretization%itype == STRUCTURED_GRID_MIMETIC)) then
         string2 = ', I=' // &
-                  trim(OutputFormatInt(grid%structured_grid%nx+1)) // &
+                  trim(StringFormatInt(grid%structured_grid%nx+1)) // &
                   ', J=' // &
-                  trim(OutputFormatInt(grid%structured_grid%ny+1)) // &
+                  trim(StringFormatInt(grid%structured_grid%ny+1)) // &
                   ', K=' // &
-                  trim(OutputFormatInt(grid%structured_grid%nz+1))
+                  trim(StringFormatInt(grid%structured_grid%nz+1))
       else if (grid%itype == IMPLICIT_UNSTRUCTURED_GRID) then
         string2 = ', N=' // &
-                  trim(OutputFormatInt(grid%unstructured_grid%num_vertices_global)) // &
+                  trim(StringFormatInt(grid%unstructured_grid%num_vertices_global)) // &
                   ', ELEMENTS=' // &
-                  trim(OutputFormatInt(grid%unstructured_grid%nmax))
+                  trim(StringFormatInt(grid%unstructured_grid%nmax))
         string2 = trim(string2) // ', ZONETYPE=FEBRICK'
       else
         string2 = ', N=' // &
-                  trim(OutputFormatInt(grid%unstructured_grid%nmax)) // &
+                  trim(StringFormatInt(grid%unstructured_grid%nmax)) // &
                   ', ELEMENTS=' // &
-                  trim(OutputFormatInt(grid%unstructured_grid%explicit_grid%num_elems))
+                  trim(StringFormatInt(grid%unstructured_grid%explicit_grid%num_elems))
         string2 = trim(string2) // ', ZONETYPE=FEBRICK'
       endif  
       
       else
         if (variable_count > 4) then
           string3 = ', VARLOCATION=([4-' // &
-                    trim(OutputFormatInt(variable_count)) // &
+                    trim(StringFormatInt(variable_count)) // &
                     ']=CELLCENTERED)'
         else
           string3 = ', VARLOCATION=([4]=CELLCENTERED)'

File src/pflotran/patch.F90

View file
  • Ignore whitespace
         vec_ptr(local_id) = vec_ptr2(grid%nL2G(local_id))
       enddo
       call GridVecRestoreArrayF90(grid,field%perm_zz_loc,vec_ptr2,ierr)
+    case(PERMEABILITY_XY)
+      call GridVecGetArrayF90(grid,field%perm_xy_loc,vec_ptr2,ierr)
+      do local_id=1,grid%nlmax
+        vec_ptr(local_id) = vec_ptr2(grid%nL2G(local_id))
+      enddo
+      call GridVecRestoreArrayF90(grid,field%perm_xy_loc,vec_ptr2,ierr)
+    case(PERMEABILITY_XZ)
+      call GridVecGetArrayF90(grid,field%perm_xz_loc,vec_ptr2,ierr)
+      do local_id=1,grid%nlmax
+        vec_ptr(local_id) = vec_ptr2(grid%nL2G(local_id))
+      enddo
+      call GridVecRestoreArrayF90(grid,field%perm_xz_loc,vec_ptr2,ierr)
+    case(PERMEABILITY_YZ)
+      call GridVecGetArrayF90(grid,field%perm_yz_loc,vec_ptr2,ierr)
+      do local_id=1,grid%nlmax
+        vec_ptr(local_id) = vec_ptr2(grid%nL2G(local_id))
+      enddo
+      call GridVecRestoreArrayF90(grid,field%perm_yz_loc,vec_ptr2,ierr)
     case(PHASE)
       call GridVecGetArrayF90(grid,field%iphas_loc,vec_ptr2,ierr)
       do local_id=1,grid%nlmax
       do local_id=1,grid%nlmax
         vec_ptr(local_id) = patch%imat(grid%nL2G(local_id))
       enddo
+    case(PROCESSOR_ID)
+      do local_id=1,grid%nlmax
+        vec_ptr(local_id) = option%myrank
+      enddo
     case default
       write(option%io_buffer, &
             '(''IVAR ('',i3,'') not found in PatchGetDataset'')') ivar

File src/pflotran/pmc_hydrogeophysics.F90

View file
  • Ignore whitespace
 
   use Realization_Base_class, only : RealizationGetDataset
   use Variables_module, only : PRIMARY_MOLALITY
+  use String_module
 
   implicit none
   
 #include "finclude/petscvec.h"
 #include "finclude/petscvec.h90"
+#include "finclude/petscviewer.h"
 
   class(pmc_base_type), pointer :: this
   class(pmc_hydrogeophysics_type), pointer :: pmc_hg
   PetscReal, pointer :: vec1_ptr(:), vec2_ptr(:)
   PetscErrorCode :: ierr
   PetscInt, save :: num_calls = 0
+  character(len=MAXSTRINGLENGTH) :: filename
+  PetscViewer :: viewer
 
   select type(pmc => this)
     class is(pmc_hydrogeophysics_type)
                                  PRIMARY_MOLALITY,ONE_INTEGER,0)
       call VecGetArrayF90(pmc%realization%field%work,vec1_ptr,ierr)
       call VecGetArrayF90(pmc_hg%solution_mpi,vec2_ptr,ierr)
-      vec2_ptr(:) = vec1_ptr(:) + num_calls
+      vec1_ptr(:) = vec1_ptr(:) + num_calls
+      vec2_ptr(:) = vec1_ptr(:)
       print *, 'PMC update to solution', vec2_ptr(16)
       call VecRestoreArrayF90(pmc%realization%field%work,vec1_ptr,ierr)
       call VecRestoreArrayF90(pmc_hg%solution_mpi,vec2_ptr,ierr)
+
+#if 0
+filename = 'pf_solution' // trim(StringFormatInt(num_calls)) // '.txt'
+call PetscViewerASCIIOpen(this%option%mycomm,filename,viewer,ierr)
+call VecView(pmc%realization%field%work,viewer,ierr)
+call PetscViewerDestroy(viewer,ierr)
+#endif
+
   end select
   num_calls = num_calls + 1
   

File src/pflotran/pmc_surface.F90

View file
  • Ignore whitespace
                                           pmc%surf_realization, &
                                           dt)
         case (TH_MODE)
-          !call SurfaceTHUpdateSubsurfSS(pmc%subsurf_realization, &
-          !                                pmc%surf_realization, &
-          !                                dt)
+          call SurfaceTHUpdateSubsurfSS(pmc%subsurf_realization, &
+                                          pmc%surf_realization,dt)
         end select
   end select
   

File src/pflotran/process_model_rt.F90

View file
  • Ignore whitespace
     call RealizationUpdateProperties(this%realization)
   endif
   
-  call MassTransferUpdate(this%realization%mass_transfer_list, &
+  call MassTransferUpdate(this%realization%rt_mass_transfer_list, &
                           this%realization%discretization, &
                           this%realization%patch%grid, &
                           this%realization%option)

File src/pflotran/process_model_top.F90

View file
  • Ignore whitespace
 module Process_Model_module
 
   use Process_Model_Base_class
-  use Process_Model_Richards_class
-  use Process_Model_RT_class
   
   implicit none
 
   use Option_module
   use Realization_class
   
-  use Richards_module
-
   implicit none
   
 #include "finclude/petscvec.h"

File src/pflotran/reactive_transport.F90

View file
  • Ignore whitespace
   call VecRestoreArrayReadF90(field%porosity_loc,porosity_loc_p,ierr)
 
   ! Mass Transfer
-  if (associated(realization%mass_transfer_list)) then
-    cur_mass_transfer => realization%mass_transfer_list
+  if (associated(realization%rt_mass_transfer_list)) then
+    cur_mass_transfer => realization%rt_mass_transfer_list
     do
       if (.not.associated(cur_mass_transfer)) exit
       call VecStrideScatter(cur_mass_transfer%vec,cur_mass_transfer%idof-1, &
   call VecRestoreArrayReadF90(field%volume, volume_p, ierr)
   
   ! Mass Transfer
-  if (associated(realization%mass_transfer_list)) then
-    cur_mass_transfer => realization%mass_transfer_list
+  if (associated(realization%rt_mass_transfer_list)) then
+    cur_mass_transfer => realization%rt_mass_transfer_list
     do
       if (.not.associated(cur_mass_transfer)) exit
       call VecGetArrayF90(cur_mass_transfer%vec, r_p, ierr)

File src/pflotran/realization.F90

View file
  • Ignore whitespace
   if (realization%option%ntrandof > 0) then
     call RealProcessTranConditions(realization)
   endif
-  if (associated(realization%mass_transfer_list)) then
-    call MassTransferInit(realization%mass_transfer_list, &
+  if (associated(realization%flow_mass_transfer_list)) then
+    call MassTransferInit(realization%flow_mass_transfer_list, &
                           realization%discretization, &
                           realization%option)
+    call MassTransferUpdate(realization%flow_mass_transfer_list, &
+                          realization%discretization, &
+                          realization%patch%grid, &
+                          realization%option)
   endif
- 
+  if (associated(realization%rt_mass_transfer_list)) then
+    call MassTransferInit(realization%rt_mass_transfer_list, &
+                          realization%discretization, &
+                          realization%option)
+    call MassTransferUpdate(realization%rt_mass_transfer_list, &
+                          realization%discretization, &
+                          realization%patch%grid, &
+                          realization%option)
+  endif
+
+
 end subroutine RealizationProcessConditions
 
 ! ************************************************************************** !
 ! currently don't use aux_vars, just condition for src/sinks
 !  call RealizationUpdateSrcSinks(realization)
 
-  call MassTransferUpdate(realization%mass_transfer_list, &
+  call MassTransferUpdate(realization%flow_mass_transfer_list, &
+                          realization%discretization, &
+                          realization%patch%grid, &
+                          realization%option)
+
+  call MassTransferUpdate(realization%rt_mass_transfer_list, &
                           realization%discretization, &
                           realization%patch%grid, &
                           realization%option)
     endif
   endif
   
-  ! add waypoints for mass transfer
-  if (associated(realization%mass_transfer_list)) then
-    cur_mass_transfer => realization%mass_transfer_list
+  ! add waypoints for flow mass transfer
+  if (associated(realization%flow_mass_transfer_list)) then
+    cur_mass_transfer => realization%flow_mass_transfer_list
+    do
+      if (.not.associated(cur_mass_transfer)) exit
+      if (associated(cur_mass_transfer%dataset%time_storage)) then
+        do itime = 1, cur_mass_transfer%dataset%time_storage%max_time_index
+          waypoint => WaypointCreate()
+          waypoint%time = cur_mass_transfer%dataset%time_storage%times(itime)
+          waypoint%update_conditions = PETSC_TRUE
+          call WaypointInsertInList(waypoint,realization%waypoints)
+        enddo
+      endif
+      cur_mass_transfer => cur_mass_transfer%next
+    enddo
+  endif  
+
+  ! add waypoints for rt mass transfer
+  if (associated(realization%rt_mass_transfer_list)) then
+    cur_mass_transfer => realization%rt_mass_transfer_list
     do
       if (.not.associated(cur_mass_transfer)) exit
       if (associated(cur_mass_transfer%dataset%time_storage)) then
   call ReactionDestroy(realization%reaction)
   
   call TranConstraintDestroy(realization%sec_transport_constraint)
-  call MassTransferDestroy(realization%mass_transfer_list)
+  call MassTransferDestroy(realization%flow_mass_transfer_list)
+  call MassTransferDestroy(realization%rt_mass_transfer_list)
   
   call WaypointListDestroy(realization%waypoints)
   

File src/pflotran/realization_base.F90

View file
  • Ignore whitespace
     type(field_type), pointer :: field
     type(flow_debug_type), pointer :: debug
     type(output_option_type), pointer :: output_option
-    type(mass_transfer_type), pointer :: mass_transfer_list
+    type(mass_transfer_type), pointer :: flow_mass_transfer_list
+    type(mass_transfer_type), pointer :: rt_mass_transfer_list
     
     type(reaction_type), pointer :: reaction
     
   nullify(realization_base%reaction)
 
   nullify(realization_base%patch)
-  nullify(realization_base%mass_transfer_list)
-  
+  nullify(realization_base%flow_mass_transfer_list)
+  nullify(realization_base%rt_mass_transfer_list)
+
 end subroutine RealizationBaseInit
 
 ! ************************************************************************** !

File src/pflotran/richards.F90

View file
  • Ignore whitespace
   use Discretization_module
   use Option_module
   use Logging_module
+  use Mass_Transfer_module, only : mass_transfer_type
 
   implicit none
 
   type(discretization_type), pointer :: discretization
   type(field_type), pointer :: field
   type(option_type), pointer :: option
+  type(mass_transfer_type), pointer :: cur_mass_transfer
   
   call PetscLogEventBegin(logging%event_r_residual,ierr)
   
 
   call PetscLogEventEnd(logging%event_r_residual,ierr)
 
+  ! Mass Transfer
+  if (associated(realization%flow_mass_transfer_list)) then
+    cur_mass_transfer => realization%flow_mass_transfer_list
+    do
+      if (.not.associated(cur_mass_transfer)) exit
+      call VecStrideScatter(cur_mass_transfer%vec,cur_mass_transfer%idof-1, &
+                            r,ADD_VALUES,ierr)
+      cur_mass_transfer => cur_mass_transfer%next
+    enddo
+  endif
+
 end subroutine RichardsResidual
 
 ! ************************************************************************** !

File src/pflotran/richards_mfd.F90

View file
  • Ignore whitespace
   use Option_module
   use Grid_module
   use Connection_module
+  use Mass_Transfer_module, only : mass_transfer_type  
 
   implicit none
 
   Vec :: xx
   Vec :: r
   type(realization_type) :: realization
+  PetscReal, pointer :: v_p(:)
   PetscViewer :: viewer
   PetscErrorCode :: ierr
 
 #if DASVYAT
 
-  PetscInt :: i, iface
+  PetscInt :: i, iface, icell
   PetscReal, pointer :: flow_xx_loc_p(:), r_p(:)
   PetscReal :: rnorm
   
   type(field_type), pointer :: field
   type(option_type), pointer :: option
   type(connection_set_type), pointer :: conn
+  type(mass_transfer_type), pointer :: cur_mass_transfer
   
   field => realization%field
   discretization => realization%discretization
    call RichardsResidualPatchMFDLP2(snes,xx,r,realization,ierr)
   !call RichardsCheckMassBalancePatch(realization)
 
+  ! Mass Transfer
+  if (associated(realization%flow_mass_transfer_list)) then
+    cur_mass_transfer => realization%flow_mass_transfer_list
+    call VecGetArrayF90(field%flow_r_loc_faces,r_p,ierr)
+    do
+      if (.not.associated(cur_mass_transfer)) exit
+      call VecGetArrayF90(cur_mass_transfer%vec,v_p,ierr)
+      do icell=1,realization%patch%grid%nlmax
+        r_p(icell+realization%patch%grid%ngmax_faces) = &
+            r_p(icell+realization%patch%grid%ngmax_faces) + v_p(icell)
+      enddo
+      call VecRestoreArrayF90(cur_mass_transfer%vec,v_p,ierr)
+      cur_mass_transfer => cur_mass_transfer%next
+    enddo
+    call VecRestoreArrayF90(field%flow_r_loc_faces,r_p,ierr)
+  endif  
+
    call VecScatterBegin( discretization%MFD%scatter_gtol_LP, field%flow_r_loc_faces, r, &
                                 ADD_VALUES,SCATTER_REVERSE, ierr)
    call VecCopy(field%flow_xx_loc_faces, field%work_loc_faces, ierr)

File src/pflotran/string.F90

View file
  • Ignore whitespace
             StringNull, &
             StringFindEntryInList, &
             StringSplit, &
-            StringSwapChar
+            StringSwapChar, &
+            StringFormatInt, &
+            StringFormatDouble
   
   interface StringCompare
     module procedure StringCompare1
   
 end function StringSplit
 
+! ************************************************************************** !
+!
+! StringFormatInt: Writes a integer to a string
+! author: Glenn Hammond
+! date: 01/13/12
+!
+! ************************************************************************** !  
+function StringFormatInt(int_value)
+
+  implicit none
+  
+  PetscInt :: int_value
+  
+  character(len=MAXWORDLENGTH) :: StringFormatInt
+
+  write(StringFormatInt,'(1i12)') int_value
+  
+  StringFormatInt = adjustl(StringFormatInt)
+  
+end function StringFormatInt
+
+! ************************************************************************** !
+!
+! StringFormatDouble: Writes a double or real to a string
+! author: Glenn Hammond
+! date: 01/13/12
+!
+! ************************************************************************** !  
+function StringFormatDouble(real_value)
+
+  implicit none
+  
+  PetscReal :: real_value
+  
+  character(len=MAXWORDLENGTH) :: StringFormatDouble
+
+  write(StringFormatDouble,'(1es13.5)') real_value
+  
+  StringFormatDouble = adjustl(StringFormatDouble)
+  
+end function StringFormatDouble
+
 end module String_module

File src/pflotran/surface_field.F90

View file
  • Ignore whitespace
 
     Vec :: area
     
-    Vec :: flux_subsurf_2_surf    ! MPI +ve value => Flow from subsurface to surface
+    Vec :: exchange_subsurf_2_surf   ! MPI +ve value => Flow from subsurface to surface
     Vec :: press_subsurf         ! MPI
 
     Vec :: Dq                    ! MPI
     Vec :: subsurf_yy            ! MPI
     Vec :: subsurf_zz            ! MPI
     Vec :: surf2subsurf_dist_gravity ! MPI
+    Vec :: surf2subsurf_dist ! MPI
+
+    ! For TH coupling
+    Vec :: temp_subsurf          ! MPI
+    Vec :: sat_ice               ! MPI
+    Vec :: ckwet                 ! MPI
+    Vec :: ckdry                 ! MPI
+    Vec :: ckice                 ! MPI
+    Vec :: th_alpha              ! MPI
+    Vec :: th_alpha_fr           ! MPI
 
     Vec :: subsurf_temp_vec_1dof ! MPI
     Vec :: subsurf_temp_vec_ndof ! MPI
   surface_field%flow_yy = 0
   surface_field%flow_accum = 0
   
-  surface_field%flux_subsurf_2_surf = 0
+  surface_field%exchange_subsurf_2_surf = 0
   surface_field%press_subsurf = 0
 
   surface_field%Dq = 0
   surface_field%subsurf_yy = 0
   surface_field%subsurf_zz = 0
   surface_field%surf2subsurf_dist_gravity = 0
+  surface_field%surf2subsurf_dist = 0
   
   surface_field%subsurf_temp_vec_1dof = 0
   surface_field%subsurf_temp_vec_ndof = 0
   surface_field%flowrate_inst = 0
   surface_field%flowrate_aveg = 0
 
+  surface_field%temp_subsurf = 0
+  surface_field%ckwet = 0
+  surface_field%ckdry = 0
+  surface_field%ckice = 0
+  surface_field%th_alpha = 0
+  surface_field%th_alpha_fr = 0
+  surface_field%sat_ice = 0
+
   SurfaceFieldCreate => surface_field
 
 end function SurfaceFieldCreate
 
   if (surface_field%area  /= 0) call VecDestroy(surface_field%area,ierr)
   
-  if (surface_field%flux_subsurf_2_surf /= 0) &
-    call VecDestroy(surface_field%flux_subsurf_2_surf,ierr)
-  if (surface_field%flux_subsurf_2_surf /= 0) &
+  if (surface_field%exchange_subsurf_2_surf /= 0) &
+    call VecDestroy(surface_field%exchange_subsurf_2_surf,ierr)
+  if (surface_field%exchange_subsurf_2_surf /= 0) &
     call VecDestroy(surface_field%press_subsurf,ierr)
   if (surface_field%press_subsurf /= 0) &
     call VecDestroy(surface_field%press_subsurf,ierr)
   if (surface_field%subsurf_zz/=0) call VecDestroy(surface_field%subsurf_zz,ierr)
   if (surface_field%surf2subsurf_dist_gravity/=0) &
     call VecDestroy(surface_field%surf2subsurf_dist_gravity,ierr)
+  if (surface_field%surf2subsurf_dist/=0) &
+    call VecDestroy(surface_field%surf2subsurf_dist,ierr)
 
   if (surface_field%flow_r /= 0) call VecDestroy(surface_field%flow_r,ierr)
   if (surface_field%flow_xx /= 0) call VecDestroy(surface_field%flow_xx,ierr)
   if (surface_field%flowrate_inst/=0) call VecDestroy(surface_field%flowrate_inst,ierr)
   if (surface_field%flowrate_aveg/=0) call VecDestroy(surface_field%flowrate_aveg,ierr)
 
+  if (surface_field%temp_subsurf /=0 ) call VecDestroy(surface_field%temp_subsurf,ierr)
+  if (surface_field%ckwet /=0 ) call VecDestroy(surface_field%ckwet,ierr)
+  if (surface_field%ckdry /=0 ) call VecDestroy(surface_field%ckdry,ierr)
+  if (surface_field%ckice /=0 ) call VecDestroy(surface_field%ckice,ierr)
+  if (surface_field%th_alpha /=0 ) call VecDestroy(surface_field%th_alpha,ierr)
+  if (surface_field%th_alpha_fr /=0 ) call VecDestroy(surface_field%th_alpha_fr,ierr)
+  if (surface_field%sat_ice /=0 ) call VecDestroy(surface_field%sat_ice,ierr)
+
 end subroutine SurfaceFieldDestroy
 
 end module Surface_Field_module

File src/pflotran/surface_flow.F90

View file
  • Ignore whitespace
                          option,vel,Res)
 
       patch%boundary_velocities(1,sum_connection) = vel
-#ifdef STORE_FLOWRATES
       patch%surf_boundary_fluxes(RICHARDS_PRESSURE_DOF,sum_connection) = Res(1)
-#endif
       if(abs(vel)>eps) max_allowable_dt = min(max_allowable_dt,dist/abs(vel)/4.d0)
       
       ff_p(local_id) = ff_p(local_id) + Res(1)/area_p(local_id)
         coupler_found = PETSC_TRUE
         
         call VecScatterBegin(dm_ptr%ugdm%scatter_bet_grids_1dof, &
-                             surf_field%flux_subsurf_2_surf, &
+                             surf_field%exchange_subsurf_2_surf, &
                              surf_field%subsurf_temp_vec_1dof, &
                              INSERT_VALUES,SCATTER_FORWARD,ierr)
         call VecScatterEnd(dm_ptr%ugdm%scatter_bet_grids_1dof, &
-                           surf_field%flux_subsurf_2_surf, &
+                           surf_field%exchange_subsurf_2_surf, &
                            surf_field%subsurf_temp_vec_1dof, &
                            INSERT_VALUES,SCATTER_FORWARD,ierr)
 
         enddo
         call VecRestoreArrayF90(surf_field%subsurf_temp_vec_1dof,vec_p,ierr)
 
-        call VecSet(surf_field%flux_subsurf_2_surf,0.d0,ierr)
+        call VecSet(surf_field%exchange_subsurf_2_surf,0.d0,ierr)
       endif
 
     endif
   PetscReal, pointer :: press_sub_p(:) ! Pressure [Pa]
   PetscReal, pointer :: icap_loc_p(:)
   PetscReal, pointer :: Dq_p(:)
-  PetscReal, pointer :: vol_p(:)
+  PetscReal, pointer :: mass_p(:)
   PetscReal, pointer :: area_p(:)
   PetscReal, pointer :: dist_p(:)
   
   call GridVecGetArrayF90(surf_grid,surf_field%flow_xx,hw_p,ierr)
   call GridVecGetArrayF90(surf_grid,surf_field%icap_loc,icap_loc_p,ierr)
   call GridVecGetArrayF90(surf_grid,surf_field%Dq,Dq_p,ierr)
-  call GridVecGetArrayF90(surf_grid,surf_field%flux_subsurf_2_surf,vol_p,ierr)
+  call GridVecGetArrayF90(surf_grid,surf_field%exchange_subsurf_2_surf,mass_p,ierr)
   call GridVecGetArrayF90(surf_grid,surf_field%surf2subsurf_dist_gravity,dist_p,ierr)
   call GridVecGetArrayF90(grid,surf_field%area,area_p,ierr)
 
           !v_darcy=0.d0
         endif
         
-        vol_p(local_id)=vol_p(local_id)+v_darcy*area_p(local_id)*option%surf_flow_dt
+        mass_p(local_id)=mass_p(local_id)+v_darcy*area_p(local_id)*option%surf_flow_dt
         !coupler%flow_aux_real_var(ONE_INTEGER,local_id)=v_darcy
         coupler%flow_aux_real_var(ONE_INTEGER,local_id)=0.d0
         hw_p(local_id) = hw_p(local_id) + v_darcy*option%surf_flow_dt
   
   call GridVecRestoreArrayF90(grid,surf_field%area,area_p,ierr)
   call GridVecRestoreArrayF90(surf_grid,surf_field%surf2subsurf_dist_gravity,dist_p,ierr)
-  call GridVecRestoreArrayF90(surf_grid,surf_field%flux_subsurf_2_surf,vol_p,ierr)
+  call GridVecRestoreArrayF90(surf_grid,surf_field%exchange_subsurf_2_surf,mass_p,ierr)
   call GridVecRestoreArrayF90(surf_grid,surf_field%Dq,Dq_p,ierr)  
   call GridVecRestoreArrayF90(surf_grid,surf_field%icap_loc,icap_loc_p,ierr)
   call GridVecRestoreArrayF90(surf_grid,surf_field%flow_xx,hw_p,ierr)

File src/pflotran/surface_init.F90

View file
  • Ignore whitespace
               call InputErrorMsg(input,option,'HDF5_WRITE_GROUP_SIZE','Group size')
             case('HYDROGRAPH')
               output_option%print_hydrograph = PETSC_TRUE
+            case('PROCESSOR_ID')
+              output_option%print_iproc = PETSC_TRUE
             case('FLOWRATES','FLOWRATE')
               mass_flowrate = PETSC_TRUE
               energy_flowrate = PETSC_TRUE

File src/pflotran/surface_realization.F90

View file
  • Ignore whitespace
                                      surf_field%work)
 
   call DiscretizationDuplicateVector(discretization,surf_field%flow_xx, &
-                                     surf_field%flux_subsurf_2_surf)
-  call DiscretizationDuplicateVector(discretization,surf_field%flow_xx, &
-                                     surf_field%press_subsurf)
+                                     surf_field%exchange_subsurf_2_surf)
 
   ! 1 degree of freedom, global
   call DiscretizationCreateVector(discretization,ONEDOF,surf_field%mannings0, &
   call DiscretizationDuplicateVector(discretization,surf_field%mannings0, &
                                      surf_field%subsurf_zz)
   call DiscretizationDuplicateVector(discretization,surf_field%mannings0, &
+                                     surf_field%surf2subsurf_dist)
+  call DiscretizationDuplicateVector(discretization,surf_field%mannings0, &
                                      surf_field%surf2subsurf_dist_gravity)
+  call DiscretizationDuplicateVector(discretization,surf_field%mannings0, &
+                                     surf_field%press_subsurf)
+  call DiscretizationDuplicateVector(discretization,surf_field%mannings0, &
+                                     surf_field%temp_subsurf)
+  call DiscretizationDuplicateVector(discretization,surf_field%mannings0, &
+                                     surf_field%sat_ice)
+  call DiscretizationDuplicateVector(discretization,surf_field%mannings0, &
+                                     surf_field%ckwet)
+  call DiscretizationDuplicateVector(discretization,surf_field%mannings0, &
+                                     surf_field%ckdry)
+  call DiscretizationDuplicateVector(discretization,surf_field%mannings0, &
+                                     surf_field%ckice)
+  call DiscretizationDuplicateVector(discretization,surf_field%mannings0, &
+                                     surf_field%th_alpha)
+  call DiscretizationDuplicateVector(discretization,surf_field%mannings0, &
+                                     surf_field%th_alpha_fr)
 
   ! n degrees of freedom, local
   call DiscretizationCreateVector(discretization,NFLOWDOF,surf_field%flow_xx_loc, &

File src/pflotran/surface_th.F90

View file
  • Ignore whitespace
   use Surface_Realization_class
   use Realization_Base_class
   use DM_Kludge_module
+  use Global_Aux_module
 
   implicit none
   
   type(surface_field_type),pointer    :: surf_field
   type(dm_ptr_type), pointer          :: dm_ptr
   type(connection_set_type), pointer  :: cur_connection_set
+  type(global_auxvar_type), pointer :: global_aux_vars(:)  
   
   Vec            :: destin_mpi_vec, source_mpi_vec
   PetscErrorCode :: ierr
-  PetscReal, pointer :: qsrc_p(:),vec_p(:)
+  PetscReal, pointer :: qsrc_p(:),press_p(:),temp_p(:)
   PetscInt :: local_id,iconn,sum_connection,ghosted_id
   PetscReal, pointer :: xx_loc_p(:)
   PetscReal, pointer :: xx_p(:)
   field      => realization%field
   surf_grid  => surf_realization%discretization%grid
   surf_field => surf_realization%surf_field
+  global_aux_vars => patch%aux%Global%aux_vars
 
-  dm_ptr => DiscretizationGetDMPtrFromIndex(realization%discretization,NFLOWDOF)
+  dm_ptr => DiscretizationGetDMPtrFromIndex(realization%discretization,ONEDOF)
   
   ! Update the surface BC
   coupler_list => patch%source_sinks
       ! Find the BC from the list of BCs
       if(StringCompare(coupler%name,'from_surface_ss')) then
 
+        ! Exchange subsurface PRESSURE
         call GridVecGetArrayF90(grid,field%flow_xx_loc,xx_loc_p, ierr)
-        call VecGetArrayF90(surf_field%subsurf_temp_vec_ndof,vec_p,ierr)
-
+        call VecGetArrayF90(surf_field%subsurf_temp_vec_1dof,press_p,ierr)
         do iconn=1,cur_connection_set%num_connections
           local_id = cur_connection_set%id_dn(iconn)
           ghosted_id = grid%nL2G(local_id)
 
           iend = ghosted_id*option%nflowdof
           istart = iend-option%nflowdof+1
-          vec_p((iconn-1)*option%nflowdof+1:iconn*option%nflowdof) = &
-            xx_loc_p(istart:iend)
+          press_p(iconn) = xx_loc_p(istart)
         enddo
-        call VecRestoreArrayF90(surf_field%subsurf_temp_vec_ndof,vec_p,ierr)
+        call VecRestoreArrayF90(surf_field%subsurf_temp_vec_1dof,press_p,ierr)
         call GridVecRestoreArrayF90(grid,field%flow_xx_loc,xx_loc_p, ierr)
 
         ! Scatter the data
-        call VecScatterBegin(dm_ptr%ugdm%scatter_bet_grids_ndof, &
-                        surf_field%subsurf_temp_vec_ndof,surf_field%press_subsurf, &
+        call VecScatterBegin(dm_ptr%ugdm%scatter_bet_grids_1dof, &
+                        surf_field%subsurf_temp_vec_1dof,surf_field%press_subsurf, &
                         INSERT_VALUES,SCATTER_FORWARD,ierr)
-        call VecScatterEnd(dm_ptr%ugdm%scatter_bet_grids_ndof, &
-                        surf_field%subsurf_temp_vec_ndof,surf_field%press_subsurf, &
+        call VecScatterEnd(dm_ptr%ugdm%scatter_bet_grids_1dof, &
+                        surf_field%subsurf_temp_vec_1dof,surf_field%press_subsurf, &
+                        INSERT_VALUES,SCATTER_FORWARD,ierr)
+
+        ! Exchange subsurface TEMPERATURE
+        call VecGetArrayF90(surf_field%subsurf_temp_vec_1dof,temp_p,ierr)
+        do iconn=1,cur_connection_set%num_connections
+          local_id = cur_connection_set%id_dn(iconn)
+          ghosted_id = grid%nL2G(local_id)
+
+          iend = ghosted_id*option%nflowdof
+          istart = iend-option%nflowdof+1
+          temp_p(iconn) = global_aux_vars(ghosted_id)%temp(1)
+        enddo
+        call VecRestoreArrayF90(surf_field%subsurf_temp_vec_1dof,temp_p,ierr)
+
+        ! Scatter the data
+        call VecScatterBegin(dm_ptr%ugdm%scatter_bet_grids_1dof, &
+                        surf_field%subsurf_temp_vec_1dof,surf_field%temp_subsurf, &
+                        INSERT_VALUES,SCATTER_FORWARD,ierr)
+        call VecScatterEnd(dm_ptr%ugdm%scatter_bet_grids_1dof, &
+                        surf_field%subsurf_temp_vec_1dof,surf_field%temp_subsurf, &
                         INSERT_VALUES,SCATTER_FORWARD,ierr)
 
       else
         coupler_found = PETSC_TRUE
         
         call VecScatterBegin(dm_ptr%ugdm%scatter_bet_grids_ndof, &
-                             surf_field%flux_subsurf_2_surf, &
+                             surf_field%exchange_subsurf_2_surf, &
                              surf_field%subsurf_temp_vec_ndof, &
                              INSERT_VALUES,SCATTER_FORWARD,ierr)
         call VecScatterEnd(dm_ptr%ugdm%scatter_bet_grids_ndof, &
-                           surf_field%flux_subsurf_2_surf, &
+                           surf_field%exchange_subsurf_2_surf, &
                            surf_field%subsurf_temp_vec_ndof, &
                            INSERT_VALUES,SCATTER_FORWARD,ierr)
 
         enddo
         call VecRestoreArrayF90(surf_field%subsurf_temp_vec_ndof,vec_p,ierr)
 
-        call VecSet(surf_field%flux_subsurf_2_surf,0.d0,ierr)
+        call VecSet(surf_field%exchange_subsurf_2_surf,0.d0,ierr)
       endif
 
     endif
       'boundary condition named from_surface_ss.'
     call printErrMsg(option)
   endif
-  call printErrMsg(option,'debugging in SurfaceTHUpdateSubsurfSS')
   
 end subroutine SurfaceTHUpdateSubsurfSS
 
   use Surface_Realization_class
   use Realization_Base_class
   use DM_Kludge_module
+  use Surface_TH_Aux_module
   
   implicit none
   
   type(surface_field_type),pointer    :: surf_field
   type(dm_ptr_type), pointer          :: dm_ptr
   type(connection_set_type), pointer  :: cur_connection_set
+  type(surface_global_auxvar_type), pointer :: surf_global_aux_vars(:)
+  type(Surface_TH_auxvar_type), pointer :: surf_aux_vars(:)
   
   Vec            :: destin_mpi_vec, source_mpi_vec
   PetscErrorCode :: ierr
-  PetscReal :: den          ! density      [kg/m^3]
-  PetscInt :: local_id,iconn
+  PetscReal :: den_surf_kg          ! density      [kg/m^3]
+  PetscReal :: den_sub_kg           ! density      [kg/m^3]
+  PetscReal :: den_aveg             ! density      [kg/m^3]
+  PetscInt :: local_id, ghosted_id, iconn
   
-  PetscReal, pointer :: hw_p(:)   ! head [m]
+  PetscReal, pointer :: xx_p(:)
   PetscReal, pointer :: press_sub_p(:) ! Pressure [Pa]
+  PetscReal, pointer :: temp_sub_p(:)  ! Temperature [C]
   PetscReal, pointer :: icap_loc_p(:)
   PetscReal, pointer :: ithrm_loc_p(:)
   PetscReal, pointer :: Dq_p(:)
-  PetscReal, pointer :: vol_p(:)
+  PetscReal, pointer :: exch_p(:)
   PetscReal, pointer :: area_p(:)
+  PetscReal, pointer :: dist_gravity_p(:)
   PetscReal, pointer :: dist_p(:)
+  PetscReal, pointer :: ckdry_p(:)
+  PetscReal, pointer :: ckwet_p(:)
+  PetscReal, pointer :: ckice_p(:)
+  PetscReal, pointer :: th_alpha_p(:)
+  PetscReal, pointer :: th_alpha_fr_p(:)
+  PetscReal, pointer :: sat_ice_p(:)
   
+  PetscReal :: hw
   PetscReal :: press_surf
   PetscReal :: dphi
   PetscReal :: press
   PetscReal :: v_darcy_max
   PetscReal :: gravity
   PetscReal :: press_up, press_dn
+  PetscReal :: k_eff_dn, k_eff_up
+  PetscReal :: Dk_eff
+  PetscReal :: Ke_up, Ke_fr
+  PetscReal :: dtemp
+  PetscReal :: Cw
+  PetscReal :: temp_half
+  PetscReal, parameter :: epsilon = 1.d-6
     
   PetscBool :: coupler_found = PETSC_FALSE
   PetscBool :: v_darcy_limit
   field      => realization%field
   surf_grid  => surf_realization%discretization%grid
   surf_field => surf_realization%surf_field
-
-  call density(option%reference_temperature,option%reference_pressure,den)
+  surf_global_aux_vars => surf_patch%surf_aux%SurfaceGlobal%aux_vars
+  surf_aux_vars => surf_patch%surf_aux%SurfaceTH%aux_vars
 
   call GridVecGetArrayF90(surf_grid,surf_field%press_subsurf,press_sub_p,ierr)
-  call GridVecGetArrayF90(surf_grid,surf_field%flow_xx_loc,hw_p,ierr)
+  call GridVecGetArrayF90(surf_grid,surf_field%temp_subsurf,temp_sub_p,ierr)
+  call GridVecGetArrayF90(surf_grid,surf_field%flow_xx_loc,xx_p,ierr)
   call GridVecGetArrayF90(surf_grid,surf_field%icap_loc,icap_loc_p,ierr)
   call GridVecGetArrayF90(surf_grid,surf_field%ithrm_loc,ithrm_loc_p,ierr)
   call GridVecGetArrayF90(surf_grid,surf_field%Dq,Dq_p,ierr)
-  call GridVecGetArrayF90(surf_grid,surf_field%flux_subsurf_2_surf,vol_p,ierr)
-  call GridVecGetArrayF90(surf_grid,surf_field%surf2subsurf_dist_gravity,dist_p,ierr)
+  call GridVecGetArrayF90(surf_grid,surf_field%exchange_subsurf_2_surf,exch_p,ierr)
+  call GridVecGetArrayF90(surf_grid,surf_field%surf2subsurf_dist_gravity,dist_gravity_p,ierr)
+  call GridVecGetArrayF90(surf_grid,surf_field%surf2subsurf_dist,dist_p,ierr)
   call GridVecGetArrayF90(grid,surf_field%area,area_p,ierr)
+  call GridVecGetArrayF90(surf_grid,surf_field%ckdry,ckdry_p,ierr)
+  call GridVecGetArrayF90(surf_grid,surf_field%ckwet,ckwet_p,ierr)
+  call GridVecGetArrayF90(surf_grid,surf_field%ckice,ckice_p,ierr)
+  call GridVecGetArrayF90(surf_grid,surf_field%th_alpha,th_alpha_p,ierr)
+  call GridVecGetArrayF90(surf_grid,surf_field%th_alpha_fr,th_alpha_fr_p,ierr)
+  call GridVecGetArrayF90(surf_grid,surf_field%sat_ice,sat_ice_p,ierr)
 
   ! Update the surface BC
   coupler_list => surf_patch%source_sinks
   coupler => coupler_list%first
+
   v_darcy_max=0.d0
   v_darcy_limit=PETSC_FALSE
   
         call printErrMsg(option)
       endif
       
-      do local_id=1,surf_grid%nlmax
+      cur_connection_set => coupler%connection_set
+
+      do iconn = 1, cur_connection_set%num_connections
+
+        local_id = cur_connection_set%id_dn(iconn)
+        ghosted_id = surf_grid%nL2G(local_id)
+
+        ! Compute densities:
+        call density(surf_global_aux_vars(ghosted_id)%temp(1), &
+                     option%reference_pressure,den_surf_kg)
+        call density(temp_sub_p(local_id),press_sub_p(local_id),den_sub_kg)
+        den_aveg = (den_surf_kg + den_sub_kg)/2.d0
 
         ! Exchange of water between surface-subsurface.
-        press_surf=hw_p((local_id-1)*option%nflowdof+1)*(abs(option%gravity(3)))*den &
-                   + option%reference_pressure
+        hw = xx_p((local_id-1)*option%nflowdof+1)
+        press_surf = hw*(abs(option%gravity(3)))*den_surf_kg + &
+                     option%reference_pressure
 
-        press_up = press_sub_p((local_id-1)*option%nflowdof+1)
+        press_up = press_sub_p(local_id)
         press_dn = press_surf
-        gravity = dist_p(local_id)*den
+        gravity = dist_gravity_p(local_id)*den_aveg
         
         dphi = press_up - press_dn + gravity
         
         v_darcy = Dq_p(local_id)*kr/visl*dphi
         if (v_darcy<=0.d0) then
           ! Flow is happening from surface to subsurface
-          if ( abs(v_darcy) > hw_p(local_id)/option%surf_flow_dt ) then
-            v_darcy = -hw_p(local_id)/option%surf_flow_dt
+          if ( abs(v_darcy) > xx_p(local_id)/option%surf_flow_dt ) then
+            v_darcy = -xx_p(local_id)/option%surf_flow_dt
             v_darcy_limit=PETSC_TRUE
           endif
+          temp_half = surf_global_aux_vars(ghosted_id)%temp(1) + 273.15d0
         else
           ! Exfiltration is occuring
+          temp_half = temp_sub_p(local_id) + 273.15d0
         endif
         
-        vol_p(local_id)=vol_p(local_id)+v_darcy*area_p(local_id)*option%surf_flow_dt
-        coupler%flow_aux_real_var(ONE_INTEGER,local_id)=v_darcy*0.d0
-        if(abs(v_darcy)>v_darcy_max) v_darcy_max=v_darcy
+        ! Mass flux
+        exch_p((local_id-1)*option%nflowdof+1) = &
+          exch_p((local_id-1)*option%nflowdof+1) + &
+          v_darcy*area_p(local_id)*option%surf_flow_dt
+        !coupler%flow_aux_real_var(ONE_INTEGER,local_id)=v_darcy
+        coupler%flow_aux_real_var(ONE_INTEGER,local_id)=0.d0
+        xx_p((local_id-1)*option%nflowdof+1) = &
+          xx_p((local_id-1)*option%nflowdof+1) + v_darcy*option%surf_flow_dt
+        if(abs(v_darcy)>v_darcy_max) v_darcy_max = v_darcy
 
-        ! Exchange of heat between surface-subsurface
+        ! Heat flux associated with mass flux
+        Cw = surf_aux_vars(ghosted_id)%Cw
+        exch_p((local_id-1)*option%nflowdof+2) = &
+          exch_p((local_id-1)*option%nflowdof+2) + &
+          den_aveg*v_darcy*temp_half*Cw*area_p(local_id)*option%surf_flow_dt
+
+        if (hw>0.d0) then
+          ! Exchange of heat between surface water--subsurface domain
+          Ke_up = (sat + epsilon)**th_alpha_p(local_id)
+          k_eff_up = ckdry_p(local_id) + &
+                      (ckwet_p(local_id) - ckdry_p(local_id))*Ke_up
+#ifdef ICE
+          Ke_fr_up = (sat_ice(local_id) + epsilon)**th_alpha_fr_p(local_id)
+          k_eff_up = ckwet_p(local_id)*Ke_up + ckice_p(local_id)*Ke_fr_up + &
+                     ckdry_p(local_id)*(1.d0 - Ke_up - Ke_fr_up)
+#endif
+          k_eff_dn = surf_aux_vars(ghosted_id)%k_therm
+
+          Dk_eff = k_eff_up*k_eff_dn/(k_eff_up*hw/2.d0 + &
+                                      k_eff_dn*dist_p(local_id))
+
+          dtemp = temp_sub_p(local_id) - surf_global_aux_vars(ghosted_id)%temp(1)
+          exch_p(local_id*option%nflowdof) = exch_p(local_id*option%nflowdof) + &
+            area_p(local_id)*Dk_eff*dTemp*option%surf_flow_dt
+        endif
 
       enddo
+    else
+
+      ! In the absence of standing water, the heat flux SS term for surface flow
+      ! is directly applied to subsurface domain.
+
+      if (associated(coupler%flow_aux_real_var) .and. &
+          coupler%flow_condition%energy_rate%itype == HET_ENERGY_RATE_SS) then
+
+        cur_connection_set => coupler%connection_set
+
+        do iconn = 1, cur_connection_set%num_connections
+
+          local_id = cur_connection_set%id_dn(iconn)
+          ghosted_id = surf_grid%nL2G(local_id)
+
+          hw = xx_p((local_id-1)*option%nflowdof+1)
+
+          if (hw==0.d0) then
+            exch_p(local_id*option%nflowdof) = exch_p(local_id*option%nflowdof) &
+              -coupler%flow_aux_real_var(TWO_INTEGER,local_id)*option%surf_flow_dt
+            ! NOTE: There is a -ve sign
+            !   coupler%flow_aux_real_var(TWO_INTEGER,:) : +ve value implies a
+            !       source of energy to surface water; while
+            !   exch_p(:) : +ve values implies transfer of energy from subsurface
+            !       to surface water.
+          endif
+
+        enddo
+
+      endif
 
     endif
-    
     coupler => coupler%next
   enddo
   
   call GridVecRestoreArrayF90(grid,surf_field%area,area_p,ierr)
-  call GridVecRestoreArrayF90(surf_grid,surf_field%surf2subsurf_dist_gravity,dist_p,ierr)
-  call GridVecRestoreArrayF90(surf_grid,surf_field%flux_subsurf_2_surf,vol_p,ierr)
+  call GridVecRestoreArrayF90(surf_grid,surf_field%surf2subsurf_dist_gravity,dist_gravity_p,ierr)
+  call GridVecRestoreArrayF90(surf_grid,surf_field%surf2subsurf_dist,dist_p,ierr)
+  call GridVecRestoreArrayF90(surf_grid,surf_field%exchange_subsurf_2_surf,exch_p,ierr)
   call GridVecRestoreArrayF90(surf_grid,surf_field%Dq,Dq_p,ierr)  
   call GridVecRestoreArrayF90(surf_grid,surf_field%icap_loc,icap_loc_p,ierr)
-  call GridVecRestoreArrayF90(surf_grid,surf_field%flow_xx_loc,hw_p,ierr)
+  call GridVecRestoreArrayF90(surf_grid,surf_field%flow_xx_loc,xx_p,ierr)
+  call GridVecRestoreArrayF90(surf_grid,surf_field%temp_subsurf,temp_sub_p,ierr)
   call GridVecRestoreArrayF90(surf_grid,surf_field%press_subsurf,press_sub_p,ierr)
+  call GridVecRestoreArrayF90(surf_grid,surf_field%ckdry,ckdry_p,ierr)
+  call GridVecRestoreArrayF90(surf_grid,surf_field%ckwet,ckwet_p,ierr)
+  call GridVecRestoreArrayF90(surf_grid,surf_field%ckice,ckice_p,ierr)
+  call GridVecRestoreArrayF90(surf_grid,surf_field%th_alpha,th_alpha_p,ierr)
+  call GridVecRestoreArrayF90(surf_grid,surf_field%th_alpha_fr,th_alpha_fr_p,ierr)
+  call GridVecRestoreArrayF90(surf_grid,surf_field%sat_ice,sat_ice_p,ierr)
 
 end subroutine SurfaceTHSurf2SubsurfFlux
 
   use Surface_Realization_class
   use Realization_Base_class
   use DM_Kludge_module
+  use TH_Aux_module
 
   implicit none
   
   type(surface_field_type),pointer    :: surf_field
   type(dm_ptr_type), pointer          :: dm_ptr
   type(connection_set_type), pointer  :: cur_connection_set
+  type(TH_parameter_type), pointer :: TH_parameter
   
   Vec            :: destin_mpi_vec, source_mpi_vec
   PetscErrorCode :: ierr
   PetscReal, pointer :: perm_yy_p(:)
   PetscReal, pointer :: perm_zz_p(:)
   PetscReal, pointer :: Dq_p(:)
+  PetscReal, pointer :: dist_gravity_p(:)
   PetscReal, pointer :: dist_p(:)
   PetscReal, pointer :: icap_loc_p(:)
   PetscReal, pointer :: ithrm_loc_p(:)
+  PetscReal, pointer :: ckwet_p(:)