Version Release 2.3.3

This is the latest Release from private Mercurial Repository
Fixes subgridding issue and implements modified integration scheme for source term function.

docs/Version_notes.txt

@@ -0,0 +1,169 @@
+
+Tue Jun 25 14:05:42 MDT 2013
+VERS_V2.3.3 (most recent version) pulled from Repository
+   hg clone https://ancho.lanl.gov/fehm/hg/plumecalc
+
+
+VERS_V2.3.3 RELEASE NOTES:
+
+The manual V2.3.2 was not updated because the updates did not change
+how the code is run or basic functionality. Updated Figures can be
+seen in Complex_subgrid_figures.pdf. See discussion below. 
+
+-------
+Response to query for differences between V2.3.2 and V2.3.3
+From Zora Dash to Ian Nienhueser May 2011
+
+The latest version of PLUMECALC corrects the subgridding issue identified 
+by Dennis Fryar.  The problem was found when using subgridding in PLUMECALC V2.3.2 
+on a simple model that used FEHM streamline particle tracking with dispersion. 
+The problem of over estimating the concentration was identified as an
+implementation issue related to when the delay applied due to diffusion was
+computed. In the original implementation of mobile concentrations, the time
+delay was calculated before any subgridding occurred. This lead to an over
+calculation of concentration no matter what level of subgridding was used.
+This problem has been corrected by doing the time delay calculations after the
+subgridding. 
+See subgrid_problem.doc
+
+-------
+
+This version also implements a modified integration scheme for the source term 
+function. We identified problems in the numerical integration scheme 
+when small decay constants were used, especially if the source term 
+wasn't constant but was increasing or decreasing. I have rerun all of my 
+test cases and found that there are only minor differences in results 
+due to the modified integration scheme. 
+
+Test cases were rerun with this version V2.3.3.  Of particular interest 
+are the test results associated with the complex grid subgrid example documented 
+in Sections 6.3.2 and 6.3.3 of the Users Manual (Dash et al, 2011)
+Figures 32 to 39 (and specifically Figures 35, 37, 38 and 39).
+With exception of the case run with both dispersion and diffusion,
+as expected, there are no visually discernable differences.
+See Complex_subgrid_figures.pdf
+
+
+SUMMARY LOG MESSAGES:
+
+Diff		Rev	Age	Author	Log Message
+
+@21:c4da1fb88f54  	20 months	zvd	 Release PLUMECALC V2.3.3
+
+@20:7af4b2014665  	20 months	zvd	 Add exxecutables for current release. Release PLUMECALC V2.3.3
+
+@19:be83f4de2fec  	20 months	zvd	 Modified code to comput time delays after subgridding when subgridding is …
+
+@18:3c35ba99e9c5  	2 years	zvd	 Added tag Release PLUMECALC V2.3.2 for changeset 7069cc9e933e
+
+@17:7069cc9e933e  	2 years	zvd	 Add executables for current release.
+
+@16:f682ce9865d0  	2 years	zvd	 Changes incorporate all updates to PLUMECALC V2.2 that were made when …
+
+@15:6fa4af72072b  	2 years	zvd	 Reorganize repository and put source code into subdirectory src.
+
+@14:406f90e402a8  	5 years	moulton	 Removed old CVSROOT directory
+
+@13:574ece0454f0  	5 years	convert-repo	 update tags
+
+@12:2dca36560150  	6 years	zvd	 Modified to account for time=0.0 in sptr2 files with XYZ data.
+
+@11:5839af6b39e3  	6 years	zvd	 Convert volumes from m3 to liters so output will be in liters.
+
+@10:81e9e6ee5e4e  	7 years	zvd	 Added option to input secondary spacing instead of computing it.
+
+@9:cffc8d4121c3  	7 years	zvd	 Eliminated check for irregular input to reduce size of error/information …
+
+@8:77bfbf757eb6  	7 years	zvd	 Added concntration threshold criteria for favg pckd output option.
+
+@7:c8bf84cd94bd  	7 years	zvd	 Fixed formatting problem on Sun.
+
+@6:a9dfe9e5bb5c  	7 years	zvd	 Added a minimum concentration threshold for output.
+
+@5:5cd5b5848a5d  	7 years	zvd	 Fixed some minor bugs.
+
+@4:66c497124bb8  	7 years	zvd	 Fixed read_mdot_info so blank lines are not counted when reading flux …
+
+@3:0887c9cff063  	7 years	zvd	 added optional velocity correction factor
+
+@2:5389a1bc5073  	7 years	zvd	 fixed binary read
+
+@1:e64988a67e42  	7 years	zvd	 Add files for PLUMECALC V2.2 to repository
+
+@0:5cce78b732d5  	7 years	zvd	 initial checkin
+
+FULL LOG MESSAGES:
+
+
+Diff		Rev	Age	Author	
+
+@21:c4da1fb88f54  	20 months	zvd	
+Release PLUMECALC V2.3.3
+
+@20:7af4b2014665  	20 months	zvd	
+Add exxecutables for current release. Release PLUMECALC V2.3.3
+
+@19:be83f4de2fec  	20 months	zvd	
+Modified code to comput time delays after subgridding when subgridding is invoked. Modified decay calculations to better handle small decay constants (large half-lifes), numerical integration was blowing up due to formulation.
+
+@18:3c35ba99e9c5  	2 years	zvd	
+Added tag Release PLUMECALC V2.3.2 for changeset 7069cc9e933e
+
+@17:7069cc9e933e  	2 years	zvd	
+Add executables for current release.
+
+@16:f682ce9865d0  	2 years	zvd	
+Changes incorporate all updates to PLUMECALC V2.2 that were made when resident mobile concentration was implemented and regular plumecalc and subgrid plumecalc were merged into a single code. The base version for this work is PLUMECALC V2.3 and final changes are tagged PLUMECALC V2.3.2.
+
+@15:6fa4af72072b  	2 years	zvd	
+Reorganize repository and put source code into subdirectory src.
+
+@14:406f90e402a8  	5 years	moulton	
+Removed old CVSROOT directory
+
+@13:574ece0454f0  	5 years	convert-repo	
+update tags
+
+@12:2dca36560150  	6 years	zvd	
+Modified to account for time=0.0 in sptr2 files with XYZ data.
+
+@11:5839af6b39e3  	6 years	zvd	
+Convert volumes from m3 to liters so output will be in liters.
+
+
+@10:81e9e6ee5e4e  	7 years	zvd	
+Added option to input secondary spacing instead of computing it.
+
+@9:cffc8d4121c3  	7 years	zvd	
+Eliminated check for irregular input to reduce size of error/information
+file. (Can be reinstated for debugging if needed.)
+
+@8:77bfbf757eb6  	7 years	zvd	
+Added concntration threshold criteria for favg pckd output option.
+
+@7:c8bf84cd94bd  	7 years	zvd	
+Fixed formatting problem on Sun.
+
+@6:a9dfe9e5bb5c  	7 years	zvd	
+Added a minimum concentration threshold for output.
+
+@5:5cd5b5848a5d  	7 years	zvd	
+Fixed some minor bugs.
+
+@4:66c497124bb8  	7 years	zvd	
+Fixed read_mdot_info so blank lines are not counted when reading flux data.
+Added packed option for output of flux averaged concentrations over the
+entire domain.
+
+@3:0887c9cff063  	7 years	zvd	
+added optional velocity correction factor
+
+@2:5389a1bc5073  	7 years	zvd	
+fixed binary read
+
+@1:e64988a67e42  	7 years	zvd	
+Add files for PLUMECALC V2.2 to repository
+
+@0:5cce78b732d5  	7 years	zvd	
+initial checkin
+

