with start and stop times for sedimentation and erosion

src/DOUAR.f90

@@ -1051,8 +1051,10 @@ do while (istep.le.params%nstep)
             if (material0.eq.0.and.params%erosion) then
                call int_to_char(ic,2,is)
                call show_time (total,step,inc,1,'Erode surface '//ic(1:2)//'$')
-               if (current_time+tiny(current_time).ge. surface(is)%activation_time) then    
-                  call erosion (surface(is),olsf,is,params)
+               if (current_time+tiny(current_time).ge. surface(is)%activation_time) then
+                  if (current_time+tiny(current_time).ge.params%er_start.and.current_time.lt.params%er_end) then    
+                     call erosion (surface(is),olsf,is,params)
+                  endif
                else
                   if (iproc.eq.0) write(8,*) 'Surface ',is,' is attached to Surface0'
                   surface(is)%z=surface(1)%z   
@@ -1060,9 +1062,18 @@ do while (istep.le.params%nstep)
             elseif (material0.eq.0.and.params%sedimentation) then
                call int_to_char(ic,2,is)
                call show_time (total,step,inc,1,'Apply sediment to surface '//ic(1:2)//'$')
-               if (current_time+tiny(current_time).ge. surface(is)%activation_time) then    
-                  if (iproc.eq.0) write(*,*) 'Calling sedimentation with surface ',is
-                  call sedimentation (surface(is),olsf,is,params,current_time)
+               if (iproc.eq.0) write(*,*) 'current time: ',current_time
+               if (iproc.eq.0) write(*,*) 'sed_start: ',params%sed_start
+               if (iproc.eq.0) write(*,*) 'sed_end: ',params%sed_end
+               if (iproc.eq.0) write(*,*) 'activation time: ',surface(is)%activation_time
+               if (current_time+tiny(current_time).ge.surface(is)%activation_time) then
+                  if (iproc.eq.0) write(*,*) 'current time: ',current_time
+                  if (iproc.eq.0) write(*,*) 'sed_start: ',params%sed_start
+                  if (iproc.eq.0) write(*,*) 'sed_end: ',params%sed_end
+                  if (current_time+tiny(current_time).ge.params%sed_start.and.current_time.lt.params%sed_end) then    
+                     if (iproc.eq.0) write(*,*) 'Calling sedimentation with surface ',is
+                     call sedimentation (surface(is),olsf,is,params,current_time, params%dt/2.d0)
+                  endif
                else
                   if (iproc.eq.0) write(8,*) 'Surface ',is,' is attached to Surface0'
                   surface(is)%z=surface(1)%z   
@@ -1365,8 +1376,10 @@ do while (istep.le.params%nstep)
       ! We apply erosion
       if (material0.eq.0.and.params%erosion) then
          call show_time (total,step,inc,1,'Erode surface '//ic(1:2)//'$')
-         if (current_time+tiny(current_time).ge. surface(is)%activation_time) then    
-            call erosion (surface(is),olsf,is,params)
+         if (current_time+tiny(current_time).ge. surface(is)%activation_time) then 
+            if (current_time+tiny(current_time).ge.params%er_start.and.current_time.lt.params%er_end) then   
+               call erosion (surface(is),olsf,is,params)
+            endif
          else
             if (iproc.eq.0) write(8,*) 'Surface ',is,' is attached to Surface0'
             surface(is)%z=surface(1)%z   
@@ -1374,12 +1387,17 @@ do while (istep.le.params%nstep)
       elseif (material0.eq.0.and.params%sedimentation) then
 !         call int_to_char(ic,2,is)
          call show_time (total,step,inc,1,'Apply sedimentation to surface '//ic(1:2)//'$')
-         if (current_time+tiny(current_time).ge. surface(is)%activation_time) then    
-            call sedimentation (surface(is),olsf,is,params,current_time)
+         if (current_time+tiny(current_time).ge. surface(is)%activation_time) then
+            if (iproc.eq.0) write(*,*) 'current time: ',current_time
+            if (iproc.eq.0) write(*,*) 'sed_start: ',params%sed_start
+            if (iproc.eq.0) write(*,*) 'sed_end: ',params%sed_end
+            if (current_time+tiny(current_time).ge.params%sed_start.and.current_time.lt.params%sed_end) then
+               call sedimentation (surface(is),olsf,is,params,current_time, 0)
+            endif
          else
             if (iproc.eq.0) write(8,*) 'Surface ',is,' is attached to Surface0'
          surface(is)%z=surface(1)%z   
-            endif   
+         endif   
       end if
       if (iproc.eq.0) write (8,*) 'Min-Max z surf ',is,':',minval(surface(is)%z),maxval(surface(is)%z)
 

src/module_definitions.f90

@@ -279,6 +279,10 @@ module definitions
       double precision minp
       double precision pressure0
       double precision damp_factor
+      double precision er_start
+      double precision er_end
+      double precision sed_start
+      double precision sed_end
       double precision zaggrade_init
       double precision x_agg_start
       double precision x_agg_end

src/read_input_file.f90

@@ -985,6 +985,22 @@ if (iproc==0) call scanfile (params%infile,'zerosion',params%zerosion,ires)
 params%zerosion=params%zerosion/params%vex
 call mpi_bcast(params%zerosion,1,mpi_double_precision,0,mpi_comm_world,ierr)
 
+params%er_start=0.d0
+if (iproc==0) call scanfile (params%infile,'er_start',params%er_start,ires)
+call mpi_bcast(params%er_start,1,mpi_double_precision,0,mpi_comm_world,ierr)
+
+params%er_end=0.d0
+if (iproc==0) call scanfile (params%infile,'er_end',params%er_end,ires)
+call mpi_bcast(params%er_end,1,mpi_double_precision,0,mpi_comm_world,ierr)
+
+params%sed_start=0.d0
+if (iproc==0) call scanfile (params%infile,'sed_start',params%sed_start,ires)
+call mpi_bcast(params%sed_start,1,mpi_double_precision,0,mpi_comm_world,ierr)
+
+params%sed_end=0.d0
+if (iproc==0) call scanfile (params%infile,'sed_end',params%sed_end,ires)
+call mpi_bcast(params%sed_end,1,mpi_double_precision,0,mpi_comm_world,ierr)
+
 params%sedimentation_type=1
 if (iproc==0) call scanfile (params%infile,'sedimentation_type',params%sedimentation_type,ires)
 call mpi_bcast(params%sedimentation_type,1,mpi_integer,0,mpi_comm_world,ierr)
@@ -1180,6 +1196,10 @@ if (params%debug.gt.0 .and. iproc.eq.0) then
   write(*,'(a,e11.4)') shift//'e2d0 ',params%e2d0
   write(*,'(a,e11.4)') shift//'minp ',params%minp
   write(*,'(a,e11.4)') shift//'pressure0 ',params%pressure0
+  write(*,'(a,e11.4)') shift//'er_start ',params%er_start
+  write(*,'(a,e11.4)') shift//'er_end ',params%er_end
+  write(*,'(a,e11.4)') shift//'sed_start ',params%sed_start
+  write(*,'(a,e11.4)') shift//'sed_end ',params%sed_end
   write(*,'(a,e11.4)') shift//'zaggrade_init ',params%zaggrade_init
   write(*,'(a,e11.4)') shift//'x_agg_start ',params%x_agg_start
   write(*,'(a,e11.4)') shift//'x_agg_end ',params%x_agg_end
@@ -1428,6 +1448,10 @@ if (params%debug.gt.1) then
   write(threadinfo%Logunit,'(a32,e11.4)') 'e2d0 ',params%e2d0
   write(threadinfo%Logunit,'(a32,e11.4)') 'minp ',params%minp
   write(threadinfo%Logunit,'(a32,e11.4)') 'pressure0 ',params%pressure0
+  write(threadinfo%Logunit,'(a32,e11.4)') 'er_start ',params%er_start
+  write(threadinfo%Logunit,'(a32,e11.4)') 'er_end ',params%er_end
+  write(threadinfo%Logunit,'(a32,e11.4)') 'sed_start ',params%sed_start
+  write(threadinfo%Logunit,'(a32,e11.4)') 'sed_end ',params%sed_end
   write(threadinfo%Logunit,'(a32,e11.4)') 'zaggrade_init ',params%zaggrade_init
   write(threadinfo%Logunit,'(a32,e11.4)') 'x_agg_start ',params%x_agg_start
   write(threadinfo%Logunit,'(a32,e11.4)') 'x_agg_end ',params%x_agg_end

src/sedimentation.f90

@@ -16,7 +16,7 @@
 !------------------------------------------------------------------------------|
 !------------------------------------------------------------------------------|
 
-subroutine sedimentation (surface,olsf,is,params, current_time)
+subroutine sedimentation (surface,olsf,is,params,current_time,time_step)
 
 !------------------------------------------------------------------------------|
 !(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
@@ -62,6 +62,7 @@ type (octreelsf) olsf
 type (parameters) params
 integer is
 double precision current_time
+double precision time_step
 
 !------------------------------------------------------------------------------|
 !(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
@@ -138,7 +139,7 @@ x_prog_start=params%x_prog_start
 x_prog_end=params%x_prog_end
 prog_rate_u=params%prog_rate_u
 x_prog_length=params%x_prog_length
-dt=params%dt
+!dt=params%dt
 
 
 if (is.ne.1) then
@@ -247,10 +248,9 @@ endif
 ! change this to have a flag (0/1) so that current time has dt/2 added for
 !minpoint call
 
-
-zaggrade_current=zaggrade_init+aggrade_rate*current_time
-x_prog1=x_prog_start+prog_rate_u*current_time
-x_prog2=x_prog_end+prog_rate_u*current_time
+zaggrade_current=zaggrade_init+aggrade_rate*(current_time+time_step)
+x_prog1=x_prog_start+prog_rate_u*(current_time+time_step)
+x_prog2=x_prog_end+prog_rate_u*(current_time+time_step)
 !zaggrade_current2=zaggrade_init2+aggrade_rate*current_time
 !if (iproc.eq.0) write(*,*) 'zaggrade_current: ', zaggrade_current