added convergence checking routine

src/atmosphere/photochem_atmosphere.f90

@@ -18,6 +18,7 @@ module photochem_atmosphere
   contains
 
     ! photochem_atmosphere_rhs.f90
+    procedure, private :: prep_atm_background_gas
     procedure :: prep_atmosphere => prep_all_background_gas
     procedure :: right_hand_side_chem
     procedure :: production_and_loss
@@ -26,6 +27,7 @@ module photochem_atmosphere
 
     ! photochem_atmosphere_integrate.f90
     procedure :: evolve
+    procedure :: check_for_convergence
     procedure :: photochemical_equilibrium
     procedure :: initialize_stepper
     procedure :: step
@@ -67,6 +69,13 @@ module function create_Atmosphere(mechanism_file, settings_file, flux_file, atmo
 
     !~~ photochem_atmosphere_rhs.f90 ~~!
 
+    module subroutine prep_atm_background_gas(self, usol_in, usol, molecules_per_particle)
+      class(Atmosphere), target, intent(inout) :: self
+      real(dp), intent(in) :: usol_in(:,:)
+      real(dp), intent(out) :: usol(:,:)
+      real(dp), intent(out) :: molecules_per_particle(:,:)
+    end subroutine
+
     !> Given usol_in, the mixing ratios of each species in the atmosphere,
     !> this subroutine calculates reaction rates, photolysis rates, etc.
     !> and puts this information into self%wrk. self%wrk contains all the
@@ -137,6 +146,13 @@ module function evolve(self, filename, tstart, usol_start, t_eval, overwrite, er
       logical :: success !! If True, then integration was successful.
       character(:), allocatable, intent(out) :: err
     end function
+
+    !> Determines if integration has converged to photochemical equilibrium or not.
+    module function check_for_convergence(self, err) result(converged)
+      class(Atmosphere), target, intent(inout) :: self
+      character(:), allocatable, intent(out) :: err
+      logical :: converged
+    end function
 
     !> Integrates to photochemical equilibrium starting from self%var%usol_init
     module subroutine photochemical_equilibrium(self, success, err)

src/atmosphere/photochem_atmosphere_integrate.f90

@@ -363,7 +363,7 @@ module subroutine photochemical_equilibrium(self, success, err)
     type(PhotochemWrk), pointer :: wrk
 
     real(dp) :: tn
-    integer :: nsteps
+    logical :: converged
 
     dat => self%dat
     var => self%var
@@ -373,17 +373,23 @@ module subroutine photochemical_equilibrium(self, success, err)
     if (allocated(err)) return
 
     success = .true.
-    tn = 0
-    nsteps = 0
-    do while(tn < var%equilibrium_time)
+    do
       tn = self%step(err)
       if (allocated(err)) then
         ! If the step fails then exit
         success = .false.
         exit
       endif
-      nsteps = nsteps + 1
-      if (nsteps > var%mxsteps) then
+
+      ! Check for convergence
+      converged = self%check_for_convergence(err)
+      if (allocated(err)) then
+        success = .false.
+        exit
+      endif
+      if (converged) exit
+
+      if (wrk%nsteps > var%mxsteps) then
         ! If we did more steps then max steps, then exit
         success = .false.
         exit
@@ -398,6 +404,71 @@ module subroutine photochemical_equilibrium(self, success, err)
 
   end subroutine
 
+  module function check_for_convergence(self, err) result(converged)
+    use, intrinsic :: iso_c_binding
+    class(Atmosphere), target, intent(inout) :: self
+    character(:), allocatable, intent(out) :: err
+    logical :: converged
+
+    integer :: i,j,ind
+    type(PhotochemData), pointer :: dat
+    type(PhotochemVars), pointer :: var
+    type(PhotochemWrk), pointer :: wrk
+
+    dat => self%dat
+    var => self%var
+    wrk => self%wrk
+
+    converged = .false.
+    if (.not.c_associated(wrk%sun%cvode_mem)) then
+      err = "You must first initialize the stepper with 'initialize_stepper'"
+      return
+    endif
+
+    ! If we reach equilibrium time, then converged
+    if (wrk%tn > var%equilibrium_time) then
+      converged = .true.
+      return
+    endif
+
+    ! Now consider step history.
+
+    ! Can't do analysis on step 0
+    if (wrk%nsteps == 0) return 
+
+    ! Find index in history closest to time of interest. Note that this will only
+    ! consider a limited step history. We can not save all history.
+    ind = minloc(abs(wrk%t_history - var%conv_hist_factor*wrk%t_history(1)),1)
+
+    ! Can't be current time, so we will check the previous step if needed.
+    if (ind == 1) ind = 2 
+
+    ! Compute difference between current mixing ratios, and mixing ratios
+    ! at our index of interest.
+    do j = 1,var%nz
+      do i = 1,dat%nq
+        if (wrk%mix_history(i,j,1) > var%conv_min_mix) then
+          wrk%dmix(i,j) = abs(wrk%mix_history(i,j,1) - wrk%mix_history(i,j,ind))
+        else 
+          ! Ignore small mixing ratios
+          wrk%dmix(i,j) = 0.0_dp
+        endif
+      enddo
+    enddo
+
+    ! Maximum normalized change
+    wrk%longdy = maxval(abs(wrk%dmix/wrk%mix_history(:,:,1)))
+    ! Also consider that change over time
+    wrk%longdydt = wrk%longdy/(wrk%t_history(1) - wrk%t_history(ind))
+
+    ! Check for convergence
+    if (wrk%longdy < var%conv_longdy .and. wrk%longdydt < var%conv_longdydt) then
+      converged = .true.
+      return
+    endif
+
+  end function
+
   module subroutine initialize_stepper(self, usol_start, err)
     use, intrinsic :: iso_c_binding
     use fcvode_mod, only: CV_BDF, CV_NORMAL, CV_ONE_STEP, FCVodeInit, FCVodeSStolerances, &
@@ -465,6 +536,17 @@ module subroutine initialize_stepper(self, usol_start, err)
         wrk%sun%yvec(i) = var%lower_fix_mr(i)
       endif
     enddo
+
+    ! Load initial conditions into history vector
+    wrk%nsteps = 0
+    wrk%t_history = -1.0_dp
+    wrk%t_history(1) = 0.0_dp
+    wrk%mix_history = -1.0_dp
+    block
+      real(c_double), pointer :: usol_tmp(:,:)
+      usol_tmp(1:dat%nq,1:var%nz) => wrk%sun%yvec(1:var%neqs)
+      wrk%mix_history(:,:,1) = usol_tmp
+    endblock
 
     ! create SUNDIALS N_Vector
     wrk%sun%sunvec_y => FN_VMake_Serial(neqs_long, wrk%sun%yvec)
@@ -549,27 +631,57 @@ module subroutine initialize_stepper(self, usol_start, err)
   end subroutine
 
   module function step(self, err) result(tn)
-    use iso_c_binding, only: c_null_ptr, c_int, c_double, c_associated
-    use fcvode_mod, only: CV_ONE_STEP, FCVode
+    use iso_c_binding, only: c_null_ptr, c_int, c_double, c_associated, c_long
+    use fcvode_mod, only: CV_ONE_STEP, FCVode, FCVodeGetNumSteps
     class(Atmosphere), target, intent(inout) :: self
     character(:), allocatable, intent(out) :: err
     real(dp) :: tn
 
     integer(c_int) :: ierr
+    integer(c_long) :: nsteps_(1)
+    integer :: i, k
     real(c_double), parameter :: dum = 0.0_dp
     real(c_double) :: tcur(1)
+    type(PhotochemData), pointer :: dat
+    type(PhotochemVars), pointer :: var
+    type(PhotochemWrk), pointer :: wrk
+
+    dat => self%dat
+    var => self%var
+    wrk => self%wrk
 
-    if (.not.c_associated(self%wrk%sun%cvode_mem)) then
+    if (.not.c_associated(wrk%sun%cvode_mem)) then
       err = "You must first initialize the stepper with 'initialize_stepper'"
       return 
     endif
 
-    ierr = FCVode(self%wrk%sun%cvode_mem, dum, self%wrk%sun%sunvec_y, tcur, CV_ONE_STEP)
+    ierr = FCVode(wrk%sun%cvode_mem, dum, wrk%sun%sunvec_y, tcur, CV_ONE_STEP)
     if (ierr /= 0) then
       err = "CVODE step failed"
       return
     endif
     tn = tcur(1)
+
+    ! Update nsteps
+    ierr = FCVodeGetNumSteps(wrk%sun%cvode_mem, nsteps_)
+    wrk%nsteps = nsteps_(1)
+
+    ! Move over t and mix history
+    k = min(wrk%nsteps+1,size(wrk%t_history))
+    do i = k,2,-1
+      wrk%t_history(i) = wrk%t_history(i-1)
+      wrk%mix_history(:,:,i) =  wrk%mix_history(:,:,i-1)
+    enddo
+
+    ! Save current t and mix
+    wrk%t_history(1) = tn
+    block
+      real(c_double), pointer :: usol_tmp(:,:)
+      usol_tmp(1:dat%nq,1:var%nz) => wrk%sun%yvec(1:var%neqs)
+      call self%prep_atm_background_gas(usol_tmp, wrk%usol, wrk%molecules_per_particle)
+      wrk%mix_history(:,:,1) = wrk%usol
+    endblock
+
   end function
 
   module subroutine destroy_stepper(self, err)

src/atmosphere/photochem_atmosphere_rhs.f90

@@ -385,7 +385,7 @@ subroutine diffusion_coefficients(dat, var, den, mubar, &
 
   end subroutine
 
-  subroutine prep_atm_background_gas(self, usol_in, usol, molecules_per_particle)
+  module subroutine prep_atm_background_gas(self, usol_in, usol, molecules_per_particle)
     use photochem_common, only: molec_per_particle
     use photochem_const, only: small_real
     use photochem_enum, only: MixingRatioBC

src/evoatmosphere/photochem_evoatmosphere.f90

@@ -50,6 +50,7 @@ function temp_dependent_albedo_fcn(T_surf) result(albedo)
 
     !!! photochem_evoatmosphere_integrate.f90 !!!
     procedure :: evolve
+    procedure :: check_for_convergence
     procedure :: initialize_stepper
     procedure :: step
     procedure :: destroy_stepper
@@ -166,6 +167,12 @@ module function evolve(self, filename, tstart, usol_start, t_eval, overwrite, re
       ! Evolve atmosphere through time, and saves output in a binary Fortran file.
     end function
 
+    module function check_for_convergence(self, err) result(converged)
+      class(EvoAtmosphere), target, intent(inout) :: self
+      character(:), allocatable, intent(out) :: err
+      logical :: converged
+    end function
+
     !> Initializes an integration starting at `usol_start`
     module subroutine initialize_stepper(self, usol_start, err)      
       class(EvoAtmosphere), target, intent(inout) :: self