src/Makefile.depends

@@ -5,7 +5,7 @@ calculate_favg.o: comgrid.o commdot.o comparttr.o comrock.o comsim.o comunits.o
 calculate_plume.o: comgrid.o commdot.o comparttr.o comrock.o comsim.o
 cell_setup.o: comparttr.o comparttr_sg.o
 compute_nxyz.o: comgrid.o comgrid_sg.o comparttr.o comparttr_sg.o comunits.o
-delayed_times.o: comgrid.o comparttr.o comrock.o
+delayed_times.o: comgrid.o comparttr.o compfrac.o comrock.o comsim.o comunits.o
 fill_id_parent_sg.o: comgrid.o comgrid_sg.o comparttr_sg.o comunits.o
 fill_ip_cell_to_cell_sg.o: comgrid.o comgrid_sg.o comparttr.o comparttr_sg.o
 fill_ip_cell_to_cell_sg.o: comunits.o
@@ -42,7 +42,7 @@ read_grid_info.o: comgrid.o comunits.o
 read_mdot_info.o: commdot.o comunits.o
 read_output_info.o: comgrid.o commdot.o comparttr.o comparttr_sg.o comsim.o
 read_output_info.o: comunits.o
-read_rock_info.o: comgrid.o comunits.o
+read_rock_info.o: comgrid.o comparttr_sg.o comunits.o
 read_sim_info.o: comgrid.o commdot.o comparttr.o comunits.o
 read_sptr_corn_dxyz_info.o: comgrid.o comparttr_sg.o comunits.o
 read_sptr_info.o: comparttr.o comparttr_sg.o comunits.o

src/Makefile.plumecalc

@@ -18,7 +18,7 @@
 # Edit FC,LIBRARY,FFLAGS for your environment and compiler 
 #***********************************************************************
 
-SRCDIR = /scratch/nts/zvd/PLUMECALC/SUBGRID/SRC_V2.3.x/
+SRCDIR = ../../
 WORKDIR = ./
 INCDIR = ./
 OBJDIR = ./
@@ -42,12 +42,12 @@ endif
 ifeq (${OPSYS}, Linux)
 OSTAG = beta_pgi${BIT}
 ifeq (${BIT},64)
-FC = /opt/pgi/linux86-64/6.2/bin/pgf90
-LIBRARY = -L/opt/pgi/linux86-64/6.2/lib
+FC = /opt/PGI/pgi_6.2-4/linux86-64/6.2/bin/pgf90
+LIBRARY = -L/opt/PGI/pgi_6.2-4/linux86-64/6.2/lib
 TARGET = k8-64e
 else
-FC = /opt/pgi/linux86/6.2/bin/pgf90
-LIBRARY = -L/opt/pgi/linux86/6.2/lib
+FC = /opt/PGI/pgi_6.2-4/linux86/6.2/bin/pgf90
+LIBRARY = -L/opt/PGI/pgi_6.2-4/linux86/6.2/lib
 TARGET = px
 endif
 DFLAG =
@@ -67,7 +67,7 @@ OBJECTS := $(patsubst ${SRCDIR}%.f,%.o,$(wildcard ${SRCDIR}*.f))
 all : plumecalc_sg
 
 plumecalc_sg : ${OBJECTS}  
-	${PURIFYIT} ${FC} ${FFLAGS} ${OBJECTS} ${LIBRARY} ${STATIC} -o plumecalc_sg_V2.3.2${OSTAG}${DATE}
+	${PURIFYIT} ${FC} ${FFLAGS} ${OBJECTS} ${LIBRARY} ${STATIC} -o plumecalc_sg_V2.3.3${OSTAG}${DATE}
 
 %.o : ${SRCDIR}%.f
 	${FC} ${DFLAG} ${FFLAGS} $< -c

src/calculate_favg.f

@@ -96,8 +96,9 @@ subroutine calculate_favg(previous_time,current_time)
          arraylen = iarrayend-iarraystart+1
 ! Integrate source function to detrmine solute that entered in 
 ! this time interval
-         if (prntvar(2))
-     2        call mdot_integral(current_time, iarraystart, iarrayend)
+         if (prntvar(4))
+     2        call mdot_integral(current_time, isource, iarraystart, 
+     3        iarrayend)
 
 !      FOR each particle in this source
          do i = istart, iend, istep

src/calculate_plume.f

@@ -96,8 +96,9 @@ subroutine calculate_plume(current_time)
 
 ! Integrate source function to detrmine solute that entered in 
 ! this time interval
-         if (prntvar(2))
-     2        call mdot_integral(current_time, iarraystart, iarrayend)
+         if (prntvar(4))
+     2        call mdot_integral(current_time, isource, iarraystart, 
+     3        iarrayend)
 
 !     FOR each particle in this source
          do i = istart, iend, istep
@@ -189,7 +190,7 @@ subroutine calculate_plume(current_time)
             else
                ps_correct = ps(cellno)
             end if
-            concentration(i) = concentration(i) / (ps(cellno)*sxliter)
+            concentration(i) = concentration(i) / (ps_correct*sxliter)
 c     Retardation now handled with mobile ratio
 c            concentration(i) = concentration(i) /
 c     2           (ps_correct * sxliter * rfac(cellno))

src/compfrac.f

@@ -73,6 +73,7 @@ module compfrac
       integer, allocatable, dimension(:,:,:) :: itf_curve
       real, allocatable, dimension(:,:,:) :: wt_curve
       real weight_factor
+      real*8 sigma_low, sigma_high, omega_low, omega_high
       real*8, allocatable, dimension(:) :: param1 
       real*8, allocatable, dimension(:) :: param2
       real*8, allocatable, dimension(:) :: param3
@@ -81,6 +82,7 @@ module compfrac
       integer, allocatable, dimension(:,:,:) :: param_density
 
       logical :: pfrac_read = .false.
+      logical :: error_flag = .false.
       integer :: ipartout = 0
 
       end module compfrac

src/comsim.f

@@ -83,6 +83,6 @@ module comsim
       integer noutnodes
       integer, allocatable :: ioutnode(:)
       integer, allocatable :: out_cell(:)
-      logical :: sparse, prntvar(3)
+      logical :: sparse = .false., prntvar(5) = .false.
 
       end module comsim

src/delayed_times.f

@@ -40,8 +40,12 @@ subroutine delayed_times
       use comparttr, only : npart, n_cells, time_packed, cell_packed,
      2     adv_packed
       use comgrid, only : n_grid_points
+      use compfrac, only : pfrac_read, error_flag, nump1, nump2,
+     2     sigma_low, sigma_high, omega_low, omega_high, param1, param2
       use comrock, only : rfac, vcf, itrc_diff, sigma_partial, 
      2     omega_partial, rd_frac, rseed, diffusion_model
+      use comsim, only : prntvar
+      use comunits, only : error_unit_number
       implicit none
 
       logical corrections
@@ -78,6 +82,7 @@ subroutine delayed_times
          if(rfac(i).ne.1..or.diffusion_model) then
 !     Set flag to denote corrections needed
             corrections = .true.
+            prntvar(5) = .true.
 !     Exit loop over each cell
             exit check_loop
          end if
@@ -108,18 +113,19 @@ subroutine delayed_times
                      sigma=sigma_partial(cell_packed(j))
                      omega=omega_partial(cell_packed(j))
                   else
-                     sigma=sigma_partial(cell_packed(j))/
-     2                    sqrt(delta_time)
-                     omega=omega_partial(cell_packed(j))*
-     2                    sqrt(delta_time)
+                     omega = omega_partial(cell_packed(j))*delta_time
+                     sigma = sigma_partial(cell_packed(j))
+c                     sigma=sigma_partial(cell_packed(j))/
+c     2                    sqrt(delta_time)
+c                     omega=omega_partial(cell_packed(j))*
+c     2                    sqrt(delta_time)
                   end if
                   par3v = 1.0
                   fm = 1
                   call time_delay(tflag,1,cell_packed(j),sigma,omega,
      2                 par3v,fm,rd_frac(itrc_diff(cell_packed(j))),
      3                 rseed,delta_time, conc_ret,time_diffused) 
-                  time_current = time_current + 
-     2                 time_diffused
+                  time_current = time_current + time_diffused
                else
                   time_current = time_current + 
      2                 rfac(cell_packed(j))*delta_time
@@ -141,7 +147,37 @@ subroutine delayed_times
          end do
 !        ENDFOR each particle
       end if
-!     ENDIF corrections are needed
-
+!     ENDIF corrections are needed      
+      if (pfrac_read) then
+         if (error_flag) then
+!        Issue warning for values outside of defined parameter space
+            write (error_unit_number, 8)
+            write (error_unit_number, 9)
+            if (sigma_low .lt.  param1(1)) then
+               write (error_unit_number, 14) param1(1)
+            end if
+            if (sigma_high .gt. param1(nump1)) then
+               write (error_unit_number, 15)
+            end if
+         else
+            write (error_unit_number, 10)
+         end if
+         write (error_unit_number, 11) param1(1), param1(nump1), 
+     2        param2(1), param2(nump2)
+         write (error_unit_number, 12)
+         write (error_unit_number, 11) sigma_low, sigma_high, 
+     2        omega_low, omega_high
+      end if
+
+ 8    format ('                ********** WARNING **********')
+ 9    format ('Diffusion model used values outside of defined ',
+     2     'parameter space:')
+ 10   format ('Diffusion model defined parameter space:')
+ 11   format ('min sigma = ', g16.9, ' max sigma = ', g16.9,
+     2     ' min omega = ', g16.9, ' max omega = ', g16.9) 
+ 12   format ('Range needed:')
+ 14   format ('Retention time set to sigma/omega for values less than',
+     2     1x, g16.9)
+ 15   format ('Retention time approximates error function solution')
       return
       end subroutine delayed_